半無限電極に挟まれたナノジャンクションの第一原理 - 日本金属学会

日本金属学会誌第 69 巻第 6 号(2005)457459
半無限電極に挟まれたナノジャンクションの
第一原理電子状態計算法1
2
佐々木　孝1,
江上喜幸1,2
谷出敦1,2
小野倫也2
後藤英和1
広瀬喜久治1
1大阪大学大学院工学研究科精密科学専攻
2大阪大学大学院工学研究科附属超精密科学研究センター
J. Japan Inst. Metals, Vol. 69, No. 6 (2005), pp. 457
459
 2005 The Japan Institute of Metals
First
Principles Calculation for Electronic Structures of Nanojunctions Suspended
between Semi
Infinite Electrodes
2, Yoshiyuki Egami1,
2, Atsushi Tanide1,
2,
Takashi Sasaki1,
2
1
1
Tomoya Ono , Hidekazu Goto and Kikuji Hirose
1Department
2Research
of Precision Science and Technology, Osaka University, Osaka 5650871
Center for UltraPrecision Science and Technology, Osaka University, Osaka 5650871
We develop a method for highspeed and highaccuracy firstprinciples calculations to derive the groundstate electronic
structure by directly minimizing the energy functional. Making efficient use of the advantages of the realspace finitedifference
method, we apply arbitrary boundary conditions and employ spatially localized orbitals. These advantages enable us to calculate
the groundstate electronic structure of a nanostructure sandwiched between crystalline electrodes. The framework of this
method and numerical examples for metallic nanowires are presented.
(Received February 15, 2005; Accepted March 28, 2005)
Keywords: firstprinciples calculation, realspace finitedifference method, orderN algorithm
算に対して周期境界条件と孤立境界条件を組み合わせて使用
はじ
1.
めに
することが可能になる．さらに，計算領域全体に広がった波
動関数の代わりに空間の限られた範囲内に局在した関数
原子，電子レベルの微視的世界では巨視的世界と異なる物
(Localized Orbital)を用いて電子状態を計算する手法を用い
理法則に支配され，電子はシュレディンガー方程式を満たす
ることが可能になり，通常は計算量がモデルサイズの 3 乗
状態でのみ存在するという特性がある．そのため，電子レベ
に比例するところをモデルサイズに比例させるアルゴリズム
ルでナノ材料の原子配列構造に基づく物性を予測するために
(オーダー N アルゴリズム)を組み込むことができる．さら
は，量子力学の第一原理に基づいて電子の従う基礎方程式を
に，この Localized Orbital の特性を利用することによって
解かなければならない．これを実現するために，密度汎関数
半無限に続く結晶電極に挟まれたナノストラクチャーの電子
理論1)の
状態計算も可能になる．
KohnSham 方程式を基底関数を用いて解く計算手
法が多数開発されてきた24) ．これらの手法は多くの成功を
本論文では，実空間差分法517,21) と Localized Orbital
収めてきたが，いくつかの欠点を持っている．例えば，原子
法22,23)を組み合わせて半無限に続く結晶に挟まれたナノスト
波動関数展開法では用いた基底関数で十分な計算精度が保証
ラクチャーの電子状態を計算するための第一原理計算手法を
されるか常に注意を払う必要があり，また，平面波展開法で
紹介する．この手法では，Mauri らにより提案されたエネル
は孤立境界条件を用いると真空領域を広く取る必要が生じる
ギー汎関数24) を Direct Minimization 法25) を用いて最小化す
ため計算量が膨大になってしまう．これらの問題点を回避す
ることにより，従来のセルフコンシステント法の反復計算を
るために，差分法や有限要素法を用いて実空間で直接
行うことなく KohnSham 方程式の解と同じ電子系のエネル
方程式を解く手法が提案されている521)．この
ギーを最小にする電子状態を得ることができる．本手法によ
手法を用いれば，境界条件を周期的に限る必要がないため，
って半無限に続く結晶電極に挟まれたナノストラクチャーの
クラスターに対して孤立境界条件を用いることや，表面の計
電子状態計算が可能であることを確かめるために，金の単原
Kohn Sham
1
Mater. Trans. 45(2004) 14191421 に掲載
2 大阪大学大学院生(Graduate Student, Osaka University)
子ワイヤーの電気伝導コンダクタンスを計算し，実験および
他の計算の結果と一致することを確認した．
458
第
日本金属学会誌(2005)
69
巻
x, y 方向には周期境界条件を課し， z 方向には半無限に結晶
2.
計算
手法
が続く境界条件を設定している．計算は次の 2 つのステッ
 3 方向に周期境界条件を課したスーパーセ
プからなる．
バルクやクラスターからなる系に対する Mauri らのエネ
ルギー汎関数は次のように与えられる．
 バルクに挟まれた散乱領域の Localized Orを計算する．
(∑Q q |－ 12 : | q ＋F [ r])
＋h N－ r(r ) dr
( f )
bital を式( 1 )のエネルギー汎関数を最小化させることによ
M
E [Q, {q}]＝2
i, j
2
i
ル法を用いてバルク領域のポテンシャルと Localized Orbital
j
り計算する．このとき，バルク領域の Localized Orbital は
i, j
(1)
 で計算されたものに固定しておき，散乱領域の Localized

Orbital のみを最適化することにより基底状態の電子状態を
ここで{q}は M 個の 1 次独立な任意の関数， F [ r]はハート
計算する．散乱領域のハートリーポテンシャルとイオンポテ
リーポテンシャル，イオンポテンシャル，交換相関ポテンシ
ンシャルは，次のポアソン方程式
ャル，外場によるポテンシャルエネルギーの和，h はケミカ
:2v(r )＝－4p( rnuc(r )＋r(r )),
(4)
ルポテンシャルであり，系内の電子数 N を正しく与えるよ
を解くことにより計算する．ここで v ( r )はハートリーポテ
うに調節する． Q は重なり積分 S を展開した 1 次の項 Q ＝
ンシャルとイオンポテンシャルの和であり，rnuc は原子核の
2I－S であり，I は単位行列である．また，電子密度は Q と
電子密度である．z 方向に半無限に結晶が続く境界条件を設
q を用いて次のように定義される．
定すると，式( 4 )は差分近似を用いて次のような連立方程
M
r(r )＝r[Q, {q}](r )＝2∑Qijqj(r )qi(r )
(2)
ij
ここでは，簡単のため{q}は実関数とし，スピンの自由度は
無視している．
最急降下法や共役勾配法を使って式( 1 )のエネルギー汎
関数を関数{q}に関して最小化することにより，基底状態に
おける電子状態を得ることができる．最小化を行うために必
要なエネルギー汎関数の{q}に関する微分は，次のように与
えられる．
M
dE [{q}, h]
＝2∑[(HKS[{q}]－h )|qj 〉Qji－|qj〉
dqi
i
×〈qj|HKS[{q}]－h |qi〉］
式に置き換えられる．
…
Axy－2Bz
 vxy(z1 ) 
Bz
0
0
Bz
Axy－2Bz S
P
vxy(z2 )
0
S
S
0
P
S
P
Bz
vxy(zN －1)
S Axy－2Bz
…

0
0
Bz
Axy－2Bz vxy(zN ) 
Pxy(z1 )－Bzvxy(z0 )

Pxy(z2 )
＝－4p P
(5)
Pxy(zN －1)
Pxy(zN )－Bzvxy(zN ＋1 )
ここで Axy は x, y 方向の 2 次微分を表す Nx×Ny 次元の行列，
Pxy(zk ) , vxy(zk )はそれぞれ r (r )＋ rnuc (r ), v (r )の z ＝ zk にお
ける xy 平面上での Nx × Ny 個の値から構成される列ベクト
ル，Bz は，

 
xy
xy


xy
xy
(3)
ここで HKS は Kohn Sham ハミルトニアン26) である．最小
化の過程においてローカルミニマムに落ち込むことを避ける
ために M は占有状態数よりも大きい数にしておく． q に波
動関数のような空間全体に広がる関数を用いるとモデルサイ
xy
B z＝
ズの 3 乗に比例する計算量が必要になるが， Localized Or-
1
I ,
h2z xy
(6)
bital を用いると重なりのない関数同士の直交化を行う必要
である．ここで， hz は z 方向のグリッド幅であり， Nx, Ny,
がないため，計算量をモデルサイズに比例させることも可能
Nz はそれぞれ x, y, z 方向の分点数である．式( 5 )の境界条
である．
件 v(z0 )と v(Nz＋1 )にはステップ( 1 )で計算したバルク領域
次に，この計算手法を Fig. 1 に示すような半無限に続く
のポテンシャルの値を用いる．
電極に挟まれたナノストラクチャーの電子状態計算に適用す
る方法について説明する．計算対象である散乱領域は，ナノ
3.
結果と考察
スケールジャンクション，走査型トンネル顕微鏡のようなチ
ップサンプルシステムにおけるトンネルジャンクション，異
この手法の有用性を確かめるために金の(001)結晶電極に
なるバルクの界面，侵入型格子欠陥や転位などで構成される．
挟まれた金の単原子ワイヤーの電気伝導コンダクタンスを計
算し，実験値および他の計算の結果と比較した．まず，ワイ
ヤーの基底状態における電子状態を 2 章で述べた計算手法
を用いて計算し，得られた電子状態を用いてコンダクタンス
を計算した．Fig. 1 に計算モデルを示す．散乱領域は，単原
子鎖，ピラミッド型の台座， 6 層の金の(001)面からなる電
極により構成されている．ワイヤーを構成する原子数 L は 1
から 3 まで変化させた．スーパーセルの xy 方向の長さは
Fig. 1 Schematic representation of the wire system for L＝1
with the transition region suspended between the left and right
bulk regions. The large and small spheres represent the atoms
on and below the cross section, respectively.
2a0 であり，a0(＝0.408 nm)は金の結晶の格子定数である．
グリッド幅は 0.036 nm であり，これは従来の平面波法では
22Ry のカットオフエネルギーに対応する．運動エネルギー
第
6
号
半無限電極に挟まれたナノジャンクションの第一原理電子状態計算法
459
晶電極に挟まれた金の単原子ワイヤーの電気伝導コンダクタ
ンスの計算結果を示した．コンダクタンスが実験および他の
理論計算の結果と矛盾しないことから，本手法によって半無
限に続く結晶電極に挟まれたナノストラクチャーの電子状態
計算が可能であることわかった．さらに，Localized Orbital
を用いれば計算量をモデルサイズに比例させるオーダー N
計算が可能になる．現在，この方法を用いた大規模計算アル
ゴリズムの開発を行っている．
本稿で紹介した研究結果の多くは文部科学省の科学研究費
補助金若手研究(B), 14750022, 2002 の支援を受けて行われ
た．計算機シミュレーションには，東京大学物性研究所，岡
Fig. 2 Conductance of the atomic wires in units of the conductance quantum (G0 ) as a function of the number of gold atoms
in the wire.
崎国立研究所，東北大学情報シナジーセンターのスーパーコ
ンピュータを用いた．
文
献
演算子から生じる微分には中心差分7)を用い，原子核からの
ポテンシャルには局所型擬ポテンシャル27) を用いた．電子
間の交換相関相互作用には密度汎関数法1)に基づく局所密度
近似26)を用いた． Localized Orbital の局在範囲は半径 0.635
nm にし， Localized Orbital の数は M ＝ 178 ＋ L とした．系
全体に広がる波動関数は overbridging boundary matching
method28) いて求め，ゼロバイアス極限でのコンダクタンス
G は次の LandauerB äuttiker の公式で計算した．
G＝
2 e2
Tij(EF )
h ∑
i, j
(7)
ここで Tij は入射電子の透過確率である．ゼロバイアス極限
での電子の伝導はフェルミレベル Ef における電子の透過に
のみ支配されるため，ここでは Tij(Ef )のみを用いて計算し
た．これまでの実験により，金の単原子鎖は約 1 G0(＝2e 2/
h )のコンダクタンスを示し，原子鎖の長さが変化するにつ
れて 2 原子周期で振動することが確認されている．さら
に，いくつかの理論計算においてもコンダクタンスが約
1 G0 に量子化されることが報告されている． Fig. 2 に原子
鎖を構成する原子数 L を変化させたときのコンダクタンス
G の変化を示す．本計算により得られたコンダクタンスは
約 1 G0 であり，偶数個の原子からなるワイヤーよりも奇数
個の原子からなるワイヤーのほうが高いコンダクタンスを示
している．コンダクタンスの計算結果が実験および他の理論
計算の結果と矛盾しないことから，本手法によって半無限に
続く結晶電極に挟まれたナノストラクチャーの電子状態計算
が可能であることがわかる．
4.
ま
と
め
実空間差分法521)と Localized Orbital 法22,23)を組み合わせ
て，半無限に続く結晶電極に挟まれたナノストラクチャーの
電子状態を計算する第一原理計算手法を開発した．本手法で
は Mauri らのエネルギー汎関数を Direct Minimization 法を
用いて最小化することにより，エネルギー汎関数の数学的に
厳密に定義された最小値として基底状態における電子状態を
求めている．本計算手法を使った例として，半無限に続く結
1) P. Hohenberg and W. Kohn: Phys. Rev. 136(1964) B864B871.
2) J. R. Chelikowsky and M. L. Cohen: Handbook on Semiconductors, ed. by P. T. Landsberg, (Elsevier, Amsterdam, 1992), Vol.
1, p. 59.
3) J. Ihm, A. Zunger and M. L. Cohen: J. Phys. C 12(1979) 4409
4422.
4) M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D.
Joannopoulos: Rev. Mod. Phys. 64(1992) 10451097.
5) Kikuji Hirose, Tomoya Ono, Yoshitaka Fujimoto and Shigeru
Tsukamoto: FirstPrinciples Calculations in RealSpace Formalism―Electronic Configurations and Transport Properties of
Nanostructures― (Imperial College Press, London, 2005).
6) T. Ono and K. Hirose: Phys. Rev. Lett. 82(1999) 50165019.
7) J. R. Chelikowsky, N. Troullier and Y. Saad: Phys. Rev. Lett.
72(1994) 12401243.
8) J. R. Chelikowsky, N. Troullier, K. Wu and Y. Saad: Phys. Rev.
B 50(1994) 1135511364.
9) S. R. White, J. W. Wilkins and M. P. Teter: Phys. Rev. B
39(1989) 58195833.
10) J. Bernholc, J. Y. Yi and D. J. Sullivan: Faraday Discuss. Chem.
Soc. 92(1991) 217.
11) X. Jing, N. Troullier, D. Dean, N. Binggeli, J. R. Chelikowsky,
K. Wu and Y. Saad: Phys. Rev. B 50(1994) 1223412237.
12) J. R. Chelikowsky, X. Jing, K. Wu and Y. Saad: Phys. Rev. B
53(1996) 1207112079.
13) A. P. Seitsonen, M. J. Puska and R. M. Nieminen: Phys. Rev. B
51(1995) 1405714061.
14) E. L. Briggs, D. J. Sullivan and J. Bernholc: Phys. Rev. B
52(1995) R5471R5474.
15) E. L. Briggs, D. J. Sullivan and J. Bernholc: Phys. Rev. B
54(1996) 1436214375.
16) T. Hoshi, M. Arai and T. Fujiwara: Phys. Rev. B 52(1995)
R5459R5462.
17) T. Hoshi and T. Fujiwara: J. Phys. Soc. Jpn. 66(1997) 3710.
18) F. Gygi and G. Galli: Phys. Rev. B 52(1995) R22292232.
19) E. Tsuchida and M. Tsukada: Phys. Rev. B 52(1995) 5573
5578.
20) E. Tsuchida and M. Tsukada: Phys. Rev. B 54(1996) 7602
7605.
21) T. L. Beck: Rev. Mod. Phys. 72(2000) 10411080.
22) J. Kim, F. Mauri and G. Galli: Phys. Rev. B 52(1995) 1640
1648.
23) F. Mauri and G. Galli: Phys. Rev. B 50(1994) 43164326.
24) F. Mauri, G. Galli and R. Car: Phys. Rev. B 47(1993) 9973
9976.
25) K. Hirose and T. Ono: Phys. Rev. B 64(2001) 085105 15.
26) W. Kohn and L. J. Sham: Phys. Rev. 140(1965) A1133A1138.
27) D. R. Hamann, M. Schluter and C. Chiang: Phys. Rev. Lett.
43(1979) 14941497.
28) Y. Fujimoto and K. Hirose: Phys. Rev. B 67(2003) 195315 1
12.

Download Report