全電子混合基底法プログラム
TOMBO を活用した材料科学
佐 原 亮 二 小 野 頌 太1) 大野かおる2)

.

は じめに

全電子混合基底法プログラム TOMBO ( TOhoku Mixed Basis Orbitals abinitio program)は，「全電子混合基底法」に基づいた第一原理計算プログラムである． TOMBO は，東北大学金属材料研究所川添研究室において，川添良幸・現東北大学名誉教授と，当時助教授であった大野かおる・現横浜国立大学教授が中心となり，理論の構築からコード化までを一貫して開発している，数少ない純国産の第一原理計算プ
ログラムのひとつである(1)．
混合基底法とは，原子に局在した原子軌道関数( atomicorbital, AO ) と 空 間 的 に 一 様 に 広 が る 平 面 波 ( plane wave,PW )の重ね合わせで一電子波動関数を記述する手法である

図1(式( 1 ))．本手法を全電子に拡張することで，内殻電子かTOMBO における一電子波動関数の記述．

ら価電子までの全電子状態を記述することが可能である．図に，本手法による一電子波動関数の記述を示す．

C l(r)＝∑cGlG

eiG･rV求めている．

l＋∑∑cJnlm
FJnlm(r－Rj)




j nlm

l＝∑cGl|G〉＋∑∑cJnlm|FJnlm〉
G

シュを用いて球ベッセル関数を含む数値積分で精度よく

(1)

j nlm

ここで右辺第一項は平面波，第二項は原子軌道関数である．本手法には以下のような特徴がある．






互いに重なり合わない原子球の中でのみ値をもつ原子軌
道関数を使用するため，隣接原子間の原子軌道関数同士の重なり行列の評価が必要なく，基底関数重ね合わせ誤差がない．また，オーバーコンプリートネスの問題も起こりにくい．

基底関数の数が大幅に減少する．例えば炭素系では，約本稿では紙面の都合上， TOMBO を用いたこれまでの数

7 Ry の平面波のカットオフ・エネルギーで分子動力学多くの材料研究の中から，時間依存密度汎関数理論( Timeをおこなえる．これに対して，従来の擬ポテンシャル法

dependent density functional theory, TDDFT)に基づいた，では40 Ry 程度必要であった．

新規エネルギー材料の理論設計に関する電子励起ダイナミッ

芯 電 子 軌 道 の 原 子 軌 道 関 数 に は 原 子 の Herman ク ス シ ミ ュ レ ー シ ョ ン 研 究 を 2 件 紹 介 す る(2)(3) ． ま た

Skillman コードなどで正確に計算された数値原子軌道

TOMBO がより広い分野に幅広く応用されるようにするたを使用するので，芯電子状態を正確に取り扱うことがでめ，現存する幾つかのバージョンを統一し，使用法の簡略化きる．

を進めている．ここでは，その新バージョンに小野らにより原子軌道関数と平面波の重なりについては，角度方向に追加された高速化手法について紹介する(4)．

ついては解析的に評価し，動径方向については対数メッ
独 物質・材料研究機構 元素戦略材料センター構造材料ユニット組織設計グループ主幹研究員(〒3050047 つくば市千現 121)
横浜国立大学大学院工学研究院知的構造の創生部門 1)研究員 2)教授
Materials Science Using All Electron Mixed Basis Program TOMBO; Ryoji Sahara, Shota Onoand Kaoru Ohno(National Institutefor Materials Science, Research Center for Strategic Materials, Tsukuba.Department of Physics, Graduate School of Engineering,
Yokohama National University, Yokohama)
Keywords: first principles calculation, all electron mixed basis approach, time dependent density functional theory, organic thin film solar cell,hydrogen storage materials, Chebyshev polynomial
2014年 5 月31日受理[doi:10.2320/materia.53.400]

ミ ニ 特 集.

Zn フタロシアニンと C60の電荷分離機構

近年，有機薄膜太陽電池が注目を集めている．これは，光を吸収する有機分子を電子給与体とし， C60を電子受容体とするものである．電子給与体にはフタロシアニンやポルフィニンなど，光捕集性デンドリマーのコア分子としての役割を果たす分子が例として挙げられる(5)(7)．ここでデンドリマ
ーとは，樹枝状の構造を持つナノメータースケールの分子である．中でも，光捕集性デンドリマーは，デンドリマー側鎖で光エネルギーを吸収して，そのエネルギーがコアに移動するという光捕集機能を有している(8)(9)．
有機薄膜太陽電池の研究において，エネルギー変換効率の向上は大きな目標のひとつである．そのためには，電子給与体の光の吸収効率がひとつの課題となるが，電子給与体に側鎖を付加して光捕集性デンドリマーとすることで光吸収バンドを広げることが可能となり，エネルギー変換向上が期待される．
本稿では， TOMBO を用いた p 共役デンドリマーの電子ダイナミックス計算(2)(10) のうち，オリゴフェニレンビニレン(OPV1)を修飾した Zn フタロシアニン(ZnPc)の光捕集機能と C60への電荷分離機構の TDDFT 電子励起ダイナミックス・シミュレーションについて紹介する(2)．

図2( a )は ZnPc と C60 から構成される系．( b )基底状態における HOMO LUMO 近傍のエネルギー準位と各準位の波動関数．( c )ZnPc LUMO の電子波束に占める各固有状態の割合の時間変化．( d )各準位の波動関数． Zn フタロシアニン(ZnPc )のエキシトンは C60 に電子を分離し，電荷分離が起こることが分かる．

なお，本研究では，電子系は時間依存密度汎関数理論により扱った．一方，原子核系については Newton 方程式により扱った．時間依存 KohnSham(TDKS)方程式は&i Cj (r,t)＝H(t)Cj(r,t)&t

り高いエネルギーレベルに移動する．
(2)

で与えられる．ここで H ( t)は電子系のハミルトニアンであ

次に ZnPc の HOMO から LUMO へ電子を選択的に励起さ せ た 状 態 か ら TDDFT 計 算 を お こ な っ た ． 図 2 ( c ) に
ZnPc LUMO の電子波束に占める各固有状態の割合の時間

り， Cj ( r,t )は j 番目の準位の電子波束である．電子波束 Cj変化を示す．時間発展により ZnPc と C60 の重ね合わせにな(r,t)はっていることが分かる．また，図 2(d )に各準位の波動関数
H(t)qk(r,t)＝ek(t)qk(r,t)(3)を示す．これらより， ZnPc の電子を選択的に励起させることで，ZnPc から C60への電子の移動が起きているといえる．

を満たす固有状態によりCk(r,t)＝∑cjk(t)qk(r,t)(4)

kZnPc を StarShaped Stilbenoid phthalocyanine (SSS1Pc)などの p 共役デンドリマーに置き換えることにより，よりのように展開することができる．ハミルトニアンがほとんど効率的な光捕集型の有機薄膜太陽電池を構成できるものと考変化しないと見なされる適当な時間幅 Dt において， TDKSえられる．

方程式をCk(r,t＋Dt)＝∑exp [－iek(t)Dt] cjk(t)qk(r,t)(5)

. Ni 触媒超微細粒子による水素分子解離反応

kのように積分することができる．

地球環境およびエネルギー供給に関する問題解決の手段と図( a )は ZnPc と C60 から構成される系を示す．二分子して，水素は理想的な二次エネルギー源として注目されてい間距離は 0.3 nm である．また図 2(b)は，基底状態におけるる(11) ．実用化可能な貯蔵性能( 6 mass以上(米国エネルギHOMOLUMO近傍のエネルギー準位と各準位の波動関数でー省目標値))達成を目的として，スピルオーバーを利用したあり，各準位がどの分子に属しているかを示している．この

水素貯蔵量向上法が注目されている(12)(14) ．ここでスピル

結果より，ZnPc と C60とのエネルギー準位が半導体の pn 接オーバーとは，水素貯蔵材料に水素解離能を持つ金属種を担合と同様の形になっている事が分かる．つまり，側鎖のエネ持することで，水素分子を原子状水素に解離し，材料中へ拡ルギーギャップがコアのエネルギーギャップより大きく，側散させる方法のことである．従来の第一原理計算による水素鎖の電子を励起させることで，電子がコアのより低いエネルのスピルオーバー機構の研究は，解離した原子状水素の最終ギーレベルに移動し，一方，側鎖に生じたホールはコアのよ的な安定吸着位置，原子状水素の拡散パスに対するものが多ま て り あ 第53巻 第 9 号(2014)
Materia Japanく，電子の基底状態による静的な研究に限られていなり，時間と共に水素分子がニッケルダイマー上へ近づきなる(15)(16)．一方，本機構のポイントは担持金属クラスター上がら原子状水素へと解離してゆく過程が示される．

における，水素分子が原子状水素へ解離するという初期過程

図 3(d)に，19番目から33番目の各エネルギーレベルの時である．そこで電子の励起状態を考慮した第一原理計算によ間発展を示す．実線が被占有軌道，点線が空軌道を示していり，スピルオーバー機構の初期段階である水素分子の解離過る．シミュレーション開始後，約 10 fs で 30番目と 31番目の程のダイナミックスを追い，水素分子の解離に対する金属ナ軌道が交差し，電子が高い軌道(originally LUMO)から低いノ粒子の効果について調査した研究を紹介する．計算手法は軌道(originally HOMO)へと移行することが分かる．この際前節と同様である．に水素分子が解離し原子状水素となる．また， 19 番目と 20

本研究では，最も単純且つ重要な系としてニッケルダイマ

番目のレベルは各々，水素の s 軌道，ニッケルの d 軌道が主ーを担持金属クラスターとして導入することで，水素解離にであるが，時間と共に両軌道が近づき，混成軌道を形成し， 電子の基及ぼす担持金属の効果を明確にした(3)．そのため

最終的にはニッケルの水素化物が形成される．ニッケルを導 ア ッ プ ス ピ ン の HOMO ( 最 高 被 占 有 軌 道 ，底状態と入しない場合，同様の一電子励起のシミュレーションをおこHighest Occupied Molecular Orbital)の電子を LUMO(最低なうと，水素分子の HOMO LUMO ギャップは 11.2 eV 程空軌道，Lowest Unoccupied Molecular Orbital)へ一電子励度と非常に大きい値を示す．これに対してニッケルダイマー起した場合について比較した．なお，本系においては，アッを導入することにより， HOMO LUMO ギャップは 0.1 eV

プスピンの HOMO は 30 番目のレベルに対応し， LUMO はのオーダーと，非常に小さい値になる事が分かる．このよう31 番目のレベルに対応する．ただし，シミュレーションをな HOMOLUMO ギャップを減少させる事がニッケルダイ遂行する際は，この約27倍である800番目のレベルまで考慮マーの効果である．このエネルギーは熱的に，あるいは分子した．これは，電子の励起状態を扱い時間発展を追う場合，の衝突のエネルギーなどで与える事が可能であると考えられ非常に高いレベルの軌道まで影響を受けるためである．

る．

原子の初期配置として，ニッケルダイマーと水素分子が互いに 90 度の角度で配置する系を導入した．ただし，テスト

このように，電子の基底状態と一電子励起状態の違いで，定性的にダイナミックスが異なる事が明確に示された．

計算として幾つかの初期配置についてシミュレーションをおこない，他の配置についても以下に述べる結果と同様の傾向

.

新バージョンに導入された高速化手法が得られたことを確認している．
図( a )( c )に，イオン位置の時間発展を示す．図( a )，新バージョンでは計算速度の高速化させるために，いくつ( b )，( c )は各々時刻 0 fs， 20 fs ， 40 fs における各イオンのかのアルゴリズムの修正をおこなっている．本節では，1 次配置のスナップショットである．この場合は基底状態とは異元球対称ポテンシャルを計算するための高速高精度 Fourier変換法と Chebyshev 多項式を利用した行列要素計算法(図の二重線) について紹介する(4)．
 次元球対称ポテンシャル計算


基底関数として平面波関数と原子軌道関数を用いる混合基

図3

一電子励起状態での，ニッケルダイマーと水素分子 の イ オ ン 位 置 の 時 間 発 展 ． ( a ) 0 fs ， ( b ) 20fs ，( c ) 40 fs のスナップショット，および( d )各エネルギーレベルの時間発展．31番目のレベルに
あ る 電 子 が ， 約 10 fs 後 に 30 番 目 の レ ベ ル に 移り，水素分子が解離する事が分かる．

図4

計算の流れ図． Kohn Sham ポテンシャル計算と行列要素計算の高速化方法については本文を参照．

ミ ニ 特 集底法においては，以下の Fourier 変換で定義される j 番目のーチンを新たに作成した．この方法を用いることで，行列要原子球内ポテンシャル Vj(rj)[rj＝|r－Rj|(Rj は j 番目の原子素を評価するための数値積分が計算コストの小さいフィッテ核位置)]を計算する必要がある．

ィング操作に置き換えられる．実際のポテンシャル関数と近Vj(rj)＝∑nã (G) exp (iG･Rj)
G

sin|G|rj| G| r j(6)似関数との誤差が極小になるようにするため(17)(18)，ポテン
シャルフィッティングに用いる関数として Chebyshev 多項

ここで， G は逆格子ベクトル， nã ( G )は Kohn Sham ポテン式を利用した(19)．具体的には，以下のようにアルゴリズム

シャルの滑らかな部分の Fourier 成分である． Vj( rj )の計算を改良した．に必要な演算数は， Fourier 空間内のグリッド数と原子球内の動径方向のグリッド数の積(～ 643 × 621 )に比例した極めて大きな数である．

あらかじめ j 次の Chebyshev 多項式に Tj(r)関する行列と〈AO|Tj(r)/r|AO〉を計算する．
要素〈PW|Tj(r)/r|AO〉



各計算ループ内で，ポテンシャル rV ( r )を Chebyshev

我々は計算時間を短縮させるために，式( 6 )の G に関す多項式で展開する V (r )は原子核位置(r ＝0)で発散するる和を G の角度に関する和と G の絶対値に関する和の 2 つので，便宜上，関数 rV ( r )を Chebyshev 多項式で展開に分け，次式で与えられる球対称ポテンシャルを導出した．

する．）．

[] sinGGr r) exp (iG′･ R j)
Vj(rj)＝∑ ∑nã (G′iG′

i j(7 )i j

ここで内側の和は Gi|G′|Gi＋1[Gi＝iDG(DG はパラメー

予め計算しておいた Chebyshev 多項式に関する行列要素を用いて，〈PW|V(r)|AO〉と〈AO|V(r)|AO〉を計算する．

タ)]を満たす G ′について，外側の和は i ＝ 1 から n ＝ 10000

このようにすると，各計算ループにおいて積分計算を実行す程度まで取られ，それぞれ G′の角度に関する和と絶対値にる必要がなくなる．以前の TOMBO のプログラムでは，原関する和に対応する．式( 7 )の演算数は 程度と見積もられ子核を原点とする動径方向の数値積分を 1 回実行するためるため，式( 7 )を用いる場合の計算時間は，式( 6 )を用いに，621回の和を実行する．一方，ポテンシャルフィッティる場合に比べて大幅に減少する．

ング法を用いると，フィッティング操作に用いる式( 7 )の近似表式の有効性を検証するため， N2 分子の全Chebyshev 多項式の数 M 回だけの和を実行する．孤立系とエネルギーを見積もった．計算には FUJITSUPRIMEHPC結晶系に対してテスト計算をおこなった結果， M  30 程度(1core/16threads)を利用した．表に示されるようで十分な計算精度が得られた．また，行列要素計算に必要なに，式( 7 )を用いることで，球対称ポテンシャル計算に必総演算数は M /621 1/ 20倍となり，行列要素計算の飛躍的要な時間が 10 分の 1 程度に減少し，計算時間が大幅に短縮な速度向上が達成された．

FX10した．さらに，式( 6 )と式( 7 )を用いて評価した全エネル

.

ギーの差は，0.001以下の誤差の範囲内で一致した．


Chebyshev 多項式を利用した行列要素計算

混合基底法を用いる場合，KohnSham ポテンシャルの行まとめ

本研究では，全電子混合基底法プログラム TOMBO を使った材料開発の一例として，時間依存密度汎関数理論に基づ列要素は，〈 PW | V ( r )| PW 〉，〈 PW | V ( r )| AO 〉，およびくエネルギー材料設計研究を紹介した．また TOMBO の新〈AO|V(r)|AO〉の 3 種類がある．ここで，|PW〉は平面波，規機能として，高速化手法の紹介をおこなった．

|AO〉は原子軌道関数，V(r) は KohnSham ポテンシャルをなお，現在開発中であり完成目前のサブルーチンとして，表す．〈 PW | V ( r )| PW 〉の計算には FFT のライブラリを用ワニエ関数を抽出して電子伝導を計算し，全電子の枠組みでいており，高速な計算が可能である．一方で，後者 2 つの初めて連続状態を取り入れた電子伝導計算がおこなえる予定

行列要素は数値積分によって計算される．このため，行列要である(20) ．また，混合基底に対して電子間クーロン相互作素を得るために必要な計算時間は平面波の数に比例して増大用を極めて効率的に計算するアルゴリズムにより RPA 分極する．用いる平面波の数が多い(エネルギーカットオフ値が関数を簡単に計算することができ，TDDFT に基づいて van大きい)場合，行列要素計算に多くの時間がかかる．

der Waals 力 の 計 算 が 非 常 に 精 密 可 能 と な る 予 定 で あ

行列要素の計算にかかる時間を短縮させるために，ポテンる(21)．また，原子変位に対する非調和項の計算(22)を可能とシャルフィッティング法を用いて行列要素を評価するサブルすることにより，熱伝導と熱電変換の計算プログラムをリンクする予定である．
さらに，TOMBO コードの LDA の部分については，ソー表1式( 6 )，( 7 )の表式を用いた場合の N2 分子の電子状態計算に必要な時間と全エネルギー．
式( 6 )式( 7 )

Vj(rj)の計算時間[s]1332133

全エネルギー[eV]－2949.310－2949.321

ま て り あ 第53巻 第 9 号(2014)
Materia Japanスコードを公開予定である．
TOMBO を使った研究に興味がある方は，筆者に問い合わせを頂きたい．
本研究の遂行にあたり，計算の一部は東北大学金属材料研
究 所 の ス ー パ ー コ ン ピ ュ ー テ ィ ン グ シ ス テ ム SR16000 /M1，東京大学物性研究所のスーパーコンピュータシステム
FUJITSU PRIMEHPCFX10 を利用した．本研究の一部は，科学研究費補助金課題番号 25289218， 24560813 および
HPCI 戦略プログラム課題番号 hp120023， hp120040の補助による．また，計算材料科学イニシアティブ CMRI からの助成を受けた．

文

献( 1 ) K. Ohno, K. Esfarjani and Y. Kawazoe: Computational
Materials Science: From ab initio to Monte Carlo Methods,Springer Series in SolidState Sciences, SpringerVerlag,Berlin, Heidelberg, 129(1999), 4246.
( 2 ) Y. Kodama and K. Ohno: Appl. Phys. Lett., 96(2010), 034101.
( 3 ) R. Sahara, H. Mizuseki, M. H. F. Sluiter, K. Ohno and Y.
Kawazoe: RSC Advances, 3(2013), 12307.
( 4 ) S. Ono and K. Ohno: Book of Abstracts of the 8th General
Meeting of Asian Consortium on Computational Materials
ScienceVirtual Organization (ACCMSVO) PS9, (2013).
( 5 ) C. Pannemann, V. Dyakonov, J. Parisi, O. Hild and D. W äohrle:
Synthetic Metals, 121(2001), 1585.
( 6 ) G. Matsunobu, Y. Oishi, M. Yokoyama and M. Hiramoto:Appl. Phys. Lett., 81(2002), 1321.
( 7 ) M. E. Milanesio, M. Gervaldo, L. A. Otero, L. Sereno, J. J.
Silber and E. N. Durantini: J. Phys. Organic. Chem., 15(2002),
844.
( 8 ) C. Devadoss, P. Bharathi and J. S. Moore: J. Am. Chem. Soc.,
118(1996), 9635.
( 9 ) M. Nakano, M. Takahata, S. Yamada, R. Kishi, T. Nitta and K.

Yamaguchi: J. Chem. Phys., 120(2004), 2359.
(10) Y. Kodama and K. Ohno: J. Chem. Phys., 125 (2006), 054501.
(11) O. M. Yaghi, M. OKeeffe,Á
N. W. Ockwig, H. K. Chae, M.
Eddaoudi and J. Kim: Science, 300(2003), 1127.
(12) Y. W. Li and R. T. Yang: J. Phys. Chem. C, 111(2007), 11086.
(13) M. A. Miller, C.Y. Wang and G. N. Merrill: J. Phys. Chem. C,113(2009), 3222.
(14) H. Nishihara, P.X. Hou, L.
X. Li, M. Ito, M. Uchiyama, T.
Kaburagi, A. Ikura, J. Katamura, T. Kawarada, K. Mizuuchi
and T. Kyotani: J. Phys. Chem. C, 113(2009), 3189.
(15) Y. Li, F. H. Yang and R. T. Yang: J. Phys. Chem. C, 113(2007), 3405.
(16) P. P. Galina, G. N. Vayssilov and N. R äosch: J. Phys. Chem. C,111(2007), 14484.
(17) 一松 信数値計算，近代数学新書，(1963).
(18) G. Dahlquist and A. Bj äorck: Numerical Methods, Englewood
Cliffs, NJ: PrenticeHall, (2003).
(19) W. M. Press et al.: Numerical Recipes in Fortran 90: The Artof Parallel Scientific Computing, Cambridge University Press,(1996).
(20) Y. Y. Liang and Y. Kawazoe: Book of Abstracts of The 20th
Anniversary of TOMBO and Russian Megagrant Opening
International Conference, (2013), 31.
(21) R. V. Belosludov, H. Mizuseki, R. Sahara, Y. Kawazoe, O. S.
Subbotin, R. K. Zhdanov and V. R. Belosludov: Computational
Materials Science and Computeraided Materials Design andProcessing, Handbook of Sustainable Engineering, J.
Kauffman. K. M. Lee (eds.), Springer Science＋BusinessMedia Dordrecht, (2013), 12151247.
(22) K. Esfarjani and H. T. Stokes: Phys. Rev. B, 77(2008),
144112.

ミ ニ 特 集