3P101
有限長ジグザグナノチューブのNBMOの組み立てとHOMOの局在化（明治薬大）溝口 則幸

Construction of NBMO and Localization of HOMO
in Finite Length Zigzag Nanotubes(Meiji Pharmaceutical University) Noriyuki Mizoguchi
1 序論
非結合性分子軌道(NBMO)は共役化合物の電子状態の研究において重要な量である。ジグザグ端を持つナノグラフェンはゼロギャップ半導体であり，端の炭素原子に局在したHOMO（エッジ状態という）を持つ。このエッジ状態は局在スピンを有することから、ナノグラフェンはスピント
ロニクスへ応用が期待されている。ジグザグカーボンナノチューブは直径が増大すると、周期２でHOMO-LUMOギャップが振動することが知られている。このHOMO-LUMOギャップの振動はHOMOがNBMOであるときとNBMOでないときの繰り返しによるものである。
HMOレベルにおいて、sum ruleを使うことにより、計算なしに、NBMOの各サイトにおける係数を組み立てられる。ここでは、sum ruleを用いて、有限長(n,0)LジグザグナノチューブのNBMO
を組み立てる。次に、STO-3G* を基底関数とするDFT B3LYP計算により、構造最適化された有限長ジグザグナノチューブのHOMOの特徴を調べる。

2

ジグザグナノチューブの NBMO の組み立てと HOMO の完全局在化(n,0)LジグザグナノチューブのNBMOをsum ruleを用いて、組み立てた。ここで、添え字のLはチ

ューブの長さ（ベンゼンベルトの数）を表す。nが偶数のときと奇数のときでは、結果が異なり、nが偶数のときは、NBMOを得られたが、nが奇数のときは、得られなかった。nが偶数のときに得られたNBMOの係数は０でないないものは、ジグザグ端上のみにあり、一方のジグザグ端における係数はa,-a,a,-a,..,a,-aであり、他方のジグザグ端における係数はb,-b,b,-b,..,b,-bであった。このNBMOは二重に縮退していて、もう一つのNBMOの係数は一方のジグザグ端における係数は-b,b,-b,b,..,-b,bであり、他方のジグザグ端における係数はa,-a,a,-a,..,a,-aであった。ただし、これらの
NBMOは規格化されていないものである。その結果の一例を図１に示した。この結果は、ジグザ
グナノチューブの長さに依らない。
ジグザグナノチューブは交互炭化水素であるので、二重に縮退しているNBMOはHOMOとLUMOになる。n=evenのジグザグナノチューブのHOMOの空間分布はジグザグ端にある炭素原子のみに完全に局在化しており、ジグザグ端上にある炭素原子に均一の存在確率を持つ。これらの結果はジグザグナノチューブの太さと長さに依らない。
STO-3G* を基底関数とするB3LYPにより、構造最適化された有限長ジグザグナノチューブのHOMOの特徴を調べた。その結果、(2n,0)のジグザグナノチューブのHOMOはチューブのサイズ（長さや直径）に依らず、完全に局在化していることを確認した。図２を参照。

図１(2,0)3 ジグザグナノチューブの二重に縮退している NBMOHMOレベルにおいては、nが奇数のときのジグザグナノチューブはNBMOを持たないことをみた。NBMOがなければ、ジグザグナノチューブは端に局在した特異な電子状態を持たないのであろうか。nが奇数のとき、ジグザグナノチューブのHOMOはサイズ（長さや直径）が小さいときは非局在化している。図２を参照。 しかし、サイズが大きくなり、HOMOのエネルギーがNBMOのエネルギーに近づくと、HOMOの局在化は進み、ついには完全に局在化する。

図２(8,0)2ジグザグナノチューブのHOMOと(７,0) 3ジグザグナノチューブのHOMO

nが偶数のときと奇数のときでは、局在化の形が異なる。nが偶数のときは、HOMOは円周上のすべてのジグザグ端に一様な大きさで局在化しているが、nが奇数のときは、円周上のすべてのジグザグ端に局在化するのではなく、円周上の一部のジグザグ端のみに完全に局在化し、この部分のジグザグ端における係数の大きさは均一ではなく、中央部分で最大であり、中央から離れるに従い減少していく。その一例を図3に示した。

図3(19,0)3ジグザグナノチューブのHOMO

横からみたものと上からみたもの

3P102

凝集系における分子の励起状態と分子間相互作用(京大福井セ 1, 分子研 TCCI2)

長谷川淳也 1, 石村和也 2

Excited states and interactions of molecules in proteins andsolutions(Kyoto Univ.1. IMS TCCI2) J. Hasegawa1, K. Ishimura2【序】周辺環境に応じて分子の励起エネルギー準位が変化する現象は、ソルバトクロミズ

ムとしてよく研究されてきた。我々はこれまで、主にレチナール蛋白質について研究し、周辺環境の静電ポテンシャルが重要な役割を果たすことを示し、他方で周辺環境の電子状態も無視できない影響を与える系が存在することも見出した[1]。分子周辺のアミノ酸残基や溶媒分子を電子状態計算に含めた場合、計算される励起エネルギーは、点電荷モデルで記述した結果と比較して、bR では約 0.3 eV、MeOH では約 0.1 eV 程度減少する。この差 0.2 eV は、ソルバトクロミズムを議論する上では無視できない寄与である[2]。周辺環境の電子状態に由来する効果について解析するために、我々は分子軌道の局在化法を開発し、これを用いて励起状態波動関数を計算することで、アミノ酸残基やペプチド結合の軌道が関与する電子励起の寄与を評価した[3]。この解析により、凝集系励起状態の波動関数の構造を示唆する結果が得られている[3]。
本研究では、同局在化軌道を基底とし、特定の励起タイプを選択的に考慮できる CI プログラムを開発した。これを蛋白質や溶液中における分子の励起状態に応用することで、凝集系の励起状態の波動関数の構造を明らかにすることを目指す。
【計算】周辺環境の電子状態が色素の励起エネルギーに及ぼす効果を記述するために、溶媒分子（以降、蛋白質も含めて周辺環境を溶媒と総称する。）を含めた分子クラスター（図 1）を構築し、CIS 計算を行った。クラスターの外側の溶媒分子の効果は静電ポテンシャルとして考慮した。分子クラスターのサイズ拡大に伴い、第一励起エネルギーは減少するが、およそ 6 Å 半径（約 1000原子系）で、励起エネルギーの減少はほぼ収束した。分子軌道の局在化は、我々が提案した方法[4]を用い、得られた局在化軌道を一電子基底を用いて CIS 計算を行なった。実験結果と比較する
際は、ONIOM 法を用いて、SAC-CI 法により計算した色素の励起エネルギーに CIS 法による環境効果を外挿した結果を用いた。
(a) Deprotonated retinal
Schiff base in bR protein(b) Deprotonated retinal Schiff basein MeOH solution(c) s-trans acroleinin water

図 1. シリンダーモデルは色素分子から 6 Å 以内の原子を含む溶媒分子であり、CIS 計算を行った。
太線の領域の原子については SAC-CI 計算を行った。

【結果と考察】水溶液中における s-trans-acrolein の n-π*励起状態と π-π*励起状態について行った解析について説明する。実験では、気相中の水溶液中では n-π*励起で 0.3 eV の青方シフ
ト、π-π*励起で 0.5 eV もの赤方シフトを示す。解析の結果、蛋白質中や溶液中のレチナール
色素の結果[3]と同様に、周辺環境－色素間の CT 型励起配置と周辺環境における局所励起配置が波動関数に寄与している[5]。前者は色素波動関数が近傍のアミノ酸の領域へと非局在化する効果であり、後者は周辺環境の波動関数の分極効果を表している。また、周辺のアミノ酸間の CT 型励起配置の寄与は非常に小さいことが明らかになった。ONIOM 法を用いて、
SAC-CI 法により計算した s-trans-acrolein の励起エネルギーに CIS 法による環境効果を補正したところ、n-π*、π-π*ともスペクトルシフトに関する実験結果を良好に再現した。
Table 1 Calculated excitation energies of truncated CIS wave function. The number in parenthesis is thechange from “all single excitations” result.
Wave functiona

DPSB
in bR (M)3 Å model

DPSB
in MeOH3 Å model

PSB
in bR (BR)6 Å model

DPSB
in MeOH
6 Å model

Acroleinn–π* state
in water
6 Å model

Acroleinπ–π* state
in water
6 Å model

All singles

3.73 ( 0.00)

4.08 ( 0.00)

3.22 ( 0.00)

4.06 (0.00)

5.03 ( 0.00)

7.01 ( 0.00)− CTenv

3.74 (+0.01)

4.09 (+0.01)

3.24 (+0.02)

4.07 (+0.01)

5.03 (+0.00)

7.03 (+0.02)− CTenv – LEx

3.86 (+0.13)

4.14 (+0.06)

3.37 (+0.15)

4.13 (+0.07)

5.03 (+0.00)

7.08 (+0.06)− CTenv – CT

4.11 (+0.39)

4.44 (+0.36)

3.65 (+0.43)

4.43 (+0.37)

5.28 (+0.25)

7.46 (+0.45)a“− CTenv,” “– LEx,” and “– CT” denote that a CIS wave function that does not contain CT excitations between theenvironmental fragments, local excitation in the environmental fragments, and CT excitations from and to chromophores,
respectively.

また、周辺環境を含んだ励起状態の波動関数
において、電子励起を制限し、算出される励起エネルギーへの影響を検討した[5]。表１に観られるように、周辺のアミノ酸間の CT 励起配置
は、数において全波動関数の大部分を占めるにも関わらず、0.02 eV 程度の誤差しか与えない。
図２に、主配置と溶媒分子の励起配置の相互作用の目安となる積分値とその頻度を示す。局在化軌道を用いているため、アミノ酸間の CT 型
励起の積分値は高々1×10-3 程度であり、主配置との相互作用は小さいこと。数値的な意味で、
溶媒間 CT 励起は first-order interacting space から外れていることを示している。
【参考文献】

図 2.主配置 I→A に相互作用する溶媒分子の励起 b→j について、励起エネルギーに重要な| の頻寄与をする積分 V 2 |度をX
log |V|として表した。
s-transacrolein の(a)n-π*状態と(b)π-π*状態。●
は b, j が異なる溶媒分子間の場合。

[1] J. Hasegawa, K. J. Fujimoto, and H. Nakatsuji ChemPhysChem 12, 3106 (2011). [2] K. Fujimoto, K. Asai,and J. Hasegawa, Phys. Chem. Chem. Phys. 12, 13107 (2010). [3] J. Hasegawa, K. J. Fujimoto, and T. Kawatsu,
J. Chem. Theo. Comp. 8, 4452-4461(2012). [4] J. Hasegawa, T. Kawatsu, K. Toyoya, and K. Matsuda, Chem.
Phys. Letters, 508, 171 (2011). [5] J. Hasegawa, Chem. Phys. Letters, 571, 77 (2013).

3P103
QM/MM 法を用いた自由エネルギー安定構造の探索法(阪府大院・理、RIMED) ○麻田俊雄
Algorithm to search stable structures on free energy
surface using QM/MM method(Osaka Prefecture Univ.)

Toshio Asada[序] 自由エネルギー面の形状は、溶液中はもとより電子デバイスや生体内等に見られる分子集合体中の化学反応に対して、極めて重要な影響を及ぼすことが知られている。これまでに、分子集合体における化学反応を取り扱うことが可能な QM/MM 法、自由エネルギー勾配(FEG1) 法 お よ び 最 小 エ ネ ル ギ ー 経 路 を 見 つ け る 方 法 の ひ と つ で あ る nudged elasticband(NEB)法を用いることで、自由エネルギー面上における反応経路を最適化する方法を提案してきた 2。FEG を使うことで自由エネルギー構造の最適化が可能になる。しかしながら、単純な最適化から得られる自由エネルギー安定構造は、初期配置からみて最近傍の安定構造に収束してしまうことになる。これは、NEB 法と組み合わせた反応経路最適化においても初期経路に依存したローカルな局所安定経路となることを意味している。そこで自由エネルギーに基づくアンサンブルを生成するために自由エネルギー面上の MD シミュレーションを適用することを計画した。自由エネルギー面における拡張 MD シミュレーションを行い、自由エネルギー面上における広範囲な自由エネルギー安定構造の探索法を検討したので報告する。
[計算方法]
自由エネルギー面上の広範囲な構造探索が目的であるため、通常のエネルギーを自由エネ
ルギーに置き換えることでハミルトニアンを構築した。ハミルトニアン自体には物理的意味は持たせていない。
p i2+ G (r1 , , rN )i =1 2miN

H ' (p, r ) = ∑(1)

ここで、 p i は運動量、 mi は質量、 G は自由エネルギーである。一方、QM/MM 法においてQM 領域の FEG は∆GiQM (r1 , , rN ) =∂E (r1 , , rN )
∂riQM
MM

(2)

で表され、QM 原子にかかる力の MM 領域のアンサンブル平均に等しい。MM 領域の構造変化から QM 原子にかかる力を得るためには QM 領域の分極を含む静電相互作用を考慮する必要がある。そこで、FEG 計算高速化のため MM 領域の時間発展に起因する QM 領域の分極を、Charge Response Kernel(CRK)法により近似した。QM 計算は RHF/6-31G(d)レベルで行い、MM 領域には Amber99 力場を用いた。計算対象とする分子は、すでに多くの研究がなされ水中の自由エネルギーマップが報告されているアラニンジペプチドとし、構造変化の自由度として Ramachandran plot に用いられる 2 つの二面角に着目した(図 1)。
a)

φ = 180.0
ψ = 180.0

φ = −131.8
ψ = −11.1b)

図 1 アラニンジペプチド分子の a) 二面角と b) 水中の自由エネルギー最適化構造。
ここで、 φ = ∠C − C A − N − C ' および φ = ∠N − C A − C − N '[結果と考察]である。

MM 領域の MD 時間(QM 原子座標は固定)と各 QM 原子にかかる FEG 平均

値の変化について図 2(a)に示す。MM 領域の時間発展につれて QM 原子にかかる平均力は収束する傾向にあり、本系では、QM 構造ごとに最低でも 3psec 以上の MM MD シミュレーシ
ョンが必要になると結論できる。MM 領域の時間発展のつど、QM 領域の分極を得るために
(a)

(b)

15

force / kcal.mol-1A

10

β構造

5
0

α 構造
R

-5
-10
-15
0

1000

2000

3000

4000

5000

time / fsec

図2(a) MM MD シミュレーションの時間に対する各 QM 原子の FEG の収束性。
(b)自由エネルギー面上の QM 領域の MD シミュレーションの軌跡。

QM 計算を行うことは時間的に困難であることから、CRK 法を用いることで実用的な時間内に平均力を得られるようにした。また、自由エネルギー面上の QM 領域の MD シミュレーシ
ョンから得られた構造探索の結果を図 2(b)に示す。自由エネルギー面上の構造探索が良好に進行していることがわかる。詳細は当日に発表する。

[参考文献] 1. M.Nagaoka, N.Okuyama-Yoshida, and T.Yamabe, J.Phys.Chem.A. 102, 8202 (1998).
2. N.Takenaka, Y.Koyano, Y.Kitamura, T.Asada, and M.Nagaoka, Theo.Chem.Acc., 130, 215-226(2011).

3P104
P104
Fragment-DFT 計算による電子カップリングの評価（安田女子大・薬，筑波大院・化）○下堂

靖代，守橋

健二

Estimation of electronic coupling values in Fragment-DFT calculation( Yasuda Women’s Univ., Univ. of Tsukuba ) ○Yasuyo Shimodo，Kenji Morihashi【序論】
電荷輸送 (Charge Transport, CT) は最も基本的な化学反応の一つであるが，生命現象の素過程や有機デバイスの基本素子を成す，きわめて重要な反応である．よって，CT 機構を詳細に理解することは，生命科学，材料科学，両分野において有用である．
分子内あるいは分子間 CT 機構は，ホッピングモデルで説明され，速度論的には Marcus
理論によって議論される．ここで，CT 速度定数を決定付ける因子の一つは，CT 前後の系の電子状態|ηa 〉と|ηb 〉との電子的な相互作用の強さを表す電子カップリング行列要素 Hab=〈ηa|Ĥ|ηb 〉であり，これは反応の頻度因子に寄与する．したがって，Hab 値の予測は，有機電子材料の改良や開発，生体分子系においては配列特異的あるいはその高効率的な CT 機構の解明につながることが期待される．
我々は，生体分子や分子集合体などの巨大分子系において Hab を効率的に評価する計算手法を確立し，CT 機構の解明に役立てることを目指している．
【方法】
Hab の評価法として最も簡便なものは，一電子近似を適用した Energy Splitting in Dimer(ESD) 法であり，Hab 値は電荷の授受に関わる分子の軌道エネルギーεのみで与えられる[1]．
しかしながら，ESD 法は，分子間配向の違いによって生じるサイトエネルギーを考慮しないと定性的にも誤った結果を与えることが報告されている[2]．また，ESD 法を分子内 CT に応用するためには，系のモデル化を必要とする．ここで，我々は，フラグメント近似を導入した解法は，ESD 法と対応していること，周囲の環境場を考慮することからサイトエネルギーの評価が可能であること，フラグメント化により分子内 CT にも対応できることに着目した．

フラグメント近似の下では，Hab はフラグメント I と J の間のカップリングとして HIJ で表し，正孔移動および電子移動についてそれぞれ次の式で書き改められる．

1
2
1
H IJelec =2
H IJhole =(εIJn IJ / 2(εIJn IJ / 2 + 2) (2− ε nIJIJ / 2 −1 − ε nI I / 2 − ε nJJ / 2

) (
2

)

2

− ε nIJIJ / 2 +1 − ε nI I / 2 +1 − ε nJJ / 2 +1)2

ただし，nI および nJ はそれぞれフラグメント I と J に属する電子の数で，nIJ = nI + nJ である．
平方根内の第二項のカッコ内がいわゆるサイトエネルギーである．
これを非経験的量子化学計算プログラム ABINIT-MP (+DFT version) へ導入し，フラグメ
ント計算の過程で HIJ を評価することを可能にした．これにより，フラグメント計算を一回実行するだけで，全てのフラグメント間においてサイトエネルギーを考慮した HIJ 値を得ることができる．本研究では，この手法を Oligoacene 結晶中における分子間 CT の HIJ 計算に応用し，その有効性を検証した．
【結果】
表には，pentacene 結晶において 1 分子に隣接する 8 分子を含む[1+8]分子系，隣接する 1分子のみを含む[1+1]分子系について，FDFT B3LYP/6-31G(d) 計算により予測された HIJ 値を
分子配向ごとにまとめた．ここで，μは電荷移動度であり，HIJ 値から Deng らの手法[3]に従って求めた．他の Oligoacene 結晶におけるμについても，実験値に対応した値が得られた．
HIJ 値によると，Oligoacene 結晶では，スタッキングした分子間 (P) よりも傾いた配向の分

子間 (T1, T2) での CT が優位であった．また，従来の ESD 法で対象となる[1+1]分子系に比べ，[1+8]分子系，さらに大きな系でサイトエネルギーの影響が小さくなることは，結晶中における等価な分子環境を再現するものである．したがって，本研究で導入した HIJ 評価法は，分子間 CT について合理的な結果を与えることが示された．
表：FDFT B3LYP/6-31G(d) による pentacene 結晶中における HIJ およびμの計算値※
正孔移動

HIJ [eV]

T1
T2
P
L

μ[cm2/Vs]
Exptl.
※

電子移動[1+8]分子系[1+1]分子系[1+8]分子系[1+1]分子系

0.140

0.184

0.156

0.196

(0.140)(0.194)(0.156)(0.204)

0.155

0.207

0.156

0.208

(0.156)(0.217)(0.156)(0.219)

0.038

0.038

0.049

0.049

(0.038)(0.038)(0.050)(0.049)

0.008

0.000

0.008

0.003

(0.065)(0.000)(0.062)(0.003)

8.60

0.156

(9.42)(0.156)3, 5－7( )内はサイトエネルギーを考慮していない値[1] M. D. Newton, et al., Chem. Rev. 91,
91 767 (1991).
[2] E. F. Valeeve, et al., J. Am. Chem. Soc. 128,
128 9882 (2006).
[3] W. Deng, et al., J. Phys. Chem. B 108,108 8614 (2004).

NA3P105
遷移金属原子を含んだ巨大系のための IMiC-MO 法の拡張
(岐阜大工) ○宇田川太郎, 酒井章吾
Development of IMiC-MO method for transition metal complexes(Gifu University) ○Taro Udagawa, Shogo Sakai【序論】

遷移金属原子は生体内において生体高分子と錯体を形成し、重要な生理機能を担っていることが知られている。また近年では、例えばハロゲン架橋錯体のような、多彩な興味深い物
性を示す無機錯体にも注目が集まっており、様々な分野で盛んに研究されている。これら遷移金
属錯体の興味深い物性の多くは、d/f 電子の存在により引き起こされるが、それらは同時に理論化
学計算による取り扱いを複雑にしている。また、生体内高分子錯体や集積型複核錯体などのような巨大系に対しては、従来の分子軌道法による直接的な取り扱いは現実的ではない。
近年我々は、系を分割することで系のエネルギーおよび各原子にかかる力を効率良く求めることが可能な Integrated Multicenter MO (IMiC-MO)法を開発してきた[1]。この手法では、全系を目的
分子(T or t)、隣接分子(A or a)、環境分子(e) (但し、大文字は IMiC 計算で取り扱う際、全量子力学
的に、小文字は点電荷として取り扱うことを意味する)の 3 つの領域に分割することで、効率良く巨大系を取り扱う。近年では様々な分割的手法が提案されているが、それらの手法を遷移金属錯体に適用した例は少なく、特に電子状態が複雑となる少数の d 電子を有する系に適用した報告は無い。そこで本研究では、簡単な単核および複核遷移金属錯体をモデルとして、B3LYP 法を用いた IMiC-MO 法により金属錯体を取り扱った。具体的には、遷移金属錯体を中心遷移金属と配位子とに分割して取り扱うことによる、中心遷移金属原子周囲の電子状態の変化を解析し、全系のエ
ネルギーの再現性、および最適化構造について検討した。

【計算方法】 遷移金属錯体のモデルとして、Ti, V, Cr, Ptの単核および複核錯体を用いた。構造はそれぞれの錯体の最も安定な電子状態において最適化した。全ての計算
は(U)B3LYP 法で行い、遷移金属錯体に LanL2DZ、その他の原子には 6-31G を用いた。

【結果・考察】 図 1 に球状電荷解析(SCA)法[2]により求めた、種々の錯体における電子密度差 Δρ0(R)を示した。
ここで Δρ0(R)は、遷移金属原子を中心とした半径 R の球
表面の電子密度の錯体形成による変化を表している。図より、遷移金属原子が異なっても酸化数の同じ錯体では良く似た挙動を示すことがわかった。これは中心金属原子が正に帯電することで、配位子の電子をより強く引きつけるためであると考えられる。また、配位子からの電子移動量に相当する Δρ0(R)は、IMiC-MO 法による全エネルギー評価の精度に直接影響すると考えられる。実際、エネルギー誤差は例えば[Ti(II)2Cl6]2-錯体では 5.1 kcal/mol と小さいのに対して、[Ti(IV)2Cl6]2+錯体では 233.9 kcal/mol と非常に大きかった。
次に分割による大きなエネルギー誤差の原因を探るため、NAO の占有数を解析した。IMiC-MO法では分割した個々の原子を目的原子とし、その全てについて(i) TA(e)計算 (ii) T(a)(e)計算 (iii)(t)A(e)計算の 3 種類を計算し、全系のエネルギーを評価する[1]。遷移金属原子が目的原子の場合、
TA(e)および T(a)(e)計算に遷移金属原子が含まれる。全系中および、TA(e)、T(a)(e)計算中における遷移金属原子の占有数を表 1 に示した。TA(e)計算では遷移金属原子は価数に関係なく全系中の電子状態を良く再現しているのに対し、T(a)(e)計算では IV 価錯体の場合に遷移金属原子の電子状態の再現性が極めて悪く、これが IV 価錯体の IMiC-MO 法によるエネルギー評価が難しい原因であることがわかった。そこで、全系における遷移金属原子の電子状態を SCA 法により解析したところ、[Ti(IV)Cl6]2-中で Ti 原子はほぼ Ti2+として存在していることが示唆された。そのため、T(a)(e)
計算において、Ti 原子を Ti(II)として取り扱い全エネルギーを評価したところ、分割による誤差が-18.1 kcal/mol と改善された。つまり金属原子単体の電子状態を改善することで、IV 価錯体の全エ
ネルギーも精度良く評価できると期待される。新しい方式については検討中である。

一方で IMiC-MO 法において各原子にかかる力は、TA(e)計算中の目的原子にかかる力として計算される。
TA(e)計算においては、遷移金属原子の電子状態は、全系中の電子状態を適切に再現していることから、各目的原子にかかる力は適切に評価されていると期待される。
図 2 に、横軸に通常の B3LYP 法による全系の構造最適化により得られた原子間距離、縦軸に IMiC-MO 法により最適化した原子間距離を示した。これらの原子間距離は非常に良く相関しており、IMiC-MO 法によって遷移金属錯体の各原子にかかる力が精度良く評価されていることがわかった。
【References】 [1] S. Morita, S. Sakai, J. Comp. Chem., 22, 1107 (2001), S. Sakai, S. Morita, J. Phys.
Chem. A, 109, 8424 (2005), T. Udagawa, S. Sakai, J. Comp. Chem., 32, 2902 (2011). [2] S. Iwata, Chem.
Phys. Lett., 69, 305 (1980).

3P106
スパースモデルによる水溶液構造の解析
水上 卓 1, 杉山 歩 2, 川崎 隆史 2, 高木 啓行 2, ホー・ツー・バオ 2, ダム・ヒョウ・チ 2(北陸先端大院・マテリアル 1、北陸先端大院・知識科学 2）
Structural analysis of protein hydration water by means of sparse model
Taku Mizukami1, Ayumu Sugiyama2, Kawasaki Takafumi2,
Takagi Hiroyuki2, Ho Tu Bao2, Dam Hieu Chi2

(JAIST Materials Sci. 1, JAIST Knowledge Sci. 2 )【序】
我々はデータマイニング法を用いた物理モデル設計支援システムの構築を目指している．
データから知識を取り出すには，物理的モデルに基づいた解釈をおこなう必要があり，それは従来人間の直観にもとづいて設計されてきた．一方，近年の計測・計算・通信技術の発展により，多次元ビッグデータを取り扱う必要性が増している．しかし，大量のデータを対象にして実際に直観を働かせるにいたるまでには様々な手続きや処理が必要であり，人間の直観のみによるモデル設計には大きな困難がともなう．
我々は分子動力学計算からのトラジェクトリデータを混合分布モデルによって特徴空間に変換し，データマイニング的手法によって，タンパク質に水和する水分子の振る舞いをクラス分けしその特徴を同定してきた．今回，水分子のダイナミックスにおける隠れた属性間の関係を明かにすることを目的に，水和・バルク 2 種類の水分子の振る舞いを定量化し，それをさらにグラフモデルによって解析した．その結果水和水・バルク水の属性の関係性を示す
特徴的なグラフ表現を得たので報告する．

図 1．混合分布モデルによる水（酸素原子）トラジェクトリのクラスタ化【方法】
分子動力学シミュレーションにおいては，力場 Amber03、TIP3P モデル水分子を約 3.3×
104 個、周期的境界条件，温度 300K，NVT アンサンブルの条件で実行した。トータルのラ
ンタイムは約 10 ns．また蛋白質(PDBID:1PSV)周辺に TIP3P 水分子を約 6×103 個配置し同条件で約 15 ns のトラジェクトリデータを作成した．
混 合 分 布 モ デ ル は 関 数 と し て ガ ウ ス 分 布 を 用 い ， EM ア ル ゴ リ ズ ム (Expectation–
maximization algorithm)を応用した Gaussian Mixture Model[1]を採用した．これを水分子の酸素原子トラジェクトリに適用し，実空間上の軌跡をさまざまな大きさ・形状の回転楕円．これらのクラスタの形状を示すパラ体に対して近似する複数のクラスタに分解した（図 1）

メータ，およびクラスタに含まれる軌跡点の個数等から特徴空間を構築し，その空間上で主成分解析(PCA)および特徴空間上でのさらなる混合分布プロセスを走らせることによって，水分子のダイナミックスをクラス分けした．バルク水の系と蛋白質+水の系を比較し，水和水のクラスを同定した．
次に Gaussian Mixture Model によってクラス分けされた水分子の振る舞い属性を，複数の時間スケールにおいて再計算し，階層的時間における特徴空間を構築した．空間上において線型モデルによる regression およびスパース性を考慮したグラフマイニング[1]を行い，属性間の関係性を表現するグラフを得た．
【結果と考察】バルク水の系および水+蛋白質の系をそれぞれデータマイニングすることによって特徴空間上の点の集合が得られる．各要素は水分子のある区間の運動の軌跡を特徴の座標軸上にマッピングされた点である．今回，蛋白質をバルク水に“投入”したことにより，大きくポピュレーションが増加したクラスを蛋白質表面に水和した水分子であると予想し，それに含ま

図 2．蛋白質 vdW 表面からの酸素原子の空間分布

図 3．第一水和水に対応すると示唆される水分子
クラスの空間分布．

れる水分子を可視化した．
蛋白質のファンデルワールス(vdW)表面と水の酸素原子との間の距離分布を図 2 に示す．2
個のピークは平均すると約 1.4Åの距離にあった．このときの水分子の空間分布を示す(図 3)．
蛋白質の疎水性，親水性のアミノ酸残基をそれぞれ赤色，青色で示し，水分子の分布はピンク色で示す．親水基の近傍に偏って水分子が分布しているのが認められ，この事から親水基に強く水素結合している第一水和水に相当するクラスであること判断される．
次に水和水，およびバルク水における属性間の関係性に関してグラフマイニングを行い，それぞれ特徴のあるグラフを得ることに成功した．これは，2 つの状態にある水分子の振る舞いの時空構造が異なることを示唆している．詳しい考察は講演でおこなう予定である．
今回おこなったグラフマイニングによる属性間の関係性の描出は，水分子のクラス分けの
次のステップに位置づけられる．この方法によって，データマイニングによってモデルの構
築を逐次的・再帰的に行い，新たな知識を獲得していくことを目指す．
[参考文献][1] Data Mining 2nd ed. I.H. Witten E. Frank, Elsevier (2005)

3P107L1 型正則化線形回帰モデルを用いたシミュレーションデータからの二元合金の物性予測○鈴木 大輔 1，川崎 隆史 1，杉山 歩 1，水上 卓 2, Dam Hieu Chi1, Ho Tu Bao1
（北陸先端大・知識科学 1，北陸先端大・マテリアル 2）
Prediction of binary alloy property from the simulation data by using L1 type regularization
linear regression modelDaisuke Suzuki1, Takahumi Kawasaki1, Ayumu Sugiyama1, Taku Mizukami2, Dam Hieu Chi1, Ho Tu Bao1

（JAIST Knowledge Sci. 1, JAIST Materials Sci. 2）【序】情報科学技術の進展とともに発展してきた計算科学は幅広い分野で成果を上げ、量子化学計算
による材料設計、分子動力学シミュレーションによるタンパク質や高分子ポリマーの解析など産業面においても必要不可欠なものとなっている。
計算機シミュレーションは実際の物質を用いた実験と比べ、多くのサンプルを多彩な条件下で簡便に実行することが出来るとの期待から時間と設備の節約になるものと考えられてきたが、近年の大規模シミュレーションでは数万コアのリソースを必要とされ、金銭的負担は看過出来ない程となっている。
一方、目的の物性と関連のあるパラメータの探索を試み、その結果を用いて予測を行う手法の研究もされており、これらの研究はいわゆる逆問題的設計法として期待が高まっている。Yousef
ら[1]は CaAg, BaPb といった二元合金の融点を、構成元素の価電子数、原子番号、電気陰性度、沸点など計 16 個の基本物性のみを用いたデータマイニングにより予測し、また、予測結果に寄与する物性を示している。このような研究は、帰納的知識発見法である機械学習やデータマイニン
グ手法の発展と共に盛んになってきているが、まだ実用的な手法として確立されたものはない。
本研究では、基礎物性データと第一原理シミュレーションから得られる物性データから成る物
性データベースを構築し、本データベースから重回帰分析により物性間の相関グラフ構造を明示化する方法を確立した。本手法の概要及び実験結果を以下の通り報告する。
【方法】
本研究は以下の手順で行う。
１．第一原理計算による合金物性データベース構築
２．予測に使用する物性特徴空間設計
３．重回帰分析による物性値（融点）予測モデル構築
４．並列した予測モデル構築からの全物性間の関係性を行列化及びグラフ化
本研究で行った全ての物性計算は、密度汎関数論（DFT）に基づいた第一原理計算であり、Dmol3
によって実行した。構造最適化及びパラメータ決定は密度勾配近似（GGA）、汎関数に PBE、基底関数に DNP を用いて計算した。尚、内殻電子処理は相対性全電子とした。
本計算によって求められた二元合金の結合距離、Mulliken 電荷など、19 種の物性値と二元合

金の構成元素の基礎物性データから物性データベースを作成し、融点予測に使用する特徴空間とした。
予測モデルとして様々な既存のモデルがあるが、本研究では予測対象とパラメータ間の関係性の解釈が容易であり、グラフ化にも適用しやすいという理由から線形回帰モデルによる予測を行う。本研究では下に示す代表的 L1 型正則化線形回帰手法である Least absolute shrinkage andselection operator（Lasso）を採用した。

𝑝

𝐿𝑎𝑠𝑠𝑜

𝑇

： 𝑚𝑖𝑛 ((𝑦 − 𝑋𝛽) (𝑦 − 𝑋𝛽) + 𝛾 ∑|𝛽𝑗 |)
𝑗=1
𝑦：実際値

𝑋：パラメータ

𝑝：パラメータ数𝛽：回帰係数𝛾：チューニングパラメータ

Lasso は正則化効果により、図１の通り、チューニングパラメータγの値が決定すると一意的にパラメータが決定されモデルが求まり、また、γの値が大きいほどモデルのパラメータ数が減少する特徴を持つ。従って、Lasso を用いることによりスパースな線形モデルを得る事が可能となる。また、汎用性のある予測モデルを得るために、データをモデル構築用とモデル評価用に分割し適用する交差検定を行った。この Lasso と交差検定の融合手法を用いて、順次γの値を変えてモデルを構築し、予測誤差が最小となるモデルを最良の予測モデルと決定した。
本手法による物性－パラメータ間の関係性の予測モデル構築を、予測対象物性毎を入れ替え並列に実行することにより、二元合金の全物性間の関係性を全て行列化した。また、この関係性行列を以って二元合金物性の関係のグラフィカルモデルを作成した。
【結果】
結果は図２の通りとなった。プロットが対角に近い
程、融点の実際値と予測値の差が小さいことを意味する。
本手法を用いて二元合金の融点予測を行ったところ、
平均二乗誤差１０１K での予測を行うことが出来た。従って、本手法を用いることにより、予測対象物性－パラメータ間の関係性を高い精度で表現することが可能であると言える。
融点予測の結果の詳細及び二元合金物性間の関係のグ
ラフ化の結果については当日報告する。
【参考】[1] Y.Saad, D.Gao, S.Bobbitt, J.R.Chelikowsky, and W.Andreoni, Phy. Rev. B 85 104104 (2012)[2]R. Tibshirani, J. R. Statist. Soc. B 58, 267 (1996)[3]B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, Annals of Statistics 32, 409 (2004)[4]N. Meinshausen and P. Buhlmann, Ann. Statist. 34 (2006)

3P108
溶液中におけるウラシル誘導体の光物理的過程に関する理論的研究(北大院総合化学1，北大院理2，弘大院理工3) ○岡井昌幸1，中山哲2，山崎祥平3，武次徹也2

Photophysical processes of uracil derivatives in solution phase(Hokkaido Univ.1, Hirosaki Univ.2) ◯Masayuki Okai1, Akira Nakayama1, Shohei Yamazaki2, andTetsuya Taketsugu1【研究背景】DNA 塩基は紫外領域に吸収帯を持つが、高効率な無輻射失活経路を有しており、この性質は DNA の光安定性として知られている。近年の精力的な研究により、この失活過程はサブピ
コ秒のオーダーで起こり[1]、ポテンシャル曲面の円錐交差(conical intersection: CI)点が重要な役割を果たしていることが分かっている。気相中の単分子塩基の失活過程に関しては実験的にも理論的にも多くの知見が得られてきたが、溶液中や DNA 骨格内での失活過程に関しては、周囲の環境が励起寿命に及ぼす影響等、未だ十分に理解はされていない。2008 年に Gustavsson らにより、様々な溶媒中におけるウラシル誘導体（ウラシル、チミ
ン、5-フルオロウラシル、図 1 参照）の時
間分解蛍光スペクトルが得られ、励起寿

図 1 ウラシル、チミン、5-フルオロウラシル[2]

命が測定された 。その結果、水溶液中ではウラシルは一つだけの緩和成分を持ち、その失活時間は約 0.1 ps と報告された。一方、チミンと 5-フルオロウラシルでは二つ以上の緩和成分が存在し、失活時間はチミンでは約 0.2と 0.6 ps、5-フルオロウラシルでは約 0.6 と 1.6 ps と測定され、5-フルオロウラシルでは失活
時間がさらに長くなることが報告された。また、様々な溶媒中で測定を行っており、緩和成分の相対比と失活時間の解析が行われたが、溶媒の極性、粘性といったマクロな物性値との明確な相関は得られなかったと報告された。

そこで、本研究では、溶媒を露わに取り扱う QM/MM−MD シミュレーションを行い、溶媒分子がウラシル誘導体（ウラシル、チミン、5-フルオロウラシル）の緩和過程に及ぼす影響を理論的に解明する。特に、置換基の効果や溶媒の性質（極性、プロトン性）による影響を吟味する。

【計算手法】
ウラシル誘導体（ウラシル、チミン、5-フルオロウラシル）に対して、溶液中での QM/MM-MD シミュレーションを行った。
QM 領域として取り扱うウラシル誘導体の量子化学計算には、基底状態に対しては MP2 法を用い、励起状態に対しては CASSCF(4,4)を参照関数とした CASPT2 法を用い た （ 以 下 CASPT2(4,4) と 記 述 ） 。 基 底 関 数 は 全 て
Sapporo‒DZP を用いた。また、MM 領域である溶媒分子には、水には SPC/F、メタノール、アセトニトリルには AMBER

図 2 QM/MM-MD 計算のスナップショット

の力場をそれぞれ用いた。一辺 16 Å（水溶液中）、24 Å（メタノール、アセトニトリル中）の立方体セル内にウラシル誘導体 1 分子と水 122 分子（メタノール 200 分子、アセトニトリル 155 分子）を配置した。MD におけるタイムステップは 0.5
fs とした。

【結果】
QM/MM−MD計算からそれぞれの分子の吸収スペクトルを求めたところ、溶媒によって大きく変化することがわかった。これは主に溶媒分子との水素結合の有無による影響が大きい。
ウラシル誘導体の光励起緩和過程には、1ππ*と1nπ*の二つの状態が関与することが報告されているが、本発表では明状態である1ππ*状態のみを経由した過程に焦点を当てる。これまでに、
C5=C6のねじれを伴う円錐交差点(1ππ*(C5=C6)/S0)CIを経由した失活過程が報告されている。
基底状態で温度一定(T = 300 K)のMDシミュレーションを行い、系を平衡化させた後、一定時間間隔で座標と運動量を取り出し、1ππ*状態でのQM(QM/MM)‒MDシミュレーションを開始した。MDシミュレーションでは、計算コストの観点から活性空間の小さいCASPT2(4,4)を用いているが、トラジェクトリに沿って一定の間隔でCASPT2(12,9)のエネルギー計算を行い、ππ*状態に関してはこれらの小さな活性空間の計算で十分な精度が得られていることを確認

1した。CASPT2(12,9)の計算結果より、1nπ*状態のエネルギーを求めることができるため、時間分解蛍光スペクトルの緩和成分の解釈として、(i)明状態である1ππ*状態のみを経由して円錐交差領域に到達する過程と(ii) 1ππ*状態から暗状態の1nπ*へ遷移する過程に着目し、それぞれの時定数と割合を求め、溶媒を変えた時の相対的変化を調べた。実験結果の解釈等、詳細は当日報告する。
【参考文献】[1] C. Canuel et al., J. Chem. Phys. 122, 074316 (2005).
[2] T. Gustavsson et al., Chem. Phys. 350, 186 (2008).
[3] S. Yamazaki and T. Taketsugu, J. Phys. Chem. A 116, 491 (2012).

3P109量子波束法に基づく非断熱遷移を考慮した ab initio 分子動力学法の開発とアゾベンゼン光異性化反応への適用
1（北大院・総合化学 , 北大院・理 2）○近藤 有輔 1, 小林 孝生 1, 原渕 祐 2, 野呂 武司 2, 武次 徹也 2

A new surface hopping AIMD approach utilizing wavepacket simulation and its application toazobenzene photoisomerization(Hokkaido Univ.) ○Yusuke Kondo, Takao Kobayashi, Yu Harabuchi, Takeshi Noro, Tetsuya Taketsugu【研究背景】近年、光化学反応に対する動的反応機構やエネルギー移動など動力学の詳細を議論するために、電子励起状態のポテンシャル曲面上での原子核の時間発展を追跡する ab initio 分子
動力学(AIMD)計算が行われるようになってきた。AIMD 計算では、電子状態計算により原子核に働く力を逐次計算しながら時間発展を追っていくので、ポテンシャル関数をあらかじめ用意しな
くても全自由度が考慮された精度の高い分子動力学計算が可能であり、計算コストは大きくなる
が任意の反応素過程のダイナミクスを調べることができる。光化学反応のシミュレーションでは励起状態に対する電子状態計算が要求され、複数の電子状態が近接する非断熱領域では状態遷移を考慮するための非断熱結合項の計算が必要となる。当研究室では、Tully の最少遷移数アルゴリ
ズム[1]に基づく非断熱遷移を考慮した AIMD プログラムを開発しているが、各ステップで非断熱結合項の計算が必要となり計算コストが高くなるため、より効率的な計算手法を模索していた。
最近、南部、石田、中村により、朱-中村理論の公式に基づき非断熱遷移を考慮する AIMD 法が提案され、いくつかの応用計算が報告された[2]。非断熱遷移はポテンシャル曲面が最も近接した点で起こると想定の下、2 状態のポテンシャル曲線の形状に基づき公式にしたがって遷移確率が見積もられ、非断熱遷移の有無が判定される。朱-中村理論は 1 次元に対する理論であるので、AIMD 計算では 2 状態のポテンシャルの最近接点に到達したときに非断熱結合ベクトルを計算し、そのベクトル方向に 1 次元のポテンシャル曲線を計算して朱-中村理論が適用される。すなわち、
朱-中村理論に基づく AIMD 計算では非断熱結合ベクトルはポテンシャルの最近接点においてのみ計算すればよいので、Tully のアルゴリズムに比べて計算時間の大幅な軽減が期待される。一方、
朱-中村理論では、2 状態のポテンシャル曲線の関係が Landau-Zener 型と非断熱トンネル型に分類され、それぞれの場合について遷移確率の公式が与えられている。当研究室で朱-中村理論を励起状態 AIMD プログラムに実装し、アゾベンゼン光異性化反応へと適用したところ、S0, S1 状態の最近接点で非断熱結合方向に計算したポテンシャル曲線の形状がこの 2 つの分類に当てはまらない場合が 30%ほど見出された。そこで本研究では、断熱状態間の最近接点で切り出した 2 状態 1
次元のポテンシャルに対し量子波束法[3]を適用することで遷移確率を見積もる新たな状態遷移
AIMD プログラムを開発することとした。
【方法】量子波束法では核の量子力学的な効果は厳密に取り込まれるが、時間発展させるにはス
テップ毎に座標空間全域のポテンシャルの情報が必要となるため、原子数が増えると全自由度を扱うことはコスト的に不可能となる。本研究では、原子核の自由度は AIMD 法に基づき古典力学
で時間発展させ、非断熱領域で遷移確率を見出すときにのみ 1 次元の量子波束計算を実行することにした（図 1 参照）。
具体的には、古典軌道に沿って断熱状態間のポテンシャルエネルギー差をモニターし、エネル
ギーの近接領域でエネルギー差が極小となった原子核配置を抜き出す。つぎに、この配置における非断熱結合ベクトルを計算し、その方向に沿って各断熱状態のポテンシャル曲線を ab initio 法
により計算する。ここで得られた 1 次元 2 状態の断熱ポテンシャル曲線を透熱ポテンシャルに変換し、量子波束計算を行う。波束の時間発展は Xin Chen らの方法[3]に従い、グリッド上の振幅値で表現された波束を Split-Operator-Fourier-Transform (SOTF) 法により時間発展させる。初期波束は古典軌道から切り出してきた核座標と運動量に応じて決定し、波束が分岐した割合に従って遷移確率を算出する。その後、遷移確率と 0~1 の間で発生させた一様乱数の大小関係で状態遷移の有無を決め、原子核をいずれかの断熱状態のポテンシャル曲面に乗せて AIMD 法による古典軌道計
算を継続する。
【結果と考察】まず Xin Chen [3]らの論文に従って 2 状態 1 次元の量子波束計算プログラムを作成
し、simple avoided crossing を示すモデルポテンシャルへと適用してテスト計算を行った。初期波束として用いるガウス波束の幅は、原子核の位置や運動量がどの程度の広がりを持っているかを表す重要なパラメーターであり、特に非断熱遷移のような量子効果を取り扱う際にはきわめて重要となる。本計算では、初期波束の幅、波束の初期位置、グリッド数、グリッドの範囲を変えて分岐比を見積もり、安定的な結果が得られる適切なパラメーターについて慎重に検討した。
つぎに、量子波束計算プログラムを MOLRPO とリンクした AIMD プログラムに実装し、実際の分子系の状態遷移ダイナミクス計算が可能なようにプログラムを整備した上で、アゾベンゼンの n* 励 起 に よ る シ ス - ト ラ ン ス 光 異 性 化 反 応 へ と 適 用 し た 。 電 子 状 態 計 算 レ ベ ル は
SA2-CASSCF(6,4)/6-31G*とし、CASSCF の活性空間には 2 つの n 軌道及び* 軌道を含めた。
AIMD 計算の初期条件としては、温度 300 K の条件で原子核の座標と速度を乱数で決定した基底状態(S0)のシスアゾベンゼンを n*状態(S1)へと垂直励起させ、S1 状態のポテンシャル曲面上を運動して S0 へと非断熱遷移するときの遷移確率を量子波束計算に基づき評価した。図 2 には、S1 から S0 への非断熱遷移が 62 fs で起こった古典軌道の例を示す。非断熱遷移に関するより詳細な解析については当日報告する。

質点 量子
波束
古典
軌道 遷移
非断熱
領域

非断熱結合
ベクトル方向に沿った
ポテンシャル

質点

図 1. 量子波束法による図 . 非断熱遷移の概念図

図 2. 基底状態(S0)と n*状態(S1)のポテンシャル変化【参考文献】[1] J. C. Tully, J. Chem. Phys., 93, 1061 (1990).
[2] S. Nanbu, T. Ishida, and H. Nakamura, Chem. Science, 1, 663 (2010).
[3] X. Chen and V. S. Batista, J. Chem. Phys., 125, 124313 (2006).

3P110
AIMD-WKB 法の開発と電子励起状態におけるトロポロン分子内プロトン移動によるトンネル分裂計算への適用（北大院総化 1，北大院理 2，物材機構 3）○佐藤 彩 1，原渕 祐 2，大谷優介 3，武次徹也 2
AIMD-WKB approach to tunneling splitting due to proton transfer in the excited state of tropolone(Hokkaido Univ.1, NIMS2) ○Aya Satoh1, Yu Harabuchi1, Yusuke Ootani2, Tetsuya Taketsugu1

エネルギー【序論】プロトン移動では核の量子性に由来する
トンネル効果が重要な役割を果たす。実験的には、
トンネル効果は速度論的同位体効果や対称二重トンネル
井戸ポテンシャルにおける零点準位の分裂として観測されるが、理論的にトンネル効果を扱う場
合には、プロトン移動以外の自由度をどのように考慮するかが一つのポイントとなる。また、図 1
Re
反応座標に示すようにプロトン移動を伴う互変異性では
図 1．励起状態におけるプロトン移動
電子励起状態において一般に活性化障壁が下が反応の模式図るため、トンネル効果の寄与が大きくなる。励起状態におけるプロトン移動過程について、分子系の全自由度を考慮したトンネルダイナミクスを調べることのできる第一原理動力学計
算手法の開発が求められている。
当研究室では最近、分子動力学計算に基づきトンネル分裂を半古典的に見積もる
Makri-Miller の方法を ab initio 分子動力学(AIMD)プログラムに実装し、アンモニア及びマロン
アルデヒドの電子基底状態におけるトンネル分裂計算に適用して実験値と良い一致を得ることに成功した[1]。本発表ではこの手法を電子励起状態に拡張し、トロポロン(C7H6O2)の励起状態におけるプロトン移動（図 2）に伴うトンネル分裂計算へと適用する。分光実験により、トロポロンの第一電子励起状態における振動基底状態のトンネル分裂値は 18.9 cm-1 と報告されており、同位体効果や他の振動モードの励起の効果など詳細な実験報告がある[2]。一方理論研究では、トンネルに直接関係する 2~3 自由度に制限した座標空間の中での量子力学計算
によりトンネル分裂値を議論した報告がある[3]。
また、トロポロンの電子基底状態に対し、
Makri-Miller 法に基づき経験的ポテンシャル関数を用いた全自由度トンネル動力学計算も報告されている[4]。本発表では、全自由度を考慮した
AIMD 計算によりトロポロンの励起状態におけるトン
図 2．トロポロンの分子内プロトン移動
ネル分裂値を計算し、トンネル経路への依存性や重要な自由度に関する知見を得ることを目的とする。
【IRC-WKB 計算】トロポロンの電子基底状態および第一励起状態(状態)に対し、
MS-CASPT2/aug-cc-pVDZ 法により平衡構造と分子内プロトン移動の遷移状態構造を決定した。続いて電子励起状態における AIMD 計算を実行するためのより低コストな計算手法を検討し、MS-CASPT2 によるエネルギー障壁と反応経路距離をよく再現する LC-BOP に基づく
TDDFT 法を用いることとした。TDDFT 法により固有反応座標(IRC)を計算し、IRC に沿って振動解析を行って各振動モードのゼロ点振動エネルギー(ZPE)を加算し、得られたポテンシャ
ル曲線に基づき WKB 近似の下でトンネル分裂値を求めたところ、トンネル分裂値が実験値

の約 1/10 と過小評価されることを確認した。この過小評価は、IRC の遷移状態前後の曲りに起因するものである。
【AIMD-WKB 計算】次に AIMD-WKB 法によりz
トンネル分裂値の計算を行う。初期条件についてはマロンアルデヒドの AIMD-WKB 計算[1]で行ったやり方を踏襲し、モード間のエネルギー移動に
yより古典的に反応が起こるのを避けるために自由度の二分を試みた。すなわち、分子内プロトン移
始点
終点x動に関係する振動自由度には量子力学的な ZPEを与え、それ以外の自由度(以下 bath mode とする)
図 3．始点と終点の関係には古典的エネルギーkT(k: ボルツマン定数、T:温度)を与え、温度パラメーターT を変えることで他のモードがトンネル分裂に及ぼす効果を調べた。古典軌道に沿った各点で x-z 平面に関する鏡像の構造を生成し(図 3)、これらの構造をつなぐ 1 次元のトンネル経路を決め、速度のトンネル方向の成分をモニターすることにより古典的転回点を見出し、古典的転回点を始点としてトンネル経路に沿って作用積分を計算し、得られたトンネル振幅を積算したものを多数の古典軌道で平均化して得られたトンネル
振幅曲線の傾きを読み取ることでトンネル分裂値を計算する。
トンネル経路の定義であるが、始点と終点をデカルト座標で直線的に線形補間した経路では障壁が高くなりトンネル分裂値を大幅に過小評価されたため、IRC に沿った内部座標変化を参考に、OH と O が近づいてから H が隣の O に移動し、O と OH が離れていく特徴を有す
る経路となる定義を工夫した（始点→経由点→中点→経由点→終点となるように経由点と中点を導入した）。
現在古典軌道の数を増やしているところであるが、現時点でトンネル分裂値は、0 K で 27.5
cm -1 (古典軌道 25 本)、10 K で 27.1 cm -1 (古典軌道 30 本)、30 K で 27.2 cm -1 (古典軌道 27 本)、
50 K で 33.5 cm -1 (古典軌道 46 本)となり、bath mode のエネルギーが高い程トンネル分裂値が増加する傾向にあることがわかった。計算結果とトンネル経路に関する議論の詳細は当日報告する。

【参考文献】[1] Y. Ootani and T. Taketsugu, J. Comp. Chem., 33, 60 (2012).
[2] H. Sekiya, Y. Nagashima, T. Tsuji, Y. Nishimura, A. Mori, and H. Takeshita, J. Phys. Chem., 95,
10311 (1991).
[3] M. J. Wojcik, L. Boda, and M. Boczar, J. Chem. Phys., 112, 6322 (2000).
[4] Y. Guo, T. D. Sewell, and D. L. Thompson, J. Phys. Chem. A, 102, 5040 (1998).

3P111
水―アセトニトリル溶液における CN 結合伸縮振動の周波数の環境依存性（京都大院・理 1，UHH, MPSD, CFEL2）○植野 正嗣 1，長谷川 太祐 2，Huse Nils2

Dependency of the CN-stretching frequency on the environment inwater-acetonitrile mixture(Kyoto Univ.1, CFEL2) ○Seiji Ueno1, Taisuke Hasegawa2, Nils Huse2【序】アセトニトリル（AcCN）はニトリル基の窒素原子上に非共有電子対を持ち、水分子の水素原子と水素結合を生成する。AcCN 分子の CN 結合伸縮振動の IR/Raman 分光における周波数は水溶液中
では青方遷移することが確かめられている[1][2]。一方で X 線ラマン分光の実験から、分子間の構造において水素結合を生成しない AcCN と水分子間の双極子相互作用も重要であると示唆された[3]。
我々は理論計算を用いて周波数のシフトにおける水素結合と電場の依存性を調べる。
古典 MD を用いてクラスタの構造を得て各クラスタにおいて電子状態計算を行い、それぞれの構造
での電場・水素結合の結合係数・伸縮振動の振動数のサンプルを得る。そのマッピングから AcCN 分子の環境に周波数がどの程度依存しているかを評価する。また古典 MD からマッピングを用いて IR
スペクトルを計算する。

【計算】水と AcCN の等モル溶液系を古典ポテンシャルのＭＤで再現した。ポテンシャルパラメータとして、水は Flexible SPC/E、AcCN には 6 点モデル[4]を用いた。300 K に設定した NVT アンサンブ
ルを用いて、そのトラジェクトリからランダムにクラスタ構造を抽出した。対象となるクラスタは水分子と AcCN 分子上の窒素原子が近接している分子ペアを中心におよそ 2000 原子程度とした。そしてその中心の AcCN 分子を含めた周囲の数分子を量子化学計算（B3LYP/6-31+G(2d,p)）で解析した。
大きな電荷を持つ水分子による影響や分子間力を適切に評価するため、基底関数は分極・分散関数を追加したものを使用した。それ以外の分子は、量子化学計算において有効な電場を形成する点電荷としてみなすことでその効果を考慮した。
クラスタ解析において、中心となる AcCN 分子について構造図計算されるクラスタの概形

1 Beth. A. Lindquist, Kristina E. Furse and Steven A. Corcelli (2009). Phys. Chem. Chem. Phys., 11, 8119-8132.
2 Francesco Muniz-Miranda, Marco Pagliai, Gianni Cardini and Roberto Righini, (2012) JCP, 137, 24450
3 Huang, N., Nordlund, D., Huang, C., Bergmann, U., Weiss, T. M., Pettersson, L. G. M., & Nilsson, A. (2011). The Journalof chemical physics, 135(16), 164509.
4 Nikitin, A. M., & Lyubartsev, A. P. (2007). Journal of Computational Chemistry, 28(12), 2020–2026.

最適化後、CN 結合伸縮振動の周波数付近の 2300 cm-1 前後の周波数を持つ基準振動の座標上でSchrödinger 方程式を解く。これらのデータから基底状態から第一励起状態への周波数  、遷移モー
メント p01 、基準振動座標による双極子微分 μ  と、中心 AcCN 分子以外の原子を点電荷としたときのニトリル基の原子上での電場 EN , EC を計算した。これらの計算を多数のクラスタ構造上で行い、電場の 2 次式としてフィッティングを行う。

  a1EN2  a2 EN  a3 EC2  a4 EC  a5μ  b1EN2  b2 EN  b3 EC2  b4 EC  b5
p01  a  b
こうして得られた   EN , EC  ， p01   ， μ  EN , EC  から、B.M.Auer 氏らの方法[5]に基づいて、
CN 伸縮振動のピーク付近の領域の IR スペクトルを計算する。
一方で、Natural Bond Orbital 解析により AcCN 上の窒素原子（N）と最も近い水分子上の水素原子（H）の間に形成される結合次数 n を求め、周波数に対する影響を調べた。

【結果と考察】

図１は結合次数 n と周波数  の散布図である。また図２は電場 EN と周波数  の散布図である。これらのデータについて、2 次の多項式で回帰させた近似式と周波数  との相関係数は、結合次数   n   an2

 bn  c は 0.5058、電場   EN , EC  と  は 0.8605 であった。それらから  との偏相関係数を調べると結合次数   n  は 0.2517、電場   EN , EC  は 0.8206 であることから、周波数においては電場のほうが影響が大きいことがわかった。
また AcCN 上の N と水分子上の H との距離 r とそれらの原子間で形成される結合次数 n の関係を図３に示す。ニトリル基上の炭素原子（C）や水分子上の酸素原子（O）が関わる角度∠CNH や∠NHOに比して、 n は r に強く依存していることが分かった。
この方法により得られた IR スペクトルなどは当日発表する予定である。
水素結合が形成される状況では電場は大きくなる。一方で水素結合は形成していないが双極子相互作用がある状況や非プロトン性極性溶媒中など、電場が大きいが水素結合がない環境も十分考えられる。そのような場でも、分子振動や光学的性質において電場を十分に考慮する必要があるかどうか、これからのデータから読み取っていく。

図1

窒素原子と水素原子の結合定数と周波数の関係

図2

窒素原子上の電場と

図３

周波数の関係

n の関係

5 Auer, B. M., & Skinner, J. L. (2008). The Journal of chemical physics, 128, 224511.

距離 r と結合次数

3P112
化学反応の量子制御における非線形電場効果（原子力機構 1, プリンストン大 2） ○黒崎 譲 1，Tak-San Ho2，Herschel Rabitz2

Nonlinear field effects on the quantum control of chemical reactions（JAEA1, Princeton Univ.2） ○Yuzuru Kurosaki1，Tak-San Ho2，Herschel Rabitz2【序】電場の存在下では、一般に分子の双極子モーメントが誘起される。この効果は特に電場が強い場合に顕著となり、理論計算の際にハミルトニアンに非線形項を加える必要が生じてくる。
本研究では、オゾン分子(O3)の open-cyclic 異性化の量子制御を例として、非線形電場効果の影響を理論的に考察する。よく知られているように、オゾン分子の open 体は global minimum、cyclic
体は local minimum であり、両者を結ぶ反応経路上の遷移状態(TS)において基底状態と第一励起状
態のポテンシャル曲線が avoided cross (AC)する。ここでは、基底および第 1 励起状態の 2 本のポ
テンシャル曲線を考慮して、open 体を cyclic 体へ効率的に異性化させるレーザー場を最適制御理論(optimal control theory, OCT)により見出し、非線形電場効果が制御機構に与える影響について考察する。
【計算方法】異性化のポテンシャル曲線は MRSDCI/aug-cc-pVTZ レベルで計算し、OCT 計算に用いるため基底及び第 1 励起状態からなる diabatic モデルポテンシャルを構築した[1]。図 1 に異性化の IRC(Intrinsic Reaction Coordinate)に沿った diabatic ポテンシャルの対角項(V11d, V22d)及び非対角項(V12d)を示す。OCT 計算においては、以下に与えられる汎関数 J：
J   (T ) 

2

 0







 dt (t )  2 Re  (T )   dt  (t ) t i( H   (t ))  (t ) 
T

2

0

T

0

0

(1)

を考え、これを最大にするレーザー場(t)を求める。(t)は直線偏光しているとし、その向きはオゾン分子の頂角の 2 等分線と平行に固定する。式(1)においては target 状態、すなわち cyclic 体の振動基底状態を表す。第 2 項はレーザーのフルエンスを表し、0 は正の整数である。第 3 項の(t)はLagrange の未定乗数である。電場存在下での分子の誘起双極子モーメントは、summation convention-224.75
Energy / au

0.0010-224.80

Energy / au

-224.850.0008dV120.00050.00030.00005.65 5.70 5.75 5.80 5.85 5.90 5.95

-224.90s / amu1/2bohr

Vd22-224.95-225.00-225.05

Vd11-225.10
-2

0

2

4

6
8
10
s / amu1/2 bohr

12

14

16

図 1．IRC に沿った diabatic ポテンシャル曲線

図 2． 2 つの反応経路

を用いて一般に次のように書ける：

 i   i0  (2!) 1  ij j  (3!) 1  ijk j  k  ,

i , j , k  x, y , z(2)

ここで、i0 は永久双極子、ij は分極率テンソル、ij は超分極率テンソルである。これらの量も
MRSDCI/aug-cc-pVTZ レベルで求め、汎関数 J（(1)式）中のとして分極率および超分極率テンソ
ルまで考慮することで非線形効果を取り込んだ。wave packet の時間発展計算は、Broeckhove ら[2]による diabatic ポテンシャルを直接用いる方法に基づき split-operator 法を用いて実行した。
【結果】図 2 に OCT 計算で検討した 2 つの反応経路(Pathways 1 and 2)を示す。Pathway 1 は、initial
状態にある wave packet がまず V22d 上の浅い well に励起され、その後 AC 点を経由して最終的にtarget 状態に至る経路である。一方 Pathway 2 は、wave packet がまず V11d 上で振動励起され、次にV22d 上の深い well の高振動励起状態に励起された後、AC 点を経由して target 状態に至る経路である。表 1 に OCT 計算の主な結果を示す。Pathway 1 に対する遷移確率(yield)は Pathway 2 に対する
それより大きく、また Pathway 1 に対する電場強度は Pathway 2 に対するそれより格段に小さい。
非線形電場効果は Pathway 1 に対してはほとんど無視できる程度であるが、Pathway 2 に対しては遷移確率を引き上げ、電場強度を下げるという結果を得た。Pathway 2 に対する最適電場のスペク
トルを図 3 および 4 に示す。非線形効果が無い場合(図 3)、V11d から V22d への電子励起に対応する
スペクトルが > 0.8 a.u.の領域でみられるが、同効果を考慮した場合(図 4)、これらのスペクトルは消失し振動数ゼロの成分の寄与が非常に大きくなっている。制御機構の詳細については発表当日議論する。
表 1．OCT 計算の主な結果

0apathway

Max. field amp./a.u.

Fluence/a.u.

Yield

No nonlinear field effectsb
10.0005
0.222
140.4
0.973
2
0.00010.774
1318.60.637
cNonlinear field effects
10.0005
0.190
134.3
0.975
2
0.00001
0.855
607.1
0.858
a
Weight in the cost function of Eq. (1). bVint = - is adopted. cVint = -( + (2!)-1-13) is adopted.
1.8
1.6
1.4

4

|ε()| / arb. unit|ε()| / arb. unit

5

3
2

1.2

1.0

1.0
0.8
0.6
0.4

1

1.5

0.5
0.0
0.000

0.005

0.010

0.2

0
0.00

0.02

0.04

0.06
/ a.u.

0.08

図 3．Pathway 2 に対する最適電場のスペクトル（非線形電場効果無し）

0.10

0.0
0.00

0.02

0.04

0.06

0.08

0.10

 / a.u.

図 4．Pathway 2 に対する最適電場のスペクトル（非線形電場効果有り）[1] Y. Kurosaki, T-S. Ho, and H. Rabitz, 第 6 回分子科学討論会(2012)東京.
[2] J. Broeckhove, B. Feyen, L. Lathouwers, F. Arickx, and P. Van Leuven, Chem. Phys. Lett. 174, 504 (1990).

3P113
イミダゾール分子間のプロトン移動における透熱ポテンシャル（金沢大院・自然） ○堀 優太，井田 朋智，水野 元博

Diabatic Potential for Proton Transfer between Imidazole Molecules(Graduate School of Natural Science and Technology, Kanazawa University)○Yuta Hori, Tomonori Ida, Motohiro Mizuno[ 序 ]プロトン移動反応のような化学反応の理論的取り扱いのひとつとして，ポテンシャルエネルギー曲面の作成が挙げられる．実際の計算においては Born-Oppenheimer 近似に基づき，ある核配置に対して，量子化学計算による全電子エネルギーを求めることによって断熱ポテンシャルエネル
ギー曲面が得られる．断熱ポテンシャルが得られれば，種々の方法により散乱断面積や反応速度定数，また反応経路などの解析が可能となる．一方，断熱ポテンシャルとは別に透熱ポテンシャ
ルから化学反応を理解しようとする研究もなされている．断熱系ではなく，透熱ポテンシャルは
量子化学計算からは一意的に決定できず，一般的にある種の近似が必要ではあるが，断熱表現では考慮しにくい電子励起状態との非断熱相互作用をあらわにポテンシャルに含めることができ，運動方程式の観点からは取り扱いが容易となる．しかし，透熱ポテンシャルの研究はモデル系に限られており，実在系に対してほとんど適用例が報告されていない．
プロトン移動反応は，化学反応の中でも基礎的でかつ重要な反応であり，有機化学や生物化学などの様々な分野において注目されており，多くの場合 1 次元の反応座標に沿った断熱ポテンシ
ャル上の力学として取り扱われてきた．しかし，化学結合の切断と生成の観点から反応を理解していく上でも透熱ポテンシャル上での議論が重要であると考えられる．そこで本研究では，イミ
ダゾール分子二量体のプロトン移動反応に注目し，量子化学計算によって得られる分子間のプロ
トン移動および振動運動を考慮した２次元の断熱ポテンシャルから，透熱ポテンシャルを作成し、プロトン移動反応に関する知見を得ることを目的とする．

[ 計算・理論 ]Fig.1 のようなイミダゾール分子二量体に対して水素結合軸の伸縮運動とプロトン移動に対する２次元の断熱ポテンシャルを計算した．計算には DFT を用い，汎関数を B3LYP，基底関数として
Aug-CC-pVDZ を選択した．

Fig.1 イ ミ ダ ゾ ー ル 二 量𝑑𝑖

透熱波動関数𝜒 の断熱波動関数𝜒
列 U を用いると，𝜒

𝑎𝑑𝑎𝑑への変換は，断熱-透熱変換行体の相互作用のモデル図

𝑑𝑖= 𝑈𝜒 となる．断熱ポテンシャルエネルギー

行列𝑉 𝑎𝑑 は透熱ポテンシャルエネルギー行列𝑉 𝑑𝑖 を対角化することによって得られる．

𝑉 𝑎𝑑 = 𝑈𝑉 𝑑𝑖 𝑈 †今回は一次元の透熱ポテンシャルとしてプロトン移動前後の結合に対する２つの valence bond
配置を仮定し，
透熱ポテンシャル行列の行列要素として，対角要素は Morse ポテンシャルを用い，非対角要素である非断熱相互作用としてガウス関数を用いた．それぞれ用いた関数を以下に示す：

(

)

(

)

𝑉 𝑑𝑖 ( ) = (𝑉 𝑑𝑖 ( ) = (

(
(

)
)

)
)

𝑉 𝑑𝑖 ( ) =量子化学計算により得られた断熱ポテンシャル曲面と，透熱ポテンシャルを比較することにより各透熱系での最適なパラメータを求めた．

[ 結果 ]
ポテンシャルエネルギー (kJ/mol)

量子化学計算により得られた断熱ポテン
シャルと，透熱系を仮定し最適化された一次
元ポテンシャルを Fig.2 に示す．また，最適化により得られたパラメータを Table 1 に示
す．イミダゾール分子間のプロトン移動における非断熱相互作用の強度（A）は，プロトンの解離エネルギー（D）に対し 5 %程度であった．
また，相互作用の半値全幅は 0.298 Å
であった．

Fig.2

60

40

最適化値

量子化学計算

20

0
-0.5

-0.3

-0.1

0.1

0.3

0.5

N-N 間の中心からの H の距離 x(Å)

一次元のポテンシャルに対するフィッティング

次に，同様に最適化された２次元系へ拡張
Table1 得られたフィッティングパラメータ

した結果を Fig. 3，4 に示す．図より透熱ポテ
ンシャルを用いても，断熱系を十分に記述可

A (kJ/mol)

D (kJ/mol)k (1010m-1)

b (1020m-2)

能であることが確認された．二次元系で得ら

31.18

667.7

1.321

31.26

れた非断熱相互作用の構造や，プロトンのダN-N 間の距離 R(Å)

N-N 間の距離 R(Å)

イナミクスに関する詳細は当日発表する．

N-N 間の中心からの H の距離 x(Å)

N-N 間の中心からの H の距離 x(Å)

Fig.3 量子化学計算による断熱ポテンシャル

Fig.4 透熱ポテンシャルを用いた断熱ポテンシャル[ 参考文献 ][1] Y. T. Chang, et al, J. Phys. Chem. 94, 5884, (1990).
[2] Y. Kim, et al, J. Chem. Phys., 112, 6, (2000).

3P114
最適化振動座標に基づく振動擬縮退摂動論（理研・杉田理論分子科学）◯八木清，大滝大樹

Vibrational Quasi-degenerate Perturbation Theory with
Optimized Vibrational Coordinates(RIKEN, Theoretical Molecular Science Lab.) ◯ Kiyoshi Yagi, Hiroki Otaki[序] 近年，我々は振動座標を変分的に最適化する一般的な方法を提案した[1]．最適
化座標では，基準座標よりも非調和ポテンシャルの座標間結合が小さいため，これを用いた振動状態計算は効率的である．これまでに，最適化座標に基づく振動 CI 計算
は CI 展開が速く収束することを示した[1]．一方，大きい分子に対する応用を目指すには，振動 CI 法よりも計算コストが低い振動擬縮退摂動（VQDPT）法が有望である．
本研究では，最適化座標に基づく VQDPT 法を開発し，従来の基準座標に基づく
VQDPT 法と比較しパフォーマンスが向上することを示す．
[方法] VQDPT 法は振動 SCF（VSCF）解を出発点とし，摂動的に振動モード間の相関効果を取り込む方法である[2]．VQDPT 法では，エネルギー的に擬縮退している配置により構成される P 空間と，その補空間である Q 空間に配置空間を分割する．P/Q 空間を構成する VSCF 配置をそれぞれ p/q とすると， 2 次の有効ハミルトニアンは以下のように与えられる,( H eff(2) )#& p Ĥ q q Ĥ p!
1
1(,= p Ĥ p! + ∑ %% (0)
+
(0)
(0) (
pp!
Ep(0)
2
! − Eq 'q∈Q $ Ep − Eq

(1)

(!)

ここで， Ĥ は振動ハミルトニアン，  𝐸𝐩 は配置 p に対するゼロ次エネルギーである．
有効ハミルトニアン行列を対角化することで，VQDPT2 エネルギーが得られる．
VQDPT2 計算を有効に実行するためには，P 空間を適切に設計することが肝要である．我々のアルゴリズムでは，VSCF 配置間の量子数差，λmn = ∑ ms − ns ,(2)sが配置間の相互作用の強さを表すことに着目し，求めようとするターゲット配置 t に対し，
A(t, k) = {t! |1 ≤ λtt! ≤ k } ,(3)と定義される空間から擬縮退配置を選択する．
VQDPT2 法は基音の励起エネルギー計算に最も頻繁に用いられるが，従来の方法では基音ごとに状態選択的に励起エネルギーが計算される．状態選択的な方法は，複数

の基音の VSCF 配置が強く混ざる場合に解が曖昧になるが，そのような状況は 2 次の結合項がゼロである基準座標ではほぼ起こらないため，これまで問題にならなかった．
しかし，最適化座標では，2 次の結合項がゼロでないため，基音の VSCF 配置が強く混ざるケースが頻繁に起こる．例えば，水分子の最適化座標は 2 つのローカル OH 伸
縮座標になるが，基音はそれらが強く混ざることで対称・逆対称伸縮振動状態を構成する．
そこで，本研究では以下のように P 空間構築のアルゴリズムを改良した．まず，各ターゲット配置 ti から従来のように P 空間（Pi とする）を構築する．次に，Pi が ti 以
外のターゲット配置を含むかどうかを調べ，以下の条件(4)を満たすならば P 空間を融合する，

Pij = Pi ∪ Pj .

(5)

この操作を繰り返し，元となる Pi に対する全ての

50

(7)

[結果] 開発した方法をエチレンと 1,3 ブタジエンに対し適用した．ポテンシャルには 4 次テイラ

ps_oc(C)/2MR
oc/2MR
ps_oc(C)/3MR
oc/3MR
ps_oc(C)/4MR
oc/4MR

20

1

2

3

4

80

最適化座標（oc）に基づく VQDPT2 計算は，特に非調和性が大きい CH 伸縮振動に対し，基準座標（nc）に基づく VQDPT2 計算よりも精度が改善されている．また，nc では 2, 3, 4MR 計算による収束が悪いが，oc ではスムーズな収束が得ら
れた．

6

7

8

nc/2MR
nc/3MR
nc/4MR(b) C4H6

ー展開ポテンシャルを用いた．結果は右図に示す．

5910 11 12

Mode

100

Abs Δν [cm-1]

t i , t j ,t k ∈ Pijk .

300置が含まれなくなるまで P 空間を拡張する,(6)

nc/2MR
nc/3MR
nc/4MR

10

ターゲット配置を含み，それ以外のターゲット配

Pijk = Pi ∪ Pj ∪ Pk ,(a) C2H440

Abs Δν [cm-1]

t j ∈ Pi ∨ t i ∈ Pj ,oc/2MR
ps_oc(C)/2MR
ps_oc(C)/3MR
oc/3MR
oc/4MR
ps_oc(C)/4MR

60
40
20
0

2

4

6

8

10 12 14 16 18 20 22 24

Mode

図．基準座標(nc)および最適化座標(oc) に 基 づ く VQDPT2 計 算 に よ る
C2H4 と C4H 6 の基本振動数の VCI 計
算からの絶対誤差．ポテンシャルは，
4 次テイラー展開ポテンシャルの n 体項までを考慮した（nMR）．

	  さらに，高精度なポテンシャルを用いることで，エチレンと 1,3 ブタジエンに対し，基本振動数を実験値と比較しそれぞれ絶対平均誤差 7.8, 7.2 cm-1 で求めることが出来た．
[参考文献]
1. K. Yagi, M. Keçeli, and S. Hirata, J. Chem. Phys. 137, 204118 (2012); 分子科学討論会
2012, 講演 2E02.
2. K. Yagi, S. Hirata, and K. Hirao, Phys. Chem. Chem. Phys. 10, 1781 (2008).

3P115
量子モンテカルロ法による水素結合性分子結晶の誘電物性における同位体効果の解析（理研 1，京大院・理 2）○大滝大樹 1，安藤耕司 2
Quantum Monte Carlo Study of Isotope Effect on Dielectric Phase Transition in
Hydrogen-Bonded Molecular Crystal(RIKEN1, Graduate School of Science, Kyoto University2)○Hiroki Otaki1, Koji Ando2【序】
近年，環境調和を目的として有機材料・有機デバイスの開発が活発に行われている。誘電体もその一つであり，強誘電性の発現の鍵として注目されているのがπ電子を有する水素結合系である[1]。水素結合性物質は，水素－重水素置換により秩序相への転移温度が倍近く上昇するな
ど大きな同位体効果が見られるのが特徴である。
5-ブロモ-9-ヒドロキシフェナレノン(BHP；図1)は分子内にO－H…O型の水素結合を有し，分子内水素移動により双極子モーメントのO…O方向の成分が反転する。この物質は水素体では相転移を起こさないが，重水素置換により約37Kで双極子モーメントが秩序化し反強誘電相に転移する[2]。我々はこれまでに，フラグメント分子軌道(FMO)法による分子間相互作用の解析から，π-π相互作

図 1：BHP の分子構造

用とC－H…O型の分子間水素結合が分子の双極子モーメントを誘起すること，その誘起効果
が水素の相対的な位置に対して大きく変動することを明らかにした。さらに，FMO法の結果を用いることで双極子の誘起効果を効率的に取り込んだモンテカルロ法を開発し，それらが重水素体の転移温度に決定的な寄与をすることを明らかにしている[3, 4]。しかし，これまで
は水素の量子効果を取り入れておらず，水素結合性分子結晶の特徴を捉えるには不十分であった。本研究では，水素のトンネル効果を導入して上記のモンテカルロ法を量子モンテカル
ロ法へ拡張した。BHPに適用して計算を行い，実験値と比較することで我々の手法が同位体効果を適切に記述できることを示す。
【方法】
系を記述するハミルトニアンは次の式で与えられる。

ここで第二項は結晶中の各分子間の双極子－双極子相互作用であり，pj はサイト j の分子の双極子モーメント，rjk はサイト j, k 間の位置ベクトル，rjk はサイト j, k 間の距離，ε0 は真空の誘電率である。第一項は横磁場 Ising モデルを参考に，各サイトの水素の量子効果を表すために取り入れた項である。σj1 はパウリ行列の x 成分，Ωj はトンネル周波数を表す。この量子ハミルトニアンを Suzuki-Trotter 分解を用いて古典系のハミルトニアンに変換して，モンテカルロ法を実行する。双極子モーメントの値および周囲の分子による誘起効果（双極子モーメントの増減率）は量子化学計算で求めたものを使用し，Ωj は実験値を用いた。
【結果】
水素体と重水素体の混合比を変えた系を想定し，実効的な Ω の値を幾つか用いてモンテカルロシミュレー
ションを行った。図 2 に本手法による結果（図中 QMC）と平均場理論
による結果（図中 MFT）を実験で得られている相図で比較したものを示す。本手法でも平均場理論でも，重水素体（xH=0）から水素体の割合を増やしていくと転移温度が低下し，水素体（xH=1）では相転移が起きな

図 2：量子モンテカルロ計算(QMC)，平均場理論

いという実験結果と同じ傾向が見ら(MFT)，実験により得られた相図（黒実線）。横軸はれた。しかし曲線の曲率に注目する

転移温度，縦軸は結晶中の水素体の割合を示す。

と，平均場理論の結果は実験結果から大きく離れるのに対し，本手法で得られた結果は実験で得られた相図の曲線と近いことが分かる。これは，（重）水素の配置に起因する双極子の秩序化の効果と，水素の量子性（トンネル効果）により秩序化を避ける効
果の競合を，本手法が適切に記述できていることを示している[5]。
【参考文献】[1] S. Horiuchi, Y. Tokura, Nat. Mater., 7, 357 (2008).
[2] T. Mochida et al., J. Chem. Phys., 101, 7971 (1994).
[3] H. Otaki and K. Ando, Phys. Chem. Chem. Phys., 13, 10719 (2011).
[4] H. Otaki and K. Ando, Int. J. Quantum Chem., 113, 386 (2013).
[5] H. Otaki and K. Ando, submitted.

3P116
ランタノイドランタノイド-シクロオクタテトラエン負イオン錯体の脱離エネルギーにおけるランタノイド依存性（慶大院理工）○中條恵理華、増田友秀、藪下聡

Ln dependence in the detachment energies of Ln(COT)2－ complexes(Keio Univ.) ○Erika Nakajo, Tomohide Masuda, Satoshi Yabushita【序】ランタノイド原子とシクロオクタテトラエン【序】（COT=1,3,5,7-cyclooctatetraene）が交互に配列す
る Ln-COT 錯体は、磁気的・光学的性質を示す新規
機能性材料として応用が期待されている。Ln を含
むこの錯体の電子状態は複雑であり、最も単純なサ
ンドイッチ型錯体 Ln(COT)2 についても、光電子ス
ペクトルの形状やその Ln 依存性の解析には多くの課題が残っている。負イオン化された Ln(COT)2- は
Ln=Eu,Yb を除いて共通の形式電荷 Ln3+(COT2-)2
を 取 る が 、 光 電 子 分 光 実 験 で は 2.4~2.7 eV と3.7~3.8 eV の領域に 2 つのピークが現れ、前者の脱離エネルギーの方が強い Ln 依存性を示す（図１）[1]。HF 計算と Koopmans の定理によると、これら

図 1. 脱離エネルギーの Ln 依存性。

は HOMO と n-HOMO の軌道エネルギーに対応し、光電子分光実験[1]（実線）と HF 計算（点線）は同様の Ln 依存性を示す

その Ln 依存性の違いもよく再現される。本研究では軌道エネルギーの Ln 依存性がどのような効果に起因するかについて解析を行った。
【計算】【計算】Ln には 4f 軌道までを内殻とする Dolg らの有効内殻ポテンシャルと(7s6p5d)/[5s4p3d] 基
底関数[2]を、COT の基底関数には 6-31+G(d)を用いて HF 法による計算を行った。本錯体はイオン結合性が強く、軌道エネルギーの解析には軌道相互作用だけでなく静電的な効果も考慮する必要があるため、Ln または COT を点電荷に置き換えるモデル計算を行い、静電効果を評価した[3,4]。
【結果】HOMOと n-HOMO に相当する HF 軌道【結果】の概形を図 2 に示した。e2u 軌道は COT の π 軌道からなる反結合性軌道であり、e2g 軌道は COT のπ 軌道と Ln の 5d 軌道が同位相で重なり合う結合性軌道である。Ln の原子番号が増加すると、価電子軌道である 5s・5p 軌道が縮小することでイオン

HOMO: e2u 軌道

n-HOMO: e2g 軌道

半径は減少するため（Ln 収縮）、錯体の最適化構

図 2. Molekel による HF 軌道の概形。

造における Ln-COT 間距離は La 錯体の 2.25Åから Lu 錯体の 1.99Åにかけて単調に減少する。
したがって本錯体の Ln 依存性には電子構造依存性（Ln の種類）と幾何構造依存性（Ln-COT 結合長）の 2 つが同時に含まれる。そこで、電子構造依存性を評価するため Ln//Gd 計算：構造最適化した Gd(COT)2-において Gd を各 Ln に置き換えた計算を行い、また幾何構造依存性を評価するため Gd//Ln 計算：構造最適化した各 Ln(COT)2-について Ln を Gd に置き換えた計算を行った。
この表記では、Ln//Ln 計算が構造最適化された Ln(COT)2-についての通常の計算を示す。
1) e2u 軌道（HOMO
軌道（HOMO）のHOMO）の大きな）の大きな不安定化
大きな不安定化

表 1. e2u 軌道のエネルギーの変化の傾き(eV)。

Ln の原子番号に対する軌道エネルギーの変化の傾き（Lu と La の差を 14 で割った値）について、
Ln//Gd 計算、Gd//Ln 計算の結果を表 1 に示した。
Gd//Ln 計算は Gd 錯体の計算であるが、Ln-COT

〈 〉

〈 〉

〈 〉

Ln//Gd

0.00360.00460.0082

Gd//Ln0.0672-0.0605

0.0068間距離として各 Ln 錯体の最適値を用いて評価し、その原子番号に対する傾きを示している。軌
道エネルギー〈 〉を運動エネルギー〈 〉とポテンシャルエネルギー〈 〉に分解すると、Ln 軌道を含まない e2u 軌道では幾何構造依存性が非常に大きいことが分かる。Ln-COT 間距離の短縮にともな
い、COT 軌道が Ln3+に接近すると静電的な効果により〈 〉は低下するものの、節面の存在による〈 〉の増加がそれを上回る。このように、e2u 軌道では静電効果を節面の効果が上回ることで、軌道エネルギーは上昇することが分かった。
2) e2g 軌道（n軌道（n-HOMO）の
HOMO）のわずかな）のわずかな安定化
わずかな安定化
COT 軌道と 5d 軌道からなる e2g 軌道のエネルギーは Fock 演算子を使って次のように表される。
=

+

+2

(1)

第 1 項、第 2 項はそれぞれ錯体環境中における COT 軌道と 5d 軌道のエネルギーの寄与を示す。
第 3 項は共有結合性に由来する軌道相互作用の寄与を示すが、Ln-COT 間距離が短縮すると COT
軌道と 5d 軌道の重なりは増加するはずであり、この項は安定化に寄与すると考えられる。図 3には 2 つの軌道の重なり積分の変化を示した。Gd//Ln 計算によると、結合長が短縮するだけなら重なりは確かに増加するが、実際は La から Lu にかけて 5d 軌道も Ln 収縮するため（Ln//Gd 計算）、2 つの効果が打ち消しあうことで重なり積分は Ln 依存性をもたない。したがって e2g 軌道の安定化には軌道相互作用項は寄与しないことになる。次に第 1 項と第 2 項の寄与を評価するた
め、Ln を+3 の点電荷に置き換えた計算

0.32と、COT を-2 の点電荷に置き換えた計0.31

果が重要となることが分かった。
Ln-COT 間距離が短縮するにつれ、Ln3+に接近することで COT 軌道のエネルギーは低下し、一方 COT2-が接近すること
で 5d 軌道のエネルギーは増加する。しかしはoverlap integral算を行ったところ、どちらも静電的な効
0.30
0.29
0.28

Gd//Ln
Ln//Ln
Ln//Gd

0.27

より 10 倍程度大きいため、

e2g 軌道の主成分である COT 軌道のエネ

La Ce PrNdPmSmEuGdTbDyHoErTmYbLu

図3. e2g軌道における重なり積分の変化。

ルギーの低下が支配的となって e2g 軌道のエネルギーは低下する。このように e2g 軌道の安定化の主要因は、錯体内の大きな電荷分離に起因する静電的な効果によることが分かった。
以上のように、e2u 軌道と e2g 軌道では軌道エネルギーの Ln 依存性の中身が大きく異なるため、
図 1 のように異なった Ln 依存性を示すといえる。
【文献】[1] T. Kurikawa et. al.,JACS.,120
120,
120 11766 (1998). [2] M. Dolg et. al.,Theor. Chim. Acta., 75,
75 173(1989). [3] R. Takegami et. al., CPL,403
403,169(2005). [4]R. Takegami et. al., JPCA, 109,403
109 2476 (2005).

3P117
ICN 分子の光解離生成物 I+CN 間の相互作用ポテンシャルに関する理論的研究(慶大院・理工) ◯ 池崎 智哉、太田 悠介、藪下 聡
Theoretical study on the interaction potentials between the photoproducts of
I+CN in the photodissociation of ICN(Keio Univ) ◯Tomoya Ikezaki, Yusuke Ohta, Satoshi Yabushita【序】ICN 分子の A-band(𝜆 ≈ 230~310 nm) 励起による光解離反応には、次の 2 種類の生成物チャンネル
ICN + ℎ𝜈 → I( 2P3⁄2 ) + CN(X 2Σ + ), high 𝑁 (𝑁は CN の回転量子数)→ I ∗ ( 2P1⁄2 ) + CN(X 2Σ + ), low 𝑁が存在する。CN は、Iチャネルで高い回転量子数に、I ∗ チャネルでは低い回転量子数に分布する。さらに CN の回転準位は、その回転運動の角運動量と電子スピンの相互作用により、𝐽 = 𝑁 ± 1/2 と表現される 2 つの微細構造準位𝐹1 , 𝐹2 に分裂するが、特に興味深い点は、Iチャネルにおける𝐹1 , 𝐹2 準位分布
が非統計的で、さらにその割合が励起波長や𝑁の大きさに応じて変化することである。これは、ICNの解離反応経路が、その途中で起こる非断熱遷移によって複数のポテンシャル面（PES）に分岐し、それぞれの PES 上での核の de Broglie 波が量子干渉することに起因すると考えられているが[1,2]、詳
細は不明である。ここで、この de Broglie 波同士の量子干渉が生成物の微細構造準位分布を通して観測されるためには、非断熱遷移領域を通過後、解離極限に至るまでの PES は断熱的に振る舞うことが必要である。Chen[3] は、比較的長距離領域における相互作用ポテンシャルは、主に I の電子四重極子モーメントと CN の双極子モーメントの間の静電的な相互作用で記述できることを指摘している。I原子の四重極子モーメントは、基底状態(𝐽 = 3⁄2)だけが有限な値を持ち、I ∗ チャネル(𝐽 = 1⁄2) では 0である。本研究では、これらの特徴を検討するために特に漸近領域の PES の振る舞いを調べた。

【理論・計算方法】
ICN の励起状態の計算には COLUMBUS プログラムを用い、 1Σ + 、 1Π𝑥 、 3Π𝑥 、 3Σ + 、 3Π𝑦 、 3Π𝑦 の
6 つの spin-free(SF)状態の CI 波動関数を基底とした縮約スピン軌道(SO)CI 法を用いた。また、漸近領域での CN の双極子モーメント𝜇CN と I の四重極子モーメント𝑄I との相互作用のモデルポテンシャ
ルは以下の(1)式で計算した。ここで𝑅 は z 軸上においた CN の重心と I の間の距離を表している。

𝑉(𝑅) = (3𝜇CN𝑄I⁄2𝑅4 ){cos 𝜃𝑎 (3 cos2 𝜃𝑏 − 1) − 2 sin 𝜃𝑎 sin 𝜃𝑏 cos 𝜃𝑏 cos(𝜙𝑎 − 𝜙𝑏 )}(1)

𝜃𝑎 , 𝜃𝑏 は、z 軸に対する、双極子と四重極子の天頂角、𝜙𝑎 , 𝜙𝑏 は、ICN が xz 平面上で屈曲型構造をとるときの x 軸からの方位角である。𝜇CN には実験値[4]を、𝑄I は最近の理論値[5]を用いた。

【結果・考察】(1) 直線構造における C-I 間の伸縮に対するモデルポテンシャルについて
𝑅CI = ∞ の解離極限で縮重したエネルギーからの変化量の絶対値を𝛥𝐸(Σ + ), 𝛥𝐸(Π) とおき、横軸に
log 𝑅、縦軸に各log 𝛥𝐸 をとると、(1)式のモデルポテンシャルのグラフは、−4の傾きを持つ直線となる。Fig.1 をみると、スピンフリー(SF)状態の CI 計算値はΣ状態・Π状態ともに、それぞれのモデルポ
テンシャルの直線の上に重なっていることがわかる。同様に Fig.2 をみると、SO 相互作用を考慮した
CI 計算値は、I チャネルに相関する 1Σ0++ , 3Π2 , 3Π0− , 1Π1 に注目すると、SO 相互作用を考慮したモデル
PES の直線に乗っていることがわかる。また、I が四重極子をもたないI ∗ チャネルに相関する 3Π0+ , 3Σ1+は、直線から大きく外れ、別の弱い遠距離相互作用で表現されることが示唆される。したがって、
15 𝑎0 < 𝑅CI < 70 𝑎0 の領域では、CN の双極子と I の四重極子との相互作用が 1Σ0++ , 3Π2 , 3Π0− , 1Π1 の PESの振る舞いを支配していることが確認できた。

(2) 分子の変角に対するモデルポテンシャルについて
続いて、角度𝜃𝑎 , 𝜃𝑏 を変化させたときのモデルポテンシャルを計算した。Σ + , Πx 状態に関しては、角度𝜃𝑎 , 𝜃𝑏 に相関関係がある。与えられた𝜃𝑎 に対して、最安定あるいは最不安定な四重極子の主軸の向き
𝜃𝑏 を求めるため、変分的に考えて、

Figure 1 𝑹𝐂𝐈 の伸縮に対する SF 状態の PES とモデルの比較(𝜽𝒂 = 𝟎)(𝜟𝑬は絶対値をとっている)∂𝑉(𝑅, 𝜃𝑎 , 𝜃𝑏 , 𝜙𝑎 , 𝜙𝑏 )∂𝜃𝑏

Figure 2

𝑹𝐂𝐈 の伸縮に対する SO を含む場合のPES とモデルの比較(𝜽𝒂 = 𝟎)（𝜟𝑬は絶対値をとっている）

3= (3𝜇CN 𝑄I⁄𝑅4 ) { cos 𝜃𝑎 sin 2𝜃𝑏 + sin 𝜃𝑎 cos 2𝜃𝑏 cos(𝜙𝑎 − 𝜙𝑏 )} = 02→ tan 2𝜃𝑏 = −(2⁄3) tan 𝜃𝑎 cos(𝜙𝑎 − 𝜙𝑏 )(0 < 𝜃𝑎 , 𝜃𝑏 < 𝜋), (cos(𝜙𝑎 − 𝜙𝑏 ) = −1)

(2)
(3)

の相関があることがわかった。ここで得た 𝜃𝑎 , 𝜃𝑏 の関係をモデルポテンシャルに組み込んで、CI 計算の値と比較した。まず Fig.3 より、Σ + , Πx , Πy の SF 状態のモデルポテンシャルはいずれも CI 計算のものと良く一致していることが分かる。さらに SO を考慮した Fig.4 より、モデルポテンシャルは、
1A′ ~4A′ の PES の変化量の絶対値において良い一致を示し、四重極子を持たない5A′ ~6A′ はほとんど変化しないことが分かる。
以上をまとめると、ICN 分子の光解離反応の PES の振る舞いは、𝑅CI が 15 𝑎0 以上の領域においては、CN の双極子と I の四重極子との間の弱いながらも有限な大きさの静電的な相互作用によって支配され、その領域における断熱的な解離過程を保証する。このため、𝑅CI = 10 𝑎0 前後における非断熱相
互作用による干渉効果の結果が生成物 CN の微細構造準位分布になめらかに「転写」されるものと考えられる。

Figure 4
Figure 3変角に対する SO を含めた場合の PES とモ変角に対する SF 状態の PES とモデルの比

デルの比較(𝑹𝐂𝐈 = 𝟐𝟎. 𝟎 𝒂𝟎 )（𝜟𝑬は絶対値を較 (𝑹𝐂𝐈 = 𝟐𝟎. 𝟎 𝒂𝟎 )とっている）【参考文献】[1] H.Joswig et al., 𝐹𝑎𝑟𝑎𝑑𝑎𝑦 𝐷𝑖𝑠𝑐𝑢𝑠𝑠. 𝐶ℎ𝑒𝑚. 𝑆𝑜𝑐., 82, 79(1986). [2] 太田悠介, 慶應義塾大学大学院 修
士 論 文 , 2009 年 度 . [3] K.M.Chen, 𝐽. 𝐶ℎ𝑒𝑚. 𝑃ℎ𝑦𝑠. , 126, 214307(2007).[4] R. Thomson 𝑒𝑡𝑎𝑙.,𝐶𝑎𝑛.𝐽.𝑃ℎ𝑦𝑠,𝟒𝟔, 2815(1968). [5] Q.K.Timerghazin 𝑒𝑡 𝑎𝑙.,𝐽.𝐶ℎ𝑒𝑚.𝑃ℎ𝑦𝑠.,𝟏𝟐𝟒, 034313 (2006).

3P118
マルチカノニカル法による RNA ポリメラーゼ C 末端領域の研究（近大先端技研）○米澤 康滋
Study on the C-terminal Domain of RNA Polymerase II by using
Multicanonical Molecular Dynamics Simuation.
（Kinki.Univ. Inst.Adv.Tech.）○Yasushige Yonezawa
In eukaryotic cells, the C-terminal domain (CTD) of the RNA polymerase II plays animportant roles to regulate the transcription of genes. The CTD consists of 26-52 tandem
heptapeptide repeats with the consensus sequence YSPTSPS. The structure of the CTD in
the RNA polymerase has not been observed by X-ray crystal diffraction experiment,suggesting the CTD has intrinsically undefined structure. The serine residues at secondand fifth position in the CTD consensus are known to be major phosphorylation sites. Thus,
the regulated phosphorylation and de-phosphorylation of the CTD switch theconformation states and carry out a vital role in the recruitment and assembly oftranscription complexes. However, so far, the nature of conformational states of the CTD
is never elucidated in experimentally and theoretically. Knowledge of the conformation
states corresponding to the different phosphorylation must be useful for understanding
the CTD mechanisms.
Here, we have tackled the issue by using an advanced computational methodology.
Recently, molecular simulations are considered to be a useful tool for supporting
experiments even in molecular biology. Indeed, molecular simulations provide atomic
detail of protein and peptide under a variety of physiological conditions. However, proteinsand peptides have the huge number of degree of freedoms. Then, normal standard
molecular simulation hardly searches the conformation space sufficiently. To address thisdifficulty, we employed an extended ensemble method, multi canonical molecular
dynamics simulation (McMD). McMD is an umbrella sampling method using thetemperature space as a reaction coordinate, significantly increasing the conformational
search ability.
To investigate phosphorylation effects on the CTD conformation, we studied foursystems: unit-repeat-heptapeptide (h-peptide) without phosphorylation, h-peptide

phosphorylated on 2Ser, h-peptide phosphorylated on 5Ser, and h-peptide phosphorylated
on simultaneously 2 and 5Ser. The phosphorylated Ser residue which is absent in thestandard force field parameter was created according to the amber general procedure.
We then conducted the McMD simulations for the unphospholylated andphosphorylated h-peptides centered at cubic box in the periodic boundary conditions with
explicit waters and added ions for neutralizing the systems as seen in Fig. 1. All of thesystems were potential energy minimized for reducing hard contacts among atoms andthen constant atom-number, temperature, and volume simulations were conducted until
reaching the equilibrium condition, followed by constant atom-number, temperature,pressure simulations.
After fixed the volumes to their averaged vales, the upper limit temperature of the
McMD simulations was set to 600K which is found to ensure efficient conformation
searches for those systems. McMD simulation needs an accurate density-of-state toperform a random walk in the temperature space. To build the density-of-state accurately,
it is necessary to iteratively carry out preparation molecular dynamics-runs until thedensity-of-state converged. After over 20-times preparation-runs, we obtained theaccurate density-of-state for all systems. Snap-shots of atomic trajectory of the systemwere taken at every 10psec interval for analysis. All of the simulations are performed byusing myPresto program package.
We then analyzed the h-peptide system from the McMD simulations. First of all, theMcMD data was re-weighted to obtain ensemble data at room-temperature. Using the reweighted data, we study the conformation space of the h-peptides. Backbone dihedral
angles, φ and ψ are well suitable to characterize a peptide conformation, so we analyzed
the thermodymanically conformational properties dependence on the phosphorylation
with respect to the dihedral angles space. The details of procedures of the analysis results
would be explained and the biological significance including cis-trans isomerase enzyme
Pin1 would be also discussed.

3P119へリックス・ストランドレプリカ交換法によるペプチドの分子動力学シミュレーション（分子研，総研大）○奥村久士，伊藤暁

Helix-strand replica-exchange molecular dynamics
simulation of a peptide（Inst. Mol. Sci.，Sokendai）○Hisashi Okumura and Satoru G. Itoh【序】生体分子など複雑な自由エネルギー地形を持つ物質の分子動力学シミュレー
ションを素朴におこなうと，自由エネルギー極小状態にトラップされてしまい広い
範囲の構造を探索できない．そこで近年，マルチカノニカル法やレプリカ交換法な
どの拡張アンサンブル法がよく用いられるようになってきた．
多次元レプリカ交換法[1]（あるいはハミルトニアンレプリカ交換法[2]）では温度の代わりにポテンシャルエネルギーのパラメーターを交換する．そこで何を交換する
パラメーターに設定するかが問題になる．
【方法】最近，榮と岡本は 2 面角ポテンシャルエネルギーに新しい項を付け加える
ことによりαへリックスまたはβストランド構造を多く再現できるように力場を変形した[3]．そこで，このエネルギーに係数をかけ，その係数を交換する新しい多変数レプリカ交換法「へリックス・ストランドレプリカ交換法」開発した（図 1）[4]．

図 1 αへリックス・βストランドレプリカ交換法の概念図．

レプリカ 1 がスケーリングパラメータi を持ち，レプリカ 2 がスケーリングパラ
メータj を持っている状態を状態 X とする．この状態から 2 つのレプリカ間のパラメーターを交換した状態 X’への遷移はメトロポリス判定を用いて以下の確率で行う．

for   0,
1,w( X  X ')  
exp(), for   0

ここで



 



    V1 ( x j )  V1 ( xi )  V2 ( x j )  V2 ( xi ) である．V はパラメーターに対応するアンブレラポテンシャルで l ( x)  0 2   l ( x)  0 2 V  x     exp 

2 2
l 1


N aminoとした．(0, 0)はアンブレラポテンシャルの中心でαへリックス構造については (0,0) = (-57◦, -47◦)，βストランド構造については (0, 0) = (-130◦, 125◦)とした．
この方法をαへリックスとβヘアピン構造の両方をもつデザインペプチド[5]の分子動力学シミュレーションに応用した．比較のために温度を交換する通常のレプリカ交換法も実行した．
【結果と考察】その結果，αへリックス・βストランドレプリカ交換法の方が広い
構造空間をサンプルし，より長いαへリックス構造やβヘアピン構造を得ることができた．また RMSD を反応座標とした自由エネルギー曲面を計算した．自由エネル
ギー曲面における極小状態からデザインペプチドのフォールディング経路を提案した（図 2）．

図 2 デザインペプチドのフォールディング経路．

【参考文献】[1] Sugita Y, Kitao A, Okamoto Y: J. Chem. Phys. 113 (2000) 6042-6051.
[2] Fukunishi F, Watanabe O, Takada S: J. Chem. Phys. 116 (2002) 9058-9067.
[3] Sakae Y, Okamoto Y: J. Phys. Soc. Jpn. 75 (2006) 054802.
[4] Okumura H, Itoh S G: Phys. Chem. Chem. Phys., in press.
[5] Araki M, Tamura A: PROTEINS 66 (2007) 860-868.

3P120水和電子の結合経路(首都大院・理工)○北山 清章，田澤 遼太郎，橋本 健朗

Bond path of hydrated electron(Tokyo Metropolitan Univ.) ○Kiyoaki Kitayama, Ryotaro Tazawa, Kenro Hashimoto【序】水和電子は古くから研究され、水クラスター負イオンやアルカリ原子の水和クラスターを材料にした電子スペクトルや理論計算など多数の報告がある。近年では、余剰電子(SOMO)の広がり、局在化と SOMO に水の水素が向けられる構造から Tsurusawa, Iwata が名付けた電子―水素結合(OH{e-}HO)で水和電子の状態が語られることが多い。一方化学結合の解析では、全電子密度の極大点(Attractor)、鞍点(Bond Critical Point, BCP)とそれらを結ぶ稜線である結合経路(Bond Path)に着目する Bader の Atoms-in-Molecules (AIM)理論が広く用いられてきた。本研
究では、AIM 理論の意味で電子―水素結合が存在するのかを中心に、水和電子の Bond Path を
幾つかの系で比較しながら解析することを目的とした。
【計算】Attractor と BCP は共に全電子密度が空間座標に関して極値をとる点で、二階微分行列の三つの固有値が前者では全て負、後者は一つが正、二つが負となる。Attractor は通常核上にあるので、核から半径を変えつつ同心球面で電子密度を計算して極値方向を辿り、Bond Path とBCP を探索した。OH や水素結合の Bond Path は、分子構造から方向を予測して部分的な球面探索に抑えて計算を効率化した。電子―水素結合の領域も SOMO から予想できるが、電子密度変化が広い範囲で小さいので、半径、球面上の点の刻み幅を小さくし、また多数の立方格子上で電子密度とその勾配を網羅的に計算する方法を組み合わせた。得られた BCP で Hessian 行列の固有値を確認した。電子状態計算には Gaussian03,09 を用い、AIM のプログラムは自作した。結果に示す分子は全て拡張基底を用い電子相関を考慮して構造計算した。図には Bond Path は SOMOと対応させるため HF 法の電子密度を用いたが、電子相関を考慮しても大きくは変わらなかった。
【結果と考察】 [H3+]電子－水素結合は核がない場所に存在する極大(NNA)のなす結合の一つで、H3+はその特徴を示す好例である。図１
は最安定構造での結果で、3 本の Bond Path が各 H 原子と中央の NNAを結び、3 つの BCP は NNA の近傍にある。一方、核を結ぶ Bond Pathは無い。正三角形を保ったまま核間距離を伸ばすと、中央の NNA が消え、中央寄りに曲がった Bond Path が水素核間に現れた。
[(H2O)6-]水六量体負イオンの SOMO がクラスター内部にある構造(Internal State)と表面にある構造(Surface State)を計算した（図 2）。
内包構造では SOMO の位置に NNA があり上下の水三量体の水素結合していない水素から Bond Path が伸びている。一方、表面構造では

図 1. H3+ の電子密度(HF)の等高線と NNA（中央）、BCP（黒丸）、Bond Path.
構造は MP2/aug-cc-pVTZ.

OH や水素結合は見られるものの、NNA も電子―水素結合の Bond
Path もなかった。SOMO の電子密度だけ見れば極大はあるものの極小さく、全密度の変化の中に埋もれている。

図 2. 水六量体負イオンの SOMO（HF/aug-cc-pVDZ による 0.04 Bohr-3/2 の等値面）と Bond Path,
BCP, NNA. 左が Internal State、右が Surface State（最安定構造）. 構造計算は MP2 による。

[Na(H2O)6] 図 3 に Na 原子が六つの水に囲まれたクラスターの SOMO と Bond Path を示した。Na
から離れた不対電子を表す SOMO に、水素結合していない水素が向けられている。しかし、その領域には NNA はなく BCP があり、異なる水分子に属する
水素原子間を繋ぐ水素―水素結合の Bond Path が描けた。また、Na(H2O)4 と Na(H2O)8 も同様で、電子―水素結合の Bond Path が描けたものはなかった。
Na+(H2O)8 の Bond Path を調べると、中性で水素―水素結合を作った水から、クラスター内部によっ

図3

Na(H2O)6

の

SOMO

(HF/

6-31++G(d,p), 0.04 Bohr-3/2 の等値面とBond Path.構造は B3LYP による.

た位置に、O-O の Bond Path が得られた。余剰電子により Bond Path の始点が水素に移り、Bond
Path も SOMO の領域に移動することは、AIM の立場では結合と呼べないものの、電子―水素の相互作用が余剰電子の局在化に効いていることを支持していると考えられる。
[H3O(H2O)6,]

図 4 は、Na(H2O)6 と等電子系のH3O(H2O)6 ラジカルの Surface state の結果である。
H3O とは逆側のクラスター表面に SOMO が広がり、NNA も見つかった。3 つの BCP は NNA の極近傍にあり、図１の H3+に似た Bond Path が三つの水分子から NNA に向かって描けた。
本研究で電子―水素結合の NNA, BCP がどのようなときに出現するかを理論的に解明できたわけ

図 4. H3O(H2O)6 の SOMO(HF/aug-cc-pVDZ,
0.04 Bohr-3/2 の等値面)と NNA, BCP, Bond
Path.構造は MP2 による.

ではない。しかし、Surface state の(H2O)6-の SOMO 密度の極大点、Na(H2O)6 の BCP、H3O(H2O)6の NNA で SOMO 密度は 0.002, 0.003, 0.007 eBohr-3 で何れも極低い値だが、全電子密度に対する割合が 22, 81, 93%であることから SOMO が全電子密度の主成分の領域で BCP や NNA が見出されることを裏付けている。また (H2O)6-の Internal State や H3O(H2O)6 で NNA、電子―水素結合が得られることは、この結合が多くの水の水素が協調的に作る結合であることを示唆する。