東京大学 藤原研究室,研究室概要
当研究室では、主に下記2テーマに沿った研究を行い、物質材料設計のための統一的な第一原理電子構造理論の確立を目指しています
詳 細
(1)大規模系における量子力学的第一原理分子動力学計算の開発
産業における有用性を念頭におきつつ、 100nmスケールサイズの原子数を扱うことの出来る第一原理分子動力学シミュレーション技術の開発と応用

半導体産業では微細加工技術の基本スケールが100nmを下回り、 ナノ構造系(nm又は10nmスケール系)が精力的に研究されています。ナノ構造では、 異なる電子状態をもつ複数の領域(例えば、バルク領域と表面領域)が、比較しうる原子数をそれぞれもち、 それらの共存・競合により多彩な構造や物性が現れます。 ナノデバイスの設計・開発にあたり、動的プロセスの理解と制御が必要不可欠であり、 計算機を活用した新しいシミュレーション法が期待されています。 しかし、通常の量子力学的第一原理分子動力学計算では、扱える原子数は数百程度であり、 産業界で理論的解析が望まれる100nmスケールサイズの原子数を扱うことは困難です。

当研究室では、100万原子あるいはそれ以上の原子からなるナノ・モデル系に適用できる第一原理分子動力学シミュレーション技術の開発と応用を目指しています。

※ELSES(Extra Large Scale Electronic Structure)研究会の立ち上げ
2007年11月、東京大学産学連携本部の産業連携創出スキーム「UCR-研究会」に基づき、 当研究室をはじめ、大学・企業・研究機関の有志によりELSES(Extra Large Scale Electronic Structure)研究会を立ち上げました。(代表:藤原毅夫教授) 当研究室で開発した大規模系第一原理分子動力学シミュレーションソフトウェアの実証・発展を進め、 ナノテク、材料分野における日本発の世界的ソフトウェアに開花させることを目指しています。
詳細はELSES研究会ホームページをご覧下さい: http://www.elses.jp/
手法:
図1に、超大規模系(104-107原子系)の計算例を示しています。 従来法以外の手法は、計算量がシステムサイズ(原子数)N に比例しています(オーダーN法)。 そこでは、一電子固有状態の代わりに一体密度行列やグリーン関数が用いられ、 量子力学の数理構造にとどまりながらも、計算量が劇的に縮小されています。

一体密度行列ρと物理量<X> は、占有された一電子固有状態{φk }の重ねあわせで形式的に定義されます:

ここで演算子X(r', r ) が短距離演算子だった場合、 密度行列ρ(r', r ) の長距離成分は(たとえ存在しても)物理量<X>に全く寄与せず、 計算する必要はありません。 この事実が、超大規模計算理論の原理的基礎となっています。

当研究室では、大規模系電子構造計算に対する線型計算手法として、 一般化Wannier状態またはKrilov部分空間に基づく手法を構築しました。 Krylov部分空間計算法については、特に2種類の手法(Krylov部分空間対角化とShifted-COCG法)を開発しました。 大規模系の電子構造計算解法には、問題を固有値方程式として解くか、 あるいは連立1次方程式として解くかの選択肢があります。

固有値方程式として解く方法の1つがLanczos法であり、 あるいは我々が開発したKrylov部分空間対角化法です。 また、連立方程式として解く方法は、これも当研究室が開発したshifted-COCG法を用いています。 shifted-COCG法では、ある1つのエネルギーを参照エネルギーとして選び、 そこで連立方程式を解けば、他のエネルギーの結果はその結果の線型計算で求めることができる点が重要な特徴の1つです。

またshifted-COCG法は、計算が安定で、密度行列の非対角成分を求めるのも容易です。 さらに任意の1つのエネルギーでの結果を、Krylov部分空間の不変性を用いて、 他のエネルギーにシフトすることができる大きな利点があり、 計算を加速することができる最も重要な特徴の1つです。

また本研究では、スレーターコスター形式にマップされたハミルトニアンを用いています。 いずれの手法でも、厳密な基礎方程式から微視的(基底)自由度ごとに残差を計算し、精度保証が行われます。 特に、ナノ構造では領域ごとで電子状態が本質的に異なり、 例えば表面部分の残差がバルク部分より大きい場合があります。 微視的精度保証は、従来は必ずしも強調されてきませんでしたが、手法確立の上では必要不可欠です。

開発した第一原理分子動力学シミュレーション手法は、数万原子規模の系については問題なく計算実行可能であり、 テスト計算では1千万原子まで行っています。

    図1:超大規模電子構造計算例

金属系(fcc Cu, 液体C)、
半導体系(バルクSi)に対する計算時間。
EIG:従来手法(固有状態を計算)
KR-SD: クリロフ部分空間法
WS-VR: 一般化ワニア状態に対する変分解法
WS-PT: 摂動解法。
KR-SD, WS-VR, WS-PTではいずれも
密度行列を計算。

研究成果:
(1)大規模線型計算手法の開発
● Krylov部分空間法とshifted-COCG法の確立

図2:Si結晶の亀裂伝播過程


図3:Auナノワイヤの構造変化

[R. Takayama, T. Hoshi, T. Sogabe, S.-L. Zhang, and T. Fujiwara, Phys. Rev. B 73, 165108, 9 (2006)]
[R. Takayama, T. Hoshi, T. Fujiwara, J. Phys. Soc. Jpn, vol. 73, No.6, pp.1519-1524 (2004)]

大規模系電子構造計算に対する線型計算手法として、 取り扱う行列を比較的低次元の部分空間に射影するKrylov部分空間法を提案しました。 更に、射影したKrylov部分空間に対し効率よく計算を行うShifted-COCG法は開発しました。 Shifted-COCG法は、問題を連立1次方程式として解くため、 固有状態を計算することなくグリーン関数の計算が可能になる。 数理的問題は、多数のエネルギー点における線形方程式系に帰着され、 実質的計算は1つのエネルギー点だけ行うことで、全エネルギー点に対する数値解を得ることができます。 本手法を、半導体物質系と金属物質系の双方に適用しました。

(2)現実の大規模系への量子力学的第一原理分子動力学計算の適用
● Si結晶中における亀裂伝播と表面再構成
[T.Hoshi, Y.Iguchi and T.Fujiwara, Phys. Rev. B 72, 075323 (2005)]
[T.Hoshi, R. Takayama, Y.Iguchi and T.Fujiwara, Physica B: Condensed Matter 376, pp. 975 (2006)]
[T. Hoshi and T. Fujiwara, J. Phys.: Condens. Matter 18 10787-10802 (2006)]

Si結晶の[001]方向に対して外力を与えた場合の亀裂伝播、形成表面の構造再構成、 容易亀裂面などの解析を行いました。 シミュレーションでは、(001)面の亀裂が、(111)面、 (110)面に沿った方向に亀裂の進行方向を変える様子が観測されました。 このことは、(111)面、(110)面がSiの容易亀裂面であるという実験事実と一致します。 (図2)

● Auヘリカルナノワイヤの多重殻構造
[Y.Iguchi, T.Hoshi and T.Fujiwara, Phys. Rev. Lett. 99, 125507 (2007)]

Auナノワイヤに対する分子動力学計算により、ナノワイヤの多重殻構造におけるマジックナンバー (隣接する2ワイヤ殻の電子数の差が7であること)、 螺旋構造などの原因を明らかにしました。 Au結晶中(面心立方格子)の(001)面上で原子列のずれが発生し、 表面全体が(111)面によって覆われた構造が発生します。 この構造変化が螺旋構造を誘起し、マジックナンバーを生み出します。 この計算結果は、当研究室で開発した第一原理分子動力学シミュレーションが金属系にも有効であることを示しています。 (図3)

(2)強相関電子系に対する第一原理電子構造計算の開発
多彩な物理現象を持ち、今後ますます重要な材料となる強相関電子系に適用できる第一原理電子構造計算手法の開発。 特に、1電子バンド抽象と多電子抽象を融合した第一原理電子構造計算手法の構築、 及び強相関電子系への応用を目指す。
電子の跳び移り積分に比べ、電子間相互作用が強く働く系は強相関電子系と呼ばれています。 強相関電子系では、金属ー絶縁体転移、巨大磁気抵抗効果(CMR)、 超伝導といった強相関電子系特有の多彩な物理現象が現れ、これらの性質を元にして、 新しい化合物やデバイスが次々と作られています。 強相関電子系の物性においては局在したdまたはf電子にはたらく強い電子間相互作用が重要な役割を担っています。

従来のLDA(局所密度近似)は、1電子バンド描像に基づいており、特に電子間相互作用が弱い場合においては、 基底状態を正しく記述し非常に良い結果を与えます。 例えば、半導体や通常の金属などの系における電子のふるまいは、LDAにより定性的に理解することができます。

しかし、強相関電子系では、バンドギャップやバンド幅などの基本的な物性量を第一原理的に定量的に決めることは困難です。 さらに、正しい電子構造や磁性を再現できない傾向にあります。 この理由は、上記物性には電子の多体効果が決定的に重要な効果を与えるからです。 一方、従来のLDAでは、電子の多体効果を平均的にしか取り扱うことができません。 そこで、強相関電子系の物理現象を定量的に理解するためには、 電子間相互作用を多体問題として具体的に取り扱うことのできる新たな基礎的方法論を開発することが必要となります。

当研究室では、1電子バンド抽象と多電子抽象を融合した第一原理電子構造計算手法の構築、 及び強相関電子系への応用を目指します。

手法:
ここでは、当研究室で開発を目指している代表的な2つの手法について簡単に説明します:

●LDA+DMFT法
LDA+DMFT法は、動的平均場理論(DMFT)をLDAと複合することにより、 現実物質の第一原理的解析において、サイト内の動的な電子相関を取り扱うことができます。

DMFT(動的平均場理論)は、無限次元での電子系を取り扱う手法として提案されました。 無限次元では、周りのサイトの数が無限に多いため、その及ぼす影響を「平均場」で置き換え、 空間の揺らぎを無視することができます。 この平均場に対応するのが電子の自己エネルギーであり、周波数のみに依存し、 波数に依存せず局所的となります。 これが「動的」平均場理論と呼ばれる所以です。 このとき、格子系の多体問題は有効媒質中に埋め込まれた1不純物問題へと射影されます。 射影された1不純物問題では、多体問題を正確に扱うのがDMFTの特徴です。 以上の手続きにより、DMFTでは、サイト内の動的な電子相関の取り扱いが可能になります。

もし、局所グリーン関数G(w)を計算する際、LDAハミルトニアンを用いると、 系のバンド構造をあらわに取り入れることが可能になり、DMFTによる第一原理的解析が可能になります。 このようなLDAとDMFTとの複合手法をLDA+DMFT法と呼んでいます。 (図4)


当研究室では、1不純物問題の解法にはIPT(逐次摂動近似法)を用いています。 IPTでは、2次摂動の自己エネルギーを用いて、2つの極限(高周波数極限、原子極限)からの内挿により、 自己エネルギーを求める方法です

○当研究室でのLDA+DMFT法の特徴
LDA:+DMFT法は近年盛んに研究がなされており、様々な種類が存在します。
本研究で提案しているLDA+DMFT法の特徴は下記3点です:

1) ハミルトニアンを射影せずLDA ハミルトニアンを直接DMFT に適用
2) 有効1不純物問題の解法への逐次摂動近似法(IPT)の利用
3) 配位子場理論(LFT)の利用によるバンド描像と原子系の多電子描像と融合

以上の特徴を全て兼ね備えるという点で、 本研究でのLDA+DMFT法は複数原子系からなる化合物系はもちろんのこと、 s,p軌道とd軌道が複雑に混成した系といった広範な強相関電子系物質への取り扱いが可能な手法になっています。


図4:LDA+DMFT法の概略
●GW近似
GW近似は、RPA(乱雑位相近似)の範囲で、動的な遮蔽効果を取り扱います。 すなわち、多体摂動論を出発点とし、 自己エネルギーをグリーン関数Gと動的に遮蔽されたクーロン相互作用Wの積という最低次の摂動項で近似します(図5)。 これが「GW」近似と呼ばれる所以です。 また、GW近似の枠組みでクーロン相互作用Uの値を第一原理的に見積もることも可能になります。

○本研究でのGW近似の特徴
通常のGW近似は、計算コストが大きいため、単位胞に数原子程度の系にしか適用することはできません。 当研究室では、アルゴリズムの開発を行いGW近似の高速化を行いました。 現在、当研究室のGW近似プログラムは、単位胞に20原子以上の系(磁気秩序を含む)に適用することが可能になっています。


図5:GW近似のファイマン図形による表現。
研究成果:
GW近似によるペロブスカイト型遷移金属酸化物LaMnO3の電子構造と準粒子の寿命

図6:GW近似によるLaMnO3の部分状態密度。


図7:LaMnO3のΓ点での
LDAの各固有値に対するスペクトル。

[Y. Nohara, A. Yamasaki, S. Kobayashi, T. Fujiwara, Phys. Rev. B 74, 064417 (2006)]

GW近似を用いて、ペロブスカイト型遷移金属酸化物LaMnO3の解析を行いました。 バンドギャップは、LDAでは過小評価されていたものの(0.44eV)、GW近似では、 実験結果とよく一致した値(1.0eV)を得ました。 Mnのt2g軌道のバンド幅がeg軌道のそれより大きく収縮しており、 このことはt2g軌道よりeg軌道の方がより強い遮蔽効果が効いていることを示唆しています。
サイト内の動的遮蔽されたクーロン相互作用Wは低エネルギー領域で強く遮蔽効果が効いています。 また、正孔の寿命が電子の寿命に比べてはるかに短いことも見出しました。 これらは、Mnのeg軌道の電子が動きやすいことに起因しています。(図6、図7)

●LDA+DMFT法による強磁性Fe,Ni、反強磁性NiOの電子構造
[O,Miura and T. Fujiwara, in preparation]

LDA+DMFT法を用いて、強磁性金属Fe,Niの解析を行いました。 占有状態の3dバンド幅はLDAでは課題評価されていたが、LDA+DMFT法では、実験結果とよく一致しました。 また、Niの2正孔束縛状態サテライト(フェルミレベルより下約6eV)も観測されました。 これらの描像は、DMFTによりサイト内の動的な電子相関を取り入れたことに起因しています。

また、LDA+DMFT法を用いて、反強磁性NiOの解析を行いました。 LDAでは小さいバンドギャップ(0.2eV)、 モット=ハバード型電子構造(下部ハバードバンドが酸素の占有2pバンドより高エネルギー側に位置する) が示唆されました。 一方、LDA+DMFT法では、バンドギャップ4.3eV、電荷移動型電子構造 (下部ハバードバンドが酸素の占有2pバンドより低エネルギー側に位置する)という結果が得られ、 分光実験の結果と非常によく一致しました。 これは、Niサイト内のクーロン相互作用が、 Niの3dバンドと酸素の2pバンドとの混成が大きく変化したことに起因します。

今回の研究結果は、当研究室で開発したLDA+DMFT法が、 混成の強い系に適用可能であることを実証しており、 今後、更に大きな化合物、より複雑な軌道混成を持つ物質に対しての適用の可能性を示唆しています。 (図8)


図8:LDA+DMFT法による
反強磁性NiOのスペクトル。
●La2-xSrxNiO4における電荷・スピンストライプ秩序
[S. Yamamoto, T. Fujiwara and Y. Hatsugai, Phys. Rev. B 76, 165114 (2007)]

La2-xSrxNiO4(x<1/2)では、 電荷及びスピン秩序を示し絶縁体です。 x=1/2の状態においてはLDA+U法では金属となり、スピンストライプ秩序も再現できません。 これは、サイト間クーロン相互作用がハートリー項では充分取り扱われていないこと、 波動関数が(LDA,LDA+Uでは暗黙の了解である)単一スレーター行列式では表現できないこと等によっています。 本研究では、LDA+U法により求められた軌道エネルギーと飛び移り積分の値を用いて、 2重縮退したeg軌道の拡張ハバード模型の(有限系)の厳密対角化を行いました。 その結果、サイト間クーロン相互作用Vと飛び移り積分の異方性δの競合により、 電荷及びスピンストライプ秩序が現れることを示しました。(図9)

●Shifted-COCG法の多電子問題への応用

研究成果3に関連して、 La
2-xSrxNiO4(x=1/2)の秩序状態を調べるため、 スペクトルを計算する際に新たにShifted-COCG法を用いました。 この研究成果により、これまで大規模系で用いてきたShifted-COCG法は、 多電子問題に対しても強力な手法になるということを示しました。 この中で、新たにseed-switchingの方法を開発しました。 一般にShifted-COCG法では、ある1つのエネルギーを三章エネルギーとして選び、 そこで連立方程式を解けば、他のエネルギーの結果はその結果の線型計算で求めることができる点が重要な特徴の1つです。 参照エネルギーを選びなおさなくてはならない時でも、新しい参照エネルギーに移して更に計算を続けることができます。Shifted-COCG法は安定な計算手法で、精度を機械精度まであげることが可能です。


図9:La2-xSrxNiO4(x=1/2)に
おける励起スペクトルの
サイト間クーロン相互作用 V 依存性。