物性物理・数学・計算科学の融合

 物性物理・数学・計算科学の融合
藤原毅夫先生の東大・物工 ご退職を記念して
平成１９年２月
前書き
この度、平成１９年３月に、藤原毅夫先生が６２歳の東京大学の定年を
お迎えになられ、東京大学物理工学科におけるご研究、ご教育に一つの区
切りをつけられることとなりました。
これに先立ちまして、平成１７年６月１８日に奥様をお招きいたしまし
て文京区本郷のフォーレスト本郷にて還暦のお祝いをかねて藤原研究室の
同窓会を開かせていただきました。会場の様子につきましては本文中の写
真をご覧頂けましたら幸いです。
続きまして平成１８年１０月７日には関係者各々の近年の研究につきま
して各自の近況のご報告を兼ねて「物性物理・数学・計算科学の融合」と
題しました研究会を開かせていただきました。この研究会のタイトルなら
びに本文に掲載いたします研究会プログラムにありますように関係者が極
めて幅広いスペクトルにわたる分野において研究を行っていることは先生
の広いご見識の一端を垣間見るところであります。
また、この機会をもちまして先生にご指導ご鞭撻いただきました関係者
のこの時点での記念碑として、各々の最近の研究ならびに考えていること
等を寄稿の形であつめ、長年のご指導を感謝しつつ、先生にご覧いただこ
うと考え、この小冊子をまとめることといたしました。
最後になりますが、藤原研究室卒業生ならびに関係者一同、今日までの
数多くの辛抱強いご指導ご鞭撻に心より感謝するとともに、先生のより一
層のご健康をここに願うものであります。
初貝 安弘（藤原研究室 昭和６０年修士入学）
新井 正男（藤原研究室 昭和６１年修士入学）
山元 進 （藤原研究室 平成３年修士入学）
目次
研究会プログラム .........................................................................................iv
記念写真 ........................................................................................................v
寄稿文集
フェルミ面の異方性が超伝導体の上部臨界磁場に与える影響
物質・材料研究機構 新井正男
M. Arai, NIMS ...............................................................................1
新しい秩序相としての量子液体相
東京大学 初貝安弘
Y. Hatsugai, Univ. of Tokyo ............................................................13
電子構造理論からナノテクノロジーへ
— 藤原毅夫先生との共同研究を中心として —
鳥取大学 星健夫
T. Hoshi, Tottori Univ. ....................................................................25
準結晶の電子構造 －第一原理計算が教えること－
中央大学 石井靖
Y. Ishii, Chuo Univ. ........................................................................45
ものつくりにおける計算材料科学
(株) 東芝 伊藤聡
S. Itoh, Toshiba Corp. .....................................................................62
一化学企業における第一原理計算 —無機蛍光体を中心に—
(株)三菱化学 三上昌義
M. Mikami, Mitsubishi Chem. Corp. ...............................................73
Efficient Linear Scaling Approaches for Computing Quantum Transport in lowdimensional systems
S. Roche, CEA/DSM/DRFMC/SPSMS/GT ....................................88
ii
A method to evaluate electrostatic potential from X-ray diffraction data by using
MEM
島根大学 田中宏志
H. Tanaka, Shimane Univ. ..............................................................117
Effect of transition-metal elements on the electronic properties of quasicrystals and
complex aluminides
G. Trambly de Laissardière1 and D. Mayou2
1
Laboratoire de Physique Théorique et Modélisation, CNRS
2
Institut Néel, CNRS .......................................................................128
計算物理と計算機科学との距離
東京大学 山元進
S. Yamamoto, Univ. of Tokyo .........................................................145
共有メモリ型並列計算機向け対称密行列固有値解法の性能評価
名古屋大学 山本有作
Y. Yamamoto, Nagoya Univ. ...........................................................156
iii
研究会プログラム
物性物理・数学・計算科学の融合
～藤原毅夫先生の東大・物工 ご退職を記念して～
日時 2006年10月7日(土) 午後2時より
場所 学士会分館
スケジュール
2:00 – 2:05 開会の挨拶
2:05 – 2:40
伊藤 聡 氏　（東芝 先端機能材料ラボラトリー）
「ものつくりにおける計算材料科学」
2:40 – 3:15
時弘 哲治 氏（東大 数理科学研究科）
「超離散系に関する最近の話題」
3:15 – 3:35 Tea break
3:35 – 4:10
常次 宏一 氏（物性研）
「準結晶からフラストレーション系まで」
4:10 – 4:45
山本 有作 氏（名大 計算理工学専攻）
「大規模固有値解法の最近の進展 – 密行列向け解法を中心に –」
4:45 – 4:55 総括　(藤原毅夫先生)
4:55 – 5:00 閉会の挨拶
5:00 – 7:00 懇親会
iv
v
2005 年 6 月 18 日 フォレスト本郷にて
vi
フェルミ面の異方性が超伝導体の上部臨界磁
場に与える影響
物質・材料研究機構 計算科学センター 新井正男
藤原研究室には大学院の学生として４年、そして助手として５年近く在
籍させていただきました。その間に、固体中の電子構造に関する理論的研
究について、藤原先生から様々なご指導を戴いたことに感謝します。研究
室に在籍していた間に広範囲の研究対象に触れたことが、その後の研究の
原点となりました。ここでは、磁場中超伝導体の物性を第一原理計算から
得られる電子構造データとモデル計算を用いて定量的に解析する研究に関
して紹介します。
1
はじめに
固体物理で最初に学ぶ基本事項の一つが、フェルミ面である [1]。フェ
ルミ面の存在が電子構造から見た金属の特徴であり、電子物性に対して重
要な役割を果たす。例えば、低温で温度に比例する電子比熱や、金属的伝
導はフェルミ面の存在に起因する。フェルミ面構造の決定は古典的金属物
理学の中心課題であり、現在でも新規に発見された金属的物質群に対して
1
電子構造と物性を議論する際の出発点である。幸いにも第一原理的手法の
発展と計算機能力の進歩により、ほとんどの物質の電子構造を理論的に計
算することが可能になった。電子相関が特異な影響をおよぶす物質をのぞ
けば、実験結果と定量的に比較できるバンド構造、フェルミ面を第一原理
計算から求めることができる。
バンド構造やフェルミ面に関する情報から、電気伝導度をはじめとする
輸送現象を評価する手法としては、ボルツマン方程式による半古典的な扱
いが広く用いられている。衝突項を緩和時間 τ で近似すると、金属の電気
伝導度はフェルミ速度 v =
∂εk
∂k
を用いて、以下のように評価することが
できる。
σαβ = e2 τ hvα vβ i
(1)
ここで h· · · i はフェルミ面上での平均操作である。このような表式は単純
金属、合金等の物性評価に古くから利用されてきた。高温超伝導体が発見
されたときには、その正常金属相を電子論から評価するための試みとして
利用された。藤原先生が準結晶に対する先駆的な第一原理計算を実施した
際にも、この手法で輸送係数を評価した [2]。
最近我々は、フェルミ面構造の詳細が超伝導体の物性に対する影響を理
論的に研究している [3, 4]。いわばボルツマン理論による金属相の解析を
超伝導相へと拡張する研究である。超伝導の物性に対してフェルミ面の構
造が影響を及ぼすことは古くから指摘されてきたが、定量的な解析は進ん
でいない。本論では、超伝導体の上部臨界磁場の評価に応用した例を紹介
する。
2
2
上部臨界磁場
超伝導体に磁場を印加すると完全反磁性を示し、磁場をその内部から排
除する。この現象はマイスナー効果と呼ばれ超伝導の本質的性質である。
外部磁場がある強さを超えると、超伝導体内に磁場が侵入し、超伝導状態
が破壊される。その際、磁場の侵入方法により、超伝導体は第一種と第二
種に分類される。第一種超伝導体では磁場が侵入すると即座に超伝導が破
壊されるのに対して、第二種超伝導体ではアブリコソフ磁束と呼ばれる形
で磁場が部分的に物質内部に侵入する。磁束の中心殻部分で超伝導の秩
序変数 ∆(r) がゼロになるが、物質全体では超伝導性を維持し続ける。侵
入する磁場が増加するにつれて |∆(r)| が減少し、外部磁場が上部臨界磁
場 Hc2 に達したときに、物質中の ∆(r) が消失し、超伝導性が失われる。
Hc2 は磁場中での超伝導体の性質を考えるうえでも応用上の観点からも重
要な量であり、実験・理論の両面から盛んに研究されている。
Hc2 近傍では ∆(r) が微少量になることを考慮すると、∆(r) に対して
線形化した基礎方程式を用いて Hc2 を理論的に求めることができる。例
えば、超伝導に対する現象論であるギンツブルグ・ランダウ理論 (GL 理
論) を線形化すると、次の形の方程式が得られる。
1
2m∗
µ
e∗
−ih̄∇ − A
c
¶2
Ψ = −αΨ
(2)
α = a(T − Tc )
(3)
この方程式は磁場中のシュレディンガー方程式と同じ形をしている。上式
が非自明な解を持つという条件から、上部臨界磁場 Hc2 を決めることが
できる。ただし、GL 方程式は Tc 近傍の領域を記述する現象論であり、低
3
温での Hc2 を決定するには微視的な理論である BCS 理論に立ち返る必要
がある。Werthamer らは球状のフェルミ面に対する Hc2 を BCS 理論から
導出した [5, 6, 7]。彼らの結果によれば、不純物のない第２種超伝導体で
は規格化した上部臨界磁場 h∗ ≡ Hc2 /(−dHc2 /dt)|t=1 ,
t ≡ T /Tc の低温
極限での値が h∗ (0) ≈ 0.73 となる。この値は現在でも Tc 近傍での測定値
から低温での Hc2 を見積もる際に利用されている。
一方、典型的な第２種超伝導体である Nb については、Hc2 の温度依存
性、方向依存性が精密な測定により調べられている [9]。その結果による
と低温極限での Hc2 の値は h∗ (0) ≈ 1.06 となり、等方的フェルミ面に対す
る理論予測から大きくずれている。また、低温での Hc2 は磁場の方向に依
存して 10%程度の異方性が存在する。この温度依存性や印加方向依存性の
起源としてはフェルミ面の構造が影響していると考えられ、フェルミ速度
の異方性をパラメーターとして用いた解析が報告されている。Butler[10]
は、電子構造計算によって得られたフェルミ面データを用いて Hc2 を理論
計算し、その角度依存性について議論したが、温度依存性に関しての議論
は不充分なままであった。
我々は、第一原理計算による電子構造データを利用して、パラメーター
を用いずに Hc2 の温度依存性と方向依存性を理論的に計算することに成
功した。この計算手法を用いるとフェルミ面の異方性を完全な形で取り込
むことができるため、実験結果との定量的評価が可能になる。
計算には準古典理論から導いた代数方程式を用いた。不純物による散乱
を無視したクリーンリミットでは、次の行列の固有値がゼロになるという
4
条件から Hc2 を求めることができる [3]。
¸
∞ ·
X
T
δN N 0
−hKN N 0 i .
AN N 0 ≡ δN N 0 ln +2πT
Tc
εn
(4)
n=0
ここで N は秩序変数 ∆(r) をランダウレベルに対する関数系で展開した
際のインデックスであり、εn は松原周波数である。K はフェルミ速度と
印加磁場 H に依存する量である。第一原理計算によって得られるフェル
ミ面構造とフェルミ速度を入力として上式の hKN N 0 i を数値的に評価し、
Hc2 を印加磁場方向と温度の関数として求めた。以下、簡単なモデルと
Nb, MgB2 に対して適用した結果について紹介する。
3
タイトバインディングモデル
まず、最も簡単なタイトバインディングモデルに適用した結果を紹介す
る [3]。単純立方格子で各サイトに１種類の原子軌道を置き、最近接間の
飛び移り積分を (−t) とした場合のエネルギーバンドは次式になる。
εk = −2t {cos(kx a) + cos(ky a) + cos(kz a)}
(5)
対応するフェルミ面を図 1 に示す。以下 t = a = 1 とする。フェルミエネ
ルギー εF が-6 近傍ではバンド分散が自由電子的な性格を持ち、フェルミ
面は球状になる。εF が増加するにつれて、タイトバインディングモデル
の性格を反映したフェルミ面構造へと変化する。εF = −2 のときにフェル
ミ面がブリルアンゾーン境界に触れ、トポロジー構造が変化する。εF = 0
のときが、所謂ハーフフィルド状態である。
このバンド構造に対して、等方的な s 波超伝導を仮定して計算した Hc2
を図 2 に示す。図 2(a) は規格化した臨界磁場 h∗ の温度依存性である。実
5
図 1: タイトバインディングモデルのフェルミ面. フェルミエネルギーの
値はそれぞれ εF = -3 (a), -2 (b), -1 (c), and 0 (d).
線が ε = −3 の結果で、点線が εF = −2.02 のときの結果である。それぞ
れ３本のグラフは上から順に印加磁場方向が [111],[110],[100] の場合であ
る。一点鎖線は等方的な丸いフェルミ面に対する h∗ である。Tc 近傍では
Hc2 の磁場方向依存性が消失するが、低温になるにつれて方向依存性が増
大し、等方的モデルに対する値からずれる。特に εF = −2.02 では顕著な
方向依存性が現れ、h∗ (0) は Helfand らによる等方的フェルミ面に対する
理論値に比べて２倍近く大きな値となる。
図 2(b) に低温での h∗ を εF の関数として示す。この図から、h∗ は εF =
−2 の近傍で磁場方向依存性が増大し、低温極限での h∗ (0) が著しく増加
することが読み取れる。εF = 2 のときに、フェルミ面は X 点でブリルア
ンゾーン境界に接触する。その際に X 点でのフェルミ速度がゼロとなる
ことが、h∗ の振る舞いに特異性を与えていると考えられる。
このように単純なタイトバインディングモデルでさえも、Hc2 の挙動を
理解するためには、フェルミ面の詳細な情報を考慮することが重要であ
る。現実に存在する超伝導体の Hc2 を解析する場合にも、電子構造の詳
6
(b)
(a)
1.2
1.6
1.4
h (0.05)
0.8
0.6
0.4
1.0
0.8
0.2
0.0
0.0
1.2
*
h*(t)
1.0
0.2
0.4 0.6
t = T/Tc
0.8
1.0
0.6
-6
-4
-2
0
F
図 2: 規格化した上部臨界磁場 h∗ . (a) 温度依存性, (b) T /Tc = 0.05 での
h∗ の εF 依存性
細を取り込んだ理論が不可欠である。
4
bcc Nb
Nb は典型的な第２種超伝導体として実験、理論の両面から盛んに研究
されてきた。第一原理計算によって求められたフェルミ面を図 3 に示す。
図 3: Nb のフェルミ面
7
このフェルミ面を用いて Hc2 を評価した。図 4(a) に角度平均した h∗ の
温度依存性を示す。実験によると低温極限での h∗ は 1.06 となり、フェル
ミ面の異方性を考慮した計算によって得られた 0.96 に比べて 10%程度小
さいが、球状フェルミ面に対する値 0.73 と比較すると大幅な改善を示す。
(a)
(b)
4
0.6
10
5
-3
a6 x10
0.2
0.2
0.4
0.6
T/Tc
0.8
1.0
6
4
2
3
(c)
2
1
0
(b)
-3
0
8
0.4
0.0
0.0
15
–a10 x10
*
h (t)
0.8
(a)
-3
–a4 x10
-3
20
a8 x10
1.0
0
0.0 0.2 0.4 0.6 0.8 1.0
T/Tc
3
(d)
2
1
0
0.0 0.2 0.4 0.6 0.8 1.0
T/Tc
図 4: (a) Nb に対する Hc2 の温度依存性 (b)Hc2 の異方性パラメーター al
の温度依存性
実験と計算の不一致の原因としては、フォノンによる電子間引力の影響
の取り扱いが不充分なことが挙げられる。今回の計算は弱結合 BCS 理論
の範囲で行っているために、フォノンによる電子間引力の波数、周波数依
存性が取り込まれていない。モデル計算等の結果によるとフォノンの波
数依存性でもっとも重要な寄与はフェルミ速度の繰り込みであることが知
られている。ドハースファンアルフェン効果を用いて実験的に求められた
フェルミ速度の繰り込み係数 λ を用いて Hc2 を再計算すると、図 4(a) に
破線で示した結果が得られ、実験結果と非常に良い一致を示す。この結果
は、フォノンの効果を定量的に取り込んだ強結合理論を利用すれば、より
8
実験と一致する Hc2 が得られることを示唆している。
Nb の Hc2 は低温になるにつれて異方性が増大する。結晶の持つ立方対
称性を考慮すると、Hc2 を全対称な立方調和関数 Pl を用いて次のように
展開することができる。


Hc2 (Ω, t) = H̄c2 (t) 1 +

X


al (t)Pl (Ω)
l=4,6,8,···

(6)
ここで Ω は印加磁場の方向、H̄c2 は角度平均した値である。al (t) が、Hc2
の異方性を特徴づけるパラメーターであり、その温度依存性を図 4(b) に
示す。今回の計算によって得られた al は実験結果をよく説明しているこ
とがわかる。言い換えれば、Nb に対する Hc2 の異方性はフェルミ面構造
に起因していることがこの計算から明確となった。
5
MgB2
MgB2 は 2001 年に青山大学の秋光先生らによって発見された超伝導体
であり、その転移温度は約 40K である。ボロンによる面の間に Mg が挟ま
れた層状物質であり、電子構造から見た特徴は、フェルミ面を横切る２種
類のバンドが存在することである。図 5(a) に第一原理計算によって得ら
れたフェルミ面を示す。フェルミ面を形成する状態は主としてボロンの p
軌道から派生される。上側の２つのフェルミ面は、p 軌道の σ バンドに起
因しており２次元性が強い。それに対して下側の２つのフェルミ面は π バ
ンドに起因しており異方性はそれほど強くない。M gB2 の超伝導ではこれ
ら２種類のバンドで超伝導ギャップの大きさが異なることが実験、理論の
双方から指摘されている。２種類のギャップの存在が様々な物性に影響を
9
及ぼし、M gB2 はマルチギャップ超伝導体として盛んに研究されている。
(b)
(a)
15
Hc2 (T)
Hc2//
10
5
0
0.0
Hc2
0.2
0.4
T/Tc
0.6
0.8
1.0
図 5: MgB2 の (a) フェルミ面と (b)Hc2 の温度依存性
Hc2 を計算する際には、σ バンドと π バンドの超伝導ギャップの比 α を
パラメーターとして取り込んだ。図 5(b) に α = 0.25, 0.3, 0.35 の場合の結
果を示した。黒丸で示した点が実験値 [11] を示し、下側のグラフが磁場
をボロン面に垂直に印加した場合で、上側のグラフが平行に印加した場合
の Hc2 である。今回の計算では Hc2 の絶対値を評価することができない
ため、磁場を面に垂直に印加した場合に低温で実験結果と一致するように
計算結果をスケールしている。この結果をみると α = 0.3 程度に選んだ場
合に計算によって得られた Hc2 は実験結果を良く再現している。α = 0.3
という値は、他の理論計算や実験結果と矛盾しておらず、Hc2 の解析から
も MgB2 の２つの超伝導ギャップの存在を裏付けることになる。
10
6
まとめ
Nb,MgB2 に対して適用した結果が示すように、第一原理計算から得ら
れるフェルミ面構造を取り込んだ定量的な計算は Hc2 の温度依存性と異
方性の理解に重要な役割を果たす。今回の計算では残念ながらフォノンの
効果を単純化した形で取り込んでいるために Hc2 の絶対値を議論すると
ころまでには至っていない。現在、強結合理論の範囲で定量的な計算を実
施することを目指して研究を進めている。特に最近の第一原理計算の進展
の一つである、密度汎関数法に対する摂動論を用いるとフォノンの分散関
係や電子格子相互作用の係数が定量的に評価できることが知られており、
今回の計算手法と組み合わせることにより非経験的な手法により Hc2 の
評価が可能になると、応用研究のうえからも重要な役割を果たすことがで
きると期待される。
謝辞
ここに紹介した結果はすべて北海道大学の北孝文先生との共同研究によ
るものです。
参考文献
[1] アブリコソフ著 金属物理学の基礎 (吉岡書店)
[2] T. Fujiwara, S. yamamoto, and G. T. de Laissardière, Phys. Rev.
Lett., 714166 (1993).
11
[3] T. Kita and M. Arai, Phys. Rev. B, 70, 224522 (2004).
[4] M. Arai and T. Kita, J. Phys. Soc. Jpn. 73 2924 (2004).
[5] E. Helfand and N. R. Werthamer, Phys. Rev. Lett., 13, 686 (1964).
[6] E. Helfand and N. R. Werthamer, Phys. Rev., 147, 288 (1966).
[7] N. R. Werthamer, E. Helfand, and P. C. Hohenberg, Phys. Rev.,
147, 295 (1965).
[8] P. C. Hohenberg and N. R. Werthamer, Phys. Rev. 153, 493 (1967).
[9] F. M. Sauerzopf, E. Moser, H. W. Weber, and F. A. Schmidt, J.
Low. Tem. Phys. 66, 191 (1987).
[10] W. H. Butler, Phys. Rev. Lett. 44, 1516 (1980).
[11] M. Zehetmayer, M. Eisterer, J. Jun, S. M. Kazakov, J. Karpinski,
A. Wisniewski, and H. W. Weber, Phys. Rev. B, 66, 052505 (2002).
12
新しい秩序相としての量子液体相
東京大学大学院工学系研究科 物理工学専攻
初貝安弘
1
はじめに
この度、藤原先生が東京大学の定年をお迎えになるに当たって、大学院
入学以来の長年のご指導ご鞭撻に感謝しつつ、固体中の電子に関します私
の最近の研究の一端をご紹介させていただきたいと思います。ただし、対
象は専門家ではなく、学部の学生程度を念頭に教育的意義に主点をおいた
ものとしたいと考えます。万一誤解を招く表現等がありましたら、ご指摘
の上お許しください。
2
幾何学的位相
まず、量子力学ならび量子力学に従う物理系においては複素数が本質
的であり複素数は絶対値と位相に分けられることに注意します。そこで量
子力学での位相の効果のどうしても取り除けない「本質的」部分が「幾何
学的位相」であると言えます。 技術的には量子系における記述はある空
間での演算子並びにある基底によるその行列要素を用いてなされます。そ
こでは種々の基底変換の自由度があり、それに対応して複素量としての行
13
列要素は変更を受けることに注意しましょう。この基底変換の自由度は完
全に自由ですが、実際の量子系における物理現象はこの任意の変換の自由
度に影響を受けることは決してなく、不変な形で表現されなければなりま
せん。数学的にはここでの（基底）変換に対応してある種のゲージ変換が
引き起こされることになるのですが、物理量はこのゲージ変換に対して
不変であるというわけです。簡単な多くの場合、基底変換は行列要素の複
素数としての位相の変化をもたらします。古い量子力学の教科書等ではこ
の基底変換に伴って起こる位相変化は意味がないとの記述すらあるのです
が、波動関数の位相なら何でもゲージ変換で完全に消去できるわけではな
く、一見この勝手に変化する位相のなかにも決して無視できない物理的意
義を持つものがあり、それらを称して幾何学的位相と呼ぶのです。
物性論における重要な概念である「量子力学的揺らぎ」、「量子干渉
効果」その他、いわゆる「量子効果」は最終的にはこの幾何学的位相とし
て理解されることが多いのです。
この幾何学的位相が重要な寄与をする物理現象の典型例としては量子
ホール効果、ベリー位相、アハロノフ・ボーム効果、分数統計粒子系など
があげられます。
3
対称性の破れと秩序変数
世の中には極めてたくさんのものつまり物質があり、それらはいろいろ
な形態をとっています。たとえば、水、炭酸ガス、シリコン、鉄などの物
質が氷、水蒸気、ドライアイス、結晶、非結晶、磁石、等の形態をとって
いるわけです。これら物質の形態の違いをきちんと理解することが物性物
14
理の基本であり、物質を用いた物質科学もこの作業なくしては近年のめ
ざましい発展はあり得ないものでした。これら物質の形態は物理的には
「相」として表現されます。水蒸気と氷は水という同じ物質ではあるがそ
の「相」が異なるというわけです。 物理学は現代科学ですから「相」に
対する（原理的には必ずしも定量的であることは必要でありませんが）科
学的な記述法を必要とします。その科学的記述の作業対象が「秩序変数」
と呼ばれるものです。「秩序変数」を観察、考察することにより物質の相
を判定するわけです。水の密度、鉄の磁石としての強度（磁化）がこの秩
序変数として使われています。たとえば、温度を変化させることで水の密
度の変化を観察することが秩序変数の温度依存性を議論することとなりま
す。ここである場所、ある時間での水の密度、鉄の磁化を議論することを
局所的な秩序変数を考えるといいます。多くの場合物質の「相」を特定す
るためには、この局所的な秩序変数を考えることが重要です。（双対変換
等、非局所的な理論的読み直しがこの局所的描像を使うために必要なこと
も時々あります。）
物理学、物性科学においては物質の「相」を特定することは物理系の
対称性と密接な関係があります。現代の物理学においては対称性を用いて
物質の相を区別するといった方が適切です。気体である水蒸気では水分子
は全くランダムに運動していると考えられますので、ある水の分子の近く
を見たとき隣の水分子は特定のどこかに存在しやすいことはありません。
つまりすべての方向は同等、一様です。一方、固体の水を考えたときは、
隣の水分子は全く無関係な位置にあるのではなくある秩序を持ってある方
向に存在することは想像できると思います。つまり固体の水の相において
は方向の一様性が失われていると考えられます。空間をある水分子の周り
15
に回転させるという操作を考えたとき水蒸気はその操作で不変、つまり対
称ですが、氷はその操作で不変でなく対称でないわけです。水蒸気を球に
たとえれば、氷は角がある箱だということになります。同じ水という物質
がその温度が変化することでその対称性を変えるわけです。水分子をたく
さん集めたとき特定の方向が特別な意味を持つことはなく、どちらの方向
も当然同等なはずです。その「同等性」が温度を変えるだけで、水分子の
お互いの相互作用により外力等を必要とせず、自発的に失われ特定の方向
に隣の分子が居やすいという特別な意味が生まれるのです。これが「対称
性の自発的破れ」と呼ばれる現代物理学におけるもっとも基本的な概念の
一つです。
より定量的な科学的議論においては、今の議論において、ある特定の
水分子の周りでは次の分子がどのあたりに居るのかを考えたことが極めて
重要です。つまり局所的な密度としての秩序変数を用いて異なる場所にお
ける秩序変数同士の関係、つまり相関により物質系の「相」を特定したわ
けです。ここに秩序変数として局所的なものを考えることの重要性が典型
的に表れています。現在の物理学においては、極めて多くの場合にこのよ
うな局所的秩序変数による相転移理論が広く用いられています。物理的に
はこれらの「相」が変化する際の「相転移」が重要かつ興味深い問題とな
ります。「相転移」とは「相」が違うほど全く異なるもの同士が入れ替わ
るわけで、物質の応用、その機能性を考えたときにも、非常に大きな潜在
的可能性を含む現象が「相転移」であるわけです。この相転移を対称性の
破れの観点からみると、転移点の近傍では対称性の破れを規定する局所的
な規則が必ずしも常に満足されないという形で現れてきます。水蒸気が液
化するぎりぎりの温度ではある水分子の周りに他の水分子が完全にランダ
16
ムにあるわけではないが、水や、氷の場合ほどはっきりした方向の異方性
を持つわけではない、といったこととなるわけです。この状況は相転移の
転移点（臨界点）近傍では秩序変数並びにその相関における局所的揺らぎ
が増大すると表現できます。ここでも秩序変数の局所性の重要さは見て取
れると思います。
もちろん非局所的な概念も現代の物理学においては登場しますが、ど
うやら認識論的にいっても非局所的な概念はわかりにくいらしく、その多
くの場合、非局所的な概念は、ある種の双対変換により局所的なものへ置
き換えて理解されることが多いのです。特に場所と時間ごとに定まる「場
の量」をその作業変数とする局所場の理論の臨界現象の理論における大成
功をその背景とし、局所的秩序変数を用いた「相」の概念および「相転移
理論」は現代物理学における概念的基礎の大きな部分となっています。
それではここで説明した局所変数による「相」の分類は十分なので
しょうか？
実は近年の研究において量子系特に量子多体系においては新しい概念が
どうしても必要であることがあきらかになりつつあります。これが量子液
体相における新しい秩序概念です。これについては別な節でまたご説明し
ましょう。
4
量子液体
気体という「相」においてはその構成粒子のあいだに特定の位置の相関
がなく、固体は粒子同士が規則的に並ぶといった強い位置の相関によって
特徴付けられます。液体はその中間程度の位置の相関を持った相であると
17
いえます。このように古典的な物理としても液体は少々微妙な位置にある
ことにまず注意しましょう。 それでは現代科学がその基礎とする量子力
学的な世界観にたって現在の物性科学を考えてみましょう。通常の生活に
おいては量子力学によらず、ニュートン以来の古典力学が十分正確である
ことに疑いはありません。これは十分高温（室温は物性論的には通常十分
高温です）であれば量子論は古典論に帰着することによります。水蒸気が
氷になるように、量子的な系でも一般的な物質は高温では秩序を持たない
対称性の高い相をとり、温度を下げていくに従って自発的に対称性が破れ
た低温の秩序相へ相転移するのが基本的な振る舞いです。よく知られた例
としては（量子的な）磁性体が高温での非磁性状態から低温の磁性状態へ
転移する磁気転移があります。また量子系固有の現象として特に有名な超
伝導転移もその相転移としての性質は基本的にはここで議論した局所的な
秩序変数を用いた通常の対称性の破れとしてされます。ただし超伝導の場
合の対称性のれははゲージ対称性の破れと呼ばれる極めて特徴的なもので
はあるのですが、ここではあまり強調しません。通常の量子系は高温で対
称性の高い相をとり低温で対称性が低い相へ転移するのが一般的なシナリ
オで、その転移は自発的対称性の破れの概念により局所的秩序変数を用い
た秩序形成の物理として一般的に理解されます。
ところが最近その秩序形成がきわめて起こりにくい系がいろいろな意
味で興味を集めています。量子系においては古典的な系における熱揺らぎ
だけでなく量子系固有の量子揺らぎが存在します。標語的にいえばハイゼ
ンベルグの不確定性に対応する物理量の不確定性起因の揺らぎが存在する
と考えればよいでしょう。特に１次元、２次元系といった低次元系におい
てはその低次元性と量子性があいまって通常の秩序形成が大きく妨げられ
18
ることがしばしば起こります。その典型例が分数量子ホール系におけるラ
フリン状態とよばれる多粒子状態です。電子は電荷を持っていますのでそ
の間にはクーロン斥力が働きます。よって電子集団を考えたときの振る舞
いは電子の運動エネルギーとクーロン斥力の関係により定まることとなり
ます。簡単な考察から電子密度が高いと運動エネルギーがクーロン斥力よ
り主な寄与をすることがわかり、実際に高密度の極限では電子系は電子ガ
スとよばれる気体の相をとります。また電子密度が小さいときは運動エネ
ルギーよりもクーロン斥力が主な寄与をあたえ、電子系はウィグナー結晶
とよばれる固体の相をとることが知られています。これらの特徴的な極限
の密度以外の中間的な電子密度では複雑な相が現れる可能性があることは
容易に想像できるでしょう。
実際、電子を２次元の平面に閉じこめさらに磁場をかけたときには、
この中間の電子密度においてさらにある特定の電子密度においては固体
のような長距離の秩序は持たないが気体としての完全な一様性ももたず、
近距離の密度密度相関のみを持つ「液体」状態が現れることが知られてい
ます。これがラフリン状態とよばれるもので、標語的にはウィグナー結晶
（固体）が量子効果で融解した「量子液体」であるということができます。
この系は液体という少々わかりにくい微妙な相である一方具体的な波動関
数が書き下されているため、多くの研究がなされこの系から非常に多くの
そして驚くべきことを現代の物性科学者は学びました。
たとえば電子は電荷 e を持つフェルミ粒子なのですが、この分数量子
ホール系において活躍する粒子（準粒子）は e/奇数といった半端な電荷
をもちさらにはその統計さえもフェルミ粒子とボーズ粒子の中間のもの
であるというのです。
（この統計は Fermion と Boson の中間という意味で
19
Anyon と呼ばれます）このようにラフリン状態は極めて特異で興味深い系
であることはわかったのですが、その「相」としての特徴は「量子液体」
として固体でも液体でもないというどっちつかずの理解にとどまらざるを
得なかったのです。
このような量子液体状態は例外的なのでしょうか？
振りかえってこの１０年、２０年の物性科学をふりかえると例外どこ
ろか、実は量子液体こそが物性科学の興味の中心であったとすらいえま
す。微妙によくわからないものにこそ多くのミステリーが潜み、そして多
くの知的興味をひいてきたのです。その典型例が高温超伝導体におけるい
わゆる RVB 状態です。これは和訳すれば「共鳴結合ボンド状態」でしょ
うか。２次元反強磁性の秩序を持ったスピン系は強制的にホールをドープ
することにより磁気秩序相としての固体が融解して特異な「量子液体状
態」をつくると考えられたのです。実際の高温超伝導体におけるこの状態
の実現可能性に関しては百花総攬、人々によって意見が全く異なり、軽々
に結論はでないのですが、理論的にはその量子液体状態としての意義の大
きさは明らかです。また古い問題でありかつまた最近興味が再燃しつつあ
るフラストレートした磁性体の問題においても量子液体状態としての量子
状態が重要であることは間違いありません。
5
トポロジカル秩序と量子秩序
考察する量子系のハミルトニアンがある対称操作のもとで不変である
にも関わらず観測される物理量がその対称操作のもとで不変でないとき
（少々不正確ではありますが）自発的に対称性が破れているといいます。
20
異方性をもたないスピン同士の相互作用をもつ磁性体において低温で全体
のスピンが特定の方向を向いたり特定のスピンのパターンをもつ相が実現
する例がもっとも典型的な対称性の自発的破れです。この対称性の破れは
現代物理学において極めて基本的な役割を持っており、多様な物質の多様
な物質相の多くはこの対称性の破れの概念により基本的に区別されてきた
のです。 一方で近年の物性科学の進歩により、必ずしも対称性の破れを
伴わないが極めて特徴的な物質相が存在することが明らかとなりました。
歴史的にみて、その典型例が量子ホール効果における種々の物質相です。
２次元的に閉じこめられたお互いに斥力相互作用しあう電子群が磁場下に
置かれたとき電子密度および磁場の強度をいろいろと変化させると驚くほ
ど多様な物質相が現れることが知られるに至りました。ところがその多く
の物質相の間で古来用いられてきた局所的な秩序変数により記述されるよ
うな対称性の違い、区別はその多くの場合全くなく、対称性的には同一の
物質相として考えざるをえないのですが、明らかに一連の物質相の性質は
特徴的に異なり同一とはどうしてもいいがたいのです。そこで考えられた
のがトポロジカル秩序といわれる新しい概念です。
トポロジカル秩序とは局所秩序変数が局所場の場の理論にその起源を
もつのに対応してトポロジカルな場の理論において使われたいくつかの概
念との類似性を元に提案がなされたものです。たとえば、状態の縮退度が
物理系の大域的な構造どのように依存するか等を相の特定に使おうという
わけです。私は近年より広く、量子系固有の幾何学的位相を用いてトポロ
ジカル秩序の概念より広くは量子秩序と呼ばれるものを利用して古典的対
称性の観点からは特徴を見いだしにくい系の特徴付けを行うための具体的
手法を提案しいくつかの具体的試みを行っております。
21
振り返ってこの２０年来の物性物理を考えるといくつもの新規な物質
相が発見されてきましたが、真に興味深いのは、いわば未だによくわから
ない量子液体相であると考えられます。これら量子液体相においてこそ、
幾何学的秩序をもちいて拡張された新しい秩序概念としてのトポロジカル
秩序および量子秩序が極めて有用なものとなると考えられるのです。これ
については節を変えてまたご説明したいと思います。
6
量子液体の特徴付け
前節にて説明しましたように近年の物性科学の進歩により見いださ
れてきた量子液体とは古典的な秩序変数並びに相分類の手法、概念にうま
くマッチしないものであり、新しい物性科学における概念的進歩を要求す
るものです。それはこの１０年から２０年の間に多く興味をあつめた物性
科学における現象においてこの量子液体相が現れてきたことと無関係では
ないと考えられます。そのいくつかを列挙すれば、高温超伝導体における
RVB 状態、フラックス相等時間反転対称性を破る状態、整数、分数量子
ホール効果、分数統計による超伝導、Haldane 相等整数スピン鎖、梯子形
スピン系、フラストレートしたスピン系、カーボン２次元系とナノチュー
ブにおける磁気秩序、異方的超伝導と時間反転対称性の破れ等があげられ
ます。 量子液体はその名前にもあるように本質的に量子効果をその基盤
とするものです。そう考えると幾何学的位相が量子系の本質的側面を記述
していると見なしたとき、幾何学的位相により量子液体相を記述、特徴付
けようとするのは極めて自然であるとすらいえます。量子的状態とは近年
量子計算等の研究の進展で広く知られるに至りましたがその局所摂動に対
22
する応答ですら局所的ではなく系全体に及びます。よって量子的波動関数
を作業対象とする幾何学的位相を用いた相分類、特徴付けを考えることは
局所場の理論に基づく局所秩序変数に基づくこれまでの相分類、相転移理
論とは本質的に異なるものとなります。
すこし細かくなりますが、幾何学的位相を用いた量子液体の特徴付け
に関する私の研究をすこし紹介しましょう。量子液体では通常の秩序変数
が使えないわけと何度もいいましたが、では何を作業変数として物理を
やればよいのでしょうか？ 私は、この問いに対する答えとしてその作業
変数として幾何学的位相を用いた「トポロジカルもしくは量子的な量」を
用いることを提案し具体的なスキームを具体例とともに提示しています。
ここで物理量といわず「量」といったのは通常の古典的対応物の存在する
物理量はエルミート演算子に対応するわけですが、ここでの「量」は演算
子としての対応物を持たないからです。この理論では、エルミート演算子
ではなく、ベリー位相の議論の際に用いられたベクトルポテンシャルを一
般化したベリー接続と呼ばれるものを基本的作業変数として議論を進めま
す。その際、系の特徴付けを行うためには密度、磁化、等の様な連続量を
使うこともできますが、ここでの理論では「量子化」された「トポロジカ
ルな量」を構成することを目指します。量子化された量をもちいれば相の
分類が明確であることは明らかでしょう。これを用いると古典的秩序変数
とは異なるトポロジカルな秩序変数を構成することもできるのです。量子
液体相では本質的な局所的な秩序変数は存在しないと繰り返してきました
がこの古典的対応物を持たない量を用いることで量子液体相に対するトポ
ロジカルな局所的秩序変数を構成することができるのです。最近その一般
論をつくるとともに重要ないくつかの例であるフラストレートしたスピ
23
ン系、２量体（ダイマー）化した電子系に関してその一般論を適用しその
有効性を示すとともに、より一般的な量子液体相での有効性を探りつつあ
ります。例えば、高温超伝導体等においてはいわゆる Zhang-Rice シング
レットと呼ばれる動き回る電子対が重要な役割を果たすと考えられていま
す。シングレットと呼ぶようにこの対はスピンとしてゼロでありスピンは
運ばないのですが、実はベリー位相としてはπという非自明な値に量子化
したいわばトポロジカルなチャージを運ぶことが示せるのです。
7
終わりに
この場を借りまして少々変わった研究をご紹介させていただきました
が、この話題は、実は固体の電子論におきます大問題であるところの電子
の遍歴性と局在性に深く関わるものであることを指摘し、より詳しいこと
は今後の研究にて、ご報告させていただきたいと考えております。
24
電子構造理論からナノテクノロジーへ
–藤原毅夫先生との共同研究を中心として–
星健夫
鳥取大学工学部応用数理工学科, CREST-JST
平成 19 年 2 月 7 日
恩師である藤原毅夫先生が東京大学大学院工学系研究科物理工学専攻を
退職なさるにあたり、これまでの共同研究についてまとめる。筆者は 2006
年 10 月まで助手として藤原研究室に所属し、現在でも、藤原先生と共同
研究を行っている。[1, 2] 以下ではまず、電子構造理論の発展や藤原研究
室における最近の研究活動を概観する。その後、筆者が主に取り組んでい
る超大規模電子構造計算理論 [3, 4, 5, 6, 7, 8, 9, 10, 11, 12] について解説
する。
1
ナノテクノロジーと電子構造理論
今日、「ナノテクノロジー」という言葉は新聞にも毎日のように現れ、
1nm から数百 nm までの系がとりあげられている。中でも半導体産業の
微細加工スケールは 100nm を下回り、nm スケールから 10nm スケール系
25
図 1: 半導体表面構造と波動関数の例；Si(111)-2× 1 構造 (第 5 章参照) の
第一原理 (カーパリネロ) 計算 (左) と模式図 (右)。5 員環と 7 員環のペア
が [21̄1̄] 方向に２倍周期を与える。7 員環部分には、π 軌道型波動関数 (図
に表示) が存在する。
の重要性が指摘されている。このようなスケールでは、バルク部分と非バ
ルク (表面・界面・欠陥など) 部分の原子数は競合できる程度であり、両
者が共存・競合することで、多彩な構造や物性が現れる。理論的にこうし
た系を扱うには、電子の量子力学的性質、つまり「波」としての性質が必
須であり、電子構造理論がその役割を担っている。半導体シリコン表面に
おける特徴的な波動関数の一例を、図 1 に示した。
一方、現代的な電子構造理論は、60 年代における密度汎関数理論の構
築 [13, 14] に始まる。その究極目標は、多体波動関数理論を出発点として、
凝縮系の諸物性を記述できる定量的理論を、演繹的 (非経験的) に構成す
ることにある。[15] こうした理論は第一原理計算と呼ばれ、計算機の高
速化も伴い、発展してきた。今日では基礎科学のみならず、工学分野にお
いても確固たる地位を築いている。[16, 17] しかし、ナノテクノロジーの
理論基盤となるための道のりは、まだ始まったばかりというのが現状であ
26
る。近年の藤原研究室では、特に、
1. １電子描像と多電子描像の融合による第一原理電子構造理論の拡張
および、強相関電子系 (主に遷移金属酸化物系) への適用 [18, 19, 20,
21, 22, 23, 24, 25, 26]
2. 超大規模電子構造理論および、ナノ構造の動的プロセスと機能性
[3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 27, 28, 29]
の２テーマに取り組んできた。以下では、103 -107 原子系計算 (図 2) を可
能にした、超大規模電子構造理論について述べる。
Cu
enti
10
onal
2
1
10
4)
; fcc
103
conv
Computational time [sec]
104
0
10
fc
c
d
ui
(6
Cu
C
ng
li
ca
)
64
liq
lk
bu
de
i
Si
2)
r-s
ea
(
y
all
lin
lk
1
(5
i
S
bu
10-1
101
102
103
104
105
106
107
108
Number of atoms (N)
図 2: 超大規模電子構造計算の例; [4, 7, 10] 金属系（fcc Cu, 液体 C）、
半導体系（バルク Si）に対する計算時間。固有状態を計算する従来手法
（conventional）の他、いくつかの超大規模計算手法を用いている (第 4 章
を参照)。Intel 社製・SGI 社製 CPU を用いて計算した。並列計算おいて
は、使用した CPU(core) 数を括弧内に記した。
27
2
カーパリネロ法と大規模系への道
今日の標準的な電子構造計算では、一電子 (コーン・シャーム) ハミル
トニアン H に対して固有状態を数値的に計算する。しかし標準的な行列
固有値問題アルゴリズムを用いると、基底数 M に対し計算コストは M 3
に比例して増大し、系が大きくなると実行できなくなる。量子力学の数理
構造の枠内で、計算量を劇的に縮小させる必要がある。1985 年に登場し
たカーパリネロ法 [30] は、固有値計算を回避し、かつ平面波基底計算に
おける高速フーリエ変換 (FFT) の効能を最大限に活用することで、計算
量の劇的縮小をなしとげ、ブレークスルーを切り開いた。カーパリネロ法
は、「第一原理分子動力学法」とも呼ばれる。様々なグループによりプロ
グラムコードが開発されており、有償・無償でダウンロードできるものも
多数存在する。[31] 藤原研究室においても、基礎理論構築や高速化手法研
究 [32, 33, 34, 35, 36] のうえに、独自プログラムコードを開発した。完成
したコードは応用研究 [37, 38, 39, 40, 41, 42] に用いられてきた。
90 年代に入り、上記カーパリネロ法などが 101 -102 原子系に対して確立
されると、大規模系 (103 原子系) での電子構造計算および分子動力学計算
を目指して、さまざまな研究グループにより、種々の理論が提案されるよ
うになった。[43, 44, 45]
1996 年コーン (W. Kohn) は、数千原子系以上の電子構造計算にあるべ
き基礎原理として、一電子固有状態（コーン・シャーム軌道）の代わりに
一体密度行列の短距離成分を用いて密度汎関数（第一原理）理論を構成す
ることを提案した。[43] 一体密度行列 ρ̂ は、占有された一電子固有状態
28
{φk } の重ねあわせで形式的に定義される;
ρ̂ ≡
occ.
X
|φk ihφk |.
(1)
k
物理量 hX̂i は密度行列 ρ̂ を用いて
Z Z
hX̂i = Tr[ρ̂X̂] =
drdr 0 ρ(r, r 0 )X(r 0 , r)
(2)
と書ける。つまり、一電子固有状態 {φk } の一つ一つを計算しなくても、
密度行列 ρ̂ が与えられれば、物理量 hX̂i は得られる。密度行列の対角要素
（ρ(r, r)) は電子密度であるので、非対角要素が量子力学的自由度である。
ここで、演算子 X(r 0 , r) が短距離演算子であるなら、密度行列 ρ(r, r 0 )
の長距離成分は（たとえ存在しても）物理量 hX̂i に全く寄与せず、計算
する必要はない。この事実が、上記理論 [43] の本質となっている。例とし
て、付録に自由電子ガスでの表式を載せた。
一方、大規模系の電子構造計算におけるグリーン関数の利用は、上記の
密度行列の議論とは別に、古くから議論されてきた（例えば、連分数展開
法 [46]）。グリーン関数は逆行列 G = (z − H)−1 で定義されるが、通常の
逆行列計算は、固有値問題と同様に、計算コストのため大行列では実行で
きない。しかし、線形方程式
(z − H)|xj i = |ji
(3)
を基礎におくことで、最新の線形計算アルゴリズムが適用できる (Gij =
hi|xj i)。ここで、z は複素エネルギーであり、(z − H) は非エルミート行
列になっている。密度行列は占有分布 f (ε) を用い、以下で与えられる
Z
1 ∞
Im Gij (ε + i0) f (ε) dε.
(4)
ρij = −
π −∞
29
グリーン関数においても、目的によっては、非対角要素をすべて計算する
必要はない。例えば局所状態密度は、実空間表示における対角要素 (Gii (z))
だけから計算できる。
このように、密度行列にせよ、グリーン関数にせよ、必要な非対角要素
だけ計算できれば良いことが重要である。後で述べる手法は、このような
性質をもつ。
3
実空間表示のハミルトニアンと電子構造の普遍性
超大規模計算では、実空間基底を用いたハミルトニアン（行列）が用い
られる。本記事の超大規模計算では、密度汎関数（コーンシャーム）ハミ
ルトニアンを直接用いる代わりに、スレーターコスター形式にマップした
ハミルトニアンを用いている。問題の本質は、短距離成分のみをもつハミ
ルトニアンを構成できるか、という点にある。半経験的な試みは古くから
あるが、今日では、金属か半導体かに関わらず、こうしたハミルトニアン
の存在が第一原理から導かれる。[47] その根源は、凝縮系における多重原
子散乱により、散乱波が遮蔽されることにある。散乱波解はダイソン型方
程式の厳密解として与えられ、短距離型ハミルトニアンが導かれる。ここ
で、長距離成分を無視しているのではなく、短距離成分しか持たない厳密
解が導出されている点に注意したい。
一般に電子構造計算は断熱近似のもとで分子動力学計算を実現するが、
上述のようなスレーターコスター形式ハミルトニアンを用いた計算が広
く行われている（文献 [44, 7] の文献リストなどを参照せよ）。このよう
なハミルトニアンは、固体・液体・欠陥・表面など幅広く適用できること
30
から、しばしば ‘汎用型 (transferable)’ ハミルトニアンと呼ばれる。こう
した汎用性は、エネルギーと長さをスケールすると、電子構造が普遍で
あることに起因する。この普遍性が、「周期律表を縦に見たり横に見たり
すると、物性に系統的な傾向がある」という物性物理の基礎を与えてい
る。図 3 に、IV 族元素の表面構造における例をあげる。[7] C の (001) 面
では対称ダイマーが形成されるが、Si と Ge では非常によく似た非対称
ダイマー (図 3(a)) が形成される。一方、これら元素においては、金属度
(metallicity) と呼ばれる無次元量 αm が汎用型ハミルトニアンから一意に
決められ、αm を連続的に変化させていくことにより、上記結果が定量的
に再現される (図 3(b))。
図 3: IV 族元素における電子構造の普遍性；(a) Si,Ge の (001) 面における
非対称ダイマー。上 (真空) 側の原子を球で強調してある。(b) 金属度 αm
を用いた C,Si,Ge(001) 表面構造における１パラメーター理論 [7]。
31
4
具体的計算法
密度行列・グリーン関数の計算に関して、我々は以下の４種の方法を開
発した; 一般化ワニア状態を用いた摂動解法 (WS-PT)・変分解法 (WS-
VR)[3, 4]、クリロフ部分空間 [48] を用いた反復型線形計算解法 (部分空
間対角化法 (KR-SD) とシフト型 Conjugate-Orthogonal 共役勾配法 (KR-
sCOCG)) [6, 8, 11]。手法の詳細や比較は、文献 [10, 12] を参照せよ。
図 2 には、C, Si,Cu について、最大 107 原子までの計算例を示した。そ
れぞれの系のハミルトニアンについては、文献を参照せよ。固有状態を
もちいた従来法の他、WS-PT 法 (バルク Si)、KR-SD 法 (fcc Cu, 液体 C)
を示してある。従来法以外の手法はみな、計算量がシステムサイズ (原子
数)N に比例している。こうした性質をもつ手法を、「オーダー N 法」と
総称する。もし密度行列 ρij の非対角要素を全て求めると、このような性
質は決して得られず、短距離成分だけの計算で済むことが本質である。一
方、グリーン関数の計算において、特に微細なスペクトル構造を議論する
場合は、KR-sCOCG 法が有用である。[8, 49] 本手法は、クリロフ部分空
間の数理的研究 [50] に基づいた、汎用な線形計算アルゴリズムである。そ
のため、応用は電子構造計算に限らない。実際文献 [50] では、量子色問題
(QCD) への適用例が論じられている。
また我々は、OpenMP(http://www.openmp.org/) を用いた共有メモリ
型並列計算アルゴリズムに重点をおいており、図 2 における並列計算は、
すべて OpenMP で行われている。実測計算時間の短縮には、論理的な計算
コスト (演算回数) を減らすだけでなく、(1) 並列アルゴリズム、(2) キャッ
シュメモリの効率化、など、計算機科学的視点も重要である。特に、近年
32
のマルチコア CPU の普及により、共有メモリ型並列計算はパーソナルコ
ンピューターでも実現できる。そのため本研究の成果は、幅広い計算機環
境（スーパーコンピュータからパーソナルコンピュータまで）で利用可能
なスケーラビリティを持っている。
5
応用研究：へき開における表面生成プロセス
応用課題として、シリコン劈 (へき) 開における、表面生成プロセスを取
り上げる。[4, 7, 9] 10nm スケール（または 104 − 105 原子）系に対する、
10ps オーダーの動的プロセスに相当する。図 4(a)-(c) では、実験で得られ
る (111) 面または (110) 面へのへき開面の折れ曲がりが、他の面（この例
では (001) 面）からへき開を (人工的に) 開始しても、得られている。一方、
図 4(d)-(e) は、ステップ形成を含んだへき開プロセスを表す。表面で形成
される７員環と 5 員環の対が、実験的にも観測されている (111)-(2 × 1) 構
造（Pandey 構造、図 1）の基本構造である。得られたステップは、ステッ
プ端 (図 4(e) 中央やや右下) に６員環をもつ。また、電子物性として、グ
リーン関数を用いて局所状態密度を計算した。[9, 49] 特に、図 4(d)-(e) の
ステップ構造に対する計算は、バイアス依存の STM 像に対応する結果と
なった。へき開面におけるステップ構造は、STM 実験が存在する (文献 [7]
の文献リストを参照せよ) が、明確な原子構造の特定には至っていない。
我々の結果は理論からの予測であり、実験に対する示唆を行っている。
ナノ構造系の物理としてみると、(i) 弾性論的プロセス（亀裂形状変化
による特異性応力場の変化）および (ii) 電子論的プロセス（バルク状態か
ら表面状態への電子構造の変化）の競合が、現象の本質である。[7]
33
図 4: 超大規模電子構造計算の計算例；シリコンへき開における表面生成
プロセス。[7] (a)-(c) へき開面のベンディング（折れ曲がり）。(d)-(e) ス
テップ形成を含んだへき開プロセス; A-E は初期（結晶）構造における結
合。(d) から (e) への過程 (0.6 ps 程度) で、結合 E が破断している。関連
動画ファイルは公開されている (http://fujimac.t.u-tokyo.ac.jp/lses/)。
34
6
今後の展望
超大規模電子構造計算は、(1) 物性物理学（電子構造理論）
・(2) 数理科
学（大行列数値線形解法）・(3) 計算機科学（例えば、共有メモリ型並列
計算）の接点となっている。理論は量子力学の数理に基づく汎用な理論で
あるが、実際の応用研究はまだごく少数である。今後は、計算コードを整
備・発展させながら、応用研究を積み重ねていく必要がある。産学連携プ
ロジェクト [2] などを通じて、計算に対する需要（ニーズ）も学んでいき
たいと考えている。
ナノ構造系研究における横断的な視点として、(a) 形成プロセス (制御可
能性など)、(b) 機能性 (電気伝導など、電子系が直接関与する性質)、(c)
信頼性 (壊れにくさ)、(d) 元素を取り替えた系統的研究、があげられる。
いかなる計算においても、精度と計算コストは相反する要求である。様々
な手法を総合的・階層的に用いることで、目的にあわせてシミュレーター
をデザインすることが重要であろう。
謝辞
本記事における研究は、藤原研究室および力学教室（東京大学大学院工
学系研究科物理工学専攻理論グループ）のメンバー・元メンバーとの共同
研究である。何人かの方には、本文献リスト作成にあたって、相談した。
35
付録：自由電子ガスでの一体密度行列
自由電子ガスの固有状態は平面波である。密度行列は、２点間の距離だ
けで書かれ (ρ(r1 , r2 ) = ρdist (|r1 − r2 |))、平面波解を波数空間で積分する
ことで計算できる。具体的には、Fermi 波数 kF をもちいて
Z
dk ik·r
ρdist (r) ≡
e
3
k<kF (2π)
½
¾
2
kF
1
=
− 2 cos kF r + 3 sin kF r .
(2π)2
r
r
(5)
で与えられる。第１項は長距離成分 (フリーデル振動) である。[51] 一方、
短距離成分は、式 (5) に現れる三角関数をテイラー展開することで与え
られ、
2
ρdist (r) =
(2π)2
µ
¶
C2 2
4
C0 −
r + O(r )
2
(6)
となる。ここで、
C0 ≡
kF3
,
6
C2 ≡
kF5
.
15
(7)
は、ゼロ次と２次のテイラー係数に相当する。エネルギーは、ハミルトニ
アン H ≡ −∆/2 との対角和で計算でき、
¯
Z
¯
1
1
−∆r1
E
≡
Tr[ρH] =
dr1
ρ(r1 , r2 )¯¯
V
V
V
2
r1 =r2
= lim
ε→0
Esphere (ε)
(4πε3 /3)
(8)
と書ける。ここで、Esphere (ε) は、実空間での半径 ε の微小球におけるエ
ネルギーであり、式 (6) より、
Z
−∆r
C2 3
Esphere (ε) ≡
dr
ρdist (r) =
ε + O(ε5 ),
2
π
r<ε
36
(9)
で与えられる。以上より、単位体積あたりのエネルギーは、
Esphere (ε)
E
3C2
.
= lim
=
3
ε→0 (4πε /3)
V
4π 2
(10)
となり、エネルギーは２次のテイラー係数 C2 だけで決まっている。言い
換えれば、密度行列（式 (5)) は非対角長距離成分をもっているにも関わ
らず、エネルギー E はその短距離成分だけで決まっている。
一方、固体物理学の教科書などでは、このエネルギーは k 空間積分で計
算され
E
=
V
Z
k<kF
kF5
dk 1 2
k
=
.
(2π)3 2
20π 2
(11)
となる。C2 の定義式 (7) を用いれば、式 (11) と式 (10) は一致しているこ
とが分かる。
参考文献
[1] JST-CREST プロジェクト「複合手法を用いた電子構造計算技術の
開発」, 代表：藤原毅夫（東京大学）, 平成 16 年度-21 年度.
[2] JST シーズイノベーションプロジェクト「100nm スケール材料の量
子力学的分子動力学シミュレータの開発」, 代表：大谷泰昭（みずほ
情報総研）, 平成 18 年度-19 年度.
[3] T. Hoshi and T. Fujiwara, ‘Theory of composite-band Wannier
states and order-N electronic-structure calculations’, J. Phys. Soc.
Jpn. 69, 3773 (2000); ‘Structure of exact Wannier states of
tetrahedral-bonded crystals’ Surf. Sci. 493, 659 (2001).
37
[4] T. Hoshi and T. Fujiwara, ‘Dynamical brittle fractures of nanocrystalline silicon using large-scale electronic structure calculations’, J.
Phys. Soc. Jpn. 72, 2429 (2003).
[5] M. Geshi, T. Hoshi and T. Fujiwara, ‘Million-atom molecular dynamics simulation by order-N electronic structure theory and parallel computation’, J. Phys. Soc. Jpn. 69, 2880 (2003).
[6] R. Takayama, T. Hoshi, and T. Fujiwara, ‘Krylov subspace method
for molecular dynamics simulation based on large-scale electronic
structure theory’, J. Phys. Soc. Jpn. 73, 1519 (2004).
[7] T. Hoshi, Y. Iguchi and T. Fujiwara, ‘Nanoscale structures formed
in silicon cleavage studied with large-scale electronic structure calculations: Surface reconstruction, steps, and bending’, Phys. Rev.
B72, 075323 (2005).
[8] R. Takayama, T. Hoshi, T. Sogabe, S.-L. Zhang, and T. Fujiwara,
‘Linear algebraic calculation of the Green’s function for large-scale
electronic structure theory’, Phys. Rev. B 73, 165108 (2006).
[9] T. Hoshi, R. Takayama, Y. Iguchi and T. Fujiwara, ‘Large-scale
electronic-structure theory and nanoscale defects formed in cleavage
process of silicon’, Physica B 376-377, 975 (2006).
[10] T. Hoshi and T. Fujiwara, ‘Large-scale electronic structure theory
for simulating nanostructure process’, J. Phys.: Condens. Matter
18, 10787 (2006).
38
[11] T. Sogabe, T. Hoshi, S.-L. Zhang, and T. Fujiwara, ‘A numerical
method for calculating the Green’s function arising from electronic
structure theory’, Preprint (math.NA/0602652); to appear in Proceedings of the International Symposium on Frontiers of Computational Science 2005, Springer Verlag.
[12] 星健夫, 藤原毅夫, 「超大規模電子構造計算と 10nm スケール系の物
理」日本物理学会誌第 61 巻第 4 号, 256 (2006).
[13] P. Hohenberg and W. Kohn, ‘Inhomogeneous electron gas’, Phys.
Rev.136 B864 (1964).
[14] W. Kohn and L. S. Sham, ‘Self-consistent equations including exchange and correlation effects’, Phys. Rev.140, A1133 (1965).
[15] W. Kohn, ‘Nobel Lecture: Electronic structure of matter ―wave
functions and density functionals’, Rev. Mod. Phys. 71, 1253
(1999).
[16] 藤原毅夫, 「固体電子構造–物質設計の基礎–」, 朝倉書店 (1999).
[17] 固体物理 Vol.39 No.11, 「特集：計算機ナノマテリアルデザイン」赤
井久純・押山淳・小口多美夫・笠井秀明・常行真司・藤原毅夫・吉田
博編, アグネ社 (2004).
[18] A. Yamasaki and T. Fujiwara, ‘Electronic Structure of oxides MO
(M=Mg, Ca, Ti, V) in the GW approximation’, Phys. Rev. B 66,
245108 (2002).
39
[19] S. Yamamoto and T. Fujiwara, ‘Charge and spin order in RNiO3
(R=Nd, Y) by LSDA+U method’, J. Phys. Soc. Jpn. 71, 1226
(2002).
[20] A. Yamasaki and T. Fujiwara, ‘Electronic structure of transition
metal Fe, Ni and Cu in the GW approximation’, J. Phys. Soc. Jpn.
72, pp.607-610 (2003).
[21] S. Kobayashi, A. Yamasaki and T. Fujiwara, ‘Electronic structure
and dielectric properties of cubic zirconia’, Jpn. J. Appl. Phys. 42,
6946 (2003).
[22] T. Fujiwara, S. Yamamoto and Y. Ishii, ‘Generalization of the iterative perturbation theory and metal-insulator transition in multiorbital Hubbard bands’, J. Phys. Soc. Jpn. 72, 777 (2003).
[23] O. Miura and T. Fujiwara, ‘k-dependent spectrum and optical conductivity near metal-insulator transition in multi-orbital Hubbard
bands’, J. Phys. Soc. Jpn. 75 , 014703 (2006).
[24] Y. Nohara, A. Yamasaki, S. Kobayashi, T. Fujiwara, ‘electronic
structure of antiferromagnetic LaMnO3 and the effects of charge
polarization’, Phys. Rev. B 74, 064417 (2006).
[25] 藤原毅夫, 「多電子論に基づく LDA を超える第一原理計算手法の基
礎と開発」 , 固体物理 Vol.39 No.11, 「特集：計算機ナノマテリアル
デザイン」赤井久純・押山淳・小口多美夫・笠井秀明・常行真司・藤
原毅夫・吉田博編, 707-708 頁, アグネ社 (2004).
40
[26] 藤原毅夫, 「第一原理計算における動的平均場近似法 (DMFT) の展
開と応用」 , 固体物理 Vol.39 No.11, 「特集：計算機ナノマテリアル
デザイン」赤井久純・押山淳・小口多美夫・笠井秀明・常行真司・藤
原毅夫・吉田博編, 715-721 頁, アグネ社 (2004).
[27] Y. Iguchi, T. Hoshi and T. Fujiwara, ‘Two-stage formation model
and helicity of gold nanowires ’, Preprint (cond-mat/0611738).
[28] 品岡寛，星健夫, 藤原毅夫, 「ナノメートルスケール金属ワイヤにおけ
る非平衡電流の構造依存性」, 日本物理学会 2006 年秋季大会 26pYH-5;
H. Shinaoka, T. Hoshi and T. Fujiwara, in preparation.
[29] 藤原毅夫, 「大規模電子構造計算法の開発―マルチスケール材料物性
学と亀裂伝搬過程」 , 「知識・構造化ミッション―大学は表現する」
松本洋一郎・小宮山宏監修、藤原毅夫・丸山茂夫・伊東乾編 , 69-82
頁, 日経 BP 社 (２００５).
[30] R. Car and M. Parrinello, ‘Unified approach for molecular dynamics
and density-functional theory’, Phys. Rev. Lett. 55, 2471 (1985).
[31] 第一原理計算コード「ABINIT」(http://www.abinit.org) に関する
日本語解説として; 三上昌義, 「第一原理バンド計算プログラム
『ABINIT』」, 固体物理 Vol.38 No.3, 219-224 頁, アグネ社 (2003).
国内における計算コード開発および応用研究についての解説として
は、[17] の記事を参照せよ。
41
[32] T. Tsuge, M. Arai and T. Fujiwara, ‘Atomic-structure of Si(111)√
√
3 × 3 R30◦ studied by first principle molecular dynamics’, Jpn.
J. Appl. Phys. 30, L1583-L1585, (1991).
[33] Y. Yamamoto and T. Fujiwara, ‘Numerical stabilization of the firstprinciples molecular-dynamics method for metals’, Phys. Rev. B46,
13596 (1992).
[34] M. Mikami and T. Fujiwara, ‘Transferability of norm-conserving
pseudopotentials: Si, Al and alkali metal elements embedded in a
jellium’, J. Phys. Soc. Jpn. 64, 1629 (1995).
[35] T. Fujiwara and T. Hoshi, ‘Optimized ultrasoft pseudopotentials’,
J. Phys. Soc. Jpn. 66, 1723 (1997).
[36] S. Yamamoto, ‘Effective implementations of multi-dimensional
radix-2 FFT’, Comp. Phys. Comm. 125, 1 (2000).
[37] M. Miyata, T. Fujiwara, S. Yamamoto and T. Hoshi, ‘Dynamical
electron-ion coupling in the ionic conductor α-NaSn’, Phys. Rev.
B60, R2135, (1999).
[38] Y. Ishii and T. Takanaga, ‘Atomic and electronic structures and
dynamics in liquid alloys near eutectic point’, J. Phys. Soc. Jpn 69
3334 (2000).
[39] M. Miyata and T. Fujiwara, ‘Ab initio calculation of Peierls stress
in silicon’, Phys. Rev. B63, 045206, (2001).
42
[40] K. Nozawa and Y. Ishii, ‘First-principles structural optimization of
cubic approximant Cd6 Ca’, Mat. Res. Soc. Symp. Proc. 805, 47
(2004).
[41] Y. Iguchi, T. Fujiwara, A. Hida and K. Maeda, ‘The electronic
structure around As antisite near (110) surface of GaAs’, Phys.
Rev. B71, 125328 (2005).
[42] K. Nozawa and Y. Ishii, ‘First-principles studies on orientational
ordering of atomic clusters in Cd6 Ca’, Philos. Mag. 86, 615 (2006).
[43] W. Kohn, ‘Density functional and density matrix method scaling
linearly with the number of atoms’, Phys. Rev. Lett. 76 (1996)
3168.
[44] レビューとして；G. Galli, ‘Large-scale electronic structure calculations using linear scaling methods’, Phys. Status Solidi B217 231
(2000); S. Y. Wu and C. S. Jayamthi, ‘Order-N methodologies and
their applications’, Phys. Rep. 358 1(2002). 文献 [10, 12] の文献リ
ストも参照せよ。
[45] 筆者らによる、実空間差分を中心とした手法研究として；T. Hoshi, M.
Arai and T. Fujiwara, ‘Density-functional molecular dynamics with
real-space finite difference’, Phys. Rev.B52, R5471 (1995); T. Hoshi
and T. Fujiwara, ‘Fully-selfconsistent electronic-structure calculation using nonorthogonal localized orbitals within a finite-difference
43
real-space scheme and ultrasoft pseudopotential’, J. Phys. Soc. Jpn.
66,3710(1997).
[46] R. Haydock, in Solid State Physics 35, ed. H. Ehrenreich, F. Seitz
and D. Turnbull, Academic Press, New York, 216 (1980).
[47] 下記論文やその文献リストを参照せよ; O. K. Andersen and
O. Jepsen, ‘Explicit, first-principles tight-binding theory’, Phys.
Rev. Lett. 53, 2571 (1984); O. K. Andersen and T. Saha-Dasgupta,
‘Muffin-tin orbitals of arbitrary order’, Phys. Rev. B62, R16219
(2000).
[48] クリロフ (Krylov) 部分空間は大規模線形計算分野では確立された概
念である。レビュー記事として、張紹良, 応用数理 8, 301 (1998)。
[49] KR-sCOCG 法 [8] を用いたグリーン関数計算コードは、下記プロジェ
クトの成果物として、公開されている (http://act.jst.go.jp); JST-
ACT プロジェクト「第一原理計算によるハイブリッド分子動力学シ
ミュレーション」, 代表：蕪木英雄 (日本原子力研究所), 平成 13-15
年度.
[50] A. Frommer, ‘BiCGStab(l) for Families of Shifted Linear Systems’,
Computing 70, 87 (2003).
[51] 式 (5) は、絶対零度における表式であるが、有限温度の表式をもち
いると長距離成分は指数関数的に減衰する; S. Goedecker, ‘Decay
properties of the finite-temperature density matrix in metals’, Phys.
Rev. B58, 3501 (1998).
44
!'#
%,
&$,)
"+
0 (*
$Fs¬&š¢UtB>h%™o^pUtWJ%`P
%‘A2007 › 2 V^§%($PjUt%“e[ƒ%‹9
X«Tc!Œr/?€
.–X!&gŠB:
%_8ª367524Q„pUtkY‚Zb©¡~•!LŽ
%N¤1)y#I¬1Ž†%!1M'
9q$di1¥Œ)†
v=. Cd Sp
Ut%“e[ƒ%WJ)–X!%L—WJ
.”KŠB‰
x$
%Hœ o$…$F\ %’a&&#&ˆ£
.'%¥Œ&–X|@$R{(,/W
JŸ¦
0
)
-Hœ o$D($&+­Œ!].
¨l1E.!!)$<›$0.–X%f
˜$Cmzw
pUt)UtB#‡u}!nG#[ƒ“e[ƒ*
“ež}$"%+$;O.%!¥Œ$&pUt
45
±l0
KÂv®}/(!ž&kb$0
n‚&ˆ
jŠ%/©|q›&kb!'fjŠ%Ÿ/¢QmÄ<J;
q›izÃ!/¸2¥d
m-
1)
0/‚y1990 ¯u-ª
%,.–Ǧ%$0 Al-TMTM™NgœhˆjŠ&
fjŠ%Ÿ/©|ʦ$kb%,/"ˆjŠ*Z)QÉ&g
œYRt´!'³·¦%>5HCˆLf&©|Ž %_6E9@
•/"0kb'ˆjŠ2=FDI8
Gtg"/s¹
2)
%W$wc2À/""*%ˆjŠ&k
b%/²ˆ¦$Äʦƒº2WÅ
»¬…ˆjŠkb&nŽ2”Ä/*&!/"«€%˜¾V
&¹%&µ¿&¶M^2¨/*&!*/"s
»r!'
†Â&U—"'O~§!fjŠ&©|q›&¢QmÄiz2‹
T$1©|q›iz%[ˆjŠ+&fjŠ&©
|q›©|´’e…]q&Œ`Ê%
/"%
'&U—
3)
%
¢QmÄiz"'SˆjŠkb%{/ƒº%'
#3$*&/
#3$"1/&""2
¨x‰%fjŠ
"*&%
X£%—½p
&¢QmÄ©|q›iz%Á-0/&ƒº2‹T/
&p%ˆjŠ&©|q›&kb%0-&¹º#&,%{10
/2Cd hˆjŠ&fjŠš!/ Cd6M ŹŠ%Ÿ/i
z2È%
„\¦$jŠ&©|Ž %Ÿ
'?I9B&¤Ä“.Å&!
¼o%¡$jŠ!*?GH4JP­&°‘A7:H2Æ|‘"
©|Ž 2a)/"‡Â/„\’&$ˆjŠ%'&¹º
46
k"+#œldŒ$Yj#c§¥¶!"0º&vZx
$dŒ.„/u²^;7=$:E<?$œldŒ4[S0
$,H$®¯
2!2vZx$¡“4t¨#¤Jœ
ldŒ4 0#%"/";7=$:E<?¥¶!"0
!,$-":E<?$œldŒ$[S%1˜e31
"H®vsQ™"dŒ4ƒ0®¯4w«0!
vsQdŒ!-m0sQ€4,dŒ4ƒ0!u
·14TmZxapproximant crystal!_5
0$Tm
Zx#Ž–y$CFAdŒYj$r¯4šµ0!
œl
{$z­4 0!u·0
vZx$dŒ%š#]l4£’fn\Ol$‘°!¸N
0}#n\‚®Ol$Ol›#]l4£’n\
$vsQOl .10´l4oqg$‚ 20 °‹I ‹v
Zx
%n\¹®Ol#]l4£’0!#"0$±p}
.67@6$i¢"!2%¸N
-
}$©®b
#I1•†vsQOl—R10n\UPqU
P4¦1#~•"®b%•`UP!_&10I ‹vZx
%qUP•`UP1,UPn\
!"0fn\Ol$dŒ
4¬))]l4w«G0!#-/¸‰™"vsQdŒ
.$«W4¦!
0$!«G$qUP¨41
Ol|ž4³0!.D9BF«G•`UP¨4D87>
F«G!-(fn\Ol#5‘«W4L0!!ŸaM$Š
h!"0qUP®b!%H—$V4,‘°
„*0!qUP
fn\Ol4
$]l$ª'®«K£’ .10
$V4š#”…0!#-/sQ™"]l$ª'TmZ
x4 0!u·0qg$X
,\‡ˆ43#«
47
GyN;P!H;PZk"%$#7H
;P`NL42F&UT$op"j%$/XeC[
R 1NL42F7H;PC[&UT$DI=btp
7H;PL4e#">FL#6OeC[N;
PiC[$7H;PVrC[0W"N
;PC[&vX$
$
7H;P]$+,
*C[:EN;PfFC[lc.a&s"%$
4^%$bt8
&9hen!#7H
;P7Hg&Q!#_
`-qL4C[&UT$
Mw$u@(')`-q9:E&1YN;PS
J&5z$
$
?KN;P]$7H;P1/1 7H;PB%$`
-q 150 Adg>F&3xoP\<m$x
48
¼‡HP
‡
\zd‡
t]b€¼‡,
!Ÿ€
d‡% Al-Cu-Fe !´À›•¯*-,¡*-
,
%ˆ…$)r}hQy
_
À*ŽŒ-,.
!
Zn-Mg-REREXªËh”!̼‡
ZŸI¾
–¯*-
,-*!
-'ŸI¾ 100 jHˆ
&¸wo–/{]—TKt¹,
{]—TKt¹
–
o
h”/V
)¸wo–/
o–N*¯*-iy°¢/'¨yo
cx-,‰t',Åͦ®-o–Á`œ
¨yo–cx/qB;E
š«‹/e¬,'„Ä,1/1
Ǽ‡)+JO+\z©/ˆ 2/1 \zd‡'f¶
,
ŸI¾V%-,iy
Š 500 jHˆ+­|0
³
Š’ÂÆ/'iy42<V%-
,‰t™g;8
/'žJiŨyo–cxq.-,!%+
ƒY¦d‡£
¨yo–cx!>G=ÅÍRÇ-
,žJiÅcxk"-,½!h”
J
L¤/‘ÉyÊS
‚Ëiy²sHP
W“¼¥~5DF;1G3¼
¥~¨yo–(¹‹/cx,!¨y
T
• lv Ã
m M / §«  m A : G5 C E
ƒY¦A:G5CE£
a µ 5DF;1G3¼¥~
7:9?,
N/¯,
7:9?
!¿©±
UŠÅÍ,[†¿©\zLocal Density Approximation, LDA4)
)+¨y¿©º·*ÂmA:G5CEcx-,.1 1
u\!¨y¿©º·
^T»K
mM/¿©n°
ÂmA
apÈÂmA:G5C
}@6#˜
49
@M<IK%Z
%<JLA6M8̯…%]/´ƒÑµs
22\¬ÑµÈÆ!¸
1!1>KD;M<=@
MB#vÏ21.$´ƒ‘¤s21LDA &uŒ
%aƒ¬•-aƒ’·=F:BK#"%h­‘¤%–‡$
&¬
Ü°#֞5×21!c+/2
1À¹¢%7C
K99I?E!
Þj=F:BK%֞-šWoŸ‚[Ç$¥
Ä21mœf´ƒr%i‹$&Ô¨1!,m½†2
0Y}$§YxÛ!
,%±ÚÂS$&ªT#2'#/
#
%ŠÍ$
&´ƒ%¾·f•5±¶#h­f•
²^
1!$.0]5»1%h­f•$.Ù#ÌÍ €
1™qHD6M@6MkºLinear Muffin-Tin Orbitals, LMTOÍ
&YY%xƒ$ŔHD6M@6Mkº5h­!Ú
1.5) Cu xƒ% 3d kº!
.#ÇÛ°#]ˆØU
1!!,$_xƒ$ s, p, d kº%h­5Ú
!, 9N
zN &©RÎ$g(21xƒ•%h­f•
1%
100 zQŽ%xƒ5g*l¦©RÎ%uŒ
$‰0N!
1YÌ
¾·f•5Äy
,Ãd°ØU
Y%xƒ54$ÊR
!%7CK9Ê[-V#1uŒ~%œ£°#O¬–#"5—
w˜,1!¼
!
t³,1
h­f•$ÉÓ¾5Ú
1ÌÍ
“
ÝX
&xƒb%n%G@M<IK
%¾·f•5—Äy1+$&²^$Ú
1É
K›7CK9%sÍ®P22YÁ[2ѵ|¿n„Generalized Gradient Approximation, GGA%ҍ Ú
/21!¡
2,pe%
_³ %´ƒÑµ%Ë%)/Õ{G@M<IK5»1!
TÐ &LDA %`
«!.3
50
À´o´ŠŒ¹%},q%.!+/
Lªyªym”)fyUcp0]Í-fyO_
°Pªy†œ0¡|rg-¸Ä
4;@0°PÊH
fyO97B
·R¦'+Lªy†œ
QZ-[97B4;@ T-
¼½+.
-Lªy†œQZ-*
u+.\[97B
4;@0žIfÇ\97B4;@*".0Ã~½0žIf
Ç\97B4;@½
s_~½¶¨*,oŽ«
0(jšªy†œ0},C-\97B4;@ÆÃ-*
,G`—(W&Âa#§ÃM³.6) %X
¥VŠºÀ´0Ãbv
o–8?2»TÁk
+
'ŸF¾Š 100 ify-*
a(bvM
³-.7) +ºÀ´0ü½(
I¢ƒ
fy
+Ë0Ž«*bv-.0wD¤fyµ 0bv*,^'-)ÁhJ«¹y­ËS4:<A4=
BM³-
Cd 2000 ²¶d.D¤ Cd-MM=Ca, Yb‚c„
8)
‚c
„’Œ_¿È¼„_zc„ Cd6M ˜t-. ȼ„
Yb EN±L`—)Y¬Ìe‘(˜tm”
Ca,
™ 1/1 _z
c„®Â 20 Àš›…‹fy42@K“3>
56ŸF¾£©š‡µ .m”d-€Å‰3>56s(°•
 20 Àš›…‹0(
ˆxÀš3>56¼l0l˜t-
n+.-a
D¤É{±e‚c„0¯-$1
51
€@!I)Zƒ!\L^QeNhyGBpUd/;
12!Qj(>)'
100[K]vHQjio?Dk*
9)
1/1 HWL^
!…
"
&(nz
F}aI)
S1 (Cd 4)
S2 (Cd 20)
S4 (Cd 30)
S3 (M12)
S5 (glue Cd)
c 2Cd6M †{^!f 20 l/;12OU0,<
f 20 l/;12!bwfTl/;12;=39
P)!VY%pUQjJS
u
!&;
=3951+rt)"[ƒfl/;12
!b
"K#{P+PfTl/;12+!p
UQj+JS)LMTO |
C
10)
+c
&(s'*pU`m~q(DOS)!L
XM=Ca !LC+M)6,<8-5<.!]
_!`m~q
g„R$E.:47M'*!_A
"`m~q]n)pd ex!wx`m~qpDOS, d-DOS+$)!‚U"%]"()6,<8-
52
27/= p-DOS l($.27/
0 [Ryd]vE
d-DOS 40J
(& p-DOSd-DOS
l(ƒ;
)()
$<yax DOS $f?J&
Cd 5p _k Ca 3d 381n_kUaI
S,W'9m>I@J([(<yM=Mg ^S
5-76.27/vE DOS €T"
J&)
^S d-DOS ,J(|&%5-76.27/vE d
ax_k
)(Z
iVUaD#Q
&
Cd-Mg
dN~HCd6M %FhC>Sw
ba)$$LX_k5-76.27/EG
d _kiV Mg Ua%(9m>r&)s
:Pg,$>Sweba)*oj&Rx
AtOpXxOZK
11)
%(5-76.27/EzpX_k
d pXaxBM+)'LX_k
uc} d _k
Ua Cd-M HEYI]%!\I]D(B
#&)(
` 3‚y] Cd6MM=Ca, MgpX_k{q
53
† m‡ ( "h }m ‡ #= F C J7 H "` r! )+ a
4E9@j-,s*-,
/Zz!),¤z"
W”!)+©Ž k ©Ž k + GG #cZzA5:I"©]
ˆ
h2 2 h2
2
k =
k + G 0³©ŽA5:I".G
2m
2m
" 2 ¥­µ3;I44E9@X†m‡#¡Š G
=3 Å-1 °+!(211111), (221001) {Ž«-,e?G96ª
,‘ 20 µ›œ‰
uœ‰"&!-*" 2 ¥
­µx*-,šµ›a?HI1KPd!hj
IBd"šµ›!“,)
,>2
ž-$>2IBµ"–
Pa4E9@Xf„3;I40§,,
")
•´#½"†m‡0§,&"…¹
<GD8q
z+"®g"T¤zŽ,
|€%v
&
ž¯•§¾",bÀ!o,
*T¤zŽ(r™–
¦ Cd6M ¼±‡!Cd Qv¬0x,p—"ƒ¿!)a
4E9@"r™!ž
£
O,
|€#r™
_S[
‹n(T¤zŽMY"Q[£
¨…¹
¶\0U,
0~,Cd-M k!#M p—²i d Œ
0'…¹" d Œ Cd " p Œ"^"w!)+
a4E9@k"f„’,",")
!¤
zr™ly# Cd6Mg ")
mU"V‚!ŸN
R˜£
-$
*
m‡!œ'€t,"
mU0
%º
¸-$p—!),Q[£¨"O
0»V,"!·¸
,
pd w!)+a4E9@’kL¢Q,,
54
O3;49%C%u.oy&ZcK%X_$N‘
-
"]-"!-[%coy…~%ZcK$SX_0
s*’c&™cDJ%Š!&X_drbond order"™!Q
œ . - % E t 0 Th V %  c o y% n _ $ H } + % COHPCrystal Orbital Hamilton Population%Et!-12)WaŽ
%"/!j(,$LMTO Ž!&Zc$ˆqMLr0“
-%!%,#Et$"
+Œ–!-lf0"
LMTO Ž$,-Ac†ƒLr&
i = c i,RL RL
RL
%,$Zc?{ R "IBƒ™ L !U‹.-MLr RL !€F
.-%oy%27<3\z0 i "k"™cDJ$
-X_dr$w‚-+%&
p RL , R 'L ' ( ) = ci*, RL ci , R 'L ' ( i ) i
"~P.-%™$
hRL,R ' L ' ( ) = pRL,R' L ' () RL H R' L' %,$8:<561=%^›”v!g)0™ COHP !-
hRL,R’L’ ()%z‰!.'%27<3š@%oyu." R " R’%?{$-Zc„bSX_0`-"Y#"i•
-
p
$ Cd6Ca ˜m$-GZcK% COHP 0e.13) !&
Cd-Cd " Cd-Ca %R—.. 3.23.6 [Å] >C%Zcx$
COHP 0Wa
-p|% S2-S4 "‡e&p$e
55
c 20 |i5?78$l 2 62@Cd PY!l 4 62@Cd P
Y$PYj#j[(3)$f1MX*$
.S2-S4, S4-
S5 $ COHP ;2@=3:@41J#wS1x,c#zC
.!,;2@=3:@4$B$`k% Cd PY1N
'v.%,1;2@=3:@4$_$`k%H# Cd
E$NT1]).%,1.!0.$+
"`k1
//NTsNT`k!R(3:@44>9<1J#
NTn"`k,sNTn"`k#oY`k$bDzC.~Y
%rpi$+
"I}NTb%-L
&&O,/.
,#NT`ke}/sNT`k%te}
.!
!,Cd PYE$NT1V*K.oY`k\Q.!
.A{Ca PY,".l 3 62@!h$62@_$ Cd PY$
E$ COHP %G4>9<$_B
/*x$m1!-NTn
.!1Z.$!,Cd-M L
% Cd $ p Fq! M
$ d Fq,NTFqW,/$Fq$NTsNTy€#+G4>9<d.!
ugZU/.
a 4{^ Cd6Ca $oY`k#. COHP
56
13)
)!DKX 20 u^-4/0SXLu^-4/
0 3 <Nys^_6t (111) sF+F%kF
rBHZ+:e* Gomez Lidin14) %(Cd6M
ysQX 20 u^-4/0SoqEM6`w)(
pJ(rBHZT(a)%XLu^ 2 c'
kF+P(Iys^bg 1/2 >xYv)(;jV
]T(b)%XLu^=bgqz)>
x 1/3 Yv)(;jV(!ysQ 100 [K]nCli$
h@fG7R?\)'7RO84513
XLu^-4/0kFd8A[9(%(#
I&)
(.9) -4/0SoqEM6`w)+h
MHZDK%'a
"((.15)
T 5ysQ Cd6M X 20 u^-4/0S(mWUYvx.
,2 14)
57
(™(b)#œÊÑ 1/3 #s~:5A"³·"s~3¾¬y
+7CM8#¯
w 3f*0#-!m|"$ÁȽb±#
§SsÏdG?N;IMÅÎÌ0#nY™#x#™#
-!w y³!¾¬º.1#w $ 3 [†(111)Ä
v3v
š}È£9K<=3 2 [†(001)Äv!
‹/"
[´Âk3Xw "!/2 [†#&Ä[´#v
"-®/#DL7;JN¢z0ƒ"™"„
9K<=—#w 3c¾¬¨PƼ#'#s~"$
.Ó3f*1"s~3¸ž£O°!
w 3f*q­!w _Ôp.1#$9K<=—#}È£
#\¡# S2 ;6M#nY3™#U"„™.
$(/$/»)‰0$ŒÍ!
}È£9K<=#
Âk$Z’š}È£"j
k!/}È£9K<=\
¡# S2 ;6M#Ÿu#v+a* Zn-Sc Òđ"¤0yj#w
Z›#nY
16) ǎ!
."}È£9K<=#s~#j
"0 S2 ;6M#s~$¦\¡"R¸i«"ª
s~`
hг·"_Ô1
0Ë~20#-!s~P¬#Â
W$n‘w Z›¦!V¶T~3‚1:5A"¦
t1
‡rnYǎ#!
+#210m|.º.1
0#-!“Ãw Z›"E4BD>91É"Ò#
1%m|30ˆ$¦Â",/
#00(
µ~”¥…˜+#-!w _Ô"¦Â"À^0Al-TM lj
n‘#H@Mw "§SsÏ"bw _Ô3€l$
d8I>F$!–#7CM88I>F3+µ~”¥
º.1#og
17)
+/1¿¹£²!n‘#©{"S#
Š]/!1%¦Â"Qe•
58
p 6t 20 £~/;24`‹\z(a) S1 1.=Wc•„–o
Nw¤0-7’oN˜w¤0-7+™(b) S1S2 1.=
\z–oN\zH«ZWcA„’oNH«x
P\zWcA„+™
hUjQdUjƒA 100 X@mWcJ!('›aR
|ID^œ'$›aUj‚CW¨c\zTb
iL LMTO Ÿ$'VO
,90 ‘
f*&_
%£”Kˆ‚CW¨M86>1:=Ÿ"§¥(
'$—Œ$&cn€¢Ž#WcIU
^rGv LMTO Ÿ+\z`‹D#B'\zIyŠ
?‡rVO‚CW¨M86>1:=Ÿ+$(
(eŸ9<57+u
š[VOE“¡]
Cd ShUjQdUjy' Cd6M ©žj'Tbª+kF
)‚CW¨TblQdUj} 1/1 Qd
UjY('ƒA …Wcq 150 X‰ŽS'03hUj†+{'¡g¦')
!˜s
59
BDPS8:$#'AH
*+&)(%?7@M<I
!>3I7@K>3I?7@"5
L/2MING0J
F.8R69=1"Q
;,!
1) T.Fujiwara, Physical Properties of Quasicrystals, edited by Z.M.Stadnik
(Springer-Verlag, Berlin,1999), p.169
2) J.Friedel, Helv. Phys. Acta, 61, 538 (1988).
3) C - O 4 E 74, 24 (2004); Y.Ishii and T.Fujiwara, Physics of
Quasicrystals, edited by T.Fujiwara and Y.Ishii (Elsevier, Amsterdam,
2007), to be published.
4) R.O.Jones and O.Gunnarsson, Rev. Mod. Phys. 61, 689 (1989).
5) O.K.Andersen, O.Jepsen, and D.Glötzel, Highlight in Condensed Matter
Theory, edited by F. Bassani, F. Fumi, and M. P. Tosi (North-Holland,
New York, 1985), p.59.
6) D.Vanderbilt, Phys. Rev. B41, 7892 (1990); G. Kresse and J.Hafner, J.
Phys. Condens. Matter 6 8245 (1994).
7) M.C.Payne, M.P.Teter, D.C.Allan, T.A.Arias, and J.D.Joannopoulos, Rev.
Mod. Phys. 64, 1045 (1992).
8) A.P.Tsai, J.Q.Guo, E.Abe, H.Takakura and T.J.Sato, Nature, 408 537
(2000).
9) R. Tamura, Y. Murao, S. Takeuchi, M. Ichihara, M. Isobe, and Y. Ueda,
Jpn. J. Appl. Phys. 41, L524 (2002).
10) Y.Ishii and T.Fujiwara, Phys. Rev. Lett. 87, 206408 (2001).
60
11) R.Tamura, T.Takeuchi, C.Aoki, S.Takeuchi, T.Kiss, T.Yokoya and S.Shin,
Phys. Rev. Lett. 92, 146402 (2004).
12) R.Dronskowski and P.E.Blöchl, J. Phys. Chem. 97, 8617 (1993).
13) Y.Ishii and T.Fujiwara, J. Non-cryst. Solids 334&335 336 (2004).
14) C.S.Gómez and S.Lidin, Phys. Rev. B 68, 024203 (2003).
15) K.Nozawa and Y.Ishii, Philos. Mag. 86, 615 (2006); in preparation.
16) T.Ishimasa, Y.Kasano, A.Tachibana, S.Kashimoto and K.Osaka, to be
published in Philos. Mag.
17) M.Kraj_í and J.Hafner, Phys. Rev. B 68 165202 (2003).
61
1 Pi$m`]„¥$c%+$/#0$+$
/#%”$…ž0&!%
#BG@<5B$0+$
3‘Ž#Œg˜q"D:BG=78FH3rd0[v
1%˜qfU#Y|m[v
/Pi
0
% 60 •‰.Š
{Kœ80 •‰e—.+&!$+$/$
[v7G;4>B10-#"1%s¢$L
w#"
}˜[v$T–
$%&!#%1(
0s¢$Lw#"
+$!
#}˜XQ"1
".†j
"!
+$+‡
1%
!
#0[v3ƒ)h20
q#J"0}˜3|)u!3I1+t¡"5A
C8FH
01(
#†j"k¤ud
0!#-Z0+$/$›M%b
_$0™z3pk¤".+m`wt¡
0™£#
02
%
"
x
“z3pk¤‚k%".x
}˜šŸ3
/T!
-
$-"k¤‹l™q€a#
%an
k¤NU$O' Vo190 •‰#mWUBG
9 6 7 @! y * . 1 4 @ E ? 7 AG 9 ^ ’ ˆ 0 %
JRCAT
0!
62
%ank¤NU3š~R~,h{!!+#Kˆ#y*
Š"{O3\$!%$e$m`S
$k
´X¤$0^kj´TY3”#0;GJN;KO[r.‡^
[r'!šZ0%$“$­²0
;GJN;KO[r3‡^X¤$˜¯0!&³V„,f‹
50 ¡’$&*.11j´©q%€W$e
V„
0!3s(ª2/
0)%+%/3™«0
2 2.1 p‰W$;M:OŽi!&U$!d
M:O&guœRS¢
!
0.
;
/%©q-/+±Š—#_w3Ÿ0
01&›nAE5=X¤%²y$!&
;M:O&~c%|\z‘3®0!+§`
u©tP3bž 0!S¢
0#2¨
%*EOC
7I?F6O<D4MO8%!#
`i%;M:OAE5=%–{& MOS ]BLO<=>
&$o-#=9MO8Œƒ/²*$¥hR$
-g‚¢RQ}*.1LSI

0
&AE5=3a%”
&#;=@H!
l%
%*$&B
LO<=>m%†
#2£ˆ¦°
0
=9MO8$-£ˆ
%•¬…v#
!#".€’x
MOSFET
%=9MO8Œ
63
0&" ¼µ¡ƒ‚1'"y31
Al ¨‰9
N:COI7;N=KP,>CN>I7;N=KP"-–‰
0
·“itY$*¦² DRAM ]©†\#
CAH"» 1980 §«%#01" IBM $
" 80 §‘¤.­sY"ª[“ Al "°30¨‰u½#•
v4n/Q›‡#¶.o¸žd=L<7D" WSi2
#º¹4g£
¨‰u½0.šœlQ›‡"®
)HO?>4{nwk"]©
#0
@P;>BPd$o
¸ž4|u½—.1
0
#Q›‡$ d EPD#X«¯
3%fp‡# d EPD4 zc*
G8MJ€S d EPD#’!W(#2"S˜0.
0$¾—.1
0 @P;>BP=L<7D$
!
#2@P;>BP=L<7D"$
#bŒpdŽ
t0r+Q›!+#$ C11b m‹# WSi2 0#„µ¡
$…"~/0#fZT´
0$WSi2
$ fp- ¬f p¯Á #”
_"G8MJ€S
/Q›$0o
aŸ¥$±*!
+
30#.
WSi2 $ ¼ µ¡ # ’ !
ÀEAG5#F6#³*
q(!
"j.1¿
…WSi2 #„µ¡
Vº¹1
0
Wsi2 R^"+''!ŠUd=L<7D"`0ex"-0
—h$ 80 §‘. 90 §‘"™ˆ1¢}#HO
64
DC#~‘C=[" ‚%0ÍÜ%„(Å
/ 2000 Lj!
100nm 0£-Ÿ°!Æ,"*
-›t#cÒ-Q]DC
‚0Ì-#˜Ý
ÁÆ-AZ@4I(…Û½w«
iØ0,) NiSi ) TiSi
AZ?^æä
ÁÆ.-%#9V
MAE^C0h'!… k Žá Hf d†Ð>HdÕ
"Z;¾à0¸-'!# La d†Ð"¼eÁÆ.
-"*
AZ?^Q]DCjÊ#zŽágk!*-¥z
((1-ÑÚ"$o!ATW\AX^r™
jÊ!#Ïf{
("!
-
2.2 ³¨#AZ?^a
,//
:v"³¨d†Ð#®‹
sp Š¢qÂ"y¢6K[
U[P7]B}+.-
!PY\^) CNT "Ê})"Ãa
С(%
€(lÊ!|s'+.-S4;]6\;H]J;C"m
½+(¾à׿0±–.-fÈ¡-"!®©É
¦Íž"‰" LSI #©n¤¬"N3"Ô'‡&! CNT 0“¢
¶-r™µÙ.-—pÞγ¨!ã,†R5¨
´¨0Àâe
(8;GFHª¢,ß"`¹d†Ð­
-fÈ¡-•ŒBxC1-2xNx #`¹
cx"d†Ð0y¢
( x=0. 25 " BC 2 N #·¤§by"ËÁ¯»¡”0’PL "
O;# 2.1eV º¿ ((²©<YP24HAH:VF
Q\CËÁ¯-0ƒ
-.#u֜€š-.
#L^I:VFQvÓ"_L^IØÄ"~‘Éå!*±
65
Q^5y2#&ZO
!2v( BC2N b"(
Š'Ša‰oSœ"
!—skˆ(3
3&M%2ŽK£&
i
'kƒhs5Y,~l".
'M%2kˆ‰o
2€'kˆ
H(9D=BF:&kˆ
€BC2N 'kˆ(N=8)
I( 2eV ”˜'9D=B5.
Ÿš‹&%2#—skˆcq/1
¢!2+'
'Ša'%".#.J“%.'(
kˆI"2#.‡8?E9cq0¢2'J“kˆ&
X2eU(|cq5#.WŽ&¢2'"23(
'Rl¡(\¥el‚¡u"2'"el'52#"
'J“‚(p&f….2#Sœ%'"2j( L.Pauling
']T{&.2’133'\¥el&el8?E95V1
2#"'(Rl¡'x©&(+1*#6$L‰%
BCN &X!)C-C:3.71eVB-C:2.59eVN-C:2.83eV
B-N:4.00eVB-B:2.32eVN-N:2.11eV "2Œ –Rl¡'„ƒ8
G;EA'vg0d-'5r!€&[k
ˆ533cq2# 'eU5›2 'P'¦& HartreeFock cq&/2kƒhs1'_z8?E9'4!n
'§w'N‘(™†%04-!¨}"kˆH#kˆI"@
G>9D=B`¤'—skˆM%2#5t•3(hsž
ª'm&ZO!1kˆH&( B-N '<7Cª‰o2'&
66
ZBC2N CP;=@DM8XN(N=8)
S CP7
-C-C- -B-N-
/-.-FEH'CP'#U9\*T
')U9\>"('AIM* 32 B
'JCP:NRM
46,898,310 W&'(*
VQH';=@DKaM[a'Z$'
@<%['
5G8XCP*^'204
W&
O[`!+13,Y]_L
6+13,Y?
Z 2BC2N CP;=@DM8XN(N=32)
67
“,Ä{nšœ
,m
!X·,."~¼‰"¸
¢
~uBC2N / CVD
q’,
"lvÁ!”£ƒ½
+# CH3CN+BCl3 BC2N+3HCl !)
s†*jq
ŽBC2N(N=32)!sV¨nš
4CQ5!‡ C-N jq/€
b&²½
,½,
7OF02?
BCN h"V¨ h-BN ´À(/X·,!
,
¤–!nšY‘7OF02?<2LKT@
,)
BC2N '¯»!nšY‘
+¤–¯»§pU\"nš˜
ªZ,¡[l¿´z®Ã_9JMR9NT*´
, !Šq'V¨]_jq!/™(,!—‹n
š! BC2N " BN — C —±-znš
'!
,!)
Y,—! layer-by-layer nšnš˜ªZ°*-,"e·
~k©µr"&
¤–h]q³
I3–(¦–/¯Â±-,"‚¾!¤–]q³
*œtvÁ/º,BCN hAD=MG'
q’-
+CNT W‰!EO4>1
,*|ŸHS;:
!a¼cž-,
3. •ˆ4R6?SB6:wf/Å'!+!ix^_/
…% ‰¥
Œ+$`
"d„©
"†­©
y
¶o…%"p¹«g›9P82@(¤–h]
q³!¬znš-}"
68
+#/#~¤$,,+1%f杒10bz!
1%!.!`¹$­0ÌƶOYBA#¤ˆ$`¹$
Íg¾åw¾0AGX@‹â,‡ÖÁÀ»…#‹âx
ƒ‹âIL†´¦ÊfšBZ>…#‹âeË!
#`¹$
s*fæ¾014&&bz$3*
Þfæ!+#/$/!#]Ô
fæ$)fæ"a#ɂ"0nµ¬!/.
/!/4×."!1%+#/$¨c°ž#±"v
)
+#/#~¤½Ø0Ú¹4?QTX?UZož
4#ɂ"0WW4|1"-ڹ4]m
"bzÞ!#0#-
!+#/#Ž«
$³#$0Œ‰|½1/Št$º¾+
#/0$F@CWPJN6;DSVZ<
‹â#i¿%º¾+#/"Þ!’¿!0#
$ѓ$Íg¾r ¾!z)0
0ãäh
¾"$œÒ^Ô#ÈÄkª#²t¾§0Ð8KW:Ÿ[$z
)0$00$›l¬"-z)0$0
™Ý£$¸]!$¸Çt®}#¥q(5
z)0
FM7A¯y#~¤$™uj8KW:MZH©4Î0
1$‘s~™¾!#0#\Õp{ˆt”$3*
Þ0Ŗ"#\Õ020$€d4àb0
+)Þ0ŠtRC‹â
ј¼Ù4ß'
0ÜÀá,×á$јÛ#¬˜31·¢0
À—È#A=W›l¾!†´4•0|£#ÜÀ
á,×á4Ó_0."$™„¾"Ï"0
+030N9GJE;A{¡MZH:SEO
69
zb$ƒ”]l()^’o}'Qy
(›(&~!™GUMn[+,
j'„)(k–(
t|ŠdR“(‡I$)šhS(5F=
2i:71F3%(.0F6BF<EWŒ.0FOcY
r‰\
(]lQy)()H‹vR+ŽgVZ
`WhuT˜vR‚(e*)
()Rt|‡Iˆp•†i‘P+J
(f€${
$']mˆ–NœK~
xI…_(—
XLa_$""k–sw(e*)(#
q845E@;>-26BDF3?C79>/A
70
&(v'tœ#'Áx&
6sÛs´:
%~96Ÿ
ä°/ÅØ"
6…'uÏ"u)EJ@BD<He™"
6¢&‹/'*#
b»&
%6#"1'
:Ž"(vÍ'¹:Õ
5:yÜS6##1&‘»&ÎZ
6#:ӊ!6'3%•5°/(bX{³·ªmÝ
·b`ªm"±4}97MÔȪmYÂ&_!(7#(
Ív&§à»&}97!6#&A?F>'\×(
MÔÈYÂ"(ĝ& ;"6x&†!6ÊÌˑ„
ß#'tœ&96ÊÌ"
6in~kŒ'Dz(..%
Ë¥k¨'r¯'–á#']hO‰:CJI?#!dš7
6'"A?F>S2¹q:Ö!6#—Ø"
6
t’&(€À%#1¶a®Ûã.&94%#8
1
6ˑ„ß'pgYÂ({³·2Ý·
6(
MÔÈ%$,!'a®"
6##1&'^:
%e™"1
6
¢&( LCA(Life cycle assessment)#'âld‚7!6#&
ºK!141'
5e™(­$%!6'
#­'©Û:Ú×6#aÑ"(
67:Ú×6
Ó»(U¡^“V'P/&Ú6L"
5'e™f-.
.%Qj:|Þ&76#tƒ'1'
5e™'Å؞o&
%!6v'T[»¹q:¨µU'›Ù&9!'¦
\]h&=G@B!#(z['ё»%#"
6
4. Ɣw·¼‰ã:¾s«¦'R"[+ˆ0‡Ã¤c¬in~
kŒ'£Û»%KÒɈ.!½('ËÛ:
ÛW6#(˜Ù.'x¸NsÛ̉¿à[Иt
71
@nK*6(4!4('!}@[›•hˆžH–*™
:<8F‹%/!#\e)Œ‰{uU$*V'N
qŠR%'!#5p*)_|M—S“fS,)G8k
.QM/MM –)2'=9>?;<–0xjB–'&/Y‡
-/)"4)Iƒ3”oCH$š'%*Ug)`B
)2'…d)D~EW$*'”l*Pr(W-!#5%
%8!
4%œ‘5%$*'%b52('!
6*ŽTZOXc–Qm$*'$52(i765†
[vt(*3w3]!#5$1%]76$
5ŽTZOXc–)/)"4A™(J˜$6+%b
#5
†[L’^{„haz
€ys‚.
S.Itoh, J,Phys.:Condens.Matter 2 3747 (1990).
M.O.Watanabe et al.
Phys.Rev.Lett. 77 187 (1996).
S.Itoh, Diamond Films and Techno. 7 195 (1997).
Y.Tateyama et al. Phys.Rev. B55 R10161 (1997).
V.Heine et al. Solid State Physics 35 (1980).
72
一化学企業における第一原理計算
–無機蛍光体を中心に–
(株) 三菱化学科学技術研究センター
三上昌義
1
緒言
筆者は 1991 年 4 月から 1993 年 3 月の 2 年間、藤原研究室で修士課程
を過した。修士課程修了後、三菱化成株式会社 (現 三菱化学株式会社) に
入社したが、その後も先生と藤原研究室の皆様には学会などの機会でお目
にかかる機会も多く、現在に至るまで何かとお世話になっている。筆者が
会社の研究テーマの取り進めるに際して、先生から頂いた御助言は直接・
間接の形で影響している。また、筆者の蛍光体研究を学位論文「希土類酸
硫化物の第一原理電子構造計算」としてまとめる際に受けたご指導ご鞭撻
は、誠に有り難いものであり、感謝の念を改めて深くする次第である。
本稿では修士時代から現在に至るまでの筆者の仕事を振り返ってみた
い。同時に、藤原先生から受けた影響についても思い返してみたい。
73
2
修士時代
筆者の修士論文テーマ探しは、ノルム保存型擬ポテンシャル法及びその
応用例について調べるところからスタートした。修士１年の終わり頃に、
擬ポテンシャル法で計算された擬波動関数から、本物の波動関数を復活
させられないか、とフォーマリズムを弄くり回し、試行錯誤している内に
徒に時を過ごし、結局うまく行かなかった。(「自分はアカデミック向き
ではないなぁ」と痛感したのもこの時期だったのではないかと思う) そこ
で方向転換して、ノルム保存型擬ポテンシャルの transferability について
孤立原子モデルを超えた簡易モデルで調べることにした。当時、X.Gonze
のゴースト状態に関する解析に関する論文 [1] が出始めた処であったが、
transferability を調べるモデルは孤立原子に関するモノがもっぱらであっ
た。そこで、もう少し外部環境の影響を取り入れられるモデルで擬ポテン
シャルと本物の原子の違いを調べることが出来ないだろうか、という先生
のご発案があった。ちょうどその時、Jellium model に atom を入れたモ
デル (atom-in-jellium(AIJ)) を計算するという論文を読んでいたこともあ
り「pseudoatom-in-jellium (PIJ) はどこまで AIJ を再現するのかを見て
みるのはどうでしょう？」ということを提案し、着手することにした。
当時、Hamann の擬ポテンシャルのコードは広く配布されていて、研究
室でも既に入手済みであった。このコードを改変すれば AIJ も PIJ も取
り扱える筈なのであった。しかし、数ヶ月経っても AIJ の論文結果を再現
できず、そんな中でも就職活動やアルバイトで夏休み中に居なくなったり
するものだから、先生には御心配をかけてしまった。それでも、関連論文
を漁りまくり鋭意検討した結果、何とか夏の終わりには、アルゴリズムの
74
本質的改良に気づき、AIJ の結果を再現出来るようになり、PIJ の計算結
果も尤もらしく出るようになった。この時の心境をいま言葉で表すとした
ら「5 つの【あ】(＝あわてず、あせらず、あきらめず、あなどらず、あて
にせず) が報われた」といったところである。
この時の PIJ の計算から、擬ポテンシャルは通常の電子密度のレベルで
は本当の原子の結果を再現しうるが、電子密度が高濃度になると擬ポテン
シャルは本物の原子の結果を再現出来なくなることが見えてきた。これは
擬ポテンシャルのカットオフ半径と、jellium モデルの中の原子が電荷的
に中性になる球の半径との大小関係で決まっていると理解できた。また、
アルカリ金属元素における内殻電子密度に関する補正の効果も、この PIJ
モデルでも顕著に現れた。(結果の詳細は論文 [2] に譲る)
あとは修士論文発表までにデータを揃えて解析すれば楽勝と思いきや、
Kleinman-Bylander 近似の評価に手間取ったことや、酸素やフッ素などの
結果がなぜか思うように得られなかったことがあり、なかなか捗らなかっ
た。修士発表の前日は泊まり込みになり、発表当日の朝の 5 時に論文を打
ち出す、という始末であった。(しかも、図をはめこむのが間に合わず、図
は別冊という有様であった。先生も発表会の場で如何に言い訳するかをお
考えになられたらしい) おまけに、この修士論文の結果を論文として発表
できたのは約 2 年後であった。
このように出来の良くない学生であったが、修士論文の仕事を通じて
上記の 5 つの【あ】の心構えが自然に身に付いたことは、その後の人生
において大きな影響を与えていると思う。以前読んだ森毅先生の本の中で
「『学問は積み上げ』という信仰の故もあり、分からんなりに付き合うとい
うことが出来なくなっている。分かることを急いでいたら、研究者なんか
75
ならん方が良い。すぐには分からんことを考えて、そのうちなんとかする
のが、研究というもんなのだから。」と述べられていたが、この修士時代
に「『わからんもの』との付き合い方／なんとかする方法」の端緒を掴む
ことが出来たのだと思う。
3
化学企業における理屈屋の役割–蛍光体を中心に
博士課程には進学せず就職することに決めた筆者は、先生からの示唆
もあって、最終的には三菱化成株式会社に入社することになった。先生は
「化学会社に一人学生が行ってみるのも、また面白からずや」といった風情
であった。化学会社において、物理屋さんは圧倒的に少ない存在であり、
極論を恐れず言えば「うどんの薬味」的な存在である。(あくまでも量的
な意味であって、質的な意味はまた別であることは言うまでもないが) 物
理と化学の間の文化が違いは大きく、入社してから現在に至るまで色々と
あったが、十分に刺激的な日々を過ごしている。そんな中でも曲りなりに
研究らしきものを続けてこられたのは、一つには、先生の卒業後も変わら
ぬご支援とご激励があったからであり、改めて感謝の念を強くするのであ
る。そこで、化学会社でどのような研究をしてきたかについて、入社以来
ずっと興味を持ち続けている蛍光体の分野を例にして振り返ってみたい。
3.1
投射型テレビ用蛍光体
現在のようにカラー液晶テレビが台頭していなかった 1990 年代前半、
入社間もない筆者に与えられた課題の一つが投射型テレビ用の青色蛍光体
76
の探索テーマであった。(当時、飛行機に乗られた方は投射型テレビを御
覧になられた筈であり、その画面が普通のテレビより薄暗かったこともご
記憶されておられることだろう) 投射型テレビでは従来のテレビより注入
する電流密度が約２桁大きく、従来のテレビ用青色蛍光体 (ZnS:Ag,Al) で
は輝度飽和現象を起こしてしまうため、新しい青色蛍光体が求められてい
た。従来のテレビ用蛍光体では駄目なのであるから、全く新しい母体及び
発光中心ドーパントも視野に入れなくてはいけない。ここで経験的に投射
型テレビ用蛍光体の発光中心としては希土類系が良いのでは、という話に
なっていた。希土類元素の中で青色に発光し、かつ発光寿命がテレビ用途
向けに十分長いという要件を満たすものは Tm しか見当たらない。そこ
で筆者は「Tm に相応しい結晶母体とは何か？」という課題に取り組むこ
とになったのである。
しかし、Tm 付活蛍光体で実用になっているものは殆ど知られていない。
(現在でもその状況はあまり変わっていない) Eu や Tb の発光中心が良く
光る通常の蛍光体母体に Tm を入れても全くと言って良い程光らないの
である。そこで当時の化学組成探索は「マトリックス法による絨毯爆撃」
を繰り返しているように目に映った。(これはこれで重要なやり方であり、
この方法の究極の形が現代のコンビナトリアル合成と言えるであろう。)
とはいえ、その組成の選択方針として「ただ研究例が少ないから／未実施
だから」という歴史的な理由でなく、
「こういうメカニズムを想定したら、
こういう組成を狙うべき／狙うべきでない」というような物理屋の視点を
持ち込めないだろうか、と考えていた。
そういう問題意識を持ちつつ色々と調べてみると、電子線励起によって
生じた正孔・電子対の再結合エネルギーが発光中心に移動するための最適
77
条件が蛍光体業界ではあまり表立って議論されていないことに気付いた。
実際、この業界のバイブル「蛍光体ハンドブック」には、その辺りの議論
が詳しく出ていなかった。その点に注力して調査および考察を重ねた結
果、エネルギー移動のメカニズムに関して仮説を立てるに至り、最終的に
新しい蛍光体 M HfO3 :Tm(M =アルカリ土類) を見出すに至った。[3] Tm
系蛍光体としてはかなり明るい部類に入る母体を自分の立てた仮説に基
づいて見出したという経験は何よりも代え難いものであり、これで「研究
で何か仕事が出来そうだ」という自信に繋がった。このテーマをきっかけ
で、その後も蛍光体の世界に興味を持ち続けることになる。(興味を持続
することの重要性は、今更ながら感じる)
なお、この研究では第一原理計算は殆ど使用していないのであるが、
「コ
レコレの理由で、この部分を明らかにするためには第一原理計算が必要な
のです」と上司に説得し、藤原先生にも許可を得た上で、週末に藤原研を
出入りする時期が暫く続いた。(この時の計算の仕事をちゃんとまとめて
ないので、先生には大変申し訳なく思う) こうして上司を説得出来たこと
により、三菱化学で「バンド計算を利用した研究」のスタートをうまく切
れた訳である。
3.2
テレビ用赤色蛍光体
前述の M HfO3 :Tm が「当たった」ことから、また別の蛍光体の研究テー
マが与えられることになった。従来のテレビ赤色蛍光体 Y2 O2 S:Eu の改
良策を考えよ、ということである。理想的な輝度よりも品質が落ちるとす
るならば、それは結晶の固有欠陥 (発光のキラーセンター) が影響しそう
78
であるのだが、それに関して微視的モデルは確立していなかった。Y2 O2 S
のバンド構造も分かっていなかったのであるから、そのような欠陥モデル
の議論も不十分であったのは致し方ないところである。そこで、Y2 O2 S 母
体のバンド構造計算を実施し、そして格子欠陥の第一原理計算に進むのが
正道だろうと判断した。
ある日、藤原研に立ち寄った際、「コレコレこういう事情で蛍光体の点
欠陥の計算をしてみたいのですが、そのような欠陥のバンド計算の話題に
詳しい先生は誰でしょう？」と藤原先生にお尋ねしたところ、「例えば押
山先生はお詳しそうだから、そちらで聞いてみたら？」という示唆を受け
た。押山教授は当時ＮＥＣから筑波大に移られたばかりで、SiO2 の点欠
陥の計算をされておられたので、確かに話を聞きに行きやすい状況だと
思った。その後、JRCAT 関係のワークショップか何かで押山先生にお目
にかかり、コーヒーブレイクの時に立ち話する機会を得た。その後、発足
したばかりの押山研究室を訪問させて頂いて事情を説明したところ、共同
研究をさせて頂けることとなった。
また当時、会社の研究所で計算機の更新があり、新しい計算機の一つ
として富士通 VX-1S (VPP700 の 1CPU モデル) が納入されたことも筆
者にとって追い風となった。分子軌道計算屋さんが主に使う Gaussian や
MOPAC はディスク I/O が多い計算になるためベクトル化率は高くなく、
同時期に納入された IBM ワークステーションで流すことの方が多かった。
そこで、FFT を多用する平面波基底のバンド計算を VX-1S で計算するこ
との優位性を主張し (当時の IBM のワークステーションに比べ７倍以上
計算が速かった)、筆者が優先的に使用できる状態が整った。1996-1998 年
当時、756MB-1GB 近いメモリを利用して長時間計算出来たことは非常に
79
大きかった。
こうして計算機環境が整うと、Y2 O2 S の第一原理計算もドンドン進ん
だ。Y2 O2 S のバンド構造を初めて明らかにし [4]、Y2 O2 S の点欠陥のスー
パーセル計算より実験結果を説明する固有欠陥を提案できた。[5] 結果を
要約すると、Y2 O2 S は間接ギャップ構造であり、イオン性と共有性が共
存したような構造であり、そのような母体に生じる酸素欠損・硫黄欠損は
Y-3d に由来する深い準位を形成し発光を阻害する Killer Center となり
うることを見出したのである。また、格子間位置に酸素が入りうることを
初めて指摘した。これは浅いアクセプター準位を形成し、ホール・トラッ
プの存在を示唆する実験結果を説明する点欠陥モデルを具体的に呈示出来
た。このようにして入社以後初めての第一原理計算の結果を Phys.Rev.B
論文として出版出来たことは、素直に嬉しかった。
しかし、【原因が分かること】と【実際に改善策を見出し、それを実現
すること】との間には大きなギャップがあり、そのギャップを埋めること
に苦労した。なお、現在でもそのギャップが埋まっている訳ではない。(一
種の「死の谷 (death valley)」とも言えよう) 筆者は「Y2 O2 S の単結晶が
得られていないから、拠り所となる分析結果が手元にないのだ」と主張
し、Y2 O2 S 単結晶作りを作ってくれる外部協力者を求め奔走した。そこ
で、武居文彦先生 (当時 阪大) を始めとする様々な人との出会いがあり、
諸先生方のご厚意により、Y2 O2 S 単結晶を作製する道が拓けた。(これは
「大学の独立行政法人化」以前の話であり、今よりも大らかな風潮であっ
たことが幸いした。現在では、こう上手くは事は運ばなかったであろうと
思う)
80
3.3
ブランク中の出来事
Y2 O2 S:Eu の母体・点欠陥の計算が一段落した処で、蛍光体とは全く違
うプロジェクトに主に携わることになった。しかし本心としては、蛍光体
の仕事をもう少しキリのいいところまで進めて、結果として学位論文とし
てまとめることが出来たら、と希望していた。そのためにも、蛍光体の話
を継続するネタ探しも継続していた。そう思っているうちに、共同研究先
(東京工科大・山元明教授) から興味深い実験結果が知らされた。山元先生
の処で先述した SrHfO3 :Tm の研究が続いていたのであるが、その試料を
塗布した面に電子線を当てる際、励起する電流が大きくなると、ある電流
以上で発光強度が非線型的に増強し、しかも発光色が青色から白色に変
化したというのである。しかも、蛍光体の塗布の仕方にどうやらポイント
があるらしい。そこで現地で実際に確認してみると、見事に再現された。
またその現象は黒体放射スペクトルでは説明出来ないことも後日確認され
た。理屈がどうにも分からないのであるが、企業人としてはこれは捨て置
けないという訳で、なんとか理屈をつけて特許を書いたり、論文化のお手
伝いをしたりした。
この件と並行して、蛍光体 Y2 O2 S の単結晶作製とその評価も継続して
いた。そして Y2 O2 S に関する第一原理計算は研究所における基盤技術構
築の一環として実施させてもらった。(上司の理解に感謝する次第である)
また、縁あって 1999 年から ABINIT プロジェクト (の前身) にも従わるよ
うになり、当時は未公開だったコードを富士通 VX-1S に移植し、Y2 O2 S
の格子振動・誘電特性を計算し、実際に作製した単結晶の Raman 分光の
結果と比較した結果をまとめたり [6]、シンチレーション特性を評価する方
81
法を論じたりした。[7] この計算結果と先の非線型発光増強現象の話 [8] を
引っ提げて国際会議 ICL02(Budapest) で発表したのも良い思い出となっ
ている。
時期はやや前後するが、1999 年の年末には Y2 O2 S に関する第一原理計
算の結果がほぼまとまって、藤原研のセミナーでお話をさせて頂く機会を
得た。バルクおよび欠陥モデルの計算結果及び密度汎関数摂動法で計算さ
れたフォオン解析及び誘電特性について説明した。セミナーでは良くある
ように、質疑応答の応酬が先生と筆者の間で暫く続いた結果、それまでの
仕事を総括して学位論文をまとめて良いというゴーサインを最終的に頂
いた。それから論文推敲に取りかかったが、なかなか進捗せず、最終的に
学位を授かるまでにセミナーから約１年半かかったのであった。その間、
先生には少なからぬご心配とご面倒をおかけしたことを改めて申し訳なく
思う。
3.4
白色 LED 用蛍光体
こうして暫くの間、会社のメインの仕事としては蛍光体から離れていて
いたが、上で述べたように興味を持ち続け、研究自体は止めてなかった。
そうしているうちに、全く新しい蛍光体のテーマについて考える機会が
やってきた。その一つが白色ＬＥＤ用途の蛍光体である。この蛍光体は
InGaN 発光ダイオードによる近紫外光 (約 400nm)∼ 青色光 (約 450nm)
による光励起を用いるものであり、水銀の紫外光励起 (254nm) を利用す
る従来の蛍光体と全く性質を異にする。この励起波長の違いにより、従来
の蛍光体のほとんどが使えなくなっているのである。以前の Tm 系蛍光体
82
の研究で、そのような「常識やぶり」なアイディアを考えてきた筆者とし
ては、まさに「昔取った杵柄」的な話である。
近紫外から可視光を吸収する発光中心がドープされる母体としては、酸
化物のようなイオン性の高い物質よりも共有性の高い物質が良さそうだ
というトレンドがあり、実際に酸窒化物や窒化物の探索が世界的に盛ん
になってきている。そのような新蛍光体母体として代表的なものとして、
M2 Si5 N8 や M Si2 O2 N2 (M :アルカリ土類) などがある。[9] これらはヨー
ロッパ発の物質であり、日本としては後塵を拝しているような状況であっ
たが、最近になって日本発の窒化物系蛍光体として CaAlSiN3 :Eu が上田
氏ら (当時の所属は東京工科大学、現在は弊社) によって開発された。[10]
そこで、まず手初めにこの母体について第一原理計算を試すことにした。
この構造は AlN-wurtzite 構造が基本骨格となっており、Al の一部を Ca
に置換し、電荷バランスを合わせるように Al の一部を Si で置換したも
のであると思って、まず差し支えない。(図１) ここで興味深い特徴とし
て、Al/Si がランダム配置することにより斜方晶構造 (Cmc21 ) となるこ
とである。このことは XRD・収束電子線回折の結果から結論付けられて
いる。ここで良く考えてみると、もし Al と Si の配置を単位胞内でオー
ダーさせたら、これらは単斜晶として帰属されることになる筈である。(図
１) では実際の XRD の結果をこれらの単斜晶構造から説明出来るのであ
ろうか、という興味が湧く。また、このような Al/Si ランダム配置の結晶
を Al3+ /Si4+ に関する仮想結晶近似で計算することが可能だろうか、と
いう別の興味も湧く。このような２つの問題意識から実際の計算に着手し
た。詳しくは論文 [11] に解説したが、要点をまとめると次の通りである。
(1) 最安定構造として P 21 (No.4) と Cc(No.9) が安定であり、これらは
83
図 1: Schematic views of crystal models of CaAlSiN3 . (1 × 2 × 1)-cell
views from c-axis: (a) No.36, (b) No.8, (c) No.9, (d) No.4. Primitive
cells are shown by broken lines. The light (dark) tetrahedron denotes
AlN4 (SiN4 ) in (b)-(d); note that another model with the Al/Si atoms
switched (by exchanging the colors of the tetrahedra) is also available
for each case. In (a), Al and Si occupy the tetrahedral sites randomly.
極めて斜方晶に近い単斜晶構造になっている。また、これらのモデルの間
にエネルギー差は殆どない。よって、混合エントロピーの効果を考える
と、これらの単位胞がランダム配置することになり、この結晶位置に関す
る対称性が復活する。その結果として、XRD で斜方晶 (Cmc21 ) と判定さ
れることになる。
(2) Cc(No.8) 構造は P 21 (No.4)/Cc(No.9) 構造に比べ不安定構造であ
る。これは Pauling 2nd Crystal Rule に反する構造だからと説明できる。
この結果は、頂点共有の AlN4 四面体は二配位窒素 (N[2] ) を介して隣接す
84
ることはないことを示唆するのであるが、実際、そのような構造は知られ
ていない模様である。酸化物では頂点共有する AlO4 四面体が２配位酸素
(O[2] ) を介して隣接することはないという経験則 (Loewenstein’s rule) が
知られているが、本計算結果は窒化物でも同様な経験則が成り立ちそうだ
ということを示唆している訳である。(筆者らが知る限り、この経験則は
筆者らが初めて明言したものである)
(3) Al/Si に関する仮想結晶近似 (Virtual Crystal Approximation) は実
験構造を再現出来る。普通、異電荷間の仮想結晶近似は行わないので、こ
のようなベンチマーク結果は稀有な例だと思う。定性的には Al/Si の化学
的性質 (イオン半径や電気陰性度) が類似であるからと理解しているが、よ
り定量的な説明が待たれるところである。
窒化物・酸窒化物の開発は “まだまだこれから” と言った状況であり、
新しい多元系化合物が合成されるチャンスも多く残されているであろう。
Pauling の結晶則を窒化物・酸窒化物に特化した形で具体化することが出
来れば、結晶構造に関するイメージを具体的に膨らますことが出来る筈で
ある。そうなれば、新物質探索は効率良く進められるであろうし、実験結
果に対峙する姿勢も異なってくるであろう。「こうなるのが筋が通ってい
るし、美しい」というような大局観／審美眼が理論屋と実験屋の間で共有
出来るようになればしめたモノである。上で提言した経験則は、そんな理
想の姿に向かうための一里塚である。
今後は、計算物理の道具を上手く用いて、実験屋さんの【暗黙知】を
【形式知】として上手く表出化していきたい。[12] 新規化合物合成の
feasibility に関して実験屋が持つ独特の感性と、そして実際の実験結果を
如何に物理／化学／結晶学の言葉で言語化出来るか、が目下の課題であ
85
り、実験屋さんと議論を重ねている。「全く新しい多元系化合物が出来る
のではないか」という期待が高まっては、
「その方針ではダメだ」と落胆す
る、といった一喜一憂の日々が続いているが、この状況に特に心配せず、
むしろ楽観視している。というのも、
「科学とは (エラーの) 自己修正過程
である」という Carl Sagan 教授の言葉通りの日々を過ごしていると言え
るからである。そう言えば、藤原先生は研究室のホームページで「学生諸
君に『好奇心と意欲と楽天性』を望む」という主旨のメッセージを送って
おられた。この気構えはきっと藤原研時代に涵養されたのであろう。今後
もこの心意気を持って研究を続けていきたいと思っている。
4
結語
こうして過去十数年を振り返ってみると、節目節目で藤原先生に貴重
なご指南を頂いていることに気づく。改めて感謝の念を強くする次第であ
る。世間に役立つような／世界を明るくするような材料開発に少しでもお
役に立ちたいと思って日々の研究を進めているが、そのようにして開発さ
れた材料を手に先生にご報告することが出来れば、細やかながらもご恩返
しになるのではないか、と思っている。(先生、今しばらくお待ち下さい)
参考文献
[1] X. Gonze, R. Stumpf and M. Scheffler, Phys. Rev. B44, 8503 (1991).
[2] M. Mikami and T. Fujiwara, J. Phys. Soc. Jpn. 64, 1629 (1995).
86
[3] H. Yamamoto, M. Mikami, Y. Shimomura, Y. Oguri, J. Lumin.
87-89, 1079 (2000). 特許第 3475565.
[4] M. Mikami and A. Oshiyama, Phys. Rev. B57, 8939 (1998).
[5] M. Mikami and A. Oshiyama, Phys. Rev. B60, 1707 (1999).
[6] M. Mikami, S. Nakamura, M. Itoh, K. Nakajima, and T. Shishido,
Phys. Rev. B65, 094302 (2002).
[7] M. Mikami, S. Nakamura, M. Itoh, K. Nakajima and T. Shishido,
J. Lumin. 102-103, 7 (2003).
[8] H. Yamamoto, M. Mikami, and S. Nakamura, J. Lumin. 102-103,
782 (2003).
[9] R. Mueller-Mach et al., Phys. Stat. Sol. (a) 202, 1727 (2005).
[10] K. Uheda et al., Electrochem. Solid-State Lett. 9, H22 (2006).
[11] M. Mikami, K. Uheda, and N. Kijima, Phys. Stat. Sol. (a) 203,
2705 (2006).
[12] 芸の道の名人が「コツ」と呼んでいるものは、「言語によって形式知
化された暗黙知」に他ならない。そのような言語化により、名人の技
は継承可能な形になるのである。これは芸事に限らず、企業における
実験技術の伝承においても重要である。
87
Efficient Linear Scaling Approaches for Computing
Quantum Transport in low-dimensional systems
CEA/DSM/DRFMC/SPSMS/GT,
17 avenue des Martyrs, 38054 Grenoble, France
Stephan Roche
1
Introduction
Today’s research endeavours in molecular electronics or nanoelectronics
rely on both : advanced top-down and bottom-up techniques to fabricate and engineering functional molecular based devices and circuits, but
also advanced computational schemes that enable in-depth exploration of
their novel properties and ultimate performances. The numerical simulation strategy of material and devices at the nanoscale is twofold. On
one hand, first principles (ab initio) methods, such as Density Functional
Theory (DFT) and its extensions push the limits of realistic simulation of
the physico-chemical complexity at the atomistic scale, searching for elaborated description of electron-electron correlations, to the price of limited
implementation for large scale simulations and out-of equilibrium problems.
On the other hand, semi-empirical approaches provide a more suitable basis for simulating large scale systems, albeit limited to effective hamiltonians with adjustable parameters. Both directions rely however on approximations of the true many-body problem. The combination of both
approaches opens new perspectives to circumvent intrinsic limitations of
separated methods, enabling to explore advanced quantum transport phenomena in materials and devices with large chemical complexity. To compute current-voltage characteristics of realistic nanodevices, efficient linear
88
scaling computational methodologies of quantum transport at equilibrium
and far for equilibrium are of crucial importance.
In this contribution, we will review our developed order N scaling transport methods, that allow to efficiently compute either the Kubo-Greenwood
or the Landauer-Büttiker conductance in materials or heterojunctions of
large complexity (provided we restrict to physical situations close to thermodynamic equilibrium). These methods are based on the recursive construction of either orthogonal (Lanczos-scheme) or bi-orthogonal basis in
which hamiltonian matrices are first tridiagonalized, and continued-fraction
expansion further used to accurately compute off-diagonal Green’s function
matrix elements. They turn out to be of broad range of applicability and
well suited for real space calculations, so using expansion of wavefunctions
in localized basis sets. This includes semi-empirical or ab initio methods
based on localized functions and implemented within DFT method such
as SIESTA [1]. Following the same philosophy out-of equilibrium methods
could be further investigated.
2
2.1
Kubo-Greenwood and Landauer-Büttiker real
space order N methods
Conduction mechanisms and conductance scaling
Efficient computational recursion and order N methods have been successfully developed in solid-state physics since their introduction by R. Haydock
[2, 3, 4, 5]. The recursion methods are based on an eigenvalue approach of
Lanczos [6], and rely on the computation of Green’s functions matrix elements by continuous fraction expansion, which can be implemented either
in real or reciprocal spaces. These techniques are particularly well suited for
treating disorder and defect-related problems, and were successfully implemented to tackle with impurity-level calculations in semiconductors using
tight-binding approximation [7], or with electronic structure investigations
for amorphous semiconductors, transition metals and metallic glasses based
on the linear-muffin-tin orbitals [8].
89
On the other hand, the general electronic transport theory in the linear
response regime relies on the approach derived by R. Kubo [9]. In its zero
frequency limit, it reduces to the trace of the operator δ(E − H)V̂x δ(E −
H)V̂x , that relates the spectral measure operator δ(E − H) to the velocity operator V̂x . Electronic dc-conductivity is thus seen as a measure of
autocorrelation average of wavepackets velocities. During the past decade,
we have endeavoured the implementation of the recursion method in the
calculation of the Kubo-Greenwood transport coefficients [10, 11, 12]. This
approach, optimized over the years, is mainly based on the resolution of
the time dependent Schrödinger equation (TDSE) employing an expansion
of the spectral measure onto a basis of orthogonal (typically Chebyshev)
polynomials Qn (E). By doing so, any operator can be expanded within the
same basis as
Ht
)|Ψ(0)i
h̄
¶
N µZ
X
Et
=
n(E)e−i h̄ Qn (E)dE Qn (H)|ψ(0)i
|Ψ(t)i = exp(−i
n=0
which allows an efficient and quick resolution of TDSE. Afterwards the
diffusion coefficient as well as the related Kubo conductance can be straightforwardly computed as
D(E, t) =
G(E) =
1
h(X̂(t) − X̂(0))2 iE
t
2e2
lim Tr[δ(E − H)D(E, t)]
Lsyst t→τL
given their direct relation to the time-dependence of electronic wavepackets. The calculation of both quantities can be reduced to the computation
of on-diagonal Green’s function that are achieved by a continuous fraction
expansion. From this methodology we have been able to analyze all the
main conduction mechanisms in the quantum coherent regime in various
low-dimensional systems or disordered materials. They include the
90
• Ballistic regime: G(E) = (2e2 /h)N (E)
• Diffusive regime: G(E) = (2e2 /h)`e /Lsyst
• Weak localization: G(E) =
2e2 `e
h Lsyst
• Strong localization G(E) =
2e2
h exp(−`e /Lsyst )
− δGW L
defining N (E) the conducting channels number, `e the elastic mean free
path associated to a given disorder model, and δGW L the weak localization
correction to the quantum conductance. In Fig.1, one shows an illustration
of all those conduction mechanisms as well as the obtained conductance for
the ballistic case for a metallic (10,10) nanotube with or without Anderson
type disorder (taken as random fluctuations of onsite energies within an
interval [−W/2, W/2]).
This approach was successfully applied to the study of quasiperiodic
systems [10], quantum-Hall effect [11], carbon nanotubes [12] and semiconducting nanowires [13]. Hereafter, we will show a quantitative validation of the method by comparison with analytical results for a model
Anderson-type disorder, and further illustrate the quantitative predictability of the methods for studying quantum transport in realistic models, such
as chemically modified carbon nanotubes, effects of magnetic fields, or incommensurability in multiwalled nanotubes. Finally, the essential steps for
generalizing the approach to the Landauer-Büttiker framework, essential
for simulating nanoscale devices, will be drawn.
2.2
Elastic mean free path scaling properties in Carbon
Nanotubes
For a in-depth understanding of disorder effects in low dimensional systems,
the evaluation of the elastic mean free path (`e ) is a first fundamental step.
For sufficiently weak disorder, a perturbative treatment can be performed
within the Fermi Golden Rule (FGR), giving a direct access to the elastic
mean free path `e = vF τ . In the context of disordered metallic carbon nanotubes [14, 15], an analytical result was first derived by White and Todorov
91
Figure 1: Top: time-dependent diffusion coefficient for various conduction regimes in metallic carbon nanotubes (10,10) with or without elastic
disorder (Anderson-type). Bottom: quantized energy dependent Kubo conductance obtained for disorder-free nanotube.
92
[16]. By reducing the bandstructure to a two-band approximation, and describing the disorder by the onsite Anderson-type potential (see below),
`e was analytically derived, and found to linearly scale with diameter for a
fixed disorder strength W , while at a fixed diameter (note that a (n,m) nan√
√
otube has a diameter dt = 3acc n2 + m2 + nm/π), the expected disorder
scaling `e ∼ 1/W 2 was shown. The analytical derivation of such fundamental length scale requires the calculation of total DoS in the vicinity of Fermi
level. The DoS can be generally written as ρ(E) = Tr[δ(E − H)] where the
trace has to be developed over a complete basis set, while the application
of the FGR yields
¯
¯2
¯
1
2π ¯¯
=
hΨn1 (kF )| Û |Ψn2 (−kF )i¯¯ ρ(EF ) × Nc NRing
¯
2τe (EF )
h̄
(1)
Figure 2: Main frame: Energy dependent mean free path as a function
of diameter. Inset: 1/W 2 -scaling in agreement with Fermi golden rule.
Adapted from [17].
93
with Nc and NRing , the respective number of pair atoms along the circumference and the total number of rings taken in the unit cell used for
diagonalization, whereas the eigenstates at the Fermi level are rewritten as
X
1
|Ψn1,n2 (kF )i = p
eimkF |αn1,n2 (m)i with
NRing m=1,N
Ring
|αn1 (m)i =
|αn2 (m)i =
√
1
2Nc
1
√
2Nc
Nc
X
n=1
Nc
X
e
2iπn
Nc
e
2iπn
Nc
µ
¶
|pA
z (mn)i
+
|pB
z (mn)i
−
|pB
z (mn)i
¶
µ
|pA
z (mn)i
(2)
n=1
while the disorder considered here is an uncorrelated white noise (Andersontype) distribution given by
A
0 0
hpA
z (mn) | Û | pz (m n )i = εA (m, n)δmm0 δnn0
b
0 0
hpB
z (mn) | Û | pz (m n )i = εB (m, n)δmm0 δnn0
A
0 0
hpA
z (mn) | Û | pz (m n )i = 0
(3)
where εB (m, n) and εA (m, n) are the onsite energies of electron at atoms A
and B in position (m, n), randomly distributed within the interval [−W/2, W/2]
and following some uniform distribution with probability P = 1/W . Then
by replacing Eq.2 in Eq.1, using Eq.3, a straightforward calculation gives :
µ
1
πρ(EF )
1
p
=
τe (EF )
h̄
Nc NRing
X
Nc NRing
ε2A
1
+p
Nc NRing
X
¶
ε2B
Nc NRing
Hence, if the disorder is described by random fluctuations of onsite energies
with uniform probability 1/W (W the disorder bandwidth) the mean free
path can be finally analytically [16, 18] derived as
`e =
18acc γ02 p 2
n + m2 + nm ∼
W2
µ
γ0
W
¶2
dt
For the armchair m = n = 5 nanotube, with disorder W = 0.2γ0 , applying the above equation, one finds `e ∼ 560 nm which is much more larger
94
than the circumference length. As shown in Fig. 2, numerical studies [17]
confirm the scaling law of the mean free path with the nanotube diameter
close to the charge neutrality point. For semiconducting bands, the 1/W 2
is still satisfied, but mean free paths are seen to be much smaller and do
not scale with diameter, in full agreement with experimental estimates [19].
2.3
Weak localization and energy dependent coherence lengths
The understanding of localization effects in disordered mesoscopic systems
stands as a central issue based on the quantum interference effects (QIE)
on charge transport [20, 21]. These QIE between clockwise and counterclockwise backscattering paths develop in the so-called coherent regime,
and yield an increase of the return probability to the origin for propagating
wavepackets. The contribution of QIE is usually reduced by several inelastic
scattering sources that produce decoherence of the wavepacket phase. At
low temperature the main decoherence mechanisms are e-ph and electronelectron couplings. Within the framework of weak localization theory, it has
been possible to derive perturbatively the relation between the measured
conductance G(E), its quantum correction δGW L (E) and the coherence
length Lφ that fixes the scale beyond which QIE are destroyed. The estimation of the coherence lengths is a central issue in mesoscopic physics,
and weak localization provides an elegant framework to extract the behavior of Lφ , that mainly depends on the dimensionality of charge transport
[9]. Stojetz and coworkers [19] have recently succeeded in measuring the
energy-dependence of the coherence length scale, by using an efficient backgate electrode able to move the Fermi level position and explore the physics
through different subbands. The magnetoresistance data were well fitted
by the conventional theory as [21]
δGW L
¶
µ
e2
1
W 2 e2 B 2 −1/2
= −A
+
,
πh̄L L2φ
3h̄2
(4)
with L the tube length, W the diameter and A a normalization factor.
From this expression, the coherence length can be quantitatively extracted.
95
Therefore, by using Lφ (T ) = (GDLh̄2 /2e2 kB )1/3 T −1/3 , the diffusivity was
deduced D = 100cm2 s−1 (at zero gate voltage) as well as the corresponding
elastic mean free path `e = 2D/vF ∼ 20nm. By reproducing this calculation at several values of the gate voltage, a corresponding energy-dependent
pattern was revealed. The spectacular observation was the systematic decrease of Lφ (as well as `e ) near the onsets of each new subbands (van-Hove
singularities positions). If the energy dependence of `e has been revealed
numerically [12], the energy dependence of Lφ was to date unexplored.
Finally, by studying the temperature dependence of Lφ , the decoherence mechanism was attributed to electron-electron scattering, in agreement with conventional theory. Electron-electron and e-ph scattering are
two different sources of quantum dephasing, but within weak localization
theory the resulting decoherence phenomenon is of similar nature, and for
instance the derivation of the temperature dependent coherence time follows some general formal treatment.
The energy dependence of the coherence length scaling properties in
carbon nanotubes can be analyzed in the situation where e-ph interaction becomes the dominant source of decoherence. The physical origin of
the coherence length can be simply understood as follows. First, weak
localization corrections are evaluated when computing the conductance
for wavepacket propagation between two real space positions, let us say
2
from P to Q : G = 2eh PP →Q which can be further expanded as PP →Q =
P
P
i(αi −αj ) , where the summation over probability am2
i6=j Ai Aj e
i |Ai | +
plitudes (denote Ai ) includes all possible electronic pathways i (see Fig.3
for illustration). Here one separates the total probability between the
classical and interference parts. Then averaging over (elastic) disorder
configurations reduce the contribution of interference terms to a single
class, namely the paths that include a finite loop which return to some
initial point,
expressed analytically as
¯ further
¯
¯ let’s say O, that¯2can be
¯
¯ ¯2
¯
+iα
iα
+
−
+ A− e
PO→O = ¯A+ e
¯ = 4¯A0 ¯ , since clockwise and counter
clockwise probability amplitudes have the same phase factor in case of time
reversal invariance [22]. Weak localization thus results from the enhancement (doubling) of the return probability to the origin, with a consequent
96
increase of quantum resistance. To understand decoherence, one must realize that when inelastic mechanisms come into play on top of coherent
transport, then clockwise and counterclockwise pathways accumulate a difR
ferent superimposed random phase factor eiα+ = he+iφ i = dφP (φ)eiφ and
R
eiα− = he−iϕ i = dϕP (ϕ)e−iϕ , respectively. These expressions are derived
by considering that the coherent propagating wavepackets is coupled to
some external fluctuating potential, such as a thermal bath or fluctuating
electromagnetic field that respectively encode the e-ph and e-e interaction
processes [22]. P (φ) denotes the probability distribution function for the
phase ϕ, that results from the combination of inelastic and elastic scattering processes, both having a probabilistic nature along the electronic paths
(Fig.3 (c) and (d)). Besides, if u(r, t) defines the phonon displacement field,
these phase factors are physically related to the Lagrangian of the system,
through its action along the electronic pathways, giving [22]
Z
mC
1
φ, ϕ =
v.(v.∇)u(r, t) − v2 ∇.u(r, t)dt,
(5)
h̄
d
defining mc the atom carbon mass; and d the space dimensionality, and
r, v describing the position and related velocity of the electronic pathways
[22], while the summation is performed for an elapsed time required to
close the relevant loop. The coherence time τφ is the time beyond which a
full uncertainty of the phase difference is achieved between clockwise and
counter clockwise paths (that is when the width of the distribution P (ϕ)
becomes in order of ∼ 2π) . Accordingly, the coherence length Lφ (E) is
the length associated with the loop that is closed after an elapsed time
t = τφ . Finally, within this phenomenology, the quantum correction of conductance is given by the relevant integrated return probability to the origin:
R∞
2
δGW L ∼ 2ehD 0 PO→O (t)(e−t/τφ − e−t/τe ), thus the quantitative measure
of conductance correction is driven by independent calculation of the return
probability to the origin in the coherent regime on one side, and estimation
of the relevant coherence time on the other side. The scaling properties of
the quantum correction in the coherent regime can be investigated numerically by solving the time-dependent Schrödinger equation (TDSE), whereas
97
Figure 3: (a) Representation of electronic propagating paths between two
points P and Q (a); several paths enclosing a loop returning back to some
origin (b), with clockwise (-) and counterclockwise (+) parts; mechanism
of decoherence on those paths (arrow indicate inelastic scattering events
occurring at different points for (+) and (-) paths.
98
the coherence time will be given by the considered total elapsed time, at
which the TDSE is solved [23]. By doing so, the energy-dependence of
the weak localization pattern is correctly reproduced without however providing a quantitative information about the exact value of the coherence
length scales, neither their temperature dependences. This point should
be considered on the basis of a microscopical modeling of e-ph interaction
and rigorous quantum mechanical treatment of inelastic transport. The
resolution of the TDSE can be made by expanding the evolution operator
e−iĤt as a product of short-time evolution steps e−iĤ∆T , for a total evolution time t = n∆T . Typically ∆T is one tenth of the oscillation period of
the considered phonon acoustic mode. During each elapsed time ∆T , the
hamiltonian energetics is fixed by the static part of the Anderson potential, whereas the time-dependent part, due to long range vibrational modes,
varies in conjunction with the associated overlap integrals modulations [23].
To extract the energy-dependence of the coherence length scale, we
proceed as follows. In the weak localization regime, the quantum correction of the Drude conductance is computed by solving the so-called
Cooperon equation [21]. Assuming a quasi-1D geometry of the system,
it has been shown that decoherence either due to e-ph or electron-electron
scattering is described within
the same scheme, and
µ
¶ that the conductance
reads GKubo (E) =
2e2
h
`e (E)
N⊥ (E) L(E,t)
− δGW L (E)
where L(E, t) is the
length scale that is energy-dependent due to velocity v(E), and scales as
p
L(E, t) = v(E)`e t, in the diffusive regime, whereas the term δGW L (E)
gives the contribution of QIE beyond the scale of `e . Within the weak localization theory, and for quasi-1D systems, this contribution is shown to
be related to the coherence length Lφ (E) as δGW L (E) = Lφ (E)/L(E, t)
whereas τφ (E) = L2φ (E)/v(E)`e . Therefore, by exploring the scaling behavior of δGW L (E), one can access relevant physical information about the
fluctuations of the coherence length scales in the weak localization regime
[23].
Fig. 4(left panel) shows the conductance computed at two different evolution times, t = 3500h̄/γ0 and t = 35000h̄/γ0 ' 8ps (dashed curves).
99
Figure 4: Left panel: Conductance for the disordered (10,10) tube with
W = 0.07γ0 , taken at evolution times t = 3500h̄/γ0 and t = 35000h̄/γ0
(dashed curves). The bold curve gives the classical part Gclass (E) =
G0 N⊥ `e (E)/L(E, t), whereas the dotted curve gives the quantum correction δGW L (E) (t = 35000h̄/γ0 ). Right panel: Coherence lengths deduced
from δGW L (E) = Lφ (E)/L(E) computed at t = 35000h̄/γ0 , and for several
values of disorder potential. Inset: Corresponding τφ (E) for the same parameters. In addition to the random short range disorder, a twist acousticmode modulation is introduced.
100
One clearly sees a downscaling of the conductance with time, that comes
from the classical linear downsizing in the diffusive regime `e (E)/L(E, t),
together with the increasing QIE contribution of δGW L (E, t) with time (or
equivalently length) scale. In the same figure, one also reports the classical term Gclass (E) = N⊥ (E)G0 `e (E)/L(E, t)(bold curve), along with the
quantum interference term δGW L (E, t) = GKubo (E, t) − Gclass (E) (dotted
curve), at t = 35000h̄/γ0 . These results are obtained for W = 0.07γ0 and
no phonon dephasing, and by adding the time dependence modulations of
integral overlap associated to the TW mode does not bring any appreciable
changes [23].
Fig.4(right panel) gives Lφ (E) (main frame) and τφ (E)(inset) for W =
{0.1γ0 , 0.2γ0 , 0.5γ0 } (bold, dotted and dashed curves respectively). The
values range within [10nm, 1000nm] in the considered energy window. One
notes that for energies ≤ 0.9eV and within the considered evolution time,
the quantum correction δG(E) ' 0, so that no meaningful information
about the coherence length can be deduced since transport remains quasiballistic. In contrast, the coherence time shows reversed behavior, owing
to the strong decrease of `e .
In the experimental situation [19], the fluctuation of Lφ (E) due to
electron-electron scattering were found to scan the range [10nm, 60nm],
with systematic decrease near the onsets of new subbands. The values
given here (for the chosen evolution time) are thus physically reasonable,
since e-ph scattering is expected to lead to weaker decoherence effect [22].
2.4
Weak localization regimes and magnetic fields effects
Hereafter the weak localization phenomena in metallic carbon nanotubes
are illustrated, on the basis of numerical results [24]. We first discuss for
the metallic (9, 0) nanotube, the various transport regimes under magnetic
field. In Fig.5, the behavior of the field-dependent diffusion coefficients
are shown in regards to the value of `e , that is modulated following its
analytical relation to disorder strength W . Indeed, by using the Andersontype disorder, the value of `e can thus be tuned by the disorder strength
101
W , and the several cases of interest can be explored.
Figure 5: Diffusion coefficient D(τφ , φ/φ0 ) (in units of Å2 -γ0 /h̄) for the (9,0)
nanotube evaluated at time τφ À τe , for two disorder strengths, W/γ0 = 3
and 1, such that the mean free path (`e ∼ 0.5 and 3 nm, respectively) is
either shorter (dashed line) or larger (solid line) than the nanotube circumference (|Ch | ∼ 2.3 nm). The right y-axis is for the dashed line and
the left y-axis is for the solid line. Inset: D(τφ , φ/φ0 ) for `e = 3 nm and
L(τφ ) < 2`e .
First, the weak localization regime [21] is analyzed under the condition `e < |Ch | < L(τφ ). Fig. 5 shows that the diffusivity increases at
low fields (negative magnetoresistance) and that the periodic AharonovBohm oscillations are dominated by a φ0 /2 period, i.e., D(τφ , φ + φ0 /2) =
D(τφ , φ) in agreement with weak localization theory. In contrast, when
`e > |Ch |, L(τφ < 2`e ), the system exhibits a positive magnetoresistance associated with D(τφ , φ + φ0 ) = D(τφ , φ). For the case `e > |Ch |, L(τφ > 2`e ),
negative magnetoresistance and Aharonov-Bohm oscillations with period
102
φ0 are obtained. Note that with the analytical formula for the mean free
path and estimates of disorder values, one gets `e ' 104 × |Ch |, by using |Ch |
as the circumference of the outer nanotube in the experiment of Bachtold
and coworkers [25]. This leads to some inconsistency since the theoretical value of the mean free path is a priori too large to be consistent with
a φ0 /2 Aharonov-Bohm oscillation. Although the bandgap opening and
oscillations, as well as other bandstructure changes (van-Hove singularities
splitting and shifting) are likely to be smoothen by disorder, the magnetofingerprints will obviously result from a entangled situation, that goes much
beyond the conventional theory of weak localization [24]. Several early
experiments suggested such additional complexity in analyzing magnetotransport measurements [26, 27]. Roche and Saito [12] theoretically found
that for a fixed disorder strength, magnetotransport fingerprints present
strong fluctuations as a function of Fermi level position, CNT diameter,
and orientation of the magnetic field with respect to the tube axis. Thanks
to the engineering of efficient electrostatic gating of the nanotube, some experimental evidence of such multiple Aharonov-Bohm effects, with a φ0 /2
oscillation driven by weak localization, superimposed to φ0 -periodic resistance fluctuations related to bandstructure modulations have been reported
[28].
2.5
Quantum Transport and localization effects in chemically doped nanotubes
The possibility to incorporate chemical impurities as substitutions of carbon
atoms has been demonstrated experimentally [29] and offers novel possibilities to investigate coherent charge transport, and magnetoresistance phenomena in chemically modified carbon nanotubes. Substitutional doping
by nitrogen or boron impurities has been a very intense research topic at the
theoretical level during the recent years [30, 31]. Initial works focused on
the effect of a single isolated defect on electronic and transport properties,
while further studies have addressed the issue of mesoscopic transport in
µms long nanotubes with random distributions of impurities. These trans-
103
port methods are mainly based on Kubo or Landauer-Büttiker frameworks,
and mostly employed ab initio calculations combined with semi-empirical
π − π ∗ hamiltonian [32, 33, 34]. Such studies have allowed to explore the
fundamental elastic transport length scales (elastic mean free path `e (E)),
and to investigate quantum interferences phenomena bringing the system
from the weak to the strong localization regime.
Figure 6: Top: representation of a single adsorbed potassium atom on
a CNT connected in between two metal contacts. Bottom: long-range
variations of the on-site and hopping parameters around K and N impurities
derived from ab initio calculations [32].
To elaborate an effective tight-binding model able to describe the physics
around the Fermi level, it is sufficient to properly describe the long-range
scattering potential due to the chemical impurities (assuming low density
approximation, i.e. no possible interferences between tails of individual
impurity potentials). The long-range variations of the on-site and hopping
parameters around impurities can be directly derived from ab initio calcula-
104
tions performed using atomic-like basis [32]. The obtained tight-binding parameters allow to reproduce the position of the quasibound states in perfect
agreement with the ab initio calculation [30]. As a matter of illustration,
in Fig.6 the evolution of the ab initio on-site hamiltonian matrix elements
associated with the pz orbitals (perpendicular to the plane of atoms) as a
function of the distance to the impurity in a doped (12x12) graphene sheet
[32]. The potential well created by N in substitution is clearly much deeper
that the one associated with the partially screened K+ ion. In particular,
the ability of adsorbed K ions to trap electrons is significantly reduced as
compared to N impurities.
Figure 7: Left: Ab-initio calculation of the conductance for a single nitrogen impurity in substitution of one carbon atom in the (6,6) metallic
nanotube [32]. Right: Length-dependence of the Landauer conductance for
the disordered (10,10) nitrogen doped nanotube at several energies (doping
is fixed to 0.1%). Inset: Conductance versus energy for the perfect (dashed
line) and single-impurity (solid line) cases for a single defect. Arrows show
the considered energies for the scaling analysis (main frame) [34].
In Fig.7, the Landauer conductance computed from the ab-initio method
(left) and the tight-binding model for a single nitrogen doped (10, 10)
armchair nanotube are shown (right-inset). At selected energies (left inset:arrows), the conductance scaling properties are shown for a fixed impurity density ndoping = 0.1% (main frame). The extraction of the elastic
mean free path `e is achieved by adding the contribution of the ballistic
105
term to the diffusive one, i.e. R = 1/G = R0 /N⊥ + R0 /(N⊥ Ltube /`e ),
where R0 = h/2e2 is the resistance quantum and N⊥ the number of available transverse modes at a given energy. When `e /Ltube À 1, the statistical
distribution of T is found to be narrowed and centred around N⊥ , in agreement with a ballistic limit G = G0 N⊥ , with G0 = 1/R0 the conductance
quantum. The other asymptotic case is found when `e /Ltube ¿ 1, where
the distribution of T becomes wider with a mean value downscaling with
the tube length as T (E) = N⊥ (E)`e (E)/Ltube . The conductance downscaling at a given energy exhibits a crossover from a ballistic to a diffusive
regime.
Figure 8: Main frame: `e for several values of ndoping (from lower to upper
curve: 0.3%, 0.2%, 0.1%, 0.05%). Inset(top): T at CNP, for an average
over 200 configurations(ndoping = 0.1%). The linear fit (dashed) directly
gives `e . Inset(bottom): ln T at E = 0.69eV, with linear fit (dashed) giving
access to ξ(ndoping = 0.1%).
Fig.8 shows the full energy-dependence of the elastic mean free path (`e ),
evaluated from the Kubo formula, for several values of the doping density.
It is readily shown that over the whole spectral window `e ∼ 1/ndoping
which is expected from the FGR. At charge neutrality point, the scaling of
`e shows a linear increase with the radius of the tube r as found in a simpler
Anderson model of disorder (not shown here). The comparison between `e
extracted from both the Kubo and the Landauer formalism are in very good
agreement over the full spectrum at a quantitative level [34]. For instance,
the top inset of Fig.2 shows the averaged transmission as well as the fitting
curve 2/(1 + Ltube /`e (E)), that yields `e (E = 0.00) = 495.5 ± 17.4nm
in excellent agreement with the Kubo calculation that gives ≈ 460nm.
106
For an energy at the frontier of the first subband below the CNP, i.e.
E = −0.78eV (Fig.1-a), the conductance slowly decays with length, and
the regime remains quasiballistic. This is consistent with the calculated
mean free path `e (E = −0.78) = 8371.6 ± 69.4nm which is much larger
than the maximum length (Ltube = 3000nm) of the tube in between contact
probes. For energies close to the nitrogen quasibound states (Fig.7), the
impurity-induced backscattering becomes very strong, yielding a very small
mean free path (see also Fig.8). On the “s wave” resonance, one actually
finds that `e (E = 0.69) = 8.2±0.8nm. For length larger than the mean free
path, the conductance becomes exponentially reduced with length (Fig.6,
curves c and d), defining a localized regime with ξ the length scale that
quantifies the exponential decay of exp (ln T ) [20]. Accordingly, one gets
ln T (E) = −L/ξ(E) [34].
2.6
Incommensurability and anomalous quantum transport
in multiwalled carbon nanotubes
Thanks to the order N calculation, we were able to establish a deep connection between spectral and dynamical properties of incommensurate multiwalled nanotubes (MWNTs) [24, 35]. The MWNTs are indeed mostly intrinsic incommensurate objects since, due to registry mismatch between
neighboring shells, there are very few cases in which the respective symmetries of individual shells allow a common unit cell for the whole object.
In most situations, the unit cell length (along the nanotube axis) ratio between adjacent shells is an irrational number, and the MWNT taken as a
whole becomes an incommensurate object.
The study of energy dependent diffusion coefficients in micron long
multiwalled carbon nanotubes with intrinsic incommensurability between
neighboring shells, demonstrates that the conduction regime evolves from
sub-ballistic to diffusive motion as the Fermi level moves to subbands regions of large density of states, as a consequence of enhanced contribution
of the underlying aperiodic potential and multiple scattering phenomena.
In incommensurate systems, such as the (6, 4)@(17, 0)@(15, 15) triple
107
time dependent diffusion coefficient for incomFigure 9:
Left:
mensurate/commensurate disorder-free MWNTs (with β = γ0 /8).
Right: length dependence of conductance for two Fermi energies for
(6, 4)@(17, 0)@(15, 15). Adapted from Triozon:PRB04
108
wall nanotube, the energy-dependent anomalous conduction due to geometrical incommensurability-induced aperiodic potential is particularly strong.
The region around the charge neutrality point remains almost ballistic, as
expected from the suppression of the intershell coupling at low energies due
to helicity-determined selection rules [36]. In contrast, the rest of the electronic spectrum shows a very slow expansion of the wave packet in time. In
Fig. 9 ones shows the time-dependent evolution of the diffusion coefficient
(left) together with the length dependence of the Kubo conductance (right)
by defining the diffusion coefficient as
D(E, t) =
X 2 (E, t)
trδ(E − Ĥ)(X̂ (t) − X̂ (0))2
=
.
t
trδ(E − Ĥ) t
D(E, t) either shows anomalously slow diffusion, or saturates at long
times. Whenever the saturation limit is reached, a mean free path `e (E)
for the whole object can be meaningfully extracted. The possibility to get
anomalous diffusion and elastic mean free path intrinsically associated to
an aperiodic potential of an otherwise clean system stands a unprecedented
feature in Condensed Matter Physics. In contrast to quasiperiodic systems
as model systems for quasicrystals [10, 37, 38], the simplified tight-binding
description in MWNTs correctly reproduce the electronic bandstructure of
real nanotubes, so a quantitative comparison between theory and experiments is possible! Despite recent efforts, the unique transport properties of
MWNTs predicted theoretically remain to be experimentally confirmed.
3
Bi-orthogonalization process and Landauer-Büttiker
conductance calculation
We have also developed a similar order N method for the calculation of
Landauer-Büttiker conductance [42]. This formula has the advantages to
be general, independent of the dimensionality of the system and its eventual
geometrical and physico-chemical complications. To express an analytical
form for the conductance of the system under study, it is convenient to
dissociate it from the two metallic external electrodes that serve to inject
109
current in and out of the system, and that allow to suitably assume a wellgiven form for incoming and outgoing scattering states. These electrodes
are further connected to measurement reservoirs, in which energy dissipation finally takes place.
We have shown that the problem can be reduced to the evaluation of
hψ|Gr |ψi, where |ψi is a normalized state. To implement the recursion
approach in the Landauer framework [39], a generalization to the case of
non-symmetric matrices was performed. Indeed, the effective appearing
in the Green’s function is here non-symmetric hamiltonian H = HSys +
ΣrL + ΣrR because of the presence of self-energy matrix elements (ΣrL,R )
describing the finite system coupled to the electrodes. A generalization
of the recursion method to the computation of electronic spectra of nonsymmetric hamiltonian was unsuccessfully proposed [40], but its further
implementation received little consideration until the work of C. Benoit
et al., who applied a similar approach in a totally different context [41].
As we could demonstrate in [42], the bi-orthogonalization process can be
efficiently implemented to the Landauer- Büttiker formalism, which yields
a stable, efficient novel computing tool for nanoelectronics. The scheme has
corrected the initial errors made in [40]. The basis equations are as follows:
starting from the normalized vector |ψi and from the non-hermitian matrix
H, we construct a bi-orthogonal basis {|ψn i, hφn |} defined as follows
|ψn+1 i = H|ψn i − an+1 |ψn i − bn |ψn−1 i
(6)
hφn+1 | = hφn |H − hφn |an+1 − hφn−1 |bn
(7)
with the initial conditions |ψ−1 i = |φ−1 i = 0, |ψ0 i = |φ0 i = |ψi, and the
bi-orthogonality condition hφn |ψm i = 0 if n 6= m. This last condition is
equivalent to the following relations for an and bn :
an =
bn =
110
hφn |H|ψn i
hφn |ψn i
hφn−1 |H|ψn i
hφn |ψn i
=
.
hφn−1 |ψn−1 i
hφn−1 |ψn−1 i
(8)
(9)
The four equations (6), (7), (8), and (9) allow a recursive determination of
the bi-orthogonal basis and of the coefficients an , bn . Note that in ”ket”
notation, Eq. (7) must be understood as: |φn+1 i = H† |φn i − a∗n+1 |φn i −
b∗n |φn−1 i. One starts from |φ0 i = |ψ0 i = |ψi. At step 0, one computes
H|ψ0 i and a1 = hφ0 |H|ψ0 i/hφ0 |ψ0 i by expanding all the amplitudes within
the tight-binding localized basis. |ψ1 i and |φ1 i are then obtained by computing H|ψ0 i − a1 |ψ0 i and H† |φ0 i − a∗1 |φ0 i, while the first coefficient b1 is
subsequently deduced from Eq. (9). At step 1, H|ψ1 i is computed together
with a2 = hφ1 |Ĥ|ψ1 i/hφ1 |ψ1 i. Then |ψ2 i and |φ2 i result from the computation of vectors H|ψ1 i − a2 |ψ1 i − b1 |ψ0 i and H† |φ1 i − a∗2 |φ1 i − b∗1 |φ0 i.
Finally, the coefficient b2 is then deduced from Eq. (9). Steps n ≥ 2 are
fully similar to step 1. In the basis {|ψn i}, H has thus a tridiagonal form:

a1 b1
 1 a b
2
2


H=
1 a3 b3


1 . .
. .




.


(10)
Hence the recurrence relations (6) and (7) lead to a non-symmetric
matrix and to a non-normalized bi-orthogonal basis. By choosing a different
convention, a symmetric tridiagonal matrix and/or a normalized basis could
1
be obtained. The quantity hψ|Gr (z = E ± 0+ )|ψi = hφ0 | z−H
|ψ0 i can then
be computed by the continued fraction method. It is equal to the first
diagonal element of (z − H)−1 where H is the tridiagonal matrix (10).
Let us call G0 (z) this matrix element and define Gn (z), the first diagonal
element of the matrix (z − Hn )−1 , with Hn the matrix H without its n first
lines and columns:

an
bn
 1 a
n+1 bn+1


Hn = 
1
an+2 bn+2


1
.
.
.
.




.


(11)
111
From standard linear algebra, it can be shown that
G0 (z) =
1
,
z − a1 − b1 G1 (z)
(12)
and replicating such algorithm one gets a continued-fraction of G0 (z) :
G0 (z) =
1
z − a1 −
b1
z − a2 −
(13)
b2
...
In contrast with the standard recursion method, the recursion coefficients
an and bn do not show any simple behaviour for large n. In the type of
applications considered here, a simple truncation of the continuous fraction
at sufficiently large n gave a good convergence. This method was tested on
carbon nanotube based heterojunctions [42], with perfect agreement with
more traditional decimation techniques [43].
4
Conclusion
In conclusion, we have presented the methodologies developed for a decade,
to provide to scientific community with an efficient optimized order N approach, able to give information about basic transport length scales in complex materials and nanodevices. Current efforts are devoted to extend the
possibilities to situations far from the equilibrium, which stands as a great
challenge; this include simulating nanodevices submitted to large bias voltages, and dissipation effects within the conducting channels such as those
driven by electron-phonon interactions [44].
5
Acknowledgements
S.R deeply acknowledges Professor Takeo Fujiwara for continuous support
and encouragements during his stay in Tokyo from November 1996 till august 1999. The abovementioned computational methodologies have been
substantially improved and extended during this postdoctoral period, and
112
have strongly benefited from the scientific discussion and environment provided by the University of Tokyo and warm support and comprehension
from Fujiwara-Sensei, Yamamoto-san, Hoshi-san, Kobayashi-san and the
rest of the Fujiwara’s group.
113
References
[1] J. M. Soler, E. Artacho, J. Gale, A. Garcı́a, J. Junquera, P. Ordejón, and D.
Sánchez-Portal, J. Phys.: Condens. Matter 14, 2745-2779 (2002).
[2] R. Haydock, in Solid State Physics, edited by H. Ehrenreich, F. Seitz and D.
Turnbull, (Academic Press New-York, 1980), Vol. 35, 215.
[3] Recursion Method and its Applications, edited by D.G. Petitfor and D.L.
Weaire, Springer Series in Solid States Sciences, Vol. 58 (Springer Verlag,
Berlin, 1985). The Recursion Method: Application to Many-Body Dynamics,
edited by V.S. Viswanath and G. Müller, Lectures Notes in Physics, Vol. 23
(Springer Verlag, Berlin, 1994).
[4] S. Goedecker, Rev. of Mod. Phys. 71, 1085 (1999).
[5] T. Hoshi and Takeo Fujiwara, J. Phys. Soc. Jpn. 72, 2429 (2003). R.
Takayama, T. Hoshi, T. Sogabe, S.-L. Zhang, and Takeo Fujiwara, J. Phys.
Soc. Jpn 73, 1519 (2004). R. Takayama, T. Hoshi, T. Sogabe, S.-L. Zhang,
and Takeo Fujiwara, Phys. Rev. B 73, 165108, pp.1-9 (2006). Takeo Hoshi
and Takeo Fujiwara, J. Phys.: Condens. Matter 18 10787-10802 (2006).
[6] C. Lanczos, J. Res. Nat. Bur. Stand. 45, 255 (1950).
[7] D.J. Lohrmann et al, Phys. Rev. B 40, 8404 (1989).
[8] S.K. Bose, K. Winer, O.K. Andersen, Phys. Rev. B 37, 6262 (1988). S.K.
Bose et al, Phys. Rev. B 37, 9955 (1988). S.K. Bose et al, Phys. Rev. B 41,
7988 (1990).
[9] R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957). D. Fisher and P.A. Lee, Phys.
Rev. B 23, 6852 (1981).
[10] S. Roche and D. Mayou, Phys. Rev. Lett. 79, 2518 (1997).
[11] S. Roche, Phys. Rev. B 59, 2284 (1999).
[12] S. Roche and R. Saito, Phys. Rev. Lett. 87, 246803 (2001).
[13] A. Lherbier, M. Persson, Y.M. Niquet, F. Triozon and S. Roche, Phys. Rev.
Lett. submitted
114
[14] S. Roche t al. in : Understanding Carbon Nanotubes From Basics to Applications, Series: Lecture Notes in Physics (Springer), 677 A. Loiseau, P. Launois,
P. Petit, S. Roche, J.-P. Salvetat (Eds.) 2006.
[15] J.C. Charlier, X. Blase, and S. Roche, Rev. Mod. Phys. in print
[16] C.T. White and T.N. Todorov, Nature 393, 240 (1998).
[17] F. Triozon, S. Roche, A. Rubio, and D. Mayou, Phys. Rev. B 69, 121410(R)
(2004).
[18] S. Roche, G. Dresselhaus, M.S. Dresselhaus, and R. Saito, Phys. Rev. B 62,
16092 (2000).
[19] B. Stojetz, C. Miko, L. Forro, and Ch. Strünk, Phys. Rev. Lett. 94, 186802
(2005).
[20] D.J. Thouless, Phys. Rev. Lett. 39 1167 (1977).
[21] A. Aronov, Y. Sharvin, Rev. Mod. Phys. 59, 755 (1987).
[22] S. Chakraverty and A. Schmid, Phys. Rep. 140, 193 (1986).
[23] S. Roche, J. Jiang, F. Triozon, and R. Saito, Phys. Rev. Lett. 95, 076803
(2005).
S. Roche, J. Jiang, F. Triozon, and R. Saito, Phys. Rev. B 72, 113410 (2005).
[24] S. Roche, F. Triozon, A. Rubio and D. Mayou, Phys. Rev. B 64, 121401
(2001).
[25] A. Bachtold et al., 1999, Nature 397, 673. C. Schönenberger et al., Appl.
Phys. A 69 283 (1999).
[26] A. Fujiwara et al., Phys. Rev. B 60, 13492 (1999).
[27] G. Fedorov, B. Lassagne, M. Sagnes, B. Raquet, JM. Broto, F. Triozon, S.
Roche, and E. Flahaut, Phys. Rev. Lett. 94, 066801 (2005).
[28] Ch. Strünk, B. Stojetz and S. Roche, Semicond. Sci. Technol. 21 S38-S45
(2006).
[29] Ch. Ewels and M. Glerup, J. Nanosci. Nanotech. 5, 1354 (2005).
115
[30] H.J. Choi, J. Ihm, S.G. Louie and M. L. Cohen, Phys. Rev. Lett. 84, 2917
(2000).
[31] A.H. Nevidomskyy, G. Csanyi, M.C. Payne, Phys. Rev. Lett. 91, 105502
(2003).
[32] Ch. Adessi, S. Roche, and X. Blase, Phys. Rev. B 73, 125414 (2006).
[33] S. Latil, S. Roche, D. Mayou, and J.C. Charlier, Phys. Rev. Lett. 92, 256805
(2004).
S. Latil, F. Triozon, and S. Roche, Phys. Rev. Lett. 95, 126802 (2005).
[34] R. Avriller, S. Latil, F. Triozon, X. Blase, and S. Roche, Phys. Rev. B 74,
121406 (2006).
[35] S. Roche et al., Phys. Lett. A 285, 94 (2001). S. Wang, M. Grifoni, and S.
Roche Phys. Rev. B 74, 121407(R) (2006).
[36] S. Wang and M. Grifoni, Phys. Rev. Lett. 95, 266802 (2005).
[37] Takeo Fujiwara, in Physical Properties of Quasicrystals, Springer Series in
Solid-States Sciences Springer-Verlag, Berlin, (1998).
[38] S. Roche and Takeo Fujiwara, Phys. Rev. B 58, 11338 (1998).
[39] R. Landauer, IBM J. Res. Dev. 32, 306 (1988). M. Büttiker, IBM J. Res.
Dev. 32, 317 (1988).
[40] R. Haydock and M. J. Kelly, J. Phys. C. 8, L290 (1975).
[41] C. Benoit et al, Modelling Simul. Mater.Sci. Eng. 3, 161 (1995). C. Benoit
and G. Poussigue, Eur. Phys. J. AP 23, 117 (2003).
[42] F. Triozon and S. Roche, Eur. Phys. J. B 46, 427431 (2005).
[43] F. Triozon, Ph. Lambin, S. Roche, Nanotechnology 16, 230 (2005).
[44] L.E.F. Foa Torres and S. Roche, Phys. Rev. Lett. 97, 076804 (2006).
116
A method to evaluate electrostatic potential from
X-ray diffraction data by using MEM
Hiroshi TANAKA
Dept. of Materials Science, Shimane Univ.
Abstract
A new method is proposed evaluating the electrostatic potential
in the crystalline solids from the x-ray diffraction data. It is based
on the electron charge density analyzed by maximum entropy method
(MEM) and Ewald’s technique. The algorithm is almost parameterfree, and has little ambiguity in the practical calculation. In order to
demonstrate the efficiency of the method, it is applied to crystalline
NaCl.
1
Introduction
Electrostatic potential is a key property of condensed maters. In organic
materials, for example, we can presume active parts of molecules or proteins
from the view point of chemical reactions, which is important information
for analyzing mechanism of reaction, synthesizing new materials or drug
design. On the surface or interface of inorganic materials, we can expect
where molecules and / or atoms are absorbed. Some ferroelectric materials are stabilized by electrostatic potentials. Therefore, several approaches
117
have been proposed, which evaluate the electrostatic potential from experimental X-ray diffraction data.
For example, Bertaut [1] and Stewart [2] proposed a formalism to evaluate the potential directly from the experimental structure factors by the
Fourier transformation. The summation in reciprocal space, however, suffers from Gibbs’s oscillation, unless enough number of structure factors are
available to get the converged potential. Su and Coppens evaluated the
electrostatic potential of deuterated benzen [3] based on the multipole expansion [4], in which the electron charge density is expanded into spherical
harmonic functions. A merit of the method is that we can lift a molecule
out of the crystal [5]. Difficulties are complex procedure in evaluation and
ambiguity in how to redistribute the electron charge to each constitutive
atom.
In this letter, we propose an alternative method to evaluate the electrostatic potential by using the maximum entropy method (MEM) together
with Ewald’s technique [6]. The MEM analysis of X-ray diffraction data
was proposed by Collins [7]. The method reproduces high-resolution electronic charge density from a limited number of diffraction data, and has
been applied to wide variety of materials [8, 9]. It is now applicable to
huge systems such as proteins, coupled with the recent parallel processing
in computation and the fast Fourier transfer (FFT) technique [10]. Then
we thought of utilizing the electronic charge density obtained by the MEM
analysis for the evaluation of electrostatic potential.
118
2
Method
The electrostatic potential U (r) in the solid consists of the summation of
contributions from the nucleus charge and the electron charge,
Z
XX
Zt
ρ(r)
U (r) = Unuc (r) + Uele (r) =
−
dr 0 ,
0|
|r
−
l
−
R
|
|r
−
r
t
t
(1)
l
where Rt and Zt are position vector relative to the lattice point and atomic
number of the t-th basis atom. ρ(r) and l denote electron charge density
and lattice point vector, respectively. There are two difficulties in the practical calculation of Eq. (1). Firstly, the Coulomb potential is long range
and we need a large number of summation over lattice points for Unuc (r) (or
integration over large volume for Uele (r)). Secondly, each of Unuc (r) and
Uele (r) diverges to infinity and only the sum of them converges to a finite
value. To overcome these problems, several techniques have been proposed.
The electrostatic potential due to electron charge can be rewritten by
P
the summation Uele (r) = −4π G ρ̃(G) exp(iGr)/|G|2 in the reciprocal
lattice G, where ρ̃(G) is the Fourier compornent of the charge density ρ(r).
Since the structure factor F (G) is expressed by ρ̃(G) and the volume Ω in
real space as F (G) = Ωρ̃(G), the electrostatic potential due to electron
charge can be calculated by the equation;
Uele (r) = −4π
X F (G) exp(iGr)
G
Ω|G|2
.
(2)
It is however difficult to observe sufficient number of structure factors to
get a converged value for Eq. (2). Fortunately, this problem can be overcome by MEM analysis of the X-ray diffraction data. By using MEM, we
can extrapolate the values of structure factor up to as large value of |G|
119
as needed for the convergence of Eq. (2), which reproduce a smooth and
accurate charge density of electrons. Then, we employ these extrapolated
structure factors FM EM (G) instead of observed ones in the following calculation. Namely, we approximate
F (G) ' FM EM (G).
(3)
This is the key point of our method.
On the other hand, the electrostatic potential due to nucleus charge can
be evaluated ordinary Ewald’s method. According to the method, Unuc (r)
is rewritten as the summation in real and reciprocal lattices as,
X P Zt exp(−|G|2 /η 2 ) exp(−iGRt )
t
Unuc (r) = 4π
exp(iGr)
Ω|G|2
G
XX
Zt
+
erfc(η|r − l − Rt |),
|r − l − Rt |
t
(4)
l
where erfc(x) is the complementary error function. The first term on the
right hand side of Eq. (4) (the summation over the reciprocal lattice G)
converges rapidly to a finite value due to exp(−|G|2 /η 2 ). The second
term (the summation over the real lattice l) also converges rapidly due
to erfc(η|r − l − Rt |). Thus we can get a converged value for Eq. (4), with
a finite summation over real and reciprocal lattices. The variable η is the
only parameter employed in this method, which should be chosen so that
both summations in real and reciprocal lattices converged rapidly. The
converged value is however not so sensitive within a wide range of value of
η as shown in the later section.
From Eqs. (1)-Eq. (4), we get a formulation to evaluate electrostatic
120
potential;
X P Zt exp(−|G|2 /η 2 − iGRt ) − FM EM (G)
t
U (r) = 4π
exp(iGr)
Ω|G|2
G
XX
Zt
+
erfc(η|r − l − Rt |).
(5)
|r
−
l
−
R
|
t
t
l
Although there is a singularity at G = 0 in the first term on right hand side
of Eq. (5), it is removable because Zt exp(−|G|2 /η 2 −iGRt ) and FM EM (G)
cancells each other at G = 0
In the present formulation, the thermal vibration of nucleus is not considered. However it is straightforward to take the effect into account in
our formulation. It results in a slight modification of Ewald’s sum Eq.
(4), because the distribution of nucleus charge is given by a Gaussian type
function whether it is isotropic or not.
3
Results
Figure 1 shows (a) the electron charge density obtained by MEM analysis, and (b) the electrostatic potential analyzed by the present method for
crystalline NaCl. In this analysis, we used 19 structure factors obtained by
X-ray diffraction experiment. As it is shown in Fig. 1(a), Cl− ion is greater
than Na+ ion. This is because the former is negative and the later is positive ion, and a general aspect for ionic atoms. Corresponding to this, the
contour lines at high electrostatic potential area become circular around
each ion, and their diameter are larger around Cl− ion than around Na+
ion. However, the tail part of the potential (low potential area) spreads
wider area around Na+ ion than around Cl− ion. This is because the number of valence electrons is less than that of nucleus charge around Na+ ion,
121
Cl
Na
(a)
Cl
Na
(b)
Figure 1: Contour plots of (a) electron charge density, and (b) electrostatic
potential for crystalline NaCl on (100) surface. (a) Contour lines are plotted
from 0.1 e/Å2 to 2.0 e/Å2 with interval of 0.1 e/Å2 . (b) Contour lines are
plotted from -0.4 e/Å to 1.6 e/Å with interval of 0.2 e/Å.
and the valence electrons can not screen the potential due to nucleus charge
away from the nucleus.
Figure 2 shows the situation more explicitly. In this figure, an isosurface
of electron charge density is colored by the value of electrostatic potential.
As is shown, the electrostatic potential around the Na+ ion is higher than
that around the Cl− ion, corresponding to the fact that the former is a
positive ion.
Another advantage of our method is that we can also evaluate the elec-
122
Figure 2: Electrostatic potential of crystalline NaCl. The amplitude of
electrostatic potential is indicated by grey scale from -400 me/Å (black)
to 100 me/Å (white) on the isosurface of electron charge density at 250
me/Å2 .
123
Figure 3: Directions of electric field at each position are indicated by arrows
for crystalline NaCl on (100) surface, together with the contour plot of
electrostatic potential, which is same with Fig. 1(b).
tric field E(r) with the following analytic formulation.
X FM EM (G) − P Zt exp(−|G|2 /η 2 − iGRt ) exp(iGr) G
t
E(r) = 4πi
Ω|G|
|G|
G
µ
¶
XX
erfc(η|r − l − Rt |) 2η exp(−η 2 |r − l − Rt |2 )
√
+
Zt
+
|r − l − Rt |2
π|r − l − Rt |
t
l
×
r − l − Rt
|r − l − Rt |
(6)
The directions of electric field calculated by Eq. (6) are plotted on
the contour plot of the electrostatic potential in Fig. 3. It is clearly shown
that the electronic field is perpendicular to the contour lines of electrostatic
potential.
124
Figure 4: η dependence of the calculated potential. Changes in the electrostatic potentials at (0.25,0,0) are plotted as a function of the parameter η.
The CPU times consumed in the calculations are also plotted.
4
discussion and conclusion
In our method, the only parameter is η in Eq. (5). We investigated the
influence of the value of η on the results. Figure 4 shows the changes
in the electrostatic potential as a function of η. As shown in the figure,
the calculated potential is not so sensitive to the value of η within a wide
rage. For example, the potential has 5 digits accuracy for η > 0.4, and
furthermore 7 digits accuracy for η > 1.7. The CPU time needed for the
calculation is also plotted in the same figure. It decrease as increasing η and
becomes almost constant at η > 0.8. This is because the most of CPU time
is consumed for the summation in real lattice points in Eq. (5). The large η
however requires large number of extrapolated structure factors FM EM (G).
Then we should choose an appropriate value of η taking into account the
125
accuracy, CPU time, and number of FM EM (G) required.
In summary, we developed a new method evaluating the electrostatic
potential from the X-ray diffraction data obtained by experiments. It is
based on MEM charge density analysis, and almost parameter free. The
only parameter is η, which controls the convergence speed of calculation.
The final result is however little sensitive to the choice of η, and everybody
can evaluate the electrostatic potential without any experience. In order
to confirm the effectiveness of the method, we applied it to a typical ion
crystal NaCl. Other practical applications of the method will be published
elsewhere.
謝辞
藤原先生にはこれまで，色々とご指導をいただいてきました．最初は修士の
学生の時 で，私は LCAO の tight-binding で分子の計算をしていたのですが，
講義で 「tight-binding で適当な第一原理計算をすると，とんでもないこと
になる」と話さ れ，あわてて相談にうかがったのを覚えています．先生は少
し困った顔をされました が (いま思えば，指導教官ではなかったためでしょ
う)，それでも丁寧に色々と教えて くださいました．もう一つ私にとって印
象的だったのは，就職してから共著の論文を書いたときのこと です．先生が
校正された部分について，
「ここはどういう意味ですか」と質問する と，や
はり丁寧に説明された後，消し始めました．私が「なぜ消すんですか」と質
問 すると，
「田中君に分からなかったから (きっと他の人にも分かりづらいで
しょ)」と 答えられました．私の理解力不足のせいだったと思うのですが，ず
いぶん謙虚な考え 方をされるんだなあと思いました．どちらも藤原先生は覚
えていらっしゃらないと思いますが，こうした経験は私にとっ て今日学生を
126
指導する上で大きく影響していると感じます．教育者として先生をまだ 越え
られていないと思いますが，先生と同じように学生を指導していることに気
づく ことがあります．最後に，今後とも研究ならびに教育の両面でご活躍さ
れることを期待いたしております．
References
[1] Bertaut E F 1987 J. Phys. Chem. Solids 39 97
[2] Stewart R F 1976 Chem. Phys. Lett. 65 335
[3] Su Z and Coppens P 1992 Acta Crystallogr. A48 188
[4] Hansen N K and Coppens P 1978 Acta Crystallogr. A34 909
[5] Coppens P, Su Z and Becker P J Chapter 1999 International Tables for
Crystallography Vol. C Eds. Willson A J C and Prince E (Dordrecht:
Kluwer Academic Publishers) Second Edition, Chapter 8.7. p.713
[6] Ewald P P 1921 Ann. Physik 64 253
[7] Collins D M 1982 Nature 298 49
[8] Sakata M and Sato M 1990 Acta Crystallogr. 46 66
[9] Takata M, Nishibori E and Sakata M 2001 216 71, and its references
[10] Tanaka H, Takata M and Sakata M 2002 J. Phys. Soc. Jpn 71 2595
127
Effect of transition-metal elements on the electronic
properties of quasicrystals and complex aluminides
G. Trambly de Laissardière1 and D. Mayou2
1
Laboratoire de Physique Théorique et Modélisation, CNRS and
Université de Cergy–Pontoise, 95302 Cergy–Pontoise, France
2
Institut Néel, CNRS and Université Joseph Fourier, Bât D,
B.P. 166, 38042 Grenoble Cedex 9, France
1
Introduction
It is with great pleasure that we contribute to this book in honor of Prof.
Takeo Fujiwara. GTL enjoyed eighteen months of Prof. Fujiwara’s hospitality at the University of Tokyo during the early 1990’s. At that time the
work of Prof. Fujiwara in the field of electronic structure of quasicrystals
had already made a major contribution to the literature (see for instance
[1]). Since that time our research owes much to his work.
Prof. Fujiwara was the first who performed realistic calculations of the
electronic structure in quasicrystalline materials without adjustable parameters (ab-initio calculations) [2]. Indeed these complex alloys [3] have very
exotic physical properties (see Refs. [4, 5] and Refs therein), and it rapidly
appeared that realistic calculations on the actual quasicrystalline materials are necessary to understand the physical mechanism that govern this
128
properties. In particular, these calculations allow to analyze numerically
the role of transition-metal elements which is essential in those materials.
In this paper, we briefly present our work on the role of transition-metal
element in electronic structure and transport properties of quasicrystals
and related complex phases. Several Parts of these works have been done
or initiated in collaboration with Prof. Fujiwara.
2
2.1
Electronic structure
Ab-initio determination of the density of states
A way to study the electronic structure of quasicrystal is to consider the case
of approximants. Approximants are crystallines phases, with very large unit
cell, which reproduce the atomic order of quasicrystals locally. Experiments
indicate that approximant phases, like α-AlMnSi, α-AlCuFeSi, R-AlCuFe,
etc., have transport properties similar to those of quasicrystals [4, 6]. In
1989 and 1991, Prof. Fujiwara performed the first numerical calculations of
the electronic structure in realistic approximants of quasicrystals [2, 7, 8].
He showed that their density of states (DOS, see figure 1) is characterized by
a depletion near the Fermi energy EF , called “pseudo-gap”, in agreement
with experimental results (for review see Ref. [4, 9, 18]) and a HumeRothery stabilization [10, 11]. The electronic structure of simpler crystals
such as orthorhombic Al6 Mn, cubic Al12 Mn, present also a pseudo-gap near
EF which is less pronounced than in complex approximants phases (figure
1) [11].
129
DOS (states/eV.cell)
12
8
4
DOS (states/eV.cell)
0
-12 -10
120
-8
-6
-4
-2
-Al69.6Si13Mn17.4
0
2
4
2
4
EF
80
40
0
-12 -10
Figure 1:
EF
Al6Mn
-8
-6 -4 -2
Energy (eV)
0
Ab-initio total DOS of Al6 Mn (simple crystal) and α-
Al69.6 Si13.0 Mn17.4 (approximant of icosahedral quasicrystals) [11, 12].
2.2
Models to analyze the role of transition-metal element
sp–d hybridization model
The role of the transition-metal (TM, TM = Ti, Cr, Mn, Fe, Co, Ni)
elements in the pseudo-gap formation has been shown from experiments,
ab-initio calculations and model analysis [4,13–19,11]. Indeed the formation
of the pseudo-gap results from a strong sp–d coupling associated to an
ordered sub-lattice of TM atoms [19, 11]. Consequently, the electronic
structure, the magnetic properties and the stability, depend strongly on
130
the TM positions, as was shown from ab-initio calculations [28–33,20,21].
How an effective TM–TM interaction induces stability?
Just as for Hume-Rothery phases a description of the band energy can be
made in terms of pair interactions (figure 2) [17, 19]. Indeed, it has been
shown that an effective medium-range Mn–Mn interaction mediated by the
sp(Al)–d(Mn) hybridization plays a determinant role in the occurrence of
the pseudo-gap [19]. We have shown that this interaction, up to distances
10–20 Å, is essential in stabilizing these phases, since it can create a HumeRothery pseudo-gap close to EF . The band energy is then minimized as
shown on figure 3 [20, 11].
The effect of these effective Mn–Mn interactions has been also studied
by several groups [17, 20, 21] (see also Refs in [11]). It has also explained
the origin of large vacancies in the hexagonal β-Al9 Mn3 Si and ϕ-Al10 Mn3
phases on some sites, whereas equivalent sites are occupied by Mn in µAl4.12 Mn and λ-Al4 Mn, and by Co in Al5 Co2 [20]. On the other hand,
an spin-polarized effective Mn–Mn interaction is also determinant for the
existence (or not) of magnetic moments in AlMn quasicrystals and approximants [21, 22, 32].
The analysis can be applied to any Al(rich)-Mn phases, where a small
number of Mn atoms are embedded in the free electron like Al matrix. The
studied effects are not specific to quasicrystals and their approximants, but
they are more important for those alloys. Such a Hume-Rothery stabilization, governed by the effective medium-range Mn–Mn interaction, might
therefore be intrinsically linked to the emergence of quasi-periodicity in
Al(rich)-Mn system.
131
(eV / Mn atom)
with repulsive term
without repulsive term
0.04
−ar
repulsive term : b e
0.02
Mn-Mn
0
-0.02
r = 6.7 A
-0.04
r = 4.8 A
2
3
4
5
6
7
8
9
10 11 12 13
Mn−Mn distance r (A)
Figure 2: Effective medium-range Mn–Mn interaction between two nonmagnetic manganese atoms in a free electron matrix which models aluminum atoms. [11]
Energy (eV / Mn atom)
0
Al6Mn
-AlMnSi
-AlMnSi
-0.1
-0.2
-0.3
0
4
8
12
16
20 24
L (A)
28
32
36
40
Figure 3: Variation of the band energy due to the effective Mn–Mn interaction in o-Al6 Mn, α-AlMnSi and β-Al9 Mn3 Si. [20]
132
Cluster Virtual Bound states
One of the main results of the ab-initio calculations performed by Prof.
Fujiwara for realistic approximant phases, is the small energy dispersion
of electrons in the reciprocal space. Consequently, the density of states of
approximants is characterized by “spiky” peaks [2, 7, 8, 28]. In order to
analyze the origin of this spiky structure of the DOS, we developed a model
that show a new kind of localization by atomic cluster [23].
As for the local atomic order, one of the characteristics of the quasicrystals and approximants is the occurrence of atomic clusters on a scale
of 10–30 Å [25]. The role of clusters has been much debated in particular
by C. Janot [24] and G. Trambly de Laissardière [23]. Our model is based
on a standard description of inter-metallic alloys. Considering the cluster
embedded in a metallic medium, the variation ∆n(E) of the DOS due to
the cluster is calculated. For electrons, which have energy in the vicinity of
the Fermi level, transition atoms (such as Mn and Fe) are strong scatters
whereas Al atoms are weak scatters. In the figure 4 the variation, ∆n(E),
of the density of states due to different clusters are shown. The Mn icosahedron is the actual Mn icosahedron of the α-AlMnSi approximant. As an
example of a larger cluster, we consider one icosahedron of Mn icosahedra.
∆n(E) of clusters exhibits strong deviations from the Virtual Bound
States (1 Mn atom) [26]. Indeed several peaks and shoulders appear. The
width of the most narrow peaks (50 − 100 meV) are comparable to the fine
peaks of the calculated DOS in the approximants (figure 1). Each peak
indicates a resonance due to the scattering by the cluster. These peaks
correspond to states “localized” by the icosahedron or the icosahedron of
icosahedra. They are not eigenstate, they have finite lifetime of the order
133
n(E) (states / (eV transition atom))
8
7
1 Mn atom
6
E
5
1 Mn
icosahedron
4
3
1 icosahedron
of 12 Mn
icosahedra
2
1
0
6
7
8
9 10 11 12 13 14 15
Energy E (eV)
Figure 4: Variation ∆n(E) of the DOS due to Mn atoms. Mn atoms are
embedded in a metallic medium (Al matrix). From [23].
of ~/δE, where δE is the width of the peak. Therefore, the stronger the
effect of the localization by cluster is, the narrower is the peak. A large
lifetime is the proof of a localization, but in the real space these states have
a quite large extension on length scale of the cluster.
The physical origin of these states can be understood as follows. Electrons are scattered by the Mn atoms of a cluster. By an effect similar to
that of a Faraday cage, electrons can by confined by the cluster provided
that their wavelength λ satisfies λ & l, where l is the distance between two
Mn spheres. Consequently, we expect to observe such a confinement by the
cluster. This effect is a multiple scattering effect, and it is not due to an
overlap between d-orbitals because Mn atoms are not first neighbor.
134
3
Transport properties
Quasicrystals have many fascinating electronic properties, and in particular
quasicrystals with high structural quality, such as the icosahedral AlCuFe
and AlPdMn alloys, have unconventional conduction properties when compared with standard inter-metallic alloys. Their conductivities can be as
low as 150–200 (Ω cm)−1 (see Refs. [4, 5, 27] and Refs. therein). Furthermore the conductivity increases with disorder and with temperature,
a behavior just at the opposite of that of standard metal. In a sense the
most striking property is the so-called “inverse Mathiessen rule” according
to which the increases of conductivity due to different sources of disorder
seems to be additive. This is just the opposite that happens with normal
metals where the increases of resistivity due to several sources of scattering are additive. An important result is also that many approximants of
these quasicrystalline phases have similar conduction properties. For example the crystalline α-AlMnSi phase with a unit cell size of about 12 Å
and 138 atoms in the unit cell has a conductivity of about 300 (Ω cm)−1 at
low temperature [4].
3.1
Small Boltzmann velocity
Prof. Fujiwara et al. was the first to show that the electronic structure
of AlTM approximants and related phases is characterized by two energy
scales [2, 7, 8, 28, 29] (see previous section). The largest energy scale, of
about 0.5−1 eV, is the width of the pseudogap near the Fermi energy EF . It
is related to the Hume–Rothery stabilization via the scattering of electrons
by the TM sub-lattice because of a strong sp–d hybridization. The smallest
135
1e+5
Doped semi-conductors
Mott
"Perfect" stable quasicrystals
(i-AlCuFe and i- AlPdMn)
1000
"Imperfect" stable quasicrystals
(i-AlLiCu)
cm)
1e+4
(
Resistivity
100
Metastable quasicrystals (i-AlMn),
Amorphous alloys
10
Metallic alloys
1
4K
Temperature
300 K
Figure 5: Schematic temperature dependencies of the experimental resistivity of quasicrystals, amorphous and metallic crystals.
cm)
1e+4
-Al69.6Si13.0Mn17.4
1000
Figure 6:
(
100
Resistivity
10
trical
Al12Mn
inverse
Ab-initio elec-
resistivity
scattering
versus
time,
in cubic approximant α-
Al (f.c.c.)
Al69.6 Si13.0 Mn17.4 , pure Al
1
(f.c.c.), and cubic Al12 Mn.
1e+14
(s−1)
Temperature
136
1e+15
energy scale, less than 0.1 eV, is characteristic of the small dispersion of
the band energy E(k). This energy scale seems more specific to phases
related to the quasi-periodicity. The first consequence on transport is a
small velocity at Fermi energy, Boltzmann velocity, VB = (∂E/∂k)E=EF .
From numerical calculations, Prof. Fujiwara et al. evaluated the Bloch–
Boltzmann dc conductivity σB in the relaxation time approximation. With
a realistic value of scattering time, τ ∼ 10−14 s [27], one obtains σB ∼
10 − 150 (Ωcm)−1 for a α-AlMn model [8] and 1/1-AlFeCu model [28].
This corresponds to the measured values [4, 6], which are anomalously
low for metallic alloys. For decagonal approximant the anisotropy found
experimentally in the conductivity is also reproduced correctly [29].
3.2
Quantum transport in Quasicrystals and approximants
The semi-classical Bloch–Boltzmann description of transport gives interesting results for the intra-band conductivity in crystalline approximants,
but it is insufficient to take into account many aspects due to the special
localization of electrons by the quasi-periodicity (see Refs. [34–43] and
Refs. therein). Some specific transport mechanisms like the temperature
dependence of the conductivity (inverse Mathiessen rule, the defects influence, the proximity of a metal / insulator transition), require to go beyond
a Bloch–Boltzmann analysis. Thus, it appears that in quasicrystals and related complex metallic alloys a new type of breakdown of the semi-classical
Bloch-Boltzmann theory operates. In the literature, two different unconventional transport mechanisms have been proposed for these materials.
Transport could be dominated, for short relaxation time τ by hopping between “critical localized states”, whereas for long time τ the regime could
137
be dominated by non-ballistic propagation of wave packets between two
scattering events.
We develop a theory of quantum transport that applies to a normal ballistic law but also to these specific diffusion laws. As we show phenomenological models based on this theory describe correctly the experimental
transport properties [41, 42, 43] (compare figures 5 and 6).
3.3
Ab-initio calculations of quantum transport
According to the Einstein relation the conductivity σ depends on the diffusivity D(E) of electrons of energy E and the density of states n(E) (summing the spin up and spin down contribution). We assume that n(E) and
D(E) vary weakly on the thermal energy scale kT , which is justified here.
In that case, the Einstein formula writes
σ = e2 n(EF )D(EF )
(1)
where EF is the chemical potential and e is the electronic charge. The
temperature dependence of σ is due to the variation of the diffusivity D(EF )
with temperature. The central quantity is thus the diffusivity which is
related to quantum diffusion. Within the relaxation time approximation,
the diffusivity is written [41]
Z
1 +∞
D(E) =
C0 (E, t) e−|t|/τ dt
(2)
2 0
D
E
where C0 (E, t) = Vx (t)Vx (0) + Vx (0)Vx (t)
it the velocity correlation
E
functions without disorder, and τ is the relaxation time. Here, the effect of
defects and temperature (scattering by phonons ...) is taken into account
through the relaxation time τ . τ decreases as disorder increases. In the
138
case of crystals phases (such as approximants of quasicrystals), one obtains
[42, 43]:
σ
σB = e2 n(EF ) VB2 τ
=
σB + σNB
and σNB = e2 n(EF )
(3)
L2 (τ )
τ
(4)
where σB is actual the Bolzmann contribution to the conductivity and σNB
a non-Boltzmann contribution. L2 (τ ) is smaller than the square of the
unit cell size L0 . L2 (τ ) can be calculated numerically for the ab-initio
electronic structure [42]. From (3) and (4), it is clear that the Bolzmann
term dominates when L0 ¿ VB τ : The diffusion of electrons is then ballistic,
which is the case in normal metallic crystals. But, when L0 ' VB τ , i.e.
when the Bolzmann velocity VB is very low, the non-Bolzmann term is
essential. In the case of α-Al69.6 Si13.0 Mn17.4 approximant (figure 7) [42],
with realistic value of τ (τ equals a few 10−14 s [27]), σNB dominates and
σ increases when 1/τ increases, i.e. when defects or temperature increases,
in agreement with experimental measurement (compare figures 5 and 6).
To evaluate the effect of TM elements on the conductivity, we have
considered an hypothetical α-Al69.6 Si13.0 Cu17.4 constructed by putting Cu
atoms in place of Mn atoms in the actual α-Al69.6 Si13.0 Mn17.4 structure.
Cu atoms have almost the same number of sp electrons as Mn atoms, but
their d DOS is very small at EF . Therefore in α-Al69.6 Si13.0 Cu17.4 , the
effect of sp(Al)–d(TM) hybridization on electronic states with energy near
EF is very small. As a result, the pseudogap disappears in total DOS, and
the conductivity is now ballistic (metallic), σ ' σB , as shown on figure 8.
139
cm)
800
600
−1
1000
B
Conductivity
(
NB
400
200
0
1e+14
Figure 7:
2e+14
3e+14
4e+14
−1
(s )
Ab-initio dc-conductivity σ in cubic approximant α-
Conductivity
(
4e+05
cm)
−1
Al69.6 Si13.0 Mn17.4 versus inverse scattering time. [42]
3e+05
B
NB
2e+05
1e+05
0
1e+14
2e+14
3e+14
4e+14
−1
(s )
Figure 8: Ab-initio dc-conductivity σ in an hypothetical cubic approximant
α-Al69.6 Si13.0 Cu17.4 versus inverse scattering time. [43]
4
Conclusion
In this article we present the effect of transition-metal atoms on the physical
properties of quasicrystals and related complex phases. These studies lead
140
to consider these aluminides as spd electron phases [11], where a specific
electronic structure governs stability, magnetism and quantum transport
properties. The principal aspects of this new physics are now understood
particularly thanks to seminal work of Prof. T. Fujiwara and subsequent
developpements of his ideas.
References
[1] Fujiwara T, Tsunetsugu H. In: Di Vincenxo DP, Steinhart PJ, editors,
Quasicrystals: The states of the art, Singapore: World Scientific, 1991.
[2] Fujiwara T. Phys Rev 1989;B40:942.
[3] Shechtman D, Blech I, Gratias D, Cahn JW. Phys Rev Lett
1984;53:1951.
[4] Berger C. In: Hippert F, Gratias D, editors. Lecture on Quasicrystals.
Les Ulis: Les Editions de Physique, 1994; p. 463.
[5] Grenet T. In: Belin-Ferré E, Berger C, Quiquandon M, Sadoc A, editors. Quasicrystals: Current Topics. Singapor: World Scientific, 2000;
p. 455.
[6] Quivy A, Quiquandon M, Calvayrac Y, Faudot F, Gratias D, Berger C,
Brand RA, Simonet V, Hippert. J Phys Condens Matter 1996;8:4223.
[7] Fujiwara T, Yokokawa T. Phys Rev Lett 1991;66:333.
[8] Fujiwara T, Yamamoto S, Trambly de Laissardière G. Phys Rev Lett
1993;71:4166. Mat Sci Forum 1994;150-151:387.
141
[9] Mizutani U, Takeuchi T, Sato H. J Phys:
Condens Matter
2002;14:R767.
[10] Massalski TB, Mizutani U. Prog Mater Sci 1978;22:151.
[11] Trambly de Laissardière G, Nguyen Manh D, Mayou D, Prog Mater
Sci 2005;50:679.
[12] Zijlstra ES, Bose SK. Phys Rev 2003;B67:224204.
[13] Dankházi Z, Trambly de Laissardière G, Nguyen–Manh D, Belin E,
Mayou D. J Phys: Condens Matter 1993;5:3339.
[14] Trambly de Laissardière G, Mayou D, Nguyen Manh D. Europhys Lett
1993;21:25. J Non-Cryst Solids 1993;153-154:430. Trambly de Laissardière G, et al. Phys Rev 1995;B52:7920.
[15] Berger C, Belin E, Mayou D. Annales de Chimie-Science des Matériaux
1993;18:485.
[16] Mayou D, Cyrot–Lackmann F, Trambly de Laissardière G, Klein T. J
Non-Cryst Solids 1993;153-154:412.
[17] Zou J, Carlsson AE. Phys Rev Lett 1993;70:3748.
[18] Belin-Ferré E. J Non-Cryst Solids 2004;334-335:323.
[19] Trambly de Laissardière G, Nguyen Manh D, Mayou D. J Non-Cryst
Solids 2004;334-335:347.
[20] Trambly de Laissardière G. Phys Rev 2003;B68:045117.
[21] Trambly de Laissardière G, Mayou D. Phys Rev Lett 2000;85:3273.
142
[22] Simonet V, Hippert F, Audier M, Trambly de Laissardière G. Phys
Rev 1998;B58:R8865.
[23] Trambly de Laissardière G, Mayou M. Phys Rev 1997;B55:2890. Trambly de Laissardière G, Roche S, Mayou D. Mat Sci Eng 1997;A226228:986.
[24] Janot C, de Boissieu M. Phys Rev Lett 1994;72:1674.
[25] Gratias D, Puyraimond F, Quiquandon M, Katz A. Phys Rev
2000;B63:24202.
[26] Friedel J. Can J Phys 1956;34:1190.
Anderson PW. Phys Rev 1961;124:41.
[27] Mayou D, Berger C, Cyrot–Lackmann F, Klein T, Lanco P. Phys Rev
Lett 1993;70:3915.
[28] Trambly de Laissardière G, Fujiwara T. Phys Rev 1994;B50:5999.
[29] Trambly de Laissardière G, Fujiwara T. Phys Rev 1994;B50:9843.
Mat Sci Eng 1994;A181-182:722.
[30] Hafner J, Krajčı́ M. Phys Rev 1998;B57:2849.
[31] Krajčı́ M, Hafner J. Phys Rev 1998;B58:14110.
[32] Nguyen–Manh D, Trambly de Laissardière G. J Mag Mag Mater
2003;262:496.
[33] Zijlstra ES, Bose SK, Klanjšek M, Jeglič P, Dolinšek J. Phys Rev
2005;B72:174206.
143
[34] Tokihiro T, Fujiwara T, Arai M. Phys. Rev 1988;B38:5981.
[35] Fujiwara T, Mitsui T, Yamamoto S. Phys Rev B 1996;53,R2910.
[36] Roche S, Trambly de Laissardière G, Mayou D. J Math Phys
1996;38:1794.
[37] Roche S, Mayou D. Phys Rev Lett 1997;79:2518.
[38] Mayou D. In: Belin-Ferré E, Berger C, Quiquandon M, Sadoc A, editors. Quasicrystals: Current Topics. Singapor: World Scientific, 2000;
p. 412.
[39] Triozon F, Vidal J, Mosseri R, Mayou D. Phys Rev 2002;B65:220202.
[40] Bellissard J. In: Garbaczeski P, Olkieicz R, editors. Dynamics of Dissipation, Lecture Notes in Physics. Berlin: Springer, 2003; p. 413.
[41] Mayou D. Phys Rev Lett 2000;85:1290.
[42] Trambly de Laissardière G, Julien JP, Mayou D. Phys Rev Lett
2006;97:026601.
[43] Mayou D, Trambly de Laissardière G. In: Fujiwara T, Y. Ishii Y,
editors. Physics of Quasicrystals. Series “Handbook of Metal Physics”.
Elsevier Science, 2007. to appear
144
計算物理と計算機科学との距離
東京大学 大学院 物理工学専攻
山元進
長年にわたる藤原先生のご指導に深く感謝致しつつ, 先生の東京大学・
大学院・工学系研究科・物理工学専攻ご退職を記念してこの記事を寄稿さ
せていただきます.
概要
電子-電子相互作用が強い系の研究の中で, 大規模疎行列の対角化が必要
となった. 対角化コードを設計する段階で, 大規模疎行列の格納に必要な
メモリ領域が大きすぎることが問題となった. 設計時に, この問題の解決と
既存ライブラリの利用とは両立しないと考えたため, 独自コードを作成し
た. 最近, 反復解法の部分については既存ライブラリが利用できることを
知った. そのライブラリを捜し出すキーワードは Reverse Communication
Interface という我々には馴染のない言葉だった. この経験を踏まえ, 計算
物理と計算機科学との交流をもっと深めれば, 計算物理の専門家は計算機
科学の成果をもっと簡単に採り入れることができること, 計算機科学の専
145
門家は計算機ユーザーの要望から興味深い問題を拾い上げられることを論
じた.
1
背景
現在の我々の研究の紹介を兼ねて, 問題の背景となる系を詳しく説明す
る. 我々は電子-電子相互作用が強い (強相関) 系が示す秩序の多様さに興
味があり, 同じサイト上にいる電子同士の強いクーロン相互作用を取り入
れた LSDA+U 法という方法を用いて研究を行っている. LSDA+U 法は,
多くの強相関物質においての軌道・スピン・電荷秩序を記述することに成
功している. 研究対象となったのは, La2−x Srx NiO4 (LSNO) という物質
である. 回折実験などから, LSNO(x = 31 , 12 ) の電子配置として, 図 1 のよ
うな構造が提案されている. [1] また, これらは反強磁性絶縁体であること
が知られている. LSDA+U 法による計算では, LSNO(x = 13 ) について図
1(a) と整合する電子構造を得ることができたが, LSNO(x = 21 ) について
は電荷が秩序化していない金属の電子構造を得ることしかできなかった.
様々な検討の結果, 電荷の秩序化にはサイト上の電子間のクーロン相互作
用に加えて, 更に長距離のクーロン相互作用 (隣接サイトに存在する電子
とのクーロン相互作用) を採り入れる必要があること, 反強磁性状態が正
しく計算上の基底状態になるためには, 様々な電子配置を採り入れた計算
(複数のスレーター行列式を用いた計算) を行う必要があること, の 2 点が
明らかとなった. これを実現するために拡張ハバード模型を厳密対角化で
解くことにした. 拡張ハバードモデル・ハミルトニアンは次の様なもので
146
図 1: 層状ペロブスカイト酸化物 La2−x Srx NiO4 で見られるストライプ秩
序. (a)x =
1
3
の場合. (b)x =
1
3
の場合. 白丸はホール, 上下の矢印はそれ
ぞれの向きのスピンを表す. 点線は単位胞.
ある.
Ĥ =
X
tiαjβ ĉ†iασ ĉjβσ +
i,j,α,β,σ
+
1
2
X
εiα ĉ†iασ ĉiασ
i,α,σ
X
Uαβγδ ĉ†iασ ĉ†iβσ0 ĉiδσ0 ĉiγσ
i,α,β,γ,δ,σ,σ 0
1
+ V
2
X
ĉ†iασ ĉiασ ĉ†jβσ0 ĉjβσ0 .
(1)
<i,j>,α,β,σ,σ 0
tiαjβ はホッピングの行列要素, Uαβγδ は同じサイト上にある電子間のクー
ロン相互作用の行列要素, εiα はサイト上の各軌道の準位の分裂を表す項,
V は隣接する 2 サイトにそれぞれ局在する 2 電子間のクーロン反発を表
す項である. 採り入れるべき軌道が Ni イオンの eg 軌道 2 つであること
は, それまでの LSDA+U 法による計算から分かっていた. 実験を基に提
案されている図 1(b) の単位胞は 4 つの Ni サイトを含む. これを含む最
147
小の正方格子である
√
8×
√
8 超格子の系で計算をすることにした. 多電
子系のハミルトニアンを対角化するとき, 占有数表示を用いて電子配置と
2-進数とを 1-1 対応させると良いことが知られている. 今, (スピンの自由
度 2)× (軌道の自由度 2)× (サイト数 8) = 32 なので, 32bit 整数を使う
と丁度一杯の計算となる. この系に 12 電子を入れると x =
1
2
に対応す
る電子密度となる. スピン空間での回転に対するハミルトニアンの不変性
から, ↑ 電子と ↓ 電子数が等しい状態だけを考えれば良く, 考えるべき電
子配置の数は (16 C6 )2 = 80082 = 64128064 となる. 従って, 基底状態に
ついての知見を得るためには, 約 6400 万次元の大規模疎行列の最低固有
値とそれに属する固有ベクトルを求めれば良い. 対角化によって得られた
LSNO の電子構造に対する知見は, 別の場所で論じるのでそちらを参照し
ていただきたい.[2] ここではコード開発上の問題を我々の経験に則して論
じる.
2
ハミルトニアンの性質を用いたデータ量の圧縮
計算を始める前には, 大規模疎行列の固有値問題の解法には既存の並列
化されたライブラリを利用できると考えていた. しかし, 現実にはライブ
ラリは利用できなかった. それは, 行列のサイズがあまりにも大きいため
にハミルトニアンの格納方法を工夫せねばならず, その工夫がライブラリ
の利用を妨げたためである. 本節と次節でこのことを説明する.
6400 万次元の倍精度実数のベクトルを格納するには, およそ 0.5GB の
メモリを必要とする. 式 (1) のハミルトニアンの行列要素は状態ベクトル
1 次元あたりに 70 程度の非零要素を含む疎行列である. 従ってハミルト
148
ニアンの非零要素を格納するだけで 35 GB ものメモリ領域を必要とする.
非零要素の位置情報を含めれば更に多くのメモリ領域が必要である. この
様な巨大なメモリ領域を確保できる計算機は非常に限られる. ハミルトニ
アンを格納をあきらめ, ハミルトニアンの非零要素を生成しながら状態ベ
クトルに作用させることは可能だが, 時間が余計にかかるので好ましくな
い. ハミルトニアンの性質を用いてデータの量を圧縮するには, 状態ベク
トル空間を ↑ 電子の配置と ↓ 電子の配置との直積に分解すると良い.
1
(状態ベクトル空間) = (6 個の ↑ 電子の配置) ⊗ (6 個の ↓ 電子の配置)
(2)
これに従い, ハミルトニアンは下のような 3 つの部分に分解される.
H = h↑ ⊗ 1 + 1 ⊗ h↓ + h↑↓
(3)
状態ベクトル ~
u に対する h↑ ⊗ 1 + 1 ⊗ h↓ の作用を式で書けば,
uij =
X
h↑ii0 ui0 j +
X
h↓jj 0 uij 0
(4)
j0
i0
となる. ここで i, i0 は ↑ 電子配置のインデックス, j, j 0 は ↓ 電子配置のイ
ンデックスである. 先に示した式 (1) のハミルトニアンの 1 体部分にはス
ピンをフリップする要素がなく, 全て h↑ か h↓ に分類される. h↑ および
h↓ が密行列だとしても, 高々 80082 × 2 ≈ 1 億 3 千万個の要素にしかなら
ない. 実際には h↑ , h↓ とも疎行列で, 非零要素数はずっと少ない. h↑↓ の
項は直積で表せない残りの項で, 以下のように作用する
uij =
X
0 0
h↑↓
ij;i0 j 0 ui j
(5)
i0 ,j 0
1
この他に空間群を利用して状態ベクトルの次元を減らすことも可能だが, 実装が難し
いことと, 対称性を崩した計算を行う予定があることとの 2 つの理由から行わなかった.
149
式 (1) のハミルトニアン中では, このタイプの項は 2 体の部分からしか表
れない. 2 体部分の非零要素数は 1 体部分のものに比べて数が少ない. 従っ
て, 1 体部分の非零要素を格納するのに必要なメモリ容量が減るこの方法
は, ハミルトニアンのデータ圧縮に有効である. 実測によれば, この方法
を用いたとき (正確には, ハミルトニアンが実対称行列であることも利用
したとき) に h↑ , h↓ , h↑↓ の非零要素とその位置情報とを格納するのに必要
なメモリ容量は, およそ 1.5GB であった. 先ほどの, 非零要素の格納だけ
で 35 GB を要するという見積からすれば 20 倍以上の圧縮率である.
3
ライブラリ利用と独自データ構造のトレード-オフ
前節で述べた方法は計算に必要なメモリ領域が大きすぎるという問題を
解決した. その一方で, 既存のライブラリが利用できなくなるという問題
を引き起こした. 大規模疎行列の対角化ルーチンの論理構造は, 線形演算
子をベクトルに作用させる (ここではこれを基本演算と呼ぶことにする)
部分と, それを呼び出す (基本演算部分の上位ルーチンにあたる) 反復法部
分とに分けられる. ハミルトニアンを独自形式で格納する以上, 基本演算
部分は独自に実装する必要がある. 論理的には, 反復法部分のアルゴリズ
ムは基本演算の内部構造には依存しない (反復法部分には, 演算結果のベ
クトルのみが必要) ので, この部分にはライブラリを利用可能なように思
える. しかし, 多くのライブラリでは上位ルーチンの呼出法がユーザに公
開され, 下位ルーチンの実装が隠蔽される. このため下位ルーチンだけを
差し替えるのは困難であり, 結局全ての部分を独自に実装することとなっ
た. これによって, コーディングの時間だけではなく実行時間でも損をす
150
ることとなった. なぜなら我々の用意したルーチンは専門家の書いたもの
よりも効率が劣るからである. 我々は Lanczos 原理によってハミルトニ
アンを 3 重対角化して行き, 適当な次元で打ち切って得られた小行列を対
角化する方法 (Lanczos 法) によって元のハミルトニアンの近似的な固有
値固有ベクトルを得ている. Lanczos 原理による 3 重対角化は数値誤差
を考慮しなければ直交変換であるが, 実際には数値誤差によって基底の直
交性が容易に崩れる.[3] それを避けるために様々な手法が知られているが,
我々は, 適当な回数で計算をリスタートするという, もっとも原始的かつ
効率の低い方法を用いている.
我々は, 下位ルーチンの実装は隠蔽されるという前提のもと, このトレー
ド-オフの本質は, ライブラリのユーザー (我々) は問題とするハミルトニ
アンの性質を利用することができるのに対し, ライブラリの提供者は渡さ
れる線形演算子の性質を仮定しない, 汎用性のあるルーチンを提供する必
要があることだと考えていた. 線形演算子について汎用性のある表現形式
は行列であり, 行列の形で渡すと決めれば行列の直積の和という形は仕様
外の形式となる.
4
2
Reverse Communication Interface と ARPACK
つい最近 Reverse Communication Interface という概念に基づくライブ
ラリがあれば前節で説明したトレード-オフは解消できることを知った. こ
2
ベクトルへの作用を関数として実装し, 行列の代わりに関数を引数として受け入れる
ようなライブラリがあればこのトレード-オフは解消する. ただしこの様な実装の場合, 引
数として渡した関数がどの様に使われるかがユーザーには分からず, 並列計算をする場合
などに弊害を起こす可能性がある.
151
の概念は, 各種線形問題の反復解法のアルゴリズムから, 基本演算部分を
分離することを目指すものである.[4] 文献 [4] に分かりやすい例あるので,
これを図 2 として引用しておく. 反復法のループ構造 (Fortran なら DO
文. ループ内の処理はライブラリ内で処理される.) と基本演算部分とをラ
イブラリの外部にユーザが記述する形となっており, 基本演算部分を独自
実装しても反復法部分の処理についてはライブラリを利用できるように
なっている. 図から分かる通り, Reverse Communication とは, 線形演算
子を下位のルーチンに渡して利用するのではなく, 下位ルーチンからベク
トルを受け取って上位ルーチンで線形演算をすることを指している. 文献
[4] の筆頭著者は LAPACK の開発者として著明な研究者で, この Reverse
Communication Interface という概念も, 計算機科学の専門家には良く知
られているものなのだろう. Reverse Communication Interface を用いて
Fortran で書かれた大規模行列向けライブラリに ARPACK というものが
あり,[5] これを使えば最低固有値とその固有ベクトルが求められることも
分かった.
5
計算科学の専門家との交流を考える
実装開始以前に ARPACK を捜し出せなかった原因の一つに, 「必要
なメモリ容量のコンパクト化には, 系がスピン空間の回転に対して不変と
いう物理を利用している. 計算機科学の成果の中に物理を利用したライ
ブラリは見出せないだろう.」という先入観が挙げられる. 実際には, 実
装の工夫により, 見掛け上の上位ルーチンと下位ルーチンとの関係を入れ
替えるという, 想像だにしない形で問題が解決していた. もう一つには,
152
MAIN-PROGRAM
SUBROUTINE CG(N, B, X, WORK, LDW,
SUBROUTINE CGREVCOM(N, B, X, WORK, LDW,
ITER, RESID, MATVEC, PSOLVE,
INFO)
10 CALL CGREVCOM()
ITER, RESID, INFO, NDX1, NDX2, SCLR1,
SCLR2, IJOB)
c
revcom return
c
return from CGREVCOM
c
with job request (IJOB)
IF (IJOB .eq. 1) CALL MATVEC()
IJOB = 3; RETURN
c
c
IF (IJOB .eq. 2) CALL PSOLVE()
c
do matrix-vector multiply with X
ITERATIVE METHOD MAIN LOOP
do preconditoner solve
IJOB = 2; RETURN
do matrix-vector multiply with X
c
do matrix-vector multiply
IF (IJOB .eq. 3) CALL MATVEC()
IJOB = 1; RETURN
c
do stopping test #2
IF (IJOB .eq. 4) CALL STOPTEST2()
c
IF (IJOB .eq. -1) GO TO 20
GO TO 10
20 RETURN
END
Note 1: Diagram only represents reverse communication
interactions.
do stopping test
IJOB = 4; RETURN
c
END OF ITERATIVE METHOD MAIN LOOP
c
all done, return to MAIN-PROGRAM
IJOB = -1; RETURN
STOP
END
Note 2: CGREVCOM() logic to resume execution at
appropriate place not shown.
図 2: Reverse Communication Interface による共役勾配法ルーティン実
装の概念図. CGREVCOM がライブラリとして提供される部分であり, そ
れ以外の CG, MATVEC (線形演算子のベクトルへの作用) などはユーザ
が実装する. CG ループ内の処理をライブラリ内に保持したまま, ループ
構造と MATVEC 部分がユーザに公開されている. ループ構造はサンプ
ル・プログラムとして与えられているので, ユーザーが実際に実装する必
要はない.
153
ARPACK を探し出すキーワードが, Reverse Communication Interface
という計算機科学の専門家ではない我々には馴染のない用語であったこと
が挙げられる. 両者はともに, 我々には計算機科学の知識が不足していた
と言い替え可能だろう. 計算物理と計算機科学とは, ごく近い分野の様で
いて意外と遠い分野なのかもしれない.
この距離の遠さは (我々の立場からいえば) 計算機科学の専門家との交
流を深めることで補うことができる.
今回はその機会を逃してしまったが, 実装上の問題点を抱えているとき
に計算機科学の専門家と相談することができれば, 我々とは違う視点から
の有用な意見を聞くことができるだろう. ただし, その際には問題を物理
とは切り離すことができる部分を, 可能な限り詳細に説明する必要がある.
「大規模疎行列の最低固有値を求めるライブラリが欲しい.」という質問は,
物理とは切り離されているが, 線形演算子の作用を行列の形ではなくサブ
ルーチンとして記述して利用したい, という一番重要な要求が抜け落ちて
おり ARPACK という解答にたどり着かない. 細かい話になりすぎると感
じても, 行列の形で渡すとメモリ容量を取りすぎるが, ベクトルへの作用
をサブルーチンとして記述することでメモリ容量を劇的に減らせることを
説明する必要があるだろう.
交流というからには双方向であるべきである. では, 我々が提供でき
るのは何であろうか？それは「面白い問題」だと考えている. Reverse
Communication Interface なるものが考え出された背景には, ライブラリ
ユーザの「線形演算子のベクトルへの作用を反復解法のルーチンから分離
したい.」という要求があり, また, この要求を「面白い問題」だと考える
専門家が存在しているはずである. ライブラリユーザの要望を知ることは
154
開発者にとって害にはならないだろう.
藤原研は, 出身者の中に山本先生の様な計算機科学の専門家がいらっ
しゃったり, 現在山本先生と同じ講座にいらっしゃる張先生のグループと
以前から交流があったりと, 平均的な計算物理の研究室よりは計算機科学
の専門家との交流が多い. しかし直近の経験から, 計算物理と計算機科学
との距離はもっと近くても良いのではないかと感じた. ここでは特定の分
野間の交流について論じたが, 他の分野間についても同様のことが言える
だろう. 例え問題の解答そのものではなくとも, 視点の違う意見が刺激に
なることはまちがいない. 藤原先生を軸として, これからも OB や関係者
の交流を盛りたてて行ければ素晴らしい. 諸先輩がた, 同輩, 後輩のみな
さんのご賛同を得られれば幸いである.
参考文献
[1] H. Yoshizawa, T. Kakeshita, R. Kajimoto, T. Tanabe, T. Katsufuji, and
Y. Tokura, Phys. Rev. B61, R854 (2000).
[2] S. Yamamoto, T. Fujiwara, unpublished.
山元進, 学位論文, 東京大学, 2006 年.
[3] 例えば次の教科書を見よ. 森正武, 杉原正顕, 室田一雄, 「線形計算」岩波
講座 応用数学 8[方法 2], (岩波書店, 東京, 1994).
[4] J. Dongarra, V. Eijkhout, A. Kalhan, “Reverse Communication Interface for Linear Algebra Templates for Iterative Methods” (1995),
http://www.netlib.org/lapack/lawnspdf/lawn99.pdf (2006 年 11 月現在)
[5] http://www.caam.rice.edu/software/ARPACK/ (2006 年 11 月現在)
155
共有メモリ型並列計算機向け
対称密行列固有値解法の性能評価
名古屋大学 大学院工学研究科 計算理工学専攻 山本有作
1
はじめに
実対称密行列の固有値計算はもっとも基本的な線形計算の一つであり，
電子状態計算，統計計算をはじめとして幅広い応用を持つ [10][12]．近年，
シミュレーションの大規模化・精密化により解くべき行列のサイズが大型
化するにつれ，固有値計算においても並列計算機の活用が必須となってい
る．特に，複数のプロセッサが同一のメモリ空間を持つ共有メモリ型並列
計算機（Symmetric Multi-Processors; SMP）は，プログラミングの容易
さから広く使われており，その上で高い性能を発揮できる固有値ソルバの
開発が重要である．
SMP では，複数のプロセッサが同時にメモリアクセスを行うとアクセ
スの競合が生じ，性能低下の要因となる．これを防ぐには，アルゴリズ
ムにおいてデータの再利用性を高める工夫を行い，各プロセッサの持つ
キャッシュを有効利用することにより，主メモリへのアクセス回数を削減
することが重要である．
156
実対称密行列の固有値計算を行う場合，標準的なアルゴリズムでは，入
力行列 A をハウスホルダー法により実対称三重対角行列 T に変換し，T
の固有値・固有ベクトルを求め，最後にハウスホルダー逆変換により A
の固有ベクトルを求める [8][10][12]．これらの処理のうち，三重対角行列
の固有値・固有ベクトルを求める部分については，MR3 アルゴリズムな
ど，演算量が O(nm)（n は行列のサイズ，m は必要な固有ベクトルの本
数）で済む高速なアルゴリズムが開発されている．また，分割統治法も
広く使われており，演算量は O(n2 )∼O(n3 ) と大きいが，演算の大部分を
データ再利用性の高い行列乗算の形で行えるため，SMP 上での性能は高
い．また逆変換も，演算量は 2n2 m と大きいが，ブロック化を行うことで
演算のほとんどを行列乗算に帰着できる．一方，ハウスホルダー法による
三重対角化の部分は，演算量が約 34 n3 と大きく，かつオリジナルの形で
はデータ再利用性が低いことが知られている [6][7]．
そのため，ハウスホルダー法をキャッシュマシン向けに最適化したアル
ゴリズムとして，Dongarra のアルゴリズム [6][7] と Bischof のアルゴリズ
ム [3][4] が提案されている．このうち，Dongarra のアルゴリズムはハウス
ホルダー変換における行列の rank-2 更新と呼ばれる処理を複数ステップ
分まとめ，データ再利用性の高い行列乗算を用いて計算する方法である．
このアルゴリズムは LAPACK[1]，ScaLAPACK[5] などに実装されて広く
利用されている．しかし，このアルゴリズムでは行列乗算の形で計算でき
るのは全演算量の半分のみであり，残りの半分はデータ再利用性の低い行
列ベクトル積の形となる．そのため，プロセッサ台数の多い SMP では並
列化効率が低下することが知られている [7][11]．
一方，Bischof のアルゴリズムは行列 A をブロックハウスホルダー変換
157
と呼ばれる方法によりいったん帯行列 B に変換し，B をさらに三重対角
行列に変換する 2 段階のアルゴリズムである．この方法では，三重対角行
列への変換において O(n3 ) の演算量を持つ部分すべてを行列乗算で行え
るため，原理的には Dongarra のアルゴリズムより高い性能を実現するこ
とが可能である．しかし，本アルゴリズムでは逆変換も 2 段階となるた
め，その演算量が 4n2 m と倍増する（ただし，これらはすべて行列乗算を
用いて実行可能である）．そのため，Bischof のアルゴリズムが有効かど
うかは，対象とする共有メモリ型並列計算機の特性と，必要な固有ベクト
ルの本数 m に依存する．
そこで本論文では，Dongarra のアルゴリズムと Bischof のアルゴリズ
ムの性能を，2 種類の共有メモリ型並列計算機上で比較する．計算機とし
ては，Opteron SMP（4 プロセッサ）と富士通 PrimePower HPC2500（最
大 32 プロセッサまで使用）を用いる．この 2 つの計算機に対し，使用す
るプロセッサの個数，行列サイズ，求める固有ベクトルの本数を様々に変
化させ，各ケースにおいてどちらのアルゴリズムが高速か調べることを目
的とする．
以下では，まず２章で基本的なハウスホルダー法のアルゴリズムを示
し，その変形として Dongarra のアルゴリズムと Bischof のアルゴリズム
を説明する．次に３章で，上記の各マシン上で２つのアルゴリズムの性能
を比較し，それぞれが高速となる条件を調べる．最後に４章でまとめを述
べる．
158
2
SMP 向けの固有値計算手法
以下では，実対称行列 A の固有値問題
Ax = λx
(1)
を考える．ここで A は n × n 行列，λ はその固有値，x は固有ベクトルで
ある．対称密行列の固有値問題を解く標準的なアルゴリズムは，前述のよ
うに，行列 A の三重対角行列 T への相似変換，T の固有値・固有ベクト
ルの計算，固有ベクトルの逆変換からなる．このうち，特に SMP 上での
性能を出すのが難しいのは三重対角行列への変換であるため，本章ではこ
の部分に絞ってアルゴリズムを説明する．
2.1
基本的なハウスホルダー法
ハウスホルダー法による三重対角化のアルゴリズムをアルゴリズム１に
示す [8][12]．なお，以下で英大文字の太文字は行列，英小文字の太文字は
ベクトル，英大文字，英小文字およびギリシャ文字の細文字はスカラーを
表す．また，ベクトル a に対して，ai は a の第 i 要素を示す．
計算は k = 1 から k = n − 2 までの n − 2 段から成る．図 1 に示すとお
り，第 k 段（1 ≤ k ≤ n − 2）において，全体行列の第 k 列目の第 k + 1 行
目以降の要素からなるベクトルを d(k) ，第 (k + 1, k + 1) 要素を左上隅の
要素とする部分行列を C(k) とする．
第 k 段では，まず d(k) から鏡像変換ベクトル u(k) を作成した後，それ
を行列 C(k) に掛けることによってベクトル p(k) ，q(k) を作成する．u(k)
を第 k 段におけるピボット列，q(k) をピボット行と呼ぶ．最後に，u(k) と
159
q(k) を使って，行列 C(k) の更新を行う．この演算は，行列 C(k) にランク
が１の行列を２個加えるので，rank-2 更新と呼ばれる．
[アルゴリズム 1: ハウスホルダー法]
do k = 1, n − 2
[鏡像変換ベクトル u(k) の作成]
√
σ (k) = d(k)t d(k)
(k)
(k)
(k)
(k)
u(k) = (d1 − sgn(d1 )σ (k) , d2 , . . . , dN −k )
α(k) = 2/ k u(k) k2
[行列ベクトル積]
p(k) = α(k) C(k) u(k)
β (k) = α(k) u(k)t p(k) /2
q(k) = p(k) − β (k) u(k)
[行列の rank-2 更新]
C(k) := C(k) − u(k) q(k)t − q(k) u(k)t
end do
以上の処理は，行列 C(k) に対し，ハウスホルダー変換行列
を用いて相似変換
H(k) ≡ I − α(k) u(k) u(k)t
(2)
³
´−1
C(k) := H(k)
C(k) H(k)
(3)
を行ったことに相当する [8][12]．この処理を第１段から第 n − 2 段まで繰
り返すことにより，行列の三重対角化が完了する [8][12]．
160
図 1: ハウスホルダー法の第 k 段における行列
ハウスホルダー法では，行列ベクトル積と行列の rank-2 更新にそれぞ
れ約 32 n3 の演算量がかかり，それぞれが全演算量のほぼ半分を占める．し
かし，これらの演算では，行列 C(k) のサイズを l(= n − k) とするとき，
ともに O(l2 ) 回の演算に対して O(l2 ) 回のメモリアクセスが必要であり，
データの再利用性が低いという問題点がある．
2.2
Dongarra のアルゴリズム
そこで Dongarra らは，rank-2 更新の部分のデータ再利用性を高めるた
め，ハウスホルダー法を多段化したアルゴリズムを開発した [6][7]．これ
をアルゴリズム２に示す．ここで，(X)i は行列 X の第 i 列目からなる列
ベクトルを表し，[A|B] は行列 A と行列 B とを並べてできる行列を表す．
また，簡単のため，n は M で割り切れると仮定する．
161
[アルゴリズム 2: Dongarra のアルゴリズム]
do K = 1, n/M
U((K−1)∗M ) = φ, Q((K−1)∗M ) = φ
do k = (K − 1) ∗ M + 1, K ∗ M
[部分ハウスホルダー変換]
d(k) := d(k) − U(k−1) (Q(k−1)t )k−(K−1)∗M
− Q(k−1) (U(k−1)t )k−(K−1)∗M
[鏡像ベクトル u(k) の作成]
√
σ (k) = d(k)t d(k)
(k)
(k)
(k)
(k)
u(k) = (d1 − sgn(d1 )σ (k) , d2 , . . . , dN −k )
α(k) = 2/ k u(k) k2
[行列ベクトル積]
p(k) = α(k) (C(k) − U(k−1) Q(k−1)t − Q(k−1) U(k−1)t )u(k)
β (k) = α(k) u(k)t p(k) /2
q(k) = p(k) − β (k) u(k)
U(k) = [U(k−1) |u(k) ]
Q(k) = [Q(k−1) |q(k) ]
end do
[行列の rank-2M 更新]
C(K∗M ) := C((K−1)∗M ) − U(K∗M ) Q(K∗M )t − Q(K∗M ) U(K∗M )t
end do
このアルゴリズムでは，１段ごとに行列の rank-2 更新を行うのではな
く，ある整数 M を決め，M 段に１回，M 段分の rank-2 更新をまとめて
162
行う．これは，行列 C(K∗M ) に対し，M 段分の相似変換
C(K∗M ) :=
³
H(K∗M )
´−1
³
´−1
· · · H((K−1)∗M +1)
× C((K−1)∗M )
× H((K−1)∗M +1) · · · H(K∗M )
(4)
をまとめて行ったことに相当する．これにより，rank-2 更新は，rank-2M
更新で置き換えられる．この処理は一種の行列乗算であるため，データの
再利用性は高く，キャッシュの有効利用が可能となる．しかし，演算量の
残りの半分を占める行列ベクトル積については，データ再利用性が低いま
まで残り，この部分が性能上のネックとなるという問題点がある．
2.3
Bischof のアルゴリズム
この問題を解決するため，Bischof らは行列 A をまず帯行列 B に変換
し，次に B を三重対角行列に変換するという２段階のアルゴリズムを提
案した [3][4]．B の半帯幅を L とするとき，前半の帯行列への変換はサイ
ズ L × L の行列乗算のみを用いて約 34 n3 の演算量で実行でき，後半の三
重対角行列への変換は村田法 [9] あるいは Rutishauser 法 [12] と呼ばれる
方法により約 6n2 L の演算量で実行できる．したがって L ¿ n ならば演
算量はハウスホルダー法とほぼ等しく，かつ演算のほとんどが行列乗算で
行えるため，Bischof のアルゴリズムは SMP 上で高い性能を発揮できる
と考えられる．
Bischof のアルゴリズムのうち，前半の帯行列への変換の部分をアルゴ
リズム３に示す [3][4]．ここで，α などのギリシャ文字の太文字は行列を
163
表すと新たに定義する．また，簡単のため，n は半帯幅 L で割り切れると
仮定する．更に，行列の要素を大きさ L × L のブロックに区分けし，ブ
ロックを単位として参照する場合は大文字の添字を使って「第 K ブロッ
ク列」，
「第 (K, K) ブロック要素」のように呼ぶ．
[アルゴリズム 3: Bischof のアルゴリズム]
do K = 1, n/L − 1
[ブロックハウスホルダ変換の作成]
D(K) を第１ブロックが上三角行列で第２ブロック以下がゼロ行列
であるブロックベクトルに変換するブロックハウスホルダー変換
I − U(K) α(K) U(K)t を求める．
[行列-ブロックベクトル積]
P(K) = C(K) U(K) α(K)
β (K) = α(K)t U(K)t P(K) /2
Q(K) = P(K) − U(K) β (K)
[行列の rank-2L 更新]
C(K) := C(K) − U(K) Q(K)t − Q(K) U(K)t
end do
計算は K = 1 から K = n/L − 1 までの n/L − 1 段から成る．第 K 段に
おける行列の変形の様子とブロックベクトル D(K) ，行列 C(K) の定義を
図 2 に示す．各段における処理の流れは 2.1 節で説明したハウスホルダー
法とほぼ同一であり，ベクトルが幅 L のブロックベクトル（すなわち行
列）に，スカラーがサイズ L × L の行列に置き換わる形となる．特に，ベ
クトル d(k) ，u(k) ，p(k) ，q(k) はそれぞれ幅 L のブロックベクトル D(K) ，
164
図 2: Bischof のアルゴリズムの第 K 段における行列
U(K) ，P(K) ，Q(K) に，スカラー α(k) ，β (k) はそれぞれ L × L 行列 α(K) ，
β (K) に置き換わる．
第 K 段（1 ≤ K ≤ n/L − 1）の処理では，図 2 に示すとおり，まず全体
行列の第 K ブロック列の第 K + 1 番目以降の要素からなるブロックベク
トル D(K) に着目し，D(K) を，第１ブロックが上三角行列でそれ以外の
部分がゼロであるようなブロックベクトルに変換するブロックハウスホル
ダー変換 I − U(K) α(K) U(K)t を求める．このような U(K) ，α(K) は D(K)
の QR 分解および WY-representation と呼ばれる方法により容易に求め
ることができる [2]．次に U(K) を行列 C(K) に掛けることにより，ブロッ
クベクトル P(K) ，Q(K) を作成する．U(K) を第 K 段におけるブロックピ
ボット列，Q(K) をブロックピボット行と呼ぶ．最後に，U(K) と Q(K) を
使って行列 C(K) の更新を行う．この演算は rank-2L 更新となる．以上の
処理は，行列 C(K) に対し，ブロックハウスホルダー変換行列
H̃(K) ≡ I − U(K) α(K) U(K)t
(5)
165
を用いて相似変換
³
´−1
C(K) := H̃(K)
C(K) H̃(K)
(6)
を行ったことに相当する．これにより，図 2 に示すように第 K ブロック
列の帯行列化が完了する [3][4]．
上記の処理においては，行列-ブロックベクトル積と行列の rank-2L 更新
が，それぞれ演算量のほぼ半分を占める．これらは両方ともサイズ L × L
の行列乗算を用いて行うことができ，L が大きいほど１回の計算における
データの再利用性が高くなる．一方，村田法の演算量は L に比例して増
大する．このトレードオフを考慮して L を適切に決めることにより，三
重対角行列への変換を SMP 上で高速に実行できると期待される．
3
2 種の SMP 上での性能評価
3.1
評価条件
前章で述べた Dongarra のアルゴリズムと Bischof のアルゴリズムに対
し，2 種類の共有メモリ型並列計算機上で性能評価を行った．評価に使っ
た計算機は Opteron SMP（1.8GHz，4 プロセッサ）と富士通 PrimePower
HPC2500（最大 32 プロセッサまで使用）である．それぞれについて，環
境の詳細を表 1 に示す．
プログラムは，Dongarra のアルゴリズムについては LAPACK のルーチ
ン DSYTRD（三重対角化）と DORMTR（固有ベクトルの逆変換）を用い，
Bischof のアルゴリズムについては FORTRAN で作成した．また，三重
対角行列の固有値・固有ベクトルを求める部分は，LAPACK の DSTEDC
166
表 1: 評価に使用した計算機環境
項目
Optern SMP
HPC2500
CPU 数
1∼4
1∼32
メモリ
2GB
512GB
コンパイラ
PGI Fortran 77
富士通 最適化 Fortran
LAPACK
ソースよりコンパイル
富士通 並列化 LAPACK
BLAS
GOTO BLAS 1.07
富士通 並列化 BLAS
（分割統治法）を共通に用いた．
テスト行列としては，[0, 1] の一様乱数を要素に持つ対称密行列を用いた．
行列サイズは n=1500，3000，6000，12000, 24000 の５通りとし，Bischof
のアルゴリズムでのブロックサイズ（半帯幅）L はそれぞれに対して 25，
50，100，120，140 とした．また，Dongarra のアルゴリズムでのブロッ
クサイズ M は，LAPACK 内部で決められている値を使用した．
3.2
全固有値・固有ベクトルを求める場合の評価結果
まず，Opteron 上で全固有値・固有ベクトルを求めた場合の実行時間を
表 2（Dongarra のアルゴリズムの場合）および表 3（Bischof のアルゴリ
ズムの場合）に示す．ここではメモリの制限より，行列サイズは 6000 以
下とした．また，プロセッサ数を 1 から 4 まで変えて測定を行った．表 3
において，
「逆変換 1」，
「逆変換 2」は，それぞれ村田法の逆変換，帯行列
化の逆変換を表す．
表中では，各行列サイズ・プロセッサ数に対し，Dongarra，Bischof の
167
表 2: Opteron SMP 上での Dongarra のアルゴリズムの実行時間（秒）
n
プロセッサ数
1500
1
2
3000
4
1
2
6000
4
1
2
4
三重対角化
3.16 2.81 2.77 24.67
21.92 20.72 192.95
170.29 160.31
分割統治法
1.46 0.98 0.74
9.26
5.55
3.70
64.31
35.71
21.46
2.34 1.34 0.89 18.25
10.40
6.69
144.45
82.01
53.53
逆変換
合計
6.96 5.13 4.40 52.18 37.87 31.11 401.71 288.01 235.3
表 3: Opteron SMP 上での Bischof のアルゴリズムの実行時間（秒）
n
3000
1
2
4
1.09 0.83 12.04 6.81
4.34
89.04
47.48
27.39
0.42
0.23 0.16
4.90
2.17
0.97
37.61
20.57
12.90
分割統治法
1.45
0.97 0.70
9.24
5.53
3.70
64.44
35.7
21.43
逆変換 1
5.83
2.54 1.36 38.61 19.74 10.07 261.15 133.67 69.76
逆変換 2
2.34
1.17 0.58 16.71 8.40
1
帯行列化
1.80
村田法
2
4
6000
4
プロセッサ数
合計
168
1500
1
2
4.19
124.1
62.00
31.02
11.84 6.00 3.63 81.50 42.65 23.27 576.34 299.42 162.5
うち高速な方の実行時間を太字で示した．これより，全固有値・固有ベク
トルを求める場合は，プロセッサ数が 1 または 2 ならば Dongarra のアル
ゴリズム，プロセッサ数が 4 ならば Bischof のアルゴリズムが高速である
ことがわかる．
次に，HPC2500 上で全固有値・固有ベクトルを求めた場合の実行時間を
表 4，5（n=6000 の場合），表 6，7（n=12000 の場合），表 8，9（n = 24000
の場合）に示す．
表 4: HPC2500 上での Dongarra のアルゴリズムの実行時間（n = 6000）
n
1
2
4
8
16
32
三重対角化
195.12
143.90
98.91
79.75
57.36
42.77
分割統治法
52.43
30.63
19.24
14.42
11.70
10.43
逆変換
125.18
64.75
32.66
18.07
10.04
6.05
合計
372.73
239.28
150.81
112.24
79.10
59.25
表 5: HPC2500 上での Bischof のアルゴリズムの実行時間（n = 6000）
n
1
2
4
8
16
32
帯行列化
78.91
40.48
29.44
16.39
9.48
7.34
村田法
42.58
14.34
5.11
3.00
2.16
2.04
分割統治法
52.53
30.11
19.60
14.54
11.97
10.30
逆変換 1
204.44
103.99
54.38
29.63
19.18
15.84
逆変換 2
96.72
47.93
24.62
12.48
6.77
3.98
合計
475.18
236.85
133.15
76.04
49.56
39.5
表より，プロセッサ数が 1 ならば Dongarra のアルゴリズム，２以上な
らば Bischof のアルゴリズムが高速であることがわかる．また，Bischof の
アルゴリズムの加速率は非常に良く，n = 24000 の場合に，32 プロセッサ
169
表 6: HPC2500 上での Dongarra のアルゴリズムの実行時間（n = 12000）
n
1
2
4
8
16
32
三重対角化
1477.16
1020.43
973.55
670.76
339.05
330.93
分割統治法
375.16
201.92
128.38
77.35
54.10
44.48
逆変換
962.98
488.06
311.27
139.31
69.80
40.51
合計
2815.30
1710.41
1413.20
887.42
462.95
415.92
表 7: HPC2500 上での Bischof のアルゴリズムの実行時間（n = 12000）
170
n
1
2
4
8
16
32
帯行列化
561.20
296.00
214.41
110.51
53.09
36.30
村田法
222.51
105.59
45.60
14.98
8.29
6.56
分割統治法
375.16
202.80
119.90
77.37
54.17
44.94
逆変換 1
1572.62
794.63
415.98
223.05
125.67
88.10
逆変換 2
753.49
376.29
191.22
96.99
48.28
26.79
合計
3484.98
1775.31
987.11
522.90
289.50
202.69
で 23 倍の加速（並列化効率 73%）を達成している．これは，Bischof のア
ルゴリズムでは O(n3 ) の演算が必要な部分をすべて行列乗算の形で行っ
ているため，キャッシュの利用効率が良く，プロセッサ間のメモリ競合が
起こらないためだと考えられる．
表 8: HPC2500 上での Dongarra のアルゴリズムの実行時間（n = 24000）
n
1
2
4
8
16
32
三重対角化
13191.97
9246.84
5759.82
5141.13
2983.91
1682.02
分割統治法
2868.94
1500.28
796.12
474.22
332.65
217.60
逆変換
7896.60
4087.05
1934.82
1066.77
674.61
298.49
合計
23957.51
14834.17
8490.76
6682.12
3991.17
2198.11
表 9: HPC2500 上での Bischof のアルゴリズムの実行時間（n = 24000）
n
1
2
4
8
16
32
帯行列化
4821.11
2534.39
1480.05
808.83
514.75
231.10
村田法
1290.32
614.36
293.83
152.56
45.86
24.00
分割統治法
2946.52
1493.06
816.39
497.67
342.85
222.58
逆変換 1
12961.06
6270.79
3178.77
1636.07
971.14
535.62
逆変換 2
6264.65
3015.47
1504.78
744.71
419.60
198.92
合計
3.3
28283.66 13928.07 7273.82 3839.84 2294.20 1212.22
一部の固有ベクトルのみが必要な場合
各アルゴリズムで実行時間に占める割合が多い部分を調べると，Don-
garra のアルゴリズムでは，プロセッサ数が多い場合，三重対角化が大部
分の時間を占める．これは，この部分で使われる行列ベクトル積のキャッ
171
シュ利用効率が悪いためである．一方，Bischof のアルゴリズムでは，1 と
2 を合わせた逆変換が全体の 2/3 以上を占める．この部分の実行時間は必
要な固有ベクトルの本数 m に比例するため，m が小さい場合は，Bischof
のアルゴリズムが高速となる場面がより多くなるのではないかと考えら
れる．
そこで，各 n に対し，プロセッサ数を横軸，必要な固有ベクトル本数
の全固有ベクトルに対する割合 m/n を縦軸とする平面上の各点において，
Dongarra，Bischof のどちらのアルゴリズムが高速かを調べた．ここで，
固有ベクトルの逆変換の演算時間は，全固有ベクトルを求める場合の演算
時間の m/n として算出した．
Opteron SMP 上での結果を図 3 に示す．図より，求める固有ベクトル
の本数が少ない場合は Bischof のアルゴリズムが高速であることがわかる．
また，Bischof のアルゴリズムはプロセッサ数が増えるほど優位となり，4
プロセッサの場合は全固有ベクトルを求める場合でも高速となることがわ
かる．
Fraction of eigenvectors
to be computed
100%
100%
100%
80%
80%
80%
60%
60%
60%
40%
40%
40%
20%
20%
20%
0%
0%
1
2
4
n = 1500
Dongarra
Bischof
0%
1
2
4
n = 3000
1
2
4
Number of
processors
n = 6000
図 3: それぞれのアルゴリズムが高速な領域（Opteron SMP）
172
HPC2500 上での結果を図 4 に示す．この場合も Opteron SMP と同様
の傾向が見られ，プロセッサ数が 2 より多い場合は，常に Bischof のアル
ゴリズムが高速である．
Fraction of eigenvectors
to be computed
100%
100%
100%
80%
80%
80%
60%
60%
60%
40%
40%
40%
20%
20%
20%
0%
0%
1
2
4
8 16 32
n = 6000
Dongarra
Bischof
0%
1
2
4
8 16 32
n = 12000
1
2
4
Number of
8 16 32 processors
n = 24000
図 4: それぞれのアルゴリズムが高速な領域（HPC2500）
4
おわりに
本論文では，対称密行列の固有値計算において，特に三重対角化と逆
変換の部分に焦点を当て，共有メモリ型並列計算機に適した手法として
Dongarra のアルゴリズム Bischof のアルゴリズムを紹介した．また，4 プ
ロセッサの Opteron SMP と 32 プロセッサの HPC2500 を用いて両アルゴ
リズムの性能を比較した．その結果，求める固有ベクトルの本数が少ない
ほど，また，プロセッサ数が多いほど，Bischof のアルゴリズムが優位と
なることを明らかにした．
今後の課題としては，クワッドコアプロセッサや Cell プロセッサなど，
より多様なマシン上で両アルゴリズムの性能を比較すること，一般固有値
問題に適用した場合の性能を評価すること，本研究で開発したプログラム
173
を電子状態計算や第一原理分子動力学に適用することなどが挙げられる．
謝辞
修士課程を通じて暖かくご指導下さり，本研究に関しても多くの有益な
コメントを下さった東京大学大学院工学系研究科の藤原毅夫先生に感謝い
たします。
参考文献
[1] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J., Du
Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Ostrouchov, S. and Sorensen, D.: LAPACK Users’ Guide, Second edition,
SIAM, Philadelphia, 1995.
[2] Bischof, C. and van Loan, C. F.: The WY Representation for Products of Householder Matrices, SIAM Journal on Scientific and Statistical Computing, Vol. 8, No. 1, pp. s2–s13 (1987).
[3] Bischof, C., Marques, M. and Sun, X.: Parallel Bandreduction
and Tridiagonalization, Technical Report 8, PRISM Working Note,
1993.
[4] Bischof, C., Lang, B. and Sun, X.: Parallel Tridiagonalization
through Two-step Band Reduction, Technical Report 17, PRISM
Working Note, 1994.
174
[5] Blackford, L. S, Choi, J., Cleary, A., D’Azevedo, E., Demmel, J.,
Dhillon, I., Dongarra, J., Hammarling, S., Henry, G., Petitet, A.,
Stanley, K., Walker, D. and Whaley, R. C.: ScaLAPACK Users’
Guide, SIAM, Philadelphia, 1997.
[6] Dongarra, J. J., Hammarling, S. J. and Sorensen, D. C.: Block Reduction of Matrices to Condensed Forms for Eigenvalue Computations, Journal of Computational and Applied Mathematics, Vol. 27,
pp. 215–227 (1989).
[7] Dongarra, J. J. and van de Geijn, R. A.: Reduction to Condensed
Form for the Eigenvalue Problem on Distributed Architectures,
Parallel Computing, Vol. 18, No. 9, pp. 973–982 (1992).
[8] Golub, G. H. and van Loan, C. F.: Matrix Computations, Third
edition, Johns Hopkins University Press, 1996.
[9] 村田健郎, 小国力, 唐木幸比古: スーパーコンピュータ： 科学技術計
算への適用, 丸善, 1985.
[10] Parlett, B. N.:
The Symmetric Eigenvalue Problem, SIAM,
Philadelphia, 1998.
[11] Salvini, S. A. and Mulholland, L. S.: The NAG FORTRAN Library,
in Proceedings of the Ninth SIAM Conference on Parallel Processing
and Scientific Computing 1999, SIAM, Philadelphia (1999).
175
[12] Wilkinson, J. H. and Reinsch, C. (eds.): Linear Algebra, SpringerVerlag, 1971.
176