3P101
Favorskii 転位を経由するゼルンボンの反応に関する量子化学的研究
(埼玉医大 1，近畿大院 2) ○土田敦子 1，福島美幸 2，宇高芳美 2，北山隆 2
DFT study on the reaction of zerumbon via Favorskii rearrangent
(Saitama-Med Univ.1, Kinki Univ.2) ○Noriko Tsuchida1, Miyuki Fukushima2,
Yoshimi Utaka2, Takashi Kitayama2
【 背 景 】天然から得られる化合物は古くから医薬品として用いられてきたが、化学の発展ととも
に合成による化合物に取って代わられ、医薬品そのものとして利用されることは少なくなった。
しかしながら、完全に医薬品との関わりを失ったわけではなく、豊富な資源量と多様な炭素骨格
から医薬品の出発原料として現在でも重用されている。ハナショウガから抽出されるセスキテル
ペン・ゼルンボンも医薬品への応用が期待される天然低分子の一つである。ゼルンボンは、二重
共役ケトンを含む 11 員環のトリエン骨格という他の天然物には見られない特徴的な構造をもつ。
近年、塩基性条件下で、二重結合の一つが臭素化されたゼルンボン誘導体はアミン系求核剤の付
加反応により容易に縮環することが明らかになった[1]。反応機構は生成物の構造から Favorskii
転位反応が関与していると予想されるが、その詳細はわかっていない。また、求核剤として酸素
系の化合物を用いると、縮環反応に加え開環反応が進行する。生成物の生成比は求核剤の種類に
よって異なり、HO-では化合物２が選択的に生成するのに対し、CH3O-では化合物１と２の生成比
が約 3:5 となる（Table１）。本研究では、量子化学計算を用いてゼルンボン誘導体の縮環・開環
反応の経路を追跡するとともに、求核剤の種類により生成物が異なる要因を明らかにすることを
目的とする。
Table 1. 塩基性条件下、臭素化ゼルンボンと各種求核剤の反応により得られた化合物の構造と生成比
Nu
O
Br
Br
O
O
Nu
Nu
reactant
Br
2
1
Ratio(%) 2)
Nu
Solvent
Temp.1)
Time (h)
CH3NH
DMSO
rt
0.5
92
CH3O
CH3OH
rt
1
29
48
HO
H 2O
rt
3
0
99
1) Rt was approximately 25
1
2
0
2) Calculation by 1H NMR
【方法】臭化物ゼルンボン誘導体, 求核剤（CH3O, CH3NH, OH）および 11 個の溶媒分子からなる
構造を計算対象とし、密度汎関数法（B3LYP/6-31G*）を用いて縮環および開環の２つの経路を追
跡した。また、求核剤 CH3NH および CH3O の反応追跡に対しては、計算負荷軽減のため DFT 法
と半経験的分子軌道法を組み合わせた 2-layer ONIOM 法（ONIOM(B3LYP/6-31G*:PM3)を用いた。
溶媒効果は連続誘電体モデル（PCM 法）を用いて考慮した。
【結果と考察】はじめに、縮環および開環反応の反応機構を追跡するため臭化ゼルンボンと CH3Oの反応を対象として計算を行った。得られたエネルギープロファイルと TS 構造の模式図を Fig. 1
に示す。縮環反応は、CH3O-の二重結合への付加（1_ts1）、シクロプロパンの形成によるα位の炭
素からの臭化物脱離 （1_ts2）、ケトンへの２つ目の MeO- の付加によるヘミアセタールの生成
（1_ts3）、水酸基からの水素脱離（1_ts4）、βおよびα’位の炭素間の結合形成によるエステル基の
生成と環の縮小（1_ts5）の 5 段階で進行する。反応の律速段階は、シクロプロパンの形成とα炭
素からの臭化物脱離で、その活性化エネルギーは 25.18（35.83-110.65）kcal/mol である。一方、開
環反応はケトンへの求核剤の付加（2_ts1）、C−C 結合開裂（2_ts2）、臭化物の脱離（2_ts3）の３段
階で進行し、反応の律速段階は C-C 結合開裂である。律速段階の活性化エネルギーは 22.81 kcal/mol
であり、縮環反応に比べ低く、開環反応が進行しやすい実験事実と矛盾しない。次に、得られた
経路をもとに、臭化ゼルンボンと CH3NH-の反応経路の追跡を行った。縮環反応では、シクロプロ
パノン形成（1_int2）までは CH3O-と同様であるが、その後ケトンへの CH3NH の付加による Nメチルアセトアミドの形成、縮環と全４段階で進行する。ケトンへの求核剤の律速段階は 1_ts2
で、活性エネルギーは 7.46 kcal/mol であった。一方、開環反応の経路は CH3O-と同様であり、C−
C 結合解離が律速（活性化エネルギー:11.19 kcal/mol）で、縮環反応よりもエネルギー的に不利な
経路であった。
Fig. 1 ONIOM(B3LYP:PM3)により求められた CH3O-の求核付加による臭化ゼルンボンの縮環反応
（左）および開環反応（右）のエネルギープロファイル
酸素系よりも窒素系求核剤で縮環反応が進行しやすい要因を明らかにするため、律速段階直前
の中間体（1_int1）の電荷に着目した(Table 2 )。この中間体はβ’位の炭素に求核剤が付加した構
造をもつ。CH3O-、CH3NH-のα’炭素の mulliken 電荷はいずれも-0.569 である。一方、α位炭素の
それはいずれも正であるが、CH3NH-が 0.116 と CH3O-の約２倍大きい。窒素系求核剤で、縮環反
応が進行しやすいのは、α位炭素の電荷がより正に傾くことでα’位炭素と結合形成がし易くなっ
たためではないかと考えられる。発表では、HO-の反応についても報告する。
Table ２. 縮環反応における中間体(1_int)の mulliken 電荷およびαC-Br 結合距離
O
α'
C
Br
C
Nu
CH3O
CH3NH
Mulliken charge
αC
α'C
0.059
0.116
-0.569
-0.569
Bond distance(Å)
αC Br
2.082
2.149
[1] 福島 美幸・ 宇高 芳美・高橋 一生・井福 壮・河合 靖・北山 隆，日本化学会第９５春季年会 (2015)
3P102
オーダーＮ法第一原理分子動力学計算の生体分子系に対する安定性
（理研 QBiC1, 東京理科大・理工 2, ロンドン大学 3, 物材機構 4）
○大塚教雄 1, 有田通朗 2, David Bowler3, 宮崎剛 4
Demonstration of an order-N first-principles DFT-MD code CONQUEST
（RIKEN QBiC1, Tokyo Univ. of Science2, UCL3, NIMS）
Takao Otsuka1, Michiaki Arita2, David Bowler3, Tsuyoshi Miyazaki4
【序】 生体機能の理解に向けて分子論的な視点である分子シミュレーションが広く活用されており、
特に古典力場を用いた分子動力学（MD）計算による研究が数多くなされている。近年、通常の古典
力場を用いた計算では難しいとされる化学結合の切断・生成や光に伴う応答といった電子状態変化
が大きく関わる現象には、電子状態が大きく変化しうる特定部分（数十～数百原子）を量子力学的に
取り扱い他は古典力場とする（QM/MM）事により対応している。しかしながら、数十ナノメートルスケ
ールの構造によって作られている特異的な環境や反応場に起こる生体分子特有な機構を電子状態
から詳細な知見を得る事は依然として計算科学に強く求められている。
本研究は、我々が開発してきたオーダーＮ法を用いた全原子第一原理密度汎関数法（DFT）計算
手法（プログラム名：CONQUEST）[1]を生体分子系に適用し、ナノサイズ構造の電子状態計算から
得られる知見を生体機能を理解する 1 つの手がかりとして提供する事を目的としている。これまでに
我々は、天然型の DNA 系に対するオーダーN 法第一原理 DFT 計算の数値検証から、我々のオーダー
N 計算手法が極めて高精度・安定である事を示してきた[2]。応用例として、非天然型塩基対を含んだ
DNA 系を対象として、経験的な van der Waals 相互作用補正を含めたオーダーN 法第一原理 DFT 計算
の数値検証を行ってきた[3]。
最近、我々はオーダーＮ法第一原理分子動力学計算法の導入に成功し[4,5]、巨大系に対する第
一原理 MD 計算の準備を行っている。今回は、DNA 系（約 3,400 原子）を例に、我々のオーダーＮ
法第一原理 MD 計算の計算安定性を議論する。
【理論的背景】 CONQUEST で用いられているオーダーN 計算手法は、密度汎関数法における一
体の密度行列を最適化する手法である。密度行列の非対角項が局所的であることから、非対角項に
対する cutoff 半径 RL を導入する事でオーダーN 法を実現している：
(r, r)
(r) K
i , j
i
i , j
j (r ),
K 3LSL 2LSLSL,
Li , j 0 for | Ri R j | RL
ここで、行列 L は補助密度行列であり、電子数一定の条件下で全エネルギーを最小にする密度行
列 L が求められる。また i (r) は、support function と呼ばれる各原子に局在した関数である。今回
も擬原子軌道（PAO）を使った計算を行っている。行列 S は重なり行列である。系の全エネルギーは、
cutoff 半径に対して変分的であるという利点を持つことから、オーダーN 法を導入することによって生
じる誤差を評価することが可能である。また、我々のオーダーN 計算は、各原子に対する力計算を可
能としており、構造最適化や分子動力学計算を行えるという特徴を持つ。
CONQUEST におけるオーダーN 法第一原理計算は、Niklasson 等によって提案された拡張
Lagrangian Born-Oppenheimer MD（XL-BOMD）法[3]を導入しており、我々のオーダーN 計算法に
適用するよう修正されている[4]：
2
LXBO ( X,X,R,R ) , LBO ( R,R ) Tr[ X 2 ]
Tr[( LS X ) 2 ]
2
2
ここで、 X は LS 行列の範囲を持つ疎行列、 は仮想電子質量、 は振動数パラメータである。
CONQUEST における計算条件等による詳細は参考文献[4]を参考されたい。
【計算と結果】 水溶液中 DNA の 10 塩基対モデル系のオー
ダーN 第一原理 MD 計算の結果を示す。モデル系は、PDB
ID: 1WQZ に AMBER で水分子を加えることにより作られた全
原子数 3439 原子（DNA: 634 原子、Mg: 9 原子、H2O: 932
分子 = 2796 原子）で、AMBER による平衡状態計算後のス
ナップショットの１構造（図 1）からオーダーN 法第一原理 MD
計 算 を 開 始 し た 。 電 子 状 態 計 算 は 、 PBE/SZP 、
non-self-consistent のオーダーN 計算、MD 計算の初期温度
を 300 K、時間刻みを 0.5 fs とした。
図 2 は、CONQUEST の XL-BOMD 計算による、ポテンシャ
ルエネルギー、運動エネルギー、それらの和である全エネル
ギーのトラジェクトリを示している。全エネルギーは保存して
Figure1. Hydrated DNA system
おり、20 fs 以降は安定したトラジェクトリを示す事が分かる。
(3,439 atoms).
当日は、オーダーN 法第一原理 MD 計算の詳細と温度制御法として速度スケーリング法を導入した
結果を示す予定である。
15
Total and Potential E. [Ha]
-19955
10
Total E
Kinetic E
Potential E 5
-19965
-19970
0
10
20
30
Time [fs]
40
50
0
600
400
Kinetic E. [K]
-19960
Kinetic E. [Ha]
800
200
0
Figure2. Energy profile of XL-BOMD by CONQUEST.
【参考文献】 [1] http://www.order-n.org/, see our web site. [2] T. Otsuka, et al., J. Phys. Condens.
Matter, 20, 294201 (2008). [3] 大塚教雄 et al., 分子科学討論会 2013, 4P084. [4] A. M. N.
Niklasson, et al., Phys. Rev. Lett., 97, 123001 (2006), Phys. Rev. Lett., 100, 123004 (2008), J.
Chem. Phys., 130, 214109 (2009). [5] M. Arita, D. R. Bowler, T. Miyazaki, J. Chem. Thory Comput.,
10, 5419 (2014).
3P103
動的分極率を用いた非局所励起状態計算法の開発：
波動関数理論によるアプローチ
（1 早大先進理工、2 早大理工研、3JST-CREST、 4 京大 ESICB）
○吉原詢也 1、吉川武司 1、中井浩巳 1-4
Non-local excited-state calculation based on dynamic polarizability:
Approach to wave function theory
(1Department of Chemistry and Biochemistry Waseda Univ., 2 Research Institute for Science and
Engineering, Waseda Univ., 3CREST, JST Agency, 4Elements Strategy Initiative for Catalysts and
Batteries, Kyoto Univ.)
○Junya Yoshihara1, Takeshi Yoshikawa1, Hiromi Nakai1-4
【序論】
通常の励起状態計算は、原子数とともに計算コスト（計算時間、メモリ量）が増大する。この
問題を解決すべく基底状態計算と同様、分割型理論の適用が試みられてきた[1]。励起中心が存在
し、環境の効果を取り込む場合、このような取り扱いは有効である。しかしながら、系内に拡が
った軌道間や電荷移動などの非局所励起の記述には適用できない。一方、当研究室では分割統治
法(DC)[2,3]に基づく動的分極率計算から間接的に励起情報を見積る手法を開発してきた[4,5]。こ
の手法では基底状態の情報に基づいているため、非局所励起を取り扱うことが可能である。これ
まで Hartree–Fock (HF)および密度汎関数理論(DFT)に基づく手法(DC-TDCPHF/DFT)が開発され、
その有効性が示されてきた。本研究では、電子相関法に対して同様の拡張を目指す。具体的には、
動的分極率を与える 1・2 電子励起結合クラスター線形応答 (CCSD-LR)[6]に DC 法を組み合わせ
ることを検討した。
【理論】
動的分極率は周波数とともに変化し、特定の周波数（極）
に対して発散する図１。極近傍の波動関数は励起状態の情報を含
み、基底状態 0 から励起状態 m への励起エネルギーに相当する周
波数m0 および振動子強度 fm0 はそれぞれ次式で与えられる。
m0
( )2 ( )2
( ) ( )
(1)
1
1
Fig. 1. Frequency dependent
f m0 (2 2 ) /
(2)
polarizability.
( ) ( )
ここで-, +は励起近傍の 2 点の周波数である。
動的分極率は双極子モーメント行列 d と応答密度 Dを用いて算出できる。
() Tr [D()d]
(3)
CCSD-LR 法における応答密度 Dは、分子軌道(MO)基底を用いて次式で計算される。
Dpq ( ) 0 1 aq† a p , T (1) 0 0 (1) aq† a p 0
(4)
ここで 0 は参照関数、aq†, ap はそれぞれ生成・消滅演算子である。演算子の添え字{p, q}は MO
を指す。T, Λ はそれぞれ励起演算子、脱励起演算子である。無摂動状態に対する T, Λ は通常の
CCSD 法により得られ、それらを用いて周波数に対する 1 次摂動 T(1), Λ(1)を CCSD-LR より求め
る。
DC 法への拡張では、DC-HF 法により得られた部分系の MO と軌道係数を用いる。Ts, Λs は DC-
CCSD 法[8]により、T(1)s, Λ(1)s は DC-CCSD-LR により計算する。
s
D pq
( ) 0 s 1 s a †p aq , T (1) 0 s 0 s (1) a †p aq 0 s
s
s
(5)
s
ここで 0 は全ての部分系に共通の Fermi 準位以下の部分系 MO からなる Slater 行列式である。
DC 法の部分系(L(s)=S(s)∪B(s))には中央領域(S(s))の他に緩衝領域(B(s))があるため部分系間の重
なりがある(L(s)∩L(s’)≠φ)。そのため、全系の応答密度行列の計算には、原子軌道(AO)基底の分
割関数 Ps が用いる必要がある。
1 ( S( s) and S( s))
D ( ) P D ( ) P 1 / 2 ( S( s) and B( s) or B( s) and S( s))
s
0
(other)
s
s
(6)
つまり、部分系の応答密度行列 Dに対して MO 基底から AO 基底への逆変換が必要となる。図
2 に上記の手続きをまとめた。
【結果】
現在、図 2 のアルゴリズムに基づきプログラム実装を行って
いる。予備検討として、(1)励起状態計算(EOM-CCSD)と動的分極
率計算(CCSD-LR)の比較および(2)CCSD-LR の応答密度の局所
性について調べた。
用いた系は push-pull 型のポリエン NH2-(C2H2)n-COOH である。基
底関数は 6-31G**を用いた。第 1 励起状態は、HOMO-LUMO 励起
に相当する。HOMO および LUMO はいずれもポリエンの π 軌道で
あるが、HOMO は電子供与基 NH2 側に LUMO は電子吸引基 COOH
側に偏っている(図 3)。
表 1 に n=3 のときの EOM-CCSD と CCSD-LR の励起エネルギー
および振動子強度の結果を示す。3 つの許容励起について、励起エ
ネルギーおよび振動子強度が一致することが確認された。
図 4 に AO 基底に変換した応答密度行列(n = 5)を示す。密度行列
の主要素は対角付近に存在し、距離に対して減少して
Fig. 2. Algorithm of DC-CCSD-LR.
いくことがわかる。このことは、応答密度行列が局所
的であることを意味し、分割型理論の親和性が高いこ
とを示唆する。
Table 1. EOM-CCSD and CCSD-LR excitation energies (eV)
and oschillator strengths of push-pull polyene.
EOM-CCSD
CCSD-LR
Oscillator
Oscillator
Energy(eV)
Energy(eV)
strength
strength
HOMO→LUMO
5.81
1.255
5.80
1.255
HOMO-1→LUMO
7.18
0.128
7.18
0.128
HOMO-5→LUMO
7.55
0.071
7.55
0.071
States
[1] T. Yoshikawa, M. Kobayashi, A. Fujii and H. Nakai, J.
Phys. Chem. B 117, 5565 (2013). [2] W. Yang, Phys. Rev. Lett.
66, 1438 (1991). [3] T. Akama, M. Kobayashi and H. Nakai,
J. Comput. Chem. 28, 2003 (2007). [4] H. Nakai and T.
Yoshikawa, in preparation. [5] T. Touma, M. Kobayashi and
H. Nakai, Chem. Phys. Lett. 485, 247 (2010). [6] M. Kallay
and J. Gauss, J. Mol. Struct. (THEOCHEM) 768, 71 (2006).
[7] H. P. Roy, A. Gupta and P. K. Mukherjee, Int. J. Quant.
Chem. 4, 75 (1975). [8] M. Kobayashi and H. Nakai, J. Chem.
Phys. 129, 044103 (2008).
Fig. 3. HOMO and LUMO of push-pull
polyene corresponding to first excited state.
Fig. 4. Response density of push-pull
polyene calculated by CCSD-LR.
3A104
金属ナノ粒子上での吸着活性化への担体効果に関する理論的研究
(1 早大先進理工, 2 早大理工研, 3 京大 ESICB, 4JST-CREST)
○出牛史子 1, 石川敦之 2,3, 中井浩巳 1-4
Theoretical research on the support effect of adsorption activation on metal nanoparticles
(1Waseda Univ., 2Waseda Univ. RISE, 3ESICB, 4JST-CREST)
○Fumiko Deushi1, Atsushi Ishikawa2,3, Hiromi Nakai1-4
【 緒 言 】不均一触媒の多くは、活性部分である金属微粒子と酸化物などの担体から構成され
ている。担体は触媒の比表面積の拡大や活性サイトを分散させるという効果が主であるが、
近年担体が金属の電子状態に影響を及ぼす、いわゆる SMSI (strong metal-support interaction)が
重要と考えられている。これらの要因により、触媒反応の活性は活性部位である金属部分の
特性(金属種や粒子サイズ)だけでなく担体にも大きく依存することが知られている。
本研究では触媒反応の一例として、自動車触媒や燃料電池の過程で重要となる CO の酸化
反応を対象として、理論的研究を行った。CO 酸化反応は Langmuir-Hinshelwood (LH)機構で
あることが実験的に知られており、反応活性を議論するうえでは CO と O2 両者の吸着エネル
ギーが重要となる。ここでは、CO および O2 吸着において担体による違いや吸着サイトによ
る違いに着目した。
【 計 算 モ デ ル ・ 計 算 条 件 】 本研究では、4 種類の酸化物に担持した Pt クラスターに対する
CO と O2 それぞれの吸着現象に対して、平面波基底に基づいた密度汎関数(density functional
theory: DFT)法を用いて計算を行った。計算モデルはスラブモデルを用い、MgO (100)面、ZrO2
(111)面、α-Al2O3 (0001)面、rutile TiO2 (110)面にそれぞれ cuboctahedron 構造の Pt13 を担持した
(Figure 1)。スーパーセルのサイズは、Pt クラスターより酸化物のサイズを十分に大きくする
ため、面内方向の 1 辺を 16 Å 以上とした。また、真空層は 20 Å 程度とした。Pt13 の吸着構造
の最適化においては、Pt13 の高さ方向のみを最適化した。CO および O2 は Pt13 の on-top サイ
(a)
(b)
(c)
(d)
Figure 1. Structure of Pt13 cluster supported on (a) MgO, (b) α-Al2O3, (c) ZrO2, and (d) TiO2.
(Upper: side view, Lower: top view, where only the lowest layer of Pt13 was shown)
トに end-on 吸着すると仮定した。DFT 計算はスピン分極を考慮し、交換相関汎関数に RPBE[1]
を用い、projector augmented wave (PAW)法により行った。エネルギーカットオフは 350 eV と
し、k 点は 1×1×1 (Γ 点)とした。すべての計算は、プログラムパッケージ VASP 5.3 を使用
した。
【 結 果 と 考 察 】 Pt13 の 4 種類の酸化物への吸着に対して、Table 1
に Pt クラスターの最下層における Pt 原子と周辺の O 原子との距離
を示す。原子ラベルの定義は Figure 2 の通りである。また、Pt13 の
Figure 2. Definition of
labels on Pt and O atoms.
吸着エネルギーは、MgO、Al2O3、ZrO2、TiO2 でそれぞれ−0.81、−1.20、
−1.47、−0.46 eV であった。これらの結果から、ZrO2
Table 1. Atomic distance of Pt and O atoms.
や MgO においては酸化物の O 原子が Pt を捕捉す
Pt-O distance (angstrom)
る役割を持つことが示唆される。
MgO
Al2O3
ZrO 2
TiO2
O1
2.712
2.603
2.098
2.735
O2
4.802
3.114
4.756
5.103
O3
3.758
4.730
4.099
4.654
O1
3.277
3.538
2.837
5.255
O2
3.229
2.639
2.970
2.597
O3
3.635
4.431
4.522
3.252
O1
3.798
2.944
3.235
4.845
O2
5.083
4.267
3.353
4.361
O3
2.275
2.832
2.330
2.536
Minimum
2.275
2.603
2.098
2.536
Pt1
次に、担持された Pt13 への CO および O2 の吸着
について検討を行った。CO および O2 については、
Pt クラスターの最上層(top layer)と中間層(mid layer)
への吸着を検討した。これらの吸着サイトの構造は
Pt2
Figure 2 の通りである。Table 2 に、担持した Pt13 ク
ラスターへの CO と O2 の吸着エネルギー(Eads)を示
Pt3
す。これらの結果から以下のことが明らかとなった。 (1) CO の吸着エネルギーは−1.92 ~ −1.42 eV と担体
依存性が小さい、(2) O2 の吸着エネルギーは−0.22 ~
+0.50 eV と担体依存性が CO の場合より大きい、(3) ZrO2 では CO、O2 ともに Pt クラスターの
吸着位置に対して大きな依存性がみられるが、Al2O3 などでは依存性は小さい。
(a) CO(top)
(b) CO (mid)
Table 2. Adsorption energies of CO and O2 on Pt13 cluster,
supported on various metal oxides.
MgO
(c) O 2 (top)
(d) O2 (mid)
Al2O3
ZrO 2
TiO 2
free Pt13
Figure 2. CO and O2 adsorptions on top
and mid layers of Pt13.
Adsorption site
Eads(CO) /eV
Eads(O2) /eV
top
−1.83
−0.10
mid
−1.42
−0.17
top
−1.72
+0.44
mid
−1.61
+0.49
top
−1.92
−0.22
mid
−1.65
−0.01
top
−1.60
+0.50
mid
−1.60
Not conv.
top
−1.37
+0.31
mid
−1.33
+0.27
【参考文献】 [1] B. Hammer, L. B. Hansen, J. K. Nørskov, Physical Review B, 59, 7413 (1999)
3P105
異核二原子分子の領域エネルギーの核間距離依存性に関する研究
（京大院工）○内藤 健人，埜崎 寛雄，市川 和秀，立花 明知
Study on interatomic distance dependence of the regional energy of
heteronuclear diatomic molecules
(Kyoto University) ○ Kento Naito, Hiroo Nozaki, Kazuhide Ichikawa, Akitomo Tachibana
Rigged QED 理論において、電子ストレステンソル密度演算子 τ̂eΠkl [1] は式 (1) で定義
される。
†
i~c ˆ
Πkl
l
0 l
τ̂e (x) =
ψ(x)γ D̂ek (x)ψ̂(x) − D̂ek (x)ψ̂(x) γ γ ψ̂(x)
(1)
2
ここで、ψ̂(x) は電子の 4 成分ディラック場演算子であり、x は時空座標を表わす (x =
(ct, ~r))。
電子ストレステンソル密度演算子 τ̂eΠkl に対して Primary rigged QED を用いて非相対
論的近似 (2), (3) を適用し、静電ハミルトニアンの定常状態に関して期待値を取ること
で、電子ストレステンソル密度 τeSkl [2] が得られる。
1
i~σ k Dk ψ̂L (x)
2mc
Spin-dependent 項を無視 σ k eµk Dµ · σ l eνl Dν ≈ (eµk Dµ )2
4 成分電子場の Small 成分の近似 ψ̂S ≈ −
τeSkl
(2)
(3)
∂ 2 ψi (~r) ∂ψi∗ (~r) ∂ψi (~r) ∂ 2 ψi∗ (~r)
∂ψi∗ (~r) ∂ψi (~r)
~2 X
∗
νi ψi (~r)
−
+
ψi (~r) −
(~r) =
4m i
∂xk ∂xl
∂xk
∂xl
∂xk ∂xl
∂xl
∂xk
(4)
ここで、ψ̂S と ψ̂L はそれぞれ 4 成分電子場 ψ̂ の small 成分と large 成分である。また、ψi
と νi は自然軌道とその占有数であり、{k, l} = {1, 2, 3} である。
また、Rigged QED 理論において定義されている電子運動エネルギー密度演算子 T̂e [1]
に対して同様の手順を踏むことで、電子運動エネルギー密度 nT [2] を得ることができる。
†
~2 1
ˆ
ˆ
†
2
2
ˆ
~ (x)ψ̂(x) + D
~ (x)ψ̂(x) · ψ(x)
(5)
·
ψ̂ (x)D
T̂e (x) = −
e
e
2m 2
nT (~r) = −
~2 X
νi [ψi∗ (~r)∆ψi (~r) + ∆ψi∗ (~r)ψi (~r)]
4m i
(6)
我々は、この電子ストレステンソル密度 τeSkl や電子運動エネルギー密度 nT などの密度
量を用いて、化学結合に対する理論的研究を行ってきた。これまでの研究から、電子ス
トレステンソル密度 τeSkl が化学結合の共有結合性・金属結合性の指標となりうることが
わかっている [3]。
また、電子運動エネルギー密度 nT が正である領域を分子の領域と定義し、さらにその
領域をラグランジュ面で分割することで、原子固有の領域を定義する研究も行われてい
る。nT がゼロとなる面は electronic interface と呼ばれ、分子の表面を表わすと解釈でき
る。ここで、ラグランジュ面はテンション密度 τeSk のセパラトリクスとして定義される。
テンション密度 τeSk は電子ストレステンソル密度 τeSkl の発散から得られるベクトル場で
P3
ある (τeSk = l=1 ∂l τeSkl )。また、テンション密度 τeSk がゼロとなる点はラグランジュ点と
呼ばれる。
本研究では、いくつかの二原子分子を対象に、この原子固有の領域それぞれについて
エネルギー密度を積分し、それぞれの原子の領域エネルギーを求める。特に、異核二原
子分子について、領域エネルギーの核間距離依存性から、化学結合についての議論を
行う。ここで、エネルギー密度 εSτ は電子ストレステンソル密度 τeSkl のトレースから得
P3
S
Skk
られ
(ε
(~
r
)
=
(~r)/2)、全空間で積分を行うと系の全エネルギー E に一致する
τ
k=1 τ
R S
( ετ (~r)d~r = E)。
なお、核間距離が大きい場合の電子状態を正確に取り扱うために、MOLPRO パッケー
ジを使用し、多参照配置間相互作用法 (Multireference configulation interaction method,
MRCI 法) を用いて自然軌道 ψi を求める。この自然軌道を使用し、我々の開発したプロ
グラムコード QEDynamics [4] を用いて領域エネルギーの計算を行う。
図 1 は平衡核間距離 (1.59Å) における LiH 分子、図 2 は核間距離 9.00Å における LiH 分
子の図である。図 1 において、Li 原子は (-0.42Å, 0.00Å)、H 原子は (1.17Å, 0.00Å) に位
置している。また、図 2 において、Li 原子は (-1.14Å, 0.00Å)、H 原子は (7.86Å, 0.00Å)
に位置している。カラーバーはエネルギー密度 εSτ を、矢印はテンション密度 τeSk を、青
の点はラグランジュ点を、黒太線はラグランジュ面を、青細線は electronic interface を
示す。
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
LiH
10
0.4
0.4
5
0.2
0
−0.2
y[Å]
y[Å]
LiH
0.2
0
0
−0.2
−5
−0.4
−1−0.5 0 0.5 1 1.5 2 2.5 3
−0.4
−10
−10
x[Å]
図 1: 平衡核間距離 (1.59Å) における LiH
−5
0
5
10
x[Å]
図 2: 核間距離 9.00Å における LiH
参考文献
[1] A. Tachibana, J. Chem. Phys. 115, 3497 (2001).
[2] A. Tachibana, J. Chem. Phys. 115, 8 (2001).
[3] K. Ichikawa, H. Nozaki, N. Komazawa, and A. Tachibana, AIP ADVANCES, 2,
042195 (2012).
[4] QEDynamics, M. Senami, K. Ichikawa and A. Tachibana
http://www.tachibana.kues.kyoto-u.ac.jp/qed/index.html
3P106
電子ストレステンソル密度の原子核間距離依存性による化学結合性の理論的研究 ○谷内公紀 1 、埜崎寛雄 1 、市川和秀 1 、立花明知 1
1 京大院工
[email protected]
電子ストレステンソル密度は、Rigged QED 理論から得られる電子ストレステンソル密度演算
子 τ̂eΠkl [1] に由来するものである。τ̂eΠkl は電子の 4 成分ディラック場演算子 ψ̂e (x) を用いて以下の
ように書ける。
[
]
(
)†
iℏc ˆ
Πkl
l
0 l
ˆ
τ̂e (x) =
ψe (x)γ D̂ek (x)ψe (x) − D̂ek (x)ψ̂e (x) γ γ ψ̂e (x)
2
電子ストレステンソル密度演算子 τ̂eΠkl に対して Primary Rigged QED を用いて非相対論的近似
（4 成分電子場の small 成分 ψ̂eS (x) を large 成分 ψ̂eL (x) で ψ̂eS (x) ≈ − 2m1e c iℏσ k Dk ψ̂eL (x) のよう
に近似し、スピン依存項を無視する）を適用し、静電ハミルトニアンのみで計算した電子状態で
期待値をとることで電子ストレステンソル密度は以下のように得られる。
τ Skl (⃗r) =
[
]
∂ 2 ψi (⃗r) ∂ψi∗ (⃗r) ∂ψi (⃗r) ∂ 2 ψi∗ (⃗r)
ℏ2 ∑
∂ψi∗ (⃗r) ∂ψi (⃗r)
νi ψi∗ (⃗r) k l −
+
ψ
(⃗
r
)
−
i
4m
∂x ∂x
∂xk
∂xl
∂xk ∂xl
∂xl
∂xk
i
ここで ψi (⃗r) は自然軌道を、νi は占有数を示す。電子ストレステンソル密度を対角化すると 3 つの
固有値 τeSii (i = 1, 2, 3)（τeS33 (⃗r) ≥ τeS22 (⃗r) ≥ τeS11 (⃗r）) が得られる。また、Rigged QED 理論から
運動エネルギー密度演算子 T̂e (x) も得ることが出来るが、この演算子に電子ストレステンソル密
度演算子と同様の手順を適用させれば運動エネルギー密度 nT (⃗r) を得ることが出来る。
T̂e (x) = −
ℏ2 1
·
2m 2
nT (⃗r) = −
)
(
(
)†
⃗ˆ e2 (x)ψ̂(x) + D
⃗ˆ e2 (x)ψ̂(x) · ψ̂(x)
ψ̂ † (x)D
ℏ2 ∑
νi [ψi∗ (⃗r)∆ψi (⃗r) + ∆ψi∗ (⃗r) · ψi (⃗r)]
4m
i
我々のこれまでの研究から、電子ストレステンソル密度の固有値 τeSii と最大固有値 τeS33 に対応す
る固有ベクトルが共有結合性や金属結合性を特徴付ける指標となり得ることが分かっている [2]。
共有結合性が強いとされる二原子分子の原子核間では τeS33 が正となり、これは引張応力が働いて
いることを示す。なおかつ対応する固有ベクトルが二つの原子核を紡錘型に結びつけるような構
造を有する。この構造のことをスピンドル構造という [3]。一方の金属結合性が強いとされる二原
子分子の原子核間では τeS33 が負となり、これは圧縮応力が働いていることを示す。また、それに
加えて 3 つの固有値 τeSii (i = 1, 2, 3) が縮退している傾向にあることが分かっている。
原子核間にまたがる Lewis 電子対は共有結合的なスピンドル構造を示すことが議論されている
[4]。たとえ金属結合的な二原子分子であっても、核間距離を平衡核間距離よりも大きくすると同
様の構造を示し、共有結合的であるとされている。
これまで我々は、主として平衡核間距離において電子ストレステンソル密度の計算を行ってき
た。共有結合的な二原子分子に関しては平衡核間距離においてスピンドル構造の存在が数値的に
確認されている。しかし、平衡核間距離よりも核間距離を大きくした場合に関しては、電子ストレ
ステンソル密度の計算は共有結合的・金属結合的な分子の何れについてもあまり行われていない。
本発表では、H から Cl までの原子からなる二原子分子に対し、平衡核間距離から核間距離を変
化させたときの電子ストレステンソル密度がどのように変化するかを計算する。また、核間距離
が大きいときの電子状態を正確に計算すべく、多参照配置間相互作用（MRCI）法を使用する。そ
れによって得られた電子状態を使用し、我々の開発したプログラムコード QEDynamics[5] を用い
て電子ストレステンソル密度を計算する。
また、イオン結合性についても議論する。共有結合性や金属結合性が強いとされる二原子分子に
関しては基底状態での計算を主に行ってきたが、イオン結合性が強いとされる二原子分子に関し
ては同一対称性の基底状態と第一励起状態について電子ストレステンソル密度の計算を行う。同
一対称性の二状態において、平衡核間距離では一方の状態が共有結合的で他方がイオン結合的と
なる。二状態のポテンシャルエネルギー曲線は非交差則より交わることはないが、エネルギー値
が非常に近くなる点が存在し、その付近を境に結合性が入れ替わることが知られている。結合性
が入れ替わった後もやはり、一方は共有結合性を示し、すなわちスピンドル構造が現れるものと
考えられる。そして、他方はイオン結合性を示し、スピンドル構造は現れないと考えられる。こ
れを数値的に確認する。
下図は Li2 の平衡核間距離 2.70Å（左図）と平衡核間距離よりも大きい核間距離 6.35Å（右図）
での電子ストレステンソル密度の図である。図中の緑破線は運動エネルギー密度が 0 となる領域
であり electronic interface(S) という。これにより原子や分子の表面を定義できる。図中の青線は
最大固有値が 0 となる等高線を示す。左図ではスピンドル構造を確認できないが、右図ではスピ
ンドル構造を確認できる。
Li2
4
0
4
0
2
−2e−05
2
−2e−05
0
−4e−05
0
−4e−05
−6e−05
−2
−8e−05
−4
y[Å]
y[Å]
Li2
−6e−05
−2
−8e−05
−4
−0.0001
−4
−2
0
x[Å]
2
4
−0.0001
−4
−2
0
x[Å]
2
4
図 1: Li2 の平衡核間距離 2.70Å(左図) と平衡核間距離よりも大きい核間距離 6.35Å(右図) での電
子ストレステンソル密度の最大固有値と、それに対応する固有ベクトル。原子核の位置は左図が
x = ±1.35Å、右図が x = ±3.175Å である。
参考文献
[1] A. Tachibana, J. Chem. Phys., 115, 3497 (2001)
[2] K. Ichikawa, H. Nozaki, N. Komazawa, A. Tachibana, AIP ADVANCES 2, 042195 (2012)
[3] A. Tachibana, Int. J. Quant. Chem. 100, 981 (2004)
[4] A. Tachibana, In Concepts and Methods in Modern Theoretical Chemistry, CRC Press,
Florida (2013) pp235-251
[5] QEDynamics, M. Senami, K. Ichikawa, A. Tachibana
http://www.tachibana.kues.kyoto-u.ac.jp/qed/index.html
3P107
多成分-NEB 法の開発と応用計算
(岐阜大・工*, 横市大院・生命ナノ**)
⃝宇田川太郎*, 鈴木机倫**, 立川仁典**
Development of the multicomponent molecular orbital-nudged elastic band
method and its application
(Gifu Univ.*, Yokohama City Univ.**)
◯Taro Udagawa*, Kimichi Suzuki**, Masanori Tachikawa**
【序】 我々はこれまでに、分子軌道の概念をプロトンのような軽い原子核にまで拡張した
新しいタイプの第一原理計算手法である、多成分分子軌道(MC_MO)法および多成分密度汎関
数理論(MC_DFT 法)を開発し[1]、様々な系における H/D 同位体効果を解析してきた。しかし
ながら、これまでに MC_MO 法および MC_DFT 法を化学反応解析へと応用した例は僅かしか
報告されていない[2,3]。これは、MC_MO 法において量子粒子として取り扱った原子核を含
む振動モードは従来の基準振動解析では得られないため、反応の遷移状態構造を定義するこ
とが困難なことが原因である。そのため、これまでの MC_MO 法を用いた化学反応解析は、
遷移状態構造が対称性から決定可能な系[2]、もしくは、予想される遷移状態構造付近の有効
ポテンシャル曲面を直接計算可能な比較的小規模な系[3]に限られていた。そこで本研究では、
MC_MO 法による化学反応解析を様々な系へ適用可能とするため、MC_MO 法と、エネルギ
ーの一次微分のみで遷移状態構造を求めることが可能な Climbing image-nudged elastic band
(CI-NEB)法[4]とを組み合わせた MC_MO-CI-NEB 法を開発し、種々の分子内水素移動反応に
おける速度論的同位体効果の解析を試みた[5]。
【理論】 NEB 法では、 ! ! , !! , ! ! , … , ! ! で表される N+1 個のイメージを仮想的なバネでつ
なぎ、R1 から RN-1 までのイメージを最適化することで反応経路を得る。ここで、Ri は各イメ
ージの位置ベクトルであり、始点 R0 と終点 RN はそれぞれ反応物と生成物に対応するエネル
ギー極小値構造である。NEB 法ではあるイメージにかかる力は、次式により計算される。
!! = !!!
||
− ∇! ! !
!
(1)
ここで右辺第 1 項はイメージ間を繋ぐバネに働く力の、反応経路に沿った成分である。本研
究では、Henkelman らの評価式[4]に基づいて!!! || を評価した。
【計算方法】 本研究では、水素移動反応 XHCHCHCHY ↔ XCHCHCHYH (X, Y = O, NH, or
CH2)における H/D 同位体効果について、MC_MO-CI-NEB 法を用いて解析した。電子状態計
算には HF レベルの MC_MO 法を、電子基底関数には 6-31G**を用いた。移動する水素原子
核のみを波動関数として取り扱い、原子核基底関数には、s 型ガウス関数(exp −! ! − !
!
)
を 1 つ用いた。また、反応中の原子核分布(軌道指数α)の変化についても検討した。
【結果】図には、(X,Y) = (O,NH)について、軌
!"($%
道指数をα(H) = 24.1825, α(D) = 35.6214 に固定
して解析した結果を示した。MC_MO-CI-NEB
,-.,/%012%
,-.,/%032%
-456758459:%1;%
!"'$%
RON [Å]
法を用いることで、これまで困難であった非
対称な反応の遷移状態を求めることが可能と
なった。また、通常の HF 法と、MC_MO-CI-NEB
!"$$%
!"&$%
法で得られた反応経路は明らかに異なってお
り、移動する水素原子核の量子効果が、遷移
状態構造のみならず反応経路にまで影響を及
ぼしている様子が明らかとなった。また、表 1
!"#$%
)*"!%
)+"'%
+%
+"'%
*"!%
!r [Å]
図 従来の!"法および#$%#&'$(')*+法により求めた,
-./01,2,-&/)!1における!体/,3体の最小エネルギー経路
には次式から計算した T = 298 K における kH/kD 値を示した。
!! /!! = exp − !!! − !!! /!"
表 1 より、MC_MO 法により移動
する水素原子核の量子性を直接
考慮することで、従来の HF 法に
比べて kH/kD 値が大きくなること
がわかった。これは、MC_MO 法
により、原子核の量子性の違いが
引き起こす構造緩和の効果を直
(2)
表1 !"法および#$%#&'!"法による! ! (! ) 値
! ! (! )
$456153/45789!"
#$%#&'!"
<=>
?=@
<=B
>=@
<=C
D=B
"40G70H
B=@
B=B
I1610J1
<=K
>=D
*&,$! @ .
"40G70H
B=<
>=B
I1610J1
<=E
E<=C
*A!,$!@ .
"40G70H
B=C
EE=C
I1610J1
<=E
E@=C
FE
9739@>LM9NO=9I=9I43P,9Q=9MR5/S,9"#$%#&'()**+(,-$.+,9?>>,9@<9*E>??.T=
*+,-.
*&,&.
*A!,A!.
*$! @ ,$!@ .
*&,A!.
)/0123/45
1:;
E@=@ FE
接的に取り込んだ結果であると
考えられる。また、反応中の原子核分布の変化を解析するため、各イメージにおいて原子核
基底関数の軌道指数を最適化しながら MC_MO-CI-NEB 計算を行った。ここでは、(X,Y) =
(O,O)における活性化エネルギーと kH/kD 値を表
2 に示した。原子核基底関数中の軌道指数を最適
表2 単一の軌道指数+,'-.および最適化した軌道指数+5$6.に
より求めた活性化エネルギー";< と !! )!( 値
化することで、各状態のエネルギーが安定化し、
!*+,'-.
!*+5$6.
" 94
結果として活性化エネルギーが H 体では 2.4
kJ/mol, D 体では 1.8 kJ/mol ほど小さくなり、kH/kD
94
";< *=>?)@5AB
!"#$%&'%#
("#$%&'%#
/012
341/
/713
/817
"/17
"41:
! ! )! (
01/
218
412
*!*+5$6.による値と ! (fix)による値の差
値にも影響を与えた。詳細は当日報告する。
【参考文献】 [1] M. Tachikawa, K. Mori, H. Nakai, K. Iguchi, Chem. Phys. Lett., 290, 437 (1998), T. Udagawa,
M. Tachikawa, J. Chem. Phys., 125, 244105 (2006). [2] Y. Kikuta, T. Ishimoto, U. Nagashima, Bull. Chem. Soc.
Jpn., 81, 820 (2008). [3] T. Ishimoto, M. Tachikawa, H. Tokiwa, U. Nagashima, Chem. Phys., 314, 231 (2005).
[4] G. Henkelman, H. Jónsson, J. Chem. Phys., 113, 9489 (2000), G. Henkelman, B. P. Uberuaga, H. Jónsson, J.
Chem. Phys., 113, 9901 (2000). [5] T. Udagawa, K. Suzuki, M. Tachikawa, ChemPhysChem, accepted.
3P108
ジェミナル理論およびHFB法を用いた
ジラジカル性分子の非線形光学応答計算
（北大院・理1，京大・ESICB2，阪大院・基礎工3）
○小林 正人1,2，福田 幸太郎3，中野 雅由3，武次 徹也1,2
Evaluation of non-linear optical response properties of
diradical molecules based on geminal and HFB theories
(Hokkaido Univ.1, Kyoto Univ.2, Osaka Univ.3)
○Masato Kobayashi1,2, Kotaro Fukuda3, Masayoshi Nakano3, Tetsuya Taketsugu1,2
【緒言】
近年、開殻一重項分子系の開殻性と物性の相関が注目されており、中でも中間的なジラジ
カル性を有する分子系が従来の閉殻系に比べて顕著に大きな非線形光学応答を示すことが明
らかになっている[1]。これまでジラジカル性を持つ一重項分子の計算には、broken symmetry
の DFT などが用いられてきたが、これはもちろんスピン対称性を満足しない。本研究では、
低コストでスピン対称性を満足した一重項ジラジカルの計算が可能なジェミナル理論と
Hartree–Fock–Bogoliubov (HFB)法に注目し、これらを用いて非線形光学応答物性を算定した。
【ジェミナル理論と HFB 法】
ジェミナル理論は、一電子軌道（MO）ではなく反平行スピンの電子対波動関数（geminal）
に対する平均場近似計算法である。その一つである反対称化強直交ジェミナル積（APSG）法
[2]では、2 電子多軌道の CAS 波動関数の反対称化積として全電子波動関数を表す。
APSG a1 a 2 ...a n 1 ... m vac
(1)
ここで i は geminal の生成演算子であり、自然軌道（NO）を用いて以下のように表す。
i C i a a
(2)
i
(1)式では n 電子を通常の MO で表わし、m 電子対を geminal で表わしている。APSG 法では、
CASSCF 法と同様に、(2)式の C i からジラジカル因子 y を見積もることが可能である。
ジェミナル理論と似た構造の波動関数を持つ理論として、HFB 法がある。HFB 法では、一
重項分子のエネルギーは、空間軌道を用いて次式によって与えられる[3]。
E HFB ( ) 2 h P
2
P P
K K 2 N S P
(3)
ここで P は HF 法で用いられるものと同じ一体の密度行列で、K はペア行列と呼ばれる。ζ は
Staroverov と Scuseria [3]によって導入された静的電子相関の強さを制御するパラメータであ
る。電子数一定の制約条件のもと、(3)式のエネルギーを P と K に関して変分的に求めるのが
HFB 法である。KS の対角要素は、対応する AO のラジカル占有性を表している。例えば最小
基底 H2 の解離極限の計算では、 Tr( KS ) 1 となるため、ジラジカル性の指標になる。我々は
最近、HFB エネルギー勾配表式を導出し、これを用いた構造最適化計算を可能としている[4]。
【HFB 法による p-キノジメタンの最適化構造と三次非線形光学応答】
HFB 法のパラメータ ζ を変化させて、ジラジカル性の
共鳴構造を持つ p-キノジメタン（図 1）の構造最適化計
4
算を行った。また、得られた平衡構造に対して、同じ計
z
2
算レベルで静電場印加時のエネルギーに基づく有限場
2
3
1
法を用いて第二超分極率 γzzzz を見積もった。基底関数系
には 6-311G**を用いた。図 2 上に構造最適化計算で得ら
4
3
1
Fig. 1. Structure of p-quinodimethane
れた C–C 結合長を、下に γzzzz の値を示す。ζ =
0.82 を少し下回ると K 0 の HFB 安定解は存在
せず、RHF 計算と同じ結果を与える。ζ が増大
すると、最適化された C1–C2 および C3–C4 距離
は徐々に伸張し、逆に C2–C3 距離は短縮する。
図 1 からもわかるとおり、結合交替が小さくな
り、ジラジカル性が増すことがわかる。これに
伴い、ζ が増大すると第二超分極率 γzzzz の値が
急激に大きくなることがわかる（CASSCF(8,8)
計算の結果は γzzzz = 7.4 × 104 a.u.）。 p-キノジメ
タンの C1–C2 距離を変化させることでジラジカ
ル因子 y をコントロールできることが知られて
いるが、図 2 で見られる γzzzz の増大は構造に起
因するものではない（CASSCF(8,8)計算で得ら
れた構造で固定して γzzzz の計算を行っても図 2
下とほぼ同じ結果が得られる）。ジラジカル性
の指標と考えられる Tr( KS ) の値は、ζ = 0.82 で
は 0.6 だが、ζ = 0.9 で 4.5、ζ = 1.0 で 7.6 と増大
する。これらのことから、HFB 法はジラジカル
電子状態を有効的に取り込むことができるこ
Fig. 2. Dependence of the optimized bond lengths
(upper panel) and the second hyperpolarizability
γzzzz (lower panel) of p-quinodimethane on the
parameter ζ of the HFB method.
と、また安定的に有限場法による超分極率計算に利用できることがわかった。現在の計算レ
ベルでは動的電子相関が欠如しているため、定量的な議論は難しいが、これを摂動論や DFT
計算と組み合わせて取り入れることも検討中である。
当日の発表では、HFB 法のほかに、APSG 法やその摂動補正計算である MP-MCPT 法[5]に
よる超分極率計算の結果についても報告する。
[1] M. Nakano et al., Phys. Rev. Lett. 99, 033001 (2007); J. Phys. Chem. A 109, 885 (2005).
[2] P.R. Surján, in Correlation and Localization (Springer, 1999), pp. 63-88.
[3] V.N. Staroverov and G.E. Scuseria, J. Chem. Phys. 117, 11107 (2002).
[4] M. Kobayashi, J. Chem. Phys. 140, 084115 (2014).
[5] M. Kobayashi , Á. Szabados, H. Nakai, and P.R. Surján, J. Chem. Theory Comput. 6, 2024 (2010); M.
Tarumi, M. Kobayashi, and H. Nakai, J. Chem. Theory Comput. 8, 4330 (2012).
3P110
多成分系第一原理手法による電子-陽電子対消滅率の高精度算定
（横市大院・生命ナノ） ○小山田隆行, 立川仁典
Accurate evaluation of the electron-positron pair-annihilation rate using multi-component ab initio methods
(Yokohama City Univ.) ○ Takayuki Oyamada, Masanori Tachikawa
【序】陽電子 e+ は電子と質量は同じだが, 電荷が電子と逆符号の正電荷を持つ反粒子である. 陽電
子が原子・分子に近づくと, 陽電子は原子核からは斥力を感じるが, 電子雲との間の引力により一
時的に束縛され, ごく短い時間であるが陽電子化合物を形成する. そして, 陽電子化合物中の電子
と陽電子が衝突すると γ 線を放出して対消滅する [1-3]. 広い分野で既に実用化されている陽電子
対消滅法では, この時の放出 γ 線を検出することで, 固体の構造欠陥・表面の解析といった物性研
究や, 陽電子断層撮影法を始めとする医療診断等を可能としている. このように応用的な観点から
も重要な電子-陽電子対消滅であるが, 陽電子化合物の寿命 (陽電子が原子・分子に束縛されてから
対消滅するまでの時間) が 10−9 sec 程度と極めて短いことから, 対消滅機構の詳細を実験的に解明
することは困難である. 従って, 陽電子対消滅法のさらなる高精度・高質化のためにも, 陽電子束
縛機構や対消滅機構を理論的に解明する高精度な第一原理手法の開発が強く求められている. こ
のような背景から我々は多成分系第一原理手法の開発を行っており, 本研究では特に, 電子-陽電子
対消滅率を高精度に算定し, 対消滅機構の解明に向けて研究を行った.
【計算方法】電子-陽電子対消滅率 Γ2 は, 電子-陽電子衝突確率を hδep i として次式で求められる [4].
対消滅率 Γ2 ≡ πα
4
c a−1
0 hδep i
Ne
E
D X
e
p 電子-陽電子衝突確率 hδep i ≡ Ψ δ(ri − r ) Ψ
,
(1)
i
ここで, α は微細構造定数, c は光速, a0 は Bohr 半径であり, 対消滅率の算定には, hδep i を高精度
に評価することが重要である. 対消滅率 Γ2 は, 陽電子化合物の寿命の逆数に比例する.
hδep i を求める方法の 1 つは次式の衝突確率密度 (contact density) ρep (r) を用いることである [5]:
Z
Ne
E
D X
integrate over r
e
p ρep (r) ≡ Ψ −−−−−−−−−→
ρep (r) dr = hδep i
(2)
δ(r − ri ) δ(r − r ) Ψ
i
ρep (r) は, ある空間座標 r に電子と陽電子を同時に見出す確率分布を表し, 分子中で電子-陽電子衝
突の起こり易い場所について有益な知見を与える. ρep (r) を空間座標 r に渡って数値積分すれば,
(2) 式のように hδep i が数値的に求まり, 従って, 対消滅率 Γ2 を算定できる. 衝突確率密度を数値的
に積分して hδep i を評価することは, 基底関数の関数型によらず柔軟に対応でき, プログラムの実
装上も比較的簡単である. 特に, 原子や直線分子の場合, その対称性を考慮することで, 高速かつ高
精度な数値積分が可能である. しかし, 一般の複雑な対称性の分子の場合, 3 次元グリッド上で数
値積分を実施する必要があり, 計算精度と計算時間の両立は困難となる.
Hartree-Fock(HF) 法や配置間相互作用 (CI) 法等では, hδep i は次式で求められる [6]:
Z
X X X X (ep)
(e) (p) (p)
hδep i =
Γµν,λσ Sµν,λσ , Sµν,λσ = φ(e)
(3)
µ φν φλ φσ dr
µ
(ep)
ここで, Γµν,λσ
(e) (e)
ある. φµ , φν
ν
λ
σ
は電子-陽電子 2 粒子密度行列であり, Sµν,λσ は電子-陽電子 4 中心重なり積分で
(p)
(p)
は電子基底, φλ , φσ は陽電子基底を表し, 本研究では基底関数として Cartesian
Gaussian を用いる. 4 中心重なり積分 Sµν,λσ は, (I) s 型 GTF のみ (floating GTF を含む) を用い
る場合と, (II) p, d, f · · · 型 GTF を用いるが関数中心が 1ヶ所の場合には積分公式により評価でき
る. しかし, 多中心で p, d, f · · ·型 GTF を用いる計算では, 漸化式等を用いる必要がある. 本研究で
は, (3) 式の 4 中心重なり積分の評価に Obara-Saika 法 [7] を用いた. また, 原子と直線分子につい
ては, (2) 式と (3) 式の両方で衝突確率 hδep i を評価し, いずれも良く一致することを確認している.
-4
Contact density (bohr -6)
5.0 x 10
【結果】 図 1 に [LiH; e+ ] の分子軸上の衝突確率
Contact density
密度 ρep (r) を (i) HF, (ii) 電子-電子相関のみ考慮
-4
for [LiH; e+]
4.0 x 10
した CISDTQ(e-e), (iii) 電子-陽電子相関のみ考慮
した CISD(e-p), (iv) Full-CI 計算で比較している.
(iii) CISD(e-p)
-4
3.0
x
10
(iv) Full-CI
核間距離は R = 3.20 bohr とし, H 核を原点とし
(ii) CISDTQ(e-e)
た. 電子基底は 6-311++G**, 陽電子基底は 14s4p
(i) HF
2.0 x 10-4
GTFs を用いた. 陽電子基底の関数中心と軌道指
数は HF 計算で最適化している. いずれの計算方
1.0 x 10-4
法でも, 電子-陽電子衝突は Li 核近傍よりも H 核近
傍で起こり易いことが分かる. これは平衡核間距
H
0.0
離付近の LiH 分子では, イオン性成分 Liδ+ Hδ− が
Li
-5
0
5
10
z (bohr)
主成分であり, 陽電子は H 核の外側領域に拡がっ
図 1：HF および各種の CI 法で求めた分子軸
た分布を持つためである [3,8,9].
+
図 1 の Full-CI の ρep (r) は HF の場合と比較し (z 軸) 上の [LiH; e ] の衝突確率密度 ρep (r).
て, ピークが 8 倍程度も高い. 1 ∼ 4 電子励起を考慮
表 1： 各種の計算方法による [H− ; e+ ] の
して電子-電子相関を高精度評価した CISDTQ(e対消滅率 Γ2 と全エネルギー Etot の比較.
e) でも, ρep (r) の高さは HF とほぼ同様である. 他
計算方法 Γ2 (ns−1 ) Etot (hartree)
方, 電子-陽電子相関のみを考慮した CISD(e-p) で
HF (a)
0.29750
-0.66693
は, ρep (r) の高さは Full-CI と同程度か僅かに大
(a)
Full-CI
1.19999
-0.77342
きい値となる. 対消滅率 Γ2 についても同様であ
(b)
XCHF
2.4800
-0.73750
り, Full-CI の Γ2 は HF 値よりも 4 倍以上大きく
(c)
VMC
2.46(5)
-0.78675(6)
なる. CISDTQ(e-e) の Γ2 は HF 値と同程度であ
(d)
SVM
2.47178
-0.78920
り, CISD(e-p) の Γ2 は Full-CI の値より僅かに大
−
+
(a) 本研究: e /e = 8s6p4d2f/8s6p4d2f GTFs
きい. このことから, 衝突確率密度や対消滅率の高
(b) Explicitly correlated Hartree-Fock: Ref.[5]
精度評価には, 電子-電子相関よりも電子-陽電子相
(c) Variational Monte Carlo method: Ref.[4]
(d) Stochastic variational method: Ref.[10]
関の方が重要であることが確かめられる.
表 1 には, 本研究の HF および Full-CI 計算で求めた [H− ; e+ ] の対消滅率 Γ2 と全エネルギー Etot
を, 他の高精度計算の結果と比較している. 電子基底と陽電子基底は, e− /e+ =8s6p4d2f/8s6p4d2f
GTFs を用いて, 軌道指数は Full-CI で変分最適化している. より小さな基底関数の Full-CI 計算
による [H− ; e+ ] の対消滅率は, Adamson らにより報告されていて [6], 本研究の Full-CI 計算の結
果と良く一致することを確かめている. HF 法で著しく過小評価されていた対消滅率は, Full-CI 法
により相関を考慮することで大幅に改善し, 4 倍程度も大きな値となっている. しかし, XCHF[5],
VMC[4], SVM[10] といった 2 粒子間距離に顕に依存した計算結果と比べると, Full-CI の対消滅率
は半分程度に留まっている. これは, 1 粒子基底関数を用いた CI 法では厳密解に対する収束が遅い
ことを示す典型例である. 特に, 対消滅率のような電子-陽電子間距離の短距離領域に強く依存し
た特性の高精度評価においては, 粒子間距離に顕に依存した 2 粒子基底関数を用いて, 電子-陽電子
相関を高精度評価することが本質的に重要である. このため本研究では, 現在, 電子-陽電子間距離
に顕に依存した多成分第一原理手法の開発を行っている.
[1] 陽電子計測の科学, (日本アイソトープ協会,1993). [2] M.Charlton and J.W.Humberston, Positron Physics, (Cambridge University Press, Cambridge, 2001).
al., Chem. Lett. 39, 1136 (2010).
[3] 立川仁典, 北幸海, 日本物理学会誌, 67, 33, (2012).
[5] C. Swalina et al., J. Chem. Phys. 136, 164105 (2012).
[4] Y. Kita et
[6] P.E. Adamson
et al., J. Phys. Chem. A 112, 1346 (2008). [7] S.Obara and A.Saika, J. Chem. Phys. 84, 3963 (1986).
et al., J. Chem. Phys. 135, 054108 (2011).
[8] Y.Kita
[9] T.Oyamada and M.Tachikawa, Eur. Phys. J. D, 68, 231 (2014).
[10] J. Mitroy, Phys. Rev. A 73, 054502, (2006).
3P111
1,2-ブタジエンの超高速失活過程に関する理論的研究
(北大院総合化学 1 , 北大院理 2 , 北大院工 3 , JST CREST4)
○佐藤壮太 1 , 原渕祐 2,4 , 飯窪亮 3 , 藤原丈久 3 , 関川太郎 3 , 武次徹也 2
Theoretical study on ultrafast relaxation process of 1,2-butadiene
(Hokkaido Univ.1 , JST CREST2)
○Sota Satoh1, Yu Harabuchi1,2, Ryo Iikubo1, Takehisa Fujiwara1, Taro Sekikawa1,
Tetsuya Taketsugu1
[email protected]
【研究背景】1,2-ブタジエンは紫外光照射により光解
離を示す基礎的な系として知られている。過去の研究
により、1,2-ブタジエンは光励起後、解離することな
く超高速無輻射失活し、解離過程は光失活過程の余剰
エネルギーによって基底状態で起こると考えられてい
る。そのため、基底状態への失活後の解離機構に着目
した研究がなされてきた[1]。一方で、光励起後の超高
速失活過程について着目した理論研究はほとんど無く、
図 1. 1,2-ブタジエンの光失活過程
その機構は明らかになっていない。
の概念図
近年、関川らにより、1,3-ブタジエンの 2 光子吸収
による光励起後の時間分解イオン化スペクトルが測定された[2]。これと同様の実験を 1,2-ブ
タジエンについても行ったところ、失活後のイオン化スペクトルピーク位置に変化が見られ
なかったことから、1,2-ブタジエンが、励起状態で解離することなく高速失活する様子が確認
できた[3]。また、励起直後に励起分子に由来するスペクトルピークを観測した。さらに、基
底状態への失活後、スペクトルのピーク強度が振動する様子が観測され、これは分子の振動
励起に起因すると考えられた。しかし、励起状態における 1,2-ブタジエンの運動は未解明な
部分が多く、実験スペクトルの起源を解明するのは困難であった。そこで本研究では、1,2ブタジエンの超高速失活過程の解明を目指し、理論計算に基づき詳細に解析した[3]。
【計算手法】1,2-ブタジエンに対し、第一励起状態(S1)のポテンシャル曲面の解析と ab initio
分子動力学(AIMD)計算を行った。計算レベルは CASPT2/cc-pVDZ とし、活性空間は π/π*軌
道 2 セットを含む 4 電子 4 軌道とした。ポテンシャル曲面(PES)の解析では、フランク-コンド
ン構造から the avoiding model function (AMF) 法を用いて最急降下経路を計算した。AIMD 計
算の時間幅は 0.2 fs で最大で 400 fs とし、50 本の古典軌道を流した。また初期構造のサンプ
リングにはボルツマン分布を用いた。電子状態計算には Molpro2010 を用い、併用して、円錐
交差構造と S1 の極小構造の最適化及び最急降下経路計算には GRRM11 プログラムを、AIMD
計算には当研究室で開発した SPPR プログラムをそれぞれ用いた。
【計算結果】まず、基底状態の 1,2-ブタジエンの安定構造を求め、垂直励起エネルギーを計
算した。1,2-ブタジエンの励起状態と基底状態を結ぶ交差領域を調べた結果、3 つの最小エネ
ルギー円錐交差(MECI)構造が見つかり、それぞれの構造の形から cis, linear, trans と命名した
（図 2 参照）。交差領域の相対エネルギーを計算し比較したが、それぞれの MECI 構造にエネ
ルギー的な差はほとんどなく、エネルギー的にどれが有利な失活経路であるかは解らなかっ
た。これに対し、S1 フランクコンドンの構造からの最急降下経路計算を行ったところ linear
の構造に到達した。linear におけ
る S1 の極小構造を図 2 に示す。
続いて、第一励起状態のポテ
ンシャル曲面を調べるために、ア
レン部分の角度と末端の CH2 の
内部回転に対応する内部座標を
それぞれ固定し、
数値を変化させ
て各点で構造最適化計算を行い
第一励起状態の 2 次元ポテンシ
ャル曲面を作成した。ポテンシャ
ル曲面から、第一励起状態のフラ
ンクコンドン領域はちょうどエ
図 2. 1,2-ブタジエンの光緩和過程における主要な構造
ネルギー的に高い位置に存在し
ており、cis, linear, trans それぞれ
の交差領域付近にはエネルギー降下経路が存在することがわかった。以上より最急降下経路
は linear へ到達するが、第一励起状態における 1,2-ブタジエンの動的な効果を検証するため
には、ダイナミクスの考慮が必要であることが示唆された。
AIMD 計算より、光励起後 3 つの安定領域それぞれに古典軌道が到達する結果が得られ
た。古典軌道の解析から分岐比、最も有利な失活経路、交差到達時間について調べたところ、
分岐比は、cis : linear : trans = 4 : 88 : 8 となり、最急降下経路計算が到達した linear 方向の失
活経路が動力学的にも最も有利なであることがわかった。さらに、交差到達時間の平均を取
ったところおよそ 150 fs 程となり実験結果を定性的に再現した。
AIMD 計算において最も有利と予想された失活経路である linear の極小構造において、励起
分子のイオン化エネルギーを計算したところ、実験値と非常に良い一致を得たため、実験的
に得られたスペクトルを励起分子のスペクトルと帰属した。さらに、AIMD 計算の各古典軌
道において分子の運動を吟味したところ、アレン部分の折れ曲がり振動と末端 CH2 の内部回
転運動が励起状態からの緩和過程で活性化されることがわかった。これら 2 つのモードの基
底状態における振動数を計算したところ、アレン部分の折れ曲がりに対応する振動数が実験
で得られたスペクトルピーク強度の振動数を非常によく再現した。そのため、実験で得られ
たスペクトルピークの振動は励起状態で活性化された分子の振動励起によるものであると結
論付けた。
1,2-ブタジエンの失活経路について cis, linear, trans の 3 つの経路があることを明らかに
した。さらに、AIMD 計算と最急降下経路計算を行うことで、有利な失活経路は linear である
ことがわかった。実験で得られたデータについて、励起直後に観測されたスペクトルピーク
は励起状態の分子によることがわかった。また、励起失活後スペクトルピーク強度が振動し
ていたことについて、励起失活後の分子の振動励起に起因していることを明らかにした。
【参考文献】
[1] H. Y. Lee, V. V. Kislov, S. H. Lin, A. M. Mebel and D. M. Neumark, Chem. Eur. J., 3, 9 (2003).
[2] A. Makida, H. Igarashi, T. Fujiwara, T. Sekikawa, Y. Harabuchi, and T. Taketsugu, J. Phys. Chem. Lett., 5,
1760 (2014).
[3] R. Iikubo, T. Fujiwara, T. Sekikawa, Y. Harabuchi, S. Satoh, T. Taketsugu, and Y. Kayanuma, J. Phys. Chem.
Lett., 6, 2463 (2015).
3P112
時間依存コーンシャム法プログラムの多電子系への拡張
（弘前大院・理工）○岡崎
功
Expansion of our time dependent kohn-sham program
to solve many-electron systems
（Hirosaki Univ.）○Isao Okazaki
【序】
我々は波動関数の時間発展が必要な物理現象を研究対象とするために、実空間をグリッド
分割して有限要素法に基づき、時間依存コーンシャム(TD-KS)方程式を解くことで電子波動関
数の時間発展を求めるプログラムを開発している。昨年度、開発中のプログラムについて１
電子系のテスト計算を報告した[1]。今回、多電子系について、スピン分極、及びスピン非分
極の時間依存コーンシャム方程式を解くことができるようにプログラムを拡張した。さらに、
初期電子波動関数として系の固有関数を求めるために必要な波動関数の非束縛最小化(Eum)
法についても多電子系に拡張したので報告する。
【計算方法】
ℓ
開発中のプログラムは、次の TD-KS 方程式によって
ℓ
の時間発展を求める。
はスピ
ン のℓ番目の空間軌道である。
ℓ
i
ℓ
,
,
，
、ハートリーポテンシャル
の項は、核引力ポテンシャル
、交換相関ポテンシャル
の和からなる。
外場ポテンシャル
ℓ
この TD-KS 方程式は次のように解いた[2]。指数積展開法を用いて∆ 秒後の
ℓ
,
∆
、
exp
∆
ℓ
Δ
,
exp
∆
∆
exp
exp
は
∆
ℓ
,
と表せる。ここで
であり( V
≡
と
に、
V
V
Δ ,
Δ
,
Δ
とした)、 Δ 秒進んだ電子密度
を利用して
ℓ
∆
≡ ∑ℓ exp
ℓ
,
である。計算の流れは次のようになる。指数積展開法からまず
exp
∆
ℓ
,
を直接利用する代わり
Δ
ℓ
,
に exp
∆
を作
、 、 方向の順に有限要素法に基づき、連立方程式を解くことで求める
ことができる[1]。
ように
と
を求める。ここで
≡ ∑ℓ
用させる。これは
Δ
ℓ
∆
exp
ℓ
,
∆ を求める途中の段階で
を作用させて、
ℓ
を得たのち、
そして
を求める。この
を得ることができる。続いて
を得る。これは単なる掛算であり、グリッド上の
ℓ
に
ℓ
の
ℓ
位相が変わるだけである。最後に
る。 は全電子密度(
に再び exp
∆
を作用させて
)からポアソン方程式を解いて求めた。
ℓ
,
∆ を得
は
と
をもとに LIBXC ライブラリ[3]を使用して求めた。
昨年度１電子系について作成した Eum 法については差分法ではあるが、スピン分極、及び
スピン非分極の多電子系に拡張した。Mauri らによりスピン非分極の系については Eum 法が
与えられている[4]。スピン分極の系についても同様に以下の目的関数
の最小化により、全
電子エネルギーを最小化する規格直交の空間軌道を求めることができる。スピン のℓ番目の
ℓ
空間軌道を
としたとき、最小化すべき
∑ ∑ ∑
,
≡2
,
≡
ℓ
は実関数とした。ここで
,
≡
はスピン の電子数、
ℓ
≡
ℓ
,
,
は交換相関エネルギー汎関数である。 はフェルミエネ
ℓ
,
,
0 による
と
による汎関数微分を求めたのち、グリッド
ℓ
の表式を書き出し、共役勾配法による最小
化のプログラムコードを拡張した。電子波動関数の時間発展のときと同様に、
方程式から求め、
,
≡∑ ∑
,
ルギーよりも大きな適当な値を与える。 の
上の値
∑
,
である。ただし Mauri らと同じく
である。
は
と
ℓ
から現れる
はポアソン
は LIBXC ライブラリを使用して求めた。以上の
プログラムコードは三次元系だけでなく、一次元系、二次元系でも実行できるようにした。
プログラムは主に C 言語によりコード化し、ユーザ入力の解析部分は flex と bison でコ
ード化している。現状で、ポアソン方程式の解法と Eum 法を除いた主要部分を OpenMP により
並列化している。
拡張したプログラムコードを検証するために、
まず仮想的に電子間相互作用が無いとした場合に
ついて、２～４電子から成るスピン分極とスピン
非分極の系についてテスト計算を行った。算出さ
れる期待値が対応する１電子系のときの和になる
ことを確認した。次に電子間相互作用を含めて水
素分子やリチウム原子などについてテスト計算を
行った。一例を図１に示す。詳細は当日報告する。
E[ψα,ψβ] (a.u.)
【検証計算】
‐1.82
‐1.84
‐1.86
‐1.88
‐1.90
‐1.92
‐1.94
‐1.96
‐1.98
‐2.00
‐2.02
0
10
20
30
40
iteration
50
60
図 1 H2 分子の Eum 法における収束過程。グリッ
ド間隔は約 0.12 a.u.である。
は PBE を用い
た。初期空間軌道は H 原子の 1s を与えた。E は
核間反発を除いた全電子エネルギーに収束する。
[1] 岡崎功, 第 8 回分子科学討論会, 3P123 (東広島,2014). [2] N.Watanabe,M.Tsukada, Phys.Rev.E
65, 036705 (2002). [3] M.A.L.Marques, M.J.T.Oliveira, T.Burnus, Comput.Phys.Commun. 183, 2272
(2012); http://www.tddft.org/programs/Libxc . [4] F.Mauri, G.Galli, R.Car, Phys.Rev. B 47, 9973
(1993).
3P113
次世代有機太陽電池に関する理論的研究
（首都大院理工1, 理研AICS2, 岐阜大地域3, 理研CEMS4）
○今村穣1, 松井亨2, 神谷宗明2,3, 菅野翔平1, 尾坂格4, 瀧宮和男4, 中嶋隆人2, 波田雅彦1
Theoretical Study on Next-generation Organic Solar Cells
(Tokyo Metropolitan Univ. 1, RIKEN AICS2, Gifu Univ.3, RIKEN CEMS4)
○Yutaka Imamura1, Toru Matsui2, Muneaki Kamiya2,3, Syohei Kanno1, Itaru Osaka4, Kazuo
Takimiya4, Takahito Nakajima2, and Masahiko Hada1
【緒言】
人類の最重要課題の一つであるエネルギー問題の解決へ向けて、これまで太陽電池の開発
が精力的に行われてきた。最近では、シリコン太陽電池や化合物太陽電池などが実用化され
ている。しかし、現在主流の無機太陽電池ではコストの面から本質的なエネルギー問題の解
決は困難と考えられている。そのため、次世代太陽電池の開発が脚光を浴びており、特に、
軽量性・コストの面から有機系太陽電池である有機薄膜太陽電池・色素増感太陽電池などが
注目されている。
最近、東京大学の瀬川グループが色素増感太陽電池材料として開発を行った DX1 色素[1]
は、長波長領域におけるスピン禁制遷移により 10％を超える高い光電変換効率を示した。
我々は、以前少ない計算コストでスピン-軌道(SO)相互作用が考慮可能な Tamm-Dancoff 近似
に基づく 2 成分相対論的時間依存密度汎関数理論(SO-TDDFT/TDA)を用いて、DX1 色素のス
ピン禁制遷移の検討[2]を行った。今回は、さらに、DX1 色素よりも効率の良い色素の構築
を目指し、異なる色素骨格および配位子の役割に関して検討した。
現在高い光電変換効率を示す有機薄膜太陽電池では、電子供与体としてポリチオフェン誘
導体がよく用いられている。今後更なる有機薄膜太陽電池材料の設計・開発には、ポリチオ
フェン誘導体の励起状態に関する詳細な情報が重要である。そこで、本研究では、様々なポ
リチオフェン誘導体の励起スペクトルを計算し、検討[3]を行った。
COOH
【新奇な増感色素の理論的探索】
HOOC
以前報告された DX1 色素は、N749 色素などのターピリジル基を
N
含む色素骨格をベースに開発された。本研究では高い光電変換効率
NCS
N
Ru
を示す N3 色素(図 1)などのビピリジル基を含む色素骨格をベースに、
NCS
N
スピン禁制遷移が起こる新奇な色素を探索した。計算手法として、
N
HOOC
SO-TDDFT を用い、交換相関汎関数として、PBE1PBE 汎関数を採用
した。
COOH
今回、様々な配位子を検討したところ、N3 色素の NCS 基を重ハ
ロゲン基に置換した場合、比較的強いスピン禁制励起が長波長領域 図 1 N3 色素の構造
に現れることがわかった。詳細を検討したところ、以前の研究で報告された Ru 元素の d 軌道
の寄与とは異なる因子により、スピン禁制遷移が現れていることがわかった。更なる詳細な
検討、その他の色素の検討に関しては当日報告する。
【ポリチオフェンの光吸収スペクトル】
ポリチオフェン誘導体の一つである TTD-オリゴチオフェン共重合体 (図 2)の励起スペクト
ルに関して時間依存密度汎関数(TDDFT)を用いて検討を行った。交換相関汎関数は PBE1PBE
である。まず、TDDFT を用いて、TTD2T, TTD4T, 2TTD4T, 3TTD4T, 4TTD4T のオリゴマーの
Me
計算を行った。図 3 に示すようにオリゴマーが伸長す
O
るにつれ、HOMO から LUMO への遷移である最低励
S
X
S
起状態のピークは大きく長波長側にシフトし、実験と
S
X
S
矛盾しない吸収スペクトルが得られた。また、実験と
O
n
同様に重合度が大きくなるにつれ、吸収ピークの強度
Me
が増加することも再現できた。
図 2 n 量 体 の FFD4T, TTD4T,
次に、アクセプター部位の硫黄を酸素(FFD)および SSD4T の構造（X = O, S, Se）
セレン(SSD)に変更した類縁体を検討した。図 4 から
わかるように、FFD の第一励起状態は、TTD と比較して青方偏移することがわかった。一方、
SSD は、変化しないことがわかった。これは、酸素により軌道エネルギーが安定化するため、
起こったと考えられる。以上の検討から、酸素置換による光物性制御が可能なことが示唆さ
れた。
Intensity (au)
Intensity (au)
Simulation Experiment
Intensity (au)
Absorption coefficient (cm‐1M‐1)
Wavelength (nm)
Wavelength (nm)
図 3 TTD2T、nTTD4T の吸収スペクトル
（上段：計算値、下段：実験値）
Wavelength (nm)
図 4 （上段）FFD2T、nFFD4T および
（下段）SSD2T、nSSD4T の吸収スペクトル
【参考文献】
[1] T. Kinoshita, J. Ting Dy, S. Uchida, T. Kubo, and H. Segawa, Nat. Photonics 7 (2013) 535.
[2] Y. Imamura, M. Kamiya, and T. Nakajima, Chem. Phys. Lett. in press.
[3] T. Matsui, Y. Imamura, I. Osaka, K. Takimiya, and T. Nakajima, J. Phys. Chem. C submitted.
3P114
三次の展開項を用いた時間依存密度汎関数強束縛法(TD-DFTB3)
（京大・福井センター） 西本佳央
Time-Dependent Density-Functional Tight-Binding Method with the ThirdOrder Expansion of Electron Density (TD-DFTB3)
(Fukui Institute for Fundamental Chemistry, Kyoto University) Yoshio Nishimoto
Prediction of photochemical properties has been increasingly common due to the
advancement of time-dependent density functional theory (TD-DFT). Even with advancement
in TD-DFT, the computational demanding is still a heavy burden for researchers who wish to
investigate large systems. An alternative method is the density-functional tight-binding (DFTB)
method. Here, we develop a formalism for the calculation of excitation energies and excited
state gradients for the self-consistent-charge DFTB with the third-order contributions of a
Taylor series of DFT energy with respect to the fluctuation of electron density, called TDDFTB3.1
In computing excitation energies, the following non-Hermitian eigenvalue problem has to
be solved,2,3
𝐀 𝐁 𝐗
1 0
𝐗
)( ) = 𝜔(
)( )
0 −1 𝐘
𝐁 𝐀 𝐘
(
where is the excitation energy and matrix elements of A and B are defined as
𝐴𝑖𝑎𝜎,𝑗𝑏𝜏 = 𝛿𝜎𝜏 𝛿𝑎𝑏 𝛿𝑖𝑗 (𝜀𝑎𝜎 − 𝜀𝑖𝜎 ) + 𝐾𝑖𝑎𝜎,𝑗𝑏𝜏
,
and
𝐵𝑖𝑎𝜎,𝑗𝑏𝜏 = 𝐾𝑖𝑎𝜎,𝑏𝑗𝜏
.
where 𝜀𝑖𝜎 is the eigenvalue of i-th molecular orbital, and the element of the coupling matrix
𝐾𝑖𝑎𝜎,𝑗𝑏𝜏 is given as the second-order derivative of the total energy with respect to density
matrix elements. For singlet-singlet excitation, the coupling matrix can be written as
1
1
𝐾𝑖𝑎𝜎,𝑏𝑗𝜏 = 𝑆𝜇𝜈 𝑆𝜅𝜆 (𝛾𝐴𝐶 + 𝛾𝐵𝐶 + 𝛾𝐴𝐷 + 𝛾𝐵𝐷 ) + 𝑆𝜇𝜈 𝑆𝜅𝜆 [(𝛤𝐴𝐶 + 𝛤𝐴𝐷 )𝛥𝑞𝐴 + (𝛤𝐵𝐶 + 𝛤𝐵𝐷 )𝛥𝑞𝐵
4
6
+ (𝛤𝐶𝐴 + 𝛤𝐶𝐵 )𝛥𝑞𝐶 + (𝛤𝐷𝐴 + 𝛤𝐷𝐵 )𝛥𝑞𝐷
+ ∑{𝛤𝐴𝐸 (𝛿𝐴𝐶 + 𝛿𝐴𝐷 ) + 𝛤𝐵𝐸 (𝛿𝐵𝐷 + 𝛿𝐵𝐷 )}𝛥𝑞𝐸 ].
𝐸
Note that only the first term is necessary in the previous TD-DFTB2.3
Nuclear gradients can be obtained by solving Z-vector equations following the standard
TD-DFT approach.4 The third-order derivative of the total energy with respect to density matrix
elements, defined as
𝑔̂𝜇𝜈𝜎,𝜅𝜆𝜏,𝜙𝜒𝜐 =
1
𝑆 𝑆 𝑆 ∑{𝛤𝐴𝐸 (𝛿𝐴𝐶 + 𝛿𝐴𝐷 ) + 𝛤𝐵𝐸 (𝛿𝐵𝐶 + 𝛿𝐵𝐷 )}(𝛿 𝜙∈𝐸 + 𝛿𝜒∈𝐸 ),
12 𝜇𝜈 𝜅𝜆 𝜙𝜒
𝐸
is necessary in TD-DFTB3, and it is missing in the previous formulation for TD-DFTB2.5 A
general equation for excitation energy gradient is finally given as
0
𝜕𝐻𝜇𝜈
𝜕𝑆𝜇𝜈
𝜕𝐾𝜇𝜈𝜎,𝜅𝜆𝜏
𝜕𝜔
1
=∑
𝑃𝜇𝜈𝜎 − ∑
𝑊𝜇𝜈𝜎 +
∑
𝛤̂
𝜕𝑅𝛼𝑥
𝜕𝑅𝛼𝑥
𝜕𝑅𝛼𝑥
2
𝜕𝑅𝛼𝑥 𝜇𝜈𝜎,𝜅𝜆𝜏
𝜇𝜈𝜎
𝜇𝜈𝜎
𝜇𝜈𝜎𝜅𝜆𝜏
where 𝑃𝜇𝜈𝜎 , 𝑊𝜇𝜈𝜎 , and 𝛤̂𝜇𝜈𝜎,𝜅𝜆𝜏 are obtained as a solution of Z-vector equations. The
presented method has been implemented in GAMESS-US, and the implementation may be
available in a future release.
Because of the limit in length, only the computational efficiency is highlighted here. For
demonstration, trans-polyacetylene (C400H402) which contained 2,002 basis functions had been
chosen. Five lowest excitation energies were calculated, and the first excitation vectors were
used in the gradient calculation. With TD-DFTB3, ground state SCF took 213.82 s, and the
computation of excitation energies and gradient took 819.92 and 45.54 s with one CPU of Xeon
E5-1620 v3 (3.50GHz). The additional computational cost for TD-DFTB3 compared to TDDFTB2 is virtually negligible. These timings are probably one or two orders of magnitude faster
than those of TD-DFT. It is therefore possible to apply TD-DFTB for medium-size molecules
routinely. In the poster, calculated absorption and fluorescence energies of cresyl violet in
explicit water molecules (~400 atoms) will be briefly discussed as well as the comparison of
the performance of TD-DFTB2 and TD-DFTB3.
References
[1] Y. Nishimoto, submitted. [2] M. E. Casida, “Time-Dependent Density Functional Response
Theory for Molecules,” in Recent Advances in Density Functional Methods, edited by D. P.
Chong (World Scientific, Singapore, 1995) Chap. 5, pp. 155-192. [3] T. A. Niehaus, S. Suhai,
F. Della Sala, P. Lugli, M. Elstner, G. Seifert, and T. Frauenheim, Phys. Rev. B 63, 085108
(2001). [4] F. Furche, and R. Ahlrichs, J. Chem. Phys. 117, 7433 (2002). [5] D. Heringer, T. A.
Niehaus, M. Wanko, and T. Frauenheim, J. Comput. Chem. 28, 2589 (2007).
3P115
量子力学的振動状態計算法に対する関数型プログラミング言語の適用
（慶大院・理工）
○菅原道彦
Application of functional programming language to quantum mechanical vibrational
state calculation
(Graduate School of Science and Technology, Keio Univ.) ○M. Sugawara
【序】近年、開発効率を意識したソフトウェア資産のモジュール化・並列プログラミ
ングへの親和性などの理由により、関数型プログラミング言語が注目されている。一
方で、科学技術計算用のプログラミング言語としては手続き型言語である FORTRAN、
C などが広く用いられているが、数値計算用ライブラリなどを含むレガシーコードの
活用など考慮すると、今後もこれらの手続き型言語が科学技術計算の主流である可能
性は高い。しかし、後発の言語と比較するとライブラリのモジュール化（＝再利用に
おける汎用性）
、並列計算への適応性などにおいて不利な点も否めない。ライブラリ管
理はコメントや、ドキュメントベースのマニュアルなどに依存する部分が多く、大規
模なプロジェクトにおける生産性において非効率的である。特に、手続き型言語にお
いて、サブルーチン、関数、メソッド内で多くの状態変数をもつライブラリ群が複雑
に関与する計算プログラムにでは、取り得る可能性がある状態すべてについてテスト
を行うのは困難であろう。また、FORTRAN、C などでもマルチコア CPU やネットワ
ーク接続されたクラスター型計算機の普及にあわせて open-MP、MPI などの技術を駆
使した並列プログラミング環境にも対応しつつあるが、並列計算を実装するためには
既存のソースコードに対する大幅な追加改変を要求されがちであり、ハードウェア進
化に対応するための保守管理に費やされる労力は少なくない。
一方で、純粋関数型言語は一切の状態変数を持たず、入力に一対一対応する出力を
得る関数のみを使用して問題を解決する。このような参照透過性を有するため、本質
的に関数自体の単体テストに重きを置くことができ、開発効率･保守管理における優位
性は高い。特に、強く静的型付けされた関数型原語である Haskell などでは、キャスト
などを基本的に許さず不正なデータ処理する関数が原理的に書けないため、ロジック
に関する大部分のバグをコンパイラによる型チェックで検出する事が可能である[1]。
また、状態変数の値が実行中に変化する「副作用」が存在しないため、同期が必要と
される MPI などと比較してもマルチコア上の並列プログラミングに対する親和性が高
いことも期待される。
本研究では、純粋関数型言語として普及している Haskell を量子力学的振動状態計算
法の開発に適用し、その応用可能性（モナド概念の有効利用、ウェブサービスとして
のインターフェース構築、並列計算への展開など）を議論する。
【方法論】一般的な振動状態の量子力学計算では
d2
− dx 2 + V ( x) Ψ ( x) = E Ψ ( x)
（１）
のような、シュレディンガー方程式を解くことになる。多くの手法では有限個数の基
底関数系 {φi ( x)}(i = 1L N ) を導入し、波動関数を
N
Ψ ( x) = ∑ ciφi ( x)
（２）
i
と展開することによって、
（１）式の微分方程式を展開係数ベクトル c ≡ ( c1 , c2 ,L cN ) と
ハミルトニアン行列（但し、 ( H )ij ≡ φi ( x)φ j ( x) dx ）からなる固有値問題
∫
H ⋅ c = Ec
（３）
に帰着される。
本研究では、基底関数の選択として、
１）グリッド基底：波動関数を各グリッド点における離散値で表現し、運動エネ
ルギー演算子を差分近似で表現する手法。
２）DVR 基底：各種 DVR 基底を定義する特性多項式の０点における離散値で波
動関数を表現する手法。
３）EFG 基底：移動最小２乗近似を用いた関数補完によって波動関数を表現する
手法。基底関数として局所性と大域性を併せ持っているため、波動関数の局
所性と運動エネルギー項における微分評価に優れる。
など種々の基底を採用し、純粋関数型言語 Haskell での実装を試みた。
この際、状態依存計算（＝State モナド）などの関数型言語特有の計算の構造を抽象
化する手法を用いて、方法論（＝基底）の変更に依存する計算アルゴリズムを表現し、
基底選択に依存しない部分である行列作成･固有値問題の数値解法部分を共通コード
で処理可能にした。インターフェース部分に関しては、ウェブサービス公開との親和
性を考えて yesod プラットフォーム[2]を採用し、結果の可視化についてはブラウザ上
での SVG 表現が可能である d3.js を用いた[3]。この際、
グラフ化ライブラリとして d3.js
より上の階層で簡便に使用できる rickshaw を利用している[4]。
【参考】
[1]Haskell : An advanced purely-functional programming language) : https://www.haskell.org/
[2] Yesod Web Framework : http://www.yesodweb.com/
[3]3d.js : Data-driven
documents : http://d3js.org/
[4] rickshaw : http://code.shutterstock.com/rickshaw/
3P116
Divide and Conquer type Density Functional based Tight
Binding Molecular Dynamics (DC-DFTB-MD) Simulation of
Proton Transfer in Bulk Water System
(Department of Chemistry and Biochemistry, Waseda University1, Research
Institute for Science and Engineering, Waseda University2, Institute for
Molecular Science3, JST-CREST4, ESICB, Kyoto University5) Aditya Wibawa
Sakti1, Yoshifumi Nishimura2,3, Hiromi Nakai1,2,4,5
[Introduction] Diffusion of proton in bulk water has attracted long-term attention because it includes
not only a simple diffusion process but also a proton transfer via Grotthuss mechanism.1,2,3
Experimentally, the Grotthuss shuttling behavior of a proton transfer has been determined by
measurement of self-diffusion coefficient of oxygen in hydronium ion.2,4,5,6 The behavior of selfdiffusion coefficient is usually following the Arrhenius theory while the higher activation barrier of
diffusion is observed at low temperature.7 In this work, we theoretically study the molecular aspects of
the proton transfer at different temperatures.
[Computational Details] We performed molecular dynamic simulations by using divide-and-conquer
type density functional based tight-binding method (DC-DFTB-MD).8,9,10 The system contains 523
water molecules and one excess proton in 25.924 × 26.230 × 25.849 Å periodic box. Firstly, the
system was equilibrated by using Nose-Hoover NVT ensemble for 10 ps at 280, 300, 320, and 360 K.
Next, 15 ps NVE ensemble simulations was performed for production run at each temperature.
In the present research, the diffusion coefficient of proton (Dp) is obtained as the summation
of physical diffusion coefficient of water (Dw) and Grotthuss diffusion coefficient (DG)
Dw = lim
()
()
rO t − rO 0
t→∞
6t
2
, DG =
2
l
r ,
6 p
Dp = Dw + DG = lim
t→∞
()
()
rO t − rO 0
6t
2
+
l2
r
6 p
(1)
where 𝑙 is the proton hopping length set to 2.5 Å1 and rO ( t ) is the position vector of oxygen atom at
time t. The proton transfer rate rp was determined from the slope of the following accumulation
function
h (δ t ) = h (δ t − 1) + δ h (δ t )
(2)
where h ( 0 ) = 0 and δ h (δ t ) is 1, 0, and -1 for forward, no, and backward shuttling, respectively. The
activation barrier of diffusion ΔE was evaluated from the logarithm form of the Arrhenius equation,
ΔE
ln D = ln D0 −
(3)
RT
where D is one of Dw, DG, or Dp, D0 is pre-exponential factor, R is gas constant, and ΔE is one of ΔEw,
ΔEGi , or ΔEpi , which describe activation barriers of water, Grotthuss, and proton diffusions,
respectively.
[Results and Discussion] Figure 1 shows the time evolution of Eq. (2). The short and long vertical
lines indicate the oscillation shuttling and forward jumping behavior in proton transfer, respectively.
The obtained rp of 0.69 ps-1 at 300 K is in good agreement with the experimental value (rp= 0.67 ps1 4
). The low proton transfer rate of 0.11 ps-1 at 280 K indicates higher activation barrier of proton
transfer due to the frequent forward and backward shuttling. This behavior is caused by very stable
special pair interaction between the proton and water molecule, leading to the formation of Eigen
complex.
Table 1 shows the calculated diffusion coefficients at 300 K together with previous
experimental and theoretical values. This work agrees well the experimental value of Dp. The larger
system size and the use of Mulliken charges to identify the hydronium ion seems to improve Dw and
Dp from the previous work.
Figure 2 shows the Arrhenius plot of three different diffusion coefficients. The diffusion of
water shows a good linear relationship (black) while Grotthuss (red) and proton diffusion (blue) gives
a significant deviation at 280 K. It indicates that the mechanism of proton transfer via Grotthuss
shuttling is totally different from physical diffusion of water molecule.
The calculated diffusion barriers of water ( ΔEw ), Grotthuss ( ΔEG1 , ΔEG2 ), and proton ( ΔEp1 , ΔEp2
) are tabulated in Table 2. The obtained value of ΔEp1 (7.76 kJ/mol) is in agreement with the
experimental one (10.04 kJ/mol).4 The stabilization by hydrogen bonding in Eigen complex is
supposed to lead to the higher values of ΔEG2 and ΔEp2 than ΔEG1 and ΔEp1 , consistent with the
1.0
∆ E1p
●
−0.5 0.0
15
Water Diffusion
Grotthuss Diffusion
Proton Diffusion
●
●
●
ln(D)
∆ E1G
∆ E2G
−1.5
10
5
h(δ t)
rp = 0.1086 ps−1
rp = 0.6871 ps−1
rp = 0.8750 ps−1
rp = 1.1370 ps−1
280 K
300 K
320 K
360 K
0.5
20
observation in experiment.7
∆ E2p
∆ Ew
0
5
10
15
−2.5
0
●
Time(ps)
2.6
2.8
3.0
3.2
3.4 3.6
−1
1000 T (K
)
3.8
Figure 1. Time course changes of accumulation
Figure 2. Arrhenius plot for diffusion of water,
function ℎ(𝛿𝑡).
Grotthuss, and proton.
Table 1. Number of proton and water molecules 𝑁, simulation time 𝑡, temperature fluctuation in NVE ensemble
simulations ΔT , and calculated diffusion coefficients among different theoretical methods at 300 K.
DC-DFTB3
CPMDa
DFTB2a
DFTB3b
Exp.
(present)
+
N(H )/N(H2O)
1/128
1/128
1/128
1/523
t(NVT)/t(NVE) (ps)
8/35
60/80
100/82
10/15
ΔT (K)
± 100
± 100
± 100
± 15
𝐷! (Å2/ps)
0.10
0.65
0.38
0.19
0.23 c
2
DG (Å /ps)
0.72
0.70 d
2
Dp (Å /ps)
0.33
0.90
0.66
0.91
0.94 e
Dp /Dw
3.30
1.38
1.74
4.77
4.09
a
Car-Parrinello MD (CPMD); second order self-consistent charge DFTB (DFTB2) in the Reference 11, b DFTB
with third-order diagonal term correction (DFTB3) in the Reference 12. cReference 5, d Reference 2, e Reference
6.
Table 2. Activation barrier of diffusion (in kJ/mol) estimated from Eq. (3).
T = 300 – 360 K
T = 280 – 300 K
ΔEw
ΔEG1
ΔEp1
ΔEG2
ΔEp2
Calc.
8.88
7.41
7.76
64.42
41.52
a Exp.
10.04
a
Reference 4.
[References]
[1] Agmon, N. Chem. Phys. Lett., 1995, 244, 456-462. [2] Meiboom, S. J. Chem. Phys., 1961, 34, 375-388. [3]
Chen, H., Voth, G.A., Agmon, N. J. Phys. Chem. B, 2010, 114, 333-339. [4] Luz, Z., Meiboom, S. J. Am. Chem.
Soc., 1964, 86, 4768-4769. [5] Mills, R. J. Phys. Chem., 1973, 77, 685-688. [6] Roberts, N.K., Northey, H.L. J.
Chem. Soc. Faraday Trans., 1974, 70, 253-262. [7] Agmon, N. J. Phys. Chem., 1996, 100, 1072-1080. [8]
Kobayashi, M., Nakai, H. in Linear-Scaling Techniques in Computational Chemistry and Physics, Springer,
2011, 97-127. [9] Elstner, M. J. Phys. Chem. A., 2007, 111, 5614-5621. [10] Yang, Y., Yu, H., York, D., Cui,
Q., Elstner, M. J. Phys. Chem. A., 2007, 111, 10861-10873. [11] Maupin, C., Aradi, B., Voth, G.A. J. Phys.
Chem. B., 2010, 114, 6922-6931. [12] Goyal, P., Elstner, M., Cui, Q. J. Phys. Chem. B., 2011, 115, 6790-6805.
3P117
有機薄膜太陽電池に関する理論的研究
(理研・計算科学研究機構 1,首都大 2 ) ○田代 基慶 1，今村 穣 2，河東田 道夫 1, 中嶋 隆人 1
Theoretical Study on Organic Photovoltaics
(RIKEN AICS1, Tokyo Metro. Univ.2)
○Motomichi TASHIRO1, Yutaka Imamura2,
Michio Katouda1, Takahito Nakajima1
有機薄膜太陽電池は従来の太陽電池に比べ比較的製造コストが掛からず、軽量・折り曲
げ易い、複雑な形状の物体へも塗布的手法で貼り付けることが可能などの特徴を持っている
ため、次世代太陽電池の１つとして注目を集めている。一方で、シリコン系太陽電池と比べ
ると変換効率・耐久性などの面で改善の余地があるため、製品レベルでの開発・素過程に関
する基礎的研究などが多くの企業・研究機関で盛んに行われている。
多くの場合、有機薄膜太陽電池は 2 つの電極の間に電子供給体（ドナー）・電子受容体
（アクセプター）からなる発電層(有機分子薄膜)が挟まれた構造を持ち、発電層はバルクヘテ
ロジャンクションと呼ばれるドナー・アクセプター領域が入り混じった特徴的な構造を持っ
ている。また、ドナー材料としてはπ共役高分子が、アクセプター材料としてはフラーレン
誘導体が良く用いられている。光が照射された際は、 1. ドナー領域での励起子(電子・正孔
対)生成・拡散、 2. ドナー・アクセプター領域境界における励起子の電子・正孔への分離、 3.
分離した電子・正孔それぞれがドナー・アクセプター領域を通って電極に到達する、という
仕組みで電荷生成・電荷取り出しが実現していると考えられているが、その詳細な機構につ
いては議論が続いている[1]。
我々は有機薄膜太陽電池での電荷分離・移動等の微視的機構の理解を目的として P3HT・
PCPDTBT(下図)などのπ共役高分子からなるドナー、PCBM(下図)などのフラーレン誘導体か
らなるアクセプターを対象とした研究を行っているが、特に、励起子が正孔と電子に分離す
る際に生成される電荷移動状態や関与する分子軌道の性質などに着目している。電荷分離過
程においてはアクセプター領域での電子の非局在性を考慮することが重要だとする議論があ
るものの、対象とする分子が大きいこともあって系統的な計算の例はそれほど多くはない。
今回、分子動力学計算で取得した構造をもとにしてアクセプター領域におけるフラーレン誘
導体同士の電子カップリングを求め、軌道の非局在化の程度・軌道エネルギーの変化等を調
査した。電子カップリングに関しては、Fock 行列から求める方法と軌道の重なりから推定す
る方法の２通りについて比較している。分子としては PCBM, bis-PCBM, PC70BM を考慮した。
本発表では、これらの電子カップリングを計算する手法の比較検討・分子種・構造などが結
果に与える影響などを中心に報告する。以上に加えて、ドナー・アクセプター界面における
励起子から電荷移動状態への遷移・再結合に関しても議論する予定である。
PCBM
P3HT
PCBM LUMO-LUMO 電子カップリングと軌道の重なり（左）
MD 配置に対して得られた軌道の重なりの頻度分布(右).
PCBM LUMO~LUMO+2 軌道エネルギーの変化（左）、各状態(軌道)の非局在性の程度(右).
両方とも MD 計算(216 PCBM)で得られた構造に対する結果.
[1] Few, Frost and Nelson, Phys.Chem.Chem.Phys. 17 2133 (2015).
3P118
超並列量子化学計算プログラム SMASH
(分子研) ○石村 和也
SMASH: Massively parallel program for quantum chemistry
calculations
(IMS) ○Kazuya Ishimura
【序】スーパーコンピュータだけではなく研究室レベルの計算機クラスターでも並列計算は日常
的に行われるようになった。1 ノード当たりの CPU コア数は増え続けハードウェアがより複雑に
なる中で、並列性能及び実行性能の高いプログラムの必要性は高まっている。そこで、2012 年よ
りオープンソースライセンスの大規模並列量子化学計算プログラム SMASH の開発を始め、2014
年 9 月にバージョン 1.0 を公開した[1]。現時点では、PC から京コンピュータまでのスカラ型 CPU
を搭載した計算機を対象としている。
【方法】これまでに開発した 2 電子積分計算[2]、MPI/OpenMP ハイブリッド並列化 Hartree-Fock
エネルギー計算[3]アルゴリズムなどを基に、Hartree-Fock、DFT エネルギー及び構造最適化プログ
ラムを実装した。
例として、Hartree-Fock 計算における 2 電子積分計算の分散アルゴリズムを図 1 に示す。4 重ル
ープの最外ループを OpenMP でノード内並列し、第 3 ループを MPI でノード間並列している。
OpenMP の動的な負荷分散の導入と、IF 文を使わない MPI の負荷分散により、均等な計算負荷分
散を図りつつ、分散のためのコストを大幅に削減した。4 重ループ終了後に、ノード内は reduction
により、ノード間は mpi_allreduce により Fock 行列要素を集めている。行列対角化は LAPACK ラ
イブラリを利用しているため、ノード内のみ並列化されている。2 電子積分ルーチンは、表 1 に示
すように基底関数などデータのやり取りを引数のみで行い、ライブラリ化している。そのため、
どの計算方法でも 2 電子積分計算の呼び出しが容易であり、このルーチンのみ抜き出し他のプロ
グラムに組み込むことも可能である。
エネルギー微分計算においても、最もコストの高い 2 電子積分の微分計算はエネルギー計算と
同様の MPI/OpenMP ハイブリッド並列アルゴリズムを採用している。構造最適化計算では、結合、
角度、二面角の情報を使う redundant 座標を導入し、その力場パラメータから初期 Hessian を作成
している。
【結果】京コンピュータで性能を測定したところ、(C150H30)2 (cc-pVDZ, 4500 基底)の B3LYP エネ
ルギー計算は 98304CPU コアで 154 秒、並列加速率は 5 万倍以上であった(図 2)。コア数が増加し
たときの並列化効率低下の主な原因はノード内のみ並列化されている行列対角化計算で、どのコ
ア数でも 35 秒かかっている。今後ノード間も並列化されているライブラリの導入で、計算全体の
並列化効率向上を図る予定である。
図 3 の(C150H30) (cc-pVDZ, 2250 基底)の B3LYP エネルギー微分計算については、行列対角化計算
を含まないことから、並列化効率はほぼ 100％であった。構造最適化計算では、最適化サイクル数
は Cartesian 座標に比べて 1/5 から 1/6 程度と大幅に削減できた(表 2)。京コンピュータなどのスー
パーコンピュータでナノサイズ分子・クラスターの構造最適化計算が容易に行えるようになった。
表 2 B3LYP/cc-pVDZ 構造最適化回数(初期構造
HF/STO-3G)
Cartesian
Redundant
Luciferin(C11H8N2O3S2)
63
11
Taxol (C47H51NO14)
203
40
60000
並列加速率
!$OMP parallel do schedule(dynamic,1) &
!$OMP reduction(+:Fock)
do =nbasis, 1, -1
do =1,
=(+1)/2+
start=mod(+mpi_rank,nproc)+1
do =start, ,nproc
do =1,
AO2 電子積分()計算
+Fock 行列に足し込み
enddo
enddo
enddo
enddo
!$OMP end parallel do
call mpi_allreduce(Fock)
40000
20000
図 1 MPI/OpenMP ハイブリッド並列化 Fock 行
0
0
列計算アルゴリズム
24576 49152 73728 98304
CPUコア数
図 2 B3LYP エネルギー計算の並列性能
表 1 2 電子積分ルーチン
subroutine int2elec(twoeri, exijkl, coijkl, xyzijkl,
nprimijkl, nangijkl, nbfijkl, maxdim, mxprsh,
16384.0
threshex)
2 電子積分値 (Output)
exijkl
primitive 基底関数の指数 (Input)
coijkl
primitive 基底関数の係数
xyzijkl
xyz 座標
nprimijkl
primitive 基底関数の数
nangijkl
軌道角運動量(s=0, p=1, d=2,...)
nbfijkl
基底関数の数(s=1, p=3, d=5or6,...)
maxdim
最大 twoeri の次元数
mxprsh
最大 primitive 基底関数の数
threshex
exp(-x)計算の x の閾値
12288.0
並列加速率
twoeri
8192.0
4096.0
0.0
0
4096
8192 12288 16384
CPUコア数
図 3 B3LYP エネルギー微分計算の並列性能
[1] http://smash-qc.sourceforge.net/
[2] K. Ishimura, S. Nagase, Theoret. Chem. Acc., 120, 185 (2008).
[3] K. Ishimura, K. Kuramoto, Y. Ikuta, S. Hyodo, J. Chem. Theory Comput., 6, 1075 (2010).
3P119
鎖往復重合の機構解明に向けた高分子鎖交換反応の理論的研究
(名大院・情報科学)
○松本 健太郎, K. S. Sandhya, 高柳 昌芳, 古賀 伸明, 長岡 正隆
Theoretical study of polymer chain exchange reaction toward
understanding chain-shuttling polymerization
(Graduate School of Information Science, Nagoya University)
○K. Matsumoto, K. S. Sandhya, M. Takayanagi, N. Koga, M. Nagaoka
【序論】 エチレンと 1-オクテンの共重合体は、用いる金属触媒によって共重合比を調整すること
ができる。1-オクテンが比較的少ない場合には結晶性・高密度・高融点の硬質高分子が、比較的多
い場合には非結晶性・低密度・低ガラス転移点の軟質高分子が得られる。硬質高分子と軟質高分
子を交互に連結したオレフィンブロック共重合体（olefin block copolymer、OBC） は、高融点・低
密度・低ガラス転移点の実用上有用な高分子を与える可能性がある。
OBC の新たな合成法として報告された鎖往復重合（chain shuttling polymerization, CSP） [1]
は、それぞれ硬質、軟質高分子の重合反応を触媒する Zr 触媒、Hf 触媒に加え、ZnEt2 を含んだ
反応系である。それぞれの触媒上で重合する高分子鎖が ZnEt2 を介して交換されることにより、
OBC が生成される。CSP が起こる反応系の報告は複数あるが[2]、どのような因子が CSP にとっ
て重要であるかは不明であり、現時点で反応系の設計方針も提案されていない。
本研究の最終目標は、分子シミュレーションを用いて CSP の
機構を分子レベルで解明し、反応に重要な因子及びより効率的
な反応系の設計のための指針を提案することである。今回はそ
のための第一段階として、Hf 触媒及び ZnEt2（図 1）の分子力
場の開発と、それらを用いた分子動力学（MD）シミュレーショ
図 1: Hf 触媒の活性種と ZnEt2
ンの結果について報告する。
【方法】
Hf 触媒と ZnEt2 の分子力場は AMBER 12 の GAFF（general AMBER force field）を元
に、密度汎関数法（DFT）計算の結果に基づいた GAFF 既存のパラメータの調整と、GAFF に含ま
れないパラメータの新規開発を行った。DFT 計算の汎関数には M06 を用い、基底関数は H、C、
N、O には 6-31G(d, p)、Zn には有効内殻ポテンシャル基底関数 LanL2DZ、Hf には LanL2DZ に加
えて f 軌道を用いた。原子電荷は Merz-Kollman 法により決定した。分子動力学（MD）計算には
AMBER12 を用いた。その対象系は Hf 触媒 1 分子、ZnEt2 3 分子、エチレン 60 分子、及び溶媒ヘ
プタン 140 分子からなり、周期境界条件の下で任意に配置した。100 MPa、0.5 ns の MD 計算で圧
縮した後、1 ns の MD 計算で実験条件である 400 K、1.4 MPa に緩和させた。その後 10 ns の平衡
MD 計算を実行して解析に用いた。
【結果と考察】
■分子力場の妥当性 開発力場の妥当性を確認するため、真空中における 400K での MD 計算で
1000 点の構造をサンプルし、Hf 触媒と ZnEt2 に対して DFT 計算及び開発力場による全ポテンシ
ャルエネルギーを計算し比較した（図 2）。相関係数は ZnEt2 が 0.91、Hf 触媒が 0.74 であり、構造
と全ポテンシャルエネルギーの相関という意味において、開発
力場は DFT 計算の結果をよく再現していることが分かった。
■エチレンの Hf 触媒への接近
Hf 触媒の活性中心である Hf
原子は、その周辺の置換基によって作られたポケット状構造の
内部に位置しているため、重合反応はこのポケット内部で進行
すると考えられる。真空中での QM 計算により、図 3 に示すよ
うに、エチレンが Hf 触媒のポケット状構造に接近した安定複
合体と、エチレンがポケット内で Hf-メチル基間に挿入して重
合が進行する遷移状態[3]が見つかっている。非極性のヘプタン
溶液中においても同様の安定複合体（I）及び遷移状態が存在す
図 2: Hf 触媒（赤）と ZnEt2（青）における
QM および MM での全ポテンシャルエネルギー
の相関関係
ると予想され、反応はこれらの状態を経由して進行すると予測される。
図 3: Hf 触媒上でのエチレンの重合反応
一方、開発した力場を用いた MD 計算の結果、エチレンがポ
ケット状構造内部の Hf 原子に接近する様子が複数回確認でき
た。図 4 に Hf 原子に接近したエチレン分子の炭素原子の中点
と Hf 原子の間の距離の時間変化の一例を示す。50 ps 付近でエ
チレンがポケット内の Hf 原子に接近し、10 ps 程度の間、5Å
程度の近傍に留まって運動した。QM 計算によると、安定複合
体における Hf 原子-エチレン炭素間の距離は約 3Åであるが、
ヘプタンによるエチレンの溶媒和によって安定複合体におけ
る Hf 原子-エチレン炭素間距離は伸びると予想される。よっ
図 4: エチレン炭素原子の中点と Hf 原子
間の距離の時間変化
て、この結果は先の予測と整合しており、開発した力場の妥当性を示していると考えられる。当
日は、MD 計算の結果のより詳細な解析と考察の結果について報告する予定である。
図 5: Hf 原子（紫）にエチレンが最接近した様子のスナップショットのステレオビュー
（図 4 の 60 ps 付近）灰色はヘプタン、緑は ZnEt2 を表す。
【参考文献】
[1] D. J. Arriola, E. M. Carnahan, P. D. Hustad, R. L. Kuhlman, T. T. Wenzel, Science, 312, 714 (2006)
[2] A. Valente, A. Mortreux, M. Visseaux, and P. Zinck, Chem. Rev., 113, 3836 (2013)
[3] K. S. Sandhya, N. Koga, M. Nagaoka, 第 9 回分子科学討論会, 2015, 東京, 1P099
3P120
局所応答分散力(LRD)法に基づく分極型力場の開発
(早大先進理工 1，早大理工研 2，JST-CREST3，京大 ESICB4)
〇若山 和史 1，五十幡 康弘 1，中井 浩巳 1-4
Development of polarizable force field based on local response dispersion (LRD) method
(Advanced Science and Engineering, Waseda Univ.1, RISE, Waseda Univ.2, CREST, JST Agency3, ESICB, Kyoto Univ.4)
〇Kazufumi Wakayama1, Yasuhiro Ikabata1, Hiromi Nakai1-4
【緒言】
古典分子力場計算・分子動力学シミュレーションでは，通常，分子間相互作用を Lennard-Jones
(LJ)型ポテンシャルと Coulomb 相互作用の和などで表す非分極型のポテンシャルが用いられる。
分極エネルギーの第一近似として，原子間の電荷-誘起双極子相互作用の和として表すことができ
る。当研究室で開発した局所応答分散力(LRD)法
1,2
は，密度汎関数理論における分散力の記述を
補正する手法であるが，導出過程で系の電子状態に依存した原子分極率が得られる。本発表では，
LRD 法で決定した原子分極率を用いて分極型の分子間ポテンシャルの開発を目指した。
【理論】
a
LRD 法において，分散力補正エネルギーを計算する過程で有効原子分極率 eff
が得られる。
effa drwa (r )
(r )
(r) 2
2
0
(1)
ここで， w a は空間分割関数， (r) は電子密度である。空間積分を数値的に安定させるため，極小
周波数 (= 0.01 a.u.) を用いる。 0 (r ) は連続誘電体の分散関係の長波長極限であり，式(2)のよう
に表される。3
1
3
0 (r) k F2 (1 λs 2 ) 2
(2)
kF (3 2 )1 / 3 ， s (2k F ) であり，λ は希ガスに対する分散力係数の実験値を再現するよう決
a
定された経験的パラメータである。式(1)の有効原子分極率 eff
を用いると，電荷-誘起双極子相互
作用 Epol は次式で表される。4
Epol
1 N a 2
eff Fab
2 a ,b
(3)
N は系の原子数である。Fab は原子 a に対して原子 b から生じる静電場の大きさであり，次式 4 で
表される。
Fab
qb
4π 0 rab
2
(4)
qb は原子 b の電荷，ε0 は真空の誘電率，rab は原子 ab 間の距離である。
(3)式で与えられる電荷-誘起双極子相互作用にもとづく分極エネルギーEpol は，任意の分子間ポ
テンシャルとの組み合わせが可能である。本研究では，Universal Force Field (UFF)5 に Epol を組み
合わせることとした。
【結果と考察】
2 つのテストセット 6（アニオン-中性分子系の AHB21，カチオン-中性分子系の CHB6）に対し
て，UFF によるエネルギー最小化計算から得られた構造で相互作用エネルギーを計算した。ただ
し，AHB21 セットにて共有結合性を有する[F…H…F]–，[Cl…H…Cl]–，[HO…H…OH]–は計算対象
から除外した。分極エネルギーの式(3)には，LC-BOP 汎関数と 6-311G(d,p)基底関数による DFT 計
算から算出した Mulliken 電荷および LRD 法による原子分極率を用いた。
古典力場計算に分極項を加えた相
互作用エネルギーを参照値と比較し
た結果を Figure 1 に示す。UFF では
どの系も相互作用を過小評価する。
UFF に分極エネルギーを加えること
で相互作用エネルギーが参照値に近
づく。参照値からの UFF の平均誤差，
平均絶対誤差はともに 10.82 kcal/
mol であるのに対し，UFF に分極エ
ネルギーを加えた場合は平均誤差が
4.54 kcal/mol，平均絶対誤差が 7.22
kcal/mol と改善された。
Figure 1. Comparison of interaction energies obtained by the
UFF (+ polarization) with reference CCSD(T)/CBS values.
また，分極エネルギーを計算する際に用いた電
荷，原子分極率のうち，三原子系までの数値を
Table 1. Mulliken atomic charges and LRD
Table 1 に示す。同じ電荷をもつ系で比較すると， polarizabilities.
周期表で下に位置する元素ほど原子分極率は大き
い。また，アニオン状態での分極率は分子中での
値に比べて 2 倍程度大きくなった。
【参考文献】
[1] T. Sato and H. Nakai, J. Chem. Phys. 131, 224104
(2009).
[2] T. Sato and H. Nakai, J. Chem. Phys. 133, 194101
(2010).
[3] O. A. Vydrov and T. Van Voorhis, J. Chem. Phys.
130, 104105 (2009).
[4] A. J. Stone, THE THEORY OF
INTERMOLECULAR FORCES, Oxford Univ.
Press (1996).
[5] C. J. Casewit, K. S. Colwell, and A. K. Rappé, J.
Am. Chem. Soc. 114, 10024 (1992).
[6] K. U. Lao, R. Schäffer, G. Jansen, and J. M.
Herbert, J. Chem. Theory Comput. 11, 2473
(2015).
© Copyright 2024 ExpyDoc
