実空間計算手法に基づく第一原理分子動力学および電気伝導特性計算プログラムの開発
広瀬喜久治，後藤英和，小野倫也，稲垣耕司，江上善幸
大阪大学大学院工学研究科

1

はじめに

21 世紀の IT 産業を支える半導体デバイスや光通信デバイス製造産業では，デバイスの微細化や高機能化にともない，原子・分子サイズの細線を利用したデバイスの開発が必要とされている．通常の金属や半導体などのマクロスケールにおける固体中の電気伝導では，電子は，固体中の不純物・欠陥との衝突や格子振動などによるエネルギーの散逸を伴う非弾性散乱を繰り返しながら伝導するために，電場の大きさに比例する速度に落ち着く．その結果，オームの法則が成立し電気伝導 (electricconductance) は導線の長さに反比例し，断面積に比例する．一方，平均自由行程より短い導線を電極で結んだようなナノスケールの構造体では，電子の通路である配線が電子の平均自由行程よりもはるかに小さな構造をしているために，電子の波動性が顕著になり，デバイスの電気伝導特性は量子化される．このような電気伝導は，ナノデバイスの新しい動作原理として利用できるものと期待されている．
実験では，ナノスケール構造体の電子輸送特性は，走査型トンネル顕微鏡 (Scanning tunneling microscopy; STM) やブレークジャンクション，リソグラフィーといった技術を用いて電極間にナノスケールの間隙を作成し，その間に挟まれた原子鎖や分子を流れる電流を測定することによって得られる．金属原子鎖系では，量子接点を引き伸ばしていくと電極間に数原子からなる原子鎖が形成される．この原子鎖のコンダ クタンスは電極間距離に対して高さが G0 (=2e2 /h: e は電気素量，h はプラン ク定数) の階段状の振る舞いをすることが金の量子接点について最初に報告された．その後，多くの金属でコンダ クタンスの量子化が確認された．また，分子系に関しても，ベンゼンチオールや DNA 分子，自己組織化単分子膜 1 層のコンダクタンスなども報告されている．分子系については，結果にはまだ検討の余地があると思われるが，新たな機能を持つ材料の可能性として興味深い．
一方，計算物理の分野でも，これまでは人工的なモデルを仮定した計算が中心であったが，近年の計算機および計算手法の発展により，このようなナノスケールでの電子輸送現象を第一原理計算で予測することが可能になってきた．特に，これまで行われてきた実験の理論的裏づけや未知のナノスケール構造体の持つ電子輸送特性の予測が精力的に行われている．本稿では，われわれが開発した Overbridging
Boundary-Maching 法 [1, 2] を用いた数値計算における工夫とそれを用いたアプ リ
ケーションのひとつを紹介したい．

12電気伝導計算プログラムの高精度化

半無限に続く電極に挟まれたナノ構造体の電子状態を調べるには，図 1 に示すような計算領域を左側電極領域，ナノ構造体領域，右側電極領域の 3 つに分割した計算モデ
ルを用いる．無限の系全体に拡がる散乱解 {Ψ(zk )} , k = −∞, . . . , −1, 0, 1, . . . , ∞,࠽ࡁ᭴ㅧ૕Ꮐ஥㔚ᭂฝ஥㔚ᭂinǞ㩿zk㪀refǞ㩿zk㪀ฝะ߈ㅴⴕᵄฝะ߈ㅴⴕᵄᏀะ߈ㅴⴕᵄtraǞ㩿zk㪀ࠛࡃࡀ࠶࠮ࡦ࠻ᵄࠛࡃࡀ࠶࠮ࡦ࠻ᵄ

Z-1 Z0Z1Zm-1 Zm Zm+1 Zm+2Z2z

図 1:
半 無 限 に 続 く 2 つ の バ ル ク 電 極 に 挟 まれ た ナ ノ 構 造 体 モ デ ル ．
Φin (zk ), Φref (zk ), Φtra (zk ) はそれぞれバルク電極内の入射波，反射波，透過波を表わす．
を得るには，散乱の境界条件を満たす Kohn-Sham 方程式の解を求めればよい．電子が左側電極から j 番目のチャンネルを通って入射する場合の境界条件は
⎧N

⎪
⎪
in⎪
Φ
(z)+rij Φref⎪ki (zk ) · · · k ≤ 0⎨ ji=1
Ψj (zk ) =(1)N⎪⎪tra⎪⎪tij Φi (zk )··· k ≥ m+ 1⎩i=1reftraである．ここに，Φinj (zk ), Φi (zk ), Φi (zk ) は，図 1 に示すように，それぞれバルク電極内の入射波，反射波，透過波を表わす．また，rij は反射確率振幅，tij は透過確率振幅である．入射波としては進行波を考えれば十分であるが，正しい散乱解を得るには，反射波と透過波として進行波だけでなく指数関数的に激しく増減するエバネッセント波も考慮する必要がある．散乱解を求めるための数学的なフォーマ
リズムは省略し ，ここではその過程での数値計算の問題点を紹介する．
散乱解を得るには，バルク電極内に存在する進行波とエバネッセント波 (バルク解と呼ぶ) を計算しておく必要がある．図 2 にバルク解を計算するためのモデルをMの領域) における Kohn-Sham示す．第 M 番目のユニットセル (z が z1M ≤ z ≤ zm
方程式
M−1MMMBΦ(zm) + AM
1 Φ(z1 ) + BΦ(z2 ) = EΦ(z1 )2╙(M-1)⇟⋡࡙࠾࠶࠻࠮࡞╙M⇟⋡࡙࠾࠶࠻࠮࡞╙(M+1)⇟⋡࡙࠾࠶࠻࠮࡞

xM-1

M-1

Zm-1 Zm

MZ1MMZ2M

M+1 M+1Z2

Zm-1 Zm Z1Lzzyz

図 2: 周期的ポテンシャルをもつ結晶のモデル． zkM は第 M 番目ユニットセル中の第 k 番目グリッド の z 座標を表わす．
M
M
M
M
BΦ(zk−1) + AMk Φ(zk ) + BΦ(zk+1 ) = EΦ(zk )(k = 2, 3, . . . , m − 1)
M+1M
M
M
BΦ(zm−1) + AM) = EΦ(zm)m Φ(zm ) + BΦ(z1(2)を考えよう．このユニットセルを全体の系から一部を truncate した有限孤立系であるとみなしたときのハミルトニアン行列を用いると，(2) 式は⎡⎤⎤ ⎡M−1BΦ(zm)Φ(z1M )M⎢⎥⎥ ⎢0⎥ ⎢ Φ(z2 ) ⎥ ⎢⎢⎢⎥⎥.
.
..
..
(3)E − Ĥ M ⎢⎥⎥=⎢⎢⎥⎥ ⎢⎣ Φ(z M ) ⎦ ⎣⎦0m−1BΦ(z1M+1 )MΦ(zm)−1と表せる．そこで，行列 Ĥ M のグリーン関数行列 ĜM (E) = E − Ĥ Mを定義
すると，(3) 式は

⎡⎡⎤
M−1)
BΦ(zm⎢⎢⎥⎥0⎢⎢⎥⎥⎢⎢⎥⎥.
..
⎢⎥⎥ = ĜM (E) ⎢⎢⎢⎥⎥M⎣ Φ(zm−1⎣⎦0) ⎦MΦ(zm)BΦ(z1M+1 )Φ(z1M )
Φ(z2M )..
.

⎤(4)となる．(4) 式の最上行と最下行の式に注目すると，今度は第 M 番目のユニットセルに対する Overbridging Boundary-Maching 公式  M M M

M−1MBΦ(zm))
G (z1 , z1 ) GM (z1M , zm
Φ(z1M )
(5)=M
MMM
Φ(zmGM (zm), z1M ) GM (zm, zm)
BΦ(z1M+1 )3

を得る． GM (zkM , zlM ) は行列 ĜM (E) の対応するブロック行列要素である．
周期系の Kohn-Sham 方程式の解 Φ(z) に課されるもうひとつの重要な制限はブロッホ条件
Φ(z + Lz ) = λΦ(z)(6)である．ここに，λ はブロッホ因子 eikz Lz ，Lz は z 方向のユニットセルの長さ，kzはブロッホ波数の z 成分である．通常のバンド 計算の場合には kz は実数であるが，ここでは一般に kz が複素数であってもよいことにする．次に示すように，こうする
ことにより，進行波だけでなくエバネッセント波も含めた一般化ブロッホ状態を計算することが可能になる．
結局， ブロッホ条件
MM−1
Φ(zm) = λΦ(zm)M+1M) = λΦ(z1 )Φ(z1(7)Mを用いて (6) 式から Φ(z1M ) と Φ(zm) を消去すると，次の一般化固有値方程式が導かれる．


M−1M−1MMΦ(zmΦ(zm))
I0Wm,1Wm,m=λ. (8)MMW1,1W1,m0I
Φ(z1M+1 )Φ(z1M+1 )ここに，MW1,1= GM (z1M , z1M )BM
M= GM (z1M , zm)B
W1,mM
M
= GM (zm, z1M )B
Wm,1(9)MMM= GM (zm, zm)B.
Wm,m(8) 式を解くと，ブロッホ因子 λ = eikz Lz は固有値として，一般化ブロッホ状
M−1 M+1, z1での値は固有ベクトルとして求まる．(6) 式の λ の定義態 Φ(z) の z = zmからわかるように，|λ| = 1 の場合のブロッホ波数 kz は実数であり，固有ベクトルは通常の進行波に対応する．一方，|λ| > 1 (|λ| < 1) の λ に属する固有ベクトルは
右 (左) に行くほど 増加するエバネッセント波にあたる．これら右増加のエバネッセント波の数と左増加のそれの数は等しく，同じ性質を持った右増加と左増加のエバネッセント波の固有値 λ を掛け合わせると 1 になるという性質を持つ．
実はこの一般化固有値問題を解くことは非常に厄介であり，実空間のグ リッド
間隔を細かくして高精度計算を行おうとすると倍精度の数値計算では図 3 に示すように固有値を掛け合わせると 1 になるという性質が保証されないことが多い．そこで，この部分だけ桁落ちによる精度劣化が起こらないようにして計算を行う．

42Nlog10Π
|λn|n=11.0

double precision
quadruple precision

0.8
0.6
0.4
0.2
0.0
0.60

0.65

0.70

0.75

0.80

Grid spacing hµ (a.u.)
図 3: (8) 式の各固有値の積．ナトリウムバルクの場合における倍精度計算と 4 倍精度計算の結果を比較したものである．

3アルミニウム原子ワイヤーの電気伝導特性

電極間に挟まれた原子鎖について，これまで数多くの実験的研究が行われており，アルミニウム単原子ナノワイヤーの特徴的な点は，破断直前での電極間距離に対するコンダクタンスの階段状の変化 (コンダクタンストレースと呼ぶ) が平らなプラトーを描くのではなく，下に凸のカーブを描くことである．この結果は，破断直前ではナノワイヤーが長くなっていくと電気抵抗が下がるという非直感的な現象が生じていることを意味している．アルミニウム単原子ナノワイヤーの場合，電極とナノワイヤーの接続状態によりコンダ クタンスが変化することが知られているため，コンダ クタンスのわずかな変化を定量的に調べるには，結晶でできた電極を用いるか必要がある．
ここでは，Al(001) 結晶電極に挟まれたアルミニウム原子 3 個からなる原子ナノワイヤーの電極間距離に対するコンダクタンスの変化を調べた．モデル系を図 3 に示す．ナノワイヤーおよび電極部の原子核は，Troullier-Martins 型のノルム保存型
擬ポテンシャル [3, 4] を用いて扱い，電子の交換相関相互作用には局所密度近似 [9]を用いた．また，実空間差分法 [1, 5, 6, 7, 8] におけるグリッド の幅はすべての方向において 0.60 bohr とし ，二階微分には 9 点差分公式を用いた．両方の電極表面上に原子 4 個からなる一辺が a0 (a0 は格子定数) の正方形の台座を置き，その間に単原子ナノワイヤーを挟んだ．フェルミエネルギーをもった入射電子のゼロバイアス極限での透過確率を図 5 に示す．いずれの電極間距離においても，電気伝導に支配的なのは主に s-pz 軌道からなる第 1 チャネルであり，px -py 軌道が支配的な二重に縮退した第 2，第 3 チャネルによる寄与は，無視し うるほど 小さい．また，第 1 チャ5xyᏀ஥㔚ᭂzฝ஥㔚ᭂᢔੂ㗔ၞ

図 4: 計算モデル系．
ネルの透過確率は，電極間距離が 26.5 bohr の時に極小値をもつ．

ㅘㆊ₸1.0

1st ch2nd ch
3rd ch

0.50.024.0

25.0

26.0

27.0

28.029.0㔚ᭂ㑆〒㔌 L (bohr)

図 5: 各チャネルの透過確率．

4

おわりに

ナノ構造体における電子のバリスティック伝導についての最近の第一原理計算手法のひとつである Overbridging Boundary-Maching 法用いた数値計算における工夫とそれを用いたアプリケーションのひとつを紹介した．これらは，第一原理伝導計
算の一例に過ぎないため，最新の結果は参考文献 10-17 を参考にされたい．
現在の問題点を挙げると，第 3 章に示したような原子鎖の計算でも，ナノ構造体領域を記述するハミルトニアンに対する巨大なグリーン関数の計算や (8) 式の固有値，固有ベクトルの計算に膨大な計算が必要であり，SX-7 でも全ての計算に 1 日以上かかっていた．アルゴ リズムの改善や SX-7c の導入により計算時間は大幅に短縮されつつあるが，高精度な電気伝導計算を可能にしているのはスーパーコンピュータである．

6ナノ構造体は，ナノデバイスを構成する素子として注目を集めている．理論計算の分野では，計算モデルに原子構造を考慮した第一原理電気伝導計算がようやく可能になり，これまで実験のみではなかなか見えてこないナノ構造体中の電流経路や電子輸送を支配している要因が明らかになりつつある．今日，シリコンを用いたデバイスもナノスケールまで微細化されつつある．このようなデバイス用の素子の電気伝導特性も計算機の性能向上に伴い第一原理計算を通じて明らかになるものと期待される．

謝辞
本稿で紹介したシミュレーションの実行に当たっては，東北大学情報シナジーセン
ターのスーパーコンピュータを大いに利用させていただいた．また，本研究で用いた実空間差分法に基づく第一原理計算プログラムの開発は，情報シナジーセンターとの共同研究として進められ，アルゴ リズムのベクトル化・並列化に当たっては多くの助言をいただいた．また，大阪大学 21 世紀 COE プログラム「原子論的生産技術の創出拠点」，文部科学省科学研究費補助金・特定領域研究「 次世代量子シミュレータ・量子デザイン手法の開発」(課題番号:17064012) の支援を受けて行われた．
ここに記して感謝する．

参考文献[1] K. Hirose, T. Ono, Y. Fujimoto, and S. Tsukamoto, First-Principles Calculations in Real-Space Formalism, Electronic Configurations and Transport Properties of Nanostructures (Imperial College Press, London, 2005).
[2] Y. Fujimoto and K. Hirose, Phys. Rev. B 67, 195315 (2003).
[3] We used the norm-conserving pseudopotentials NCPS97 constructed by K.
Kobayashi. See K. Kobayashi, Comput. Mater. Sci. 14, 72 (1999).
[4] N. Troullier and J.L. Martins, Phys. Rev. B 43, 1993 (1991).
[5] J.R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. Lett. 72, 1240 (1994).
[6] T. Ono and K. Hirose, Phys. Rev. Lett. 82, 5016 (1999).
[7] T. Ono and K. Hirose, Phys. Rev. B 72, 085105 (2005).
[8] T. Ono and K. Hirose, Phys. Rev. B 72, 085115 (2005).
[9] J.P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
[10] Y. Egami, T. Sasaki, T. Ono, and K. Hirose, Nanotechnology, 16, S161
(2005).
[11] T. Ono and K. Hirose, Phys. Rev. Lett. 94, 206806 (2005).
[12] 佐々木孝，江上喜幸，谷出敦，小野倫也，後藤英和，広瀬喜久治，日本金属学会誌，69, 457 (2005).

7[13] Y. Egami, T. Ono, K. Hirose, Phys. Rev. B, 72, 125318 (2005).
[14] Y. Egami, T. Sasaki, T. Ono, H. Goto, and K. Hirose, Jpn. J. Appl. Phys 45,2132 (2006).
[15] S. Horie, T. Ono, Y. Kuwahara, K. Endo and K. Hirose, Jpn. J. Appl. Phys45, 2154 (2006).
[16] D. Nakagawa, K. Kutsuki, T. Ono, and K. Hirose, Physica B 376-377, 389(2006).
[17] T. Ono, S. Horie, K. Endo, and K. Hirose, Phys. Rev. B 73, 245314 (2006).

8