[共同研究成果]
スーパーコンピュータを用いたシミュレーションによる超精密加工プロセス・ナノ電子デバイスの解析
大阪大学

大学院工学研究科
広瀬喜久治※

物質・材料研究機構

後藤英和

稲垣耕司

ナノマテリアル研究所

小野倫也

※

塚本茂

我々は、量子力学に基づくコンピュータシミュレーションプログラムを開発するとともに、これを用いて加工プロセスやナノ電子デバイスの電気的特性を解明する研究を行っている。本稿では、主流になりつつある Fortran90 を用いた高速化で問題となった事例について述べるとともに、シミュレーションを利用した成果について述べる。
１．はじめに
我々は、近年実用的になりつつある「第一原理分子動力学法」と呼ばれる量子力学に基づいて固体内部や表面の電子状態を精度よく計算する方法 1)~6)を用いて、超精密加工という一見実験的研究中心の分野で、その中の現象を原子レベルから理解しようと、計算機シミュレーションによる解析を行っている。
本稿では、主流になりつつある Fortran90 を用いてプログラミングした場合にベクトル化の最適化で実際に問題となった事例について述べるとともに、金属ナノワイヤーの電気伝導の解析や、超純水電気化学的加工法におけるシミュレーションを利用した成果について述べる。
２．Fortran90 シミュレ−ションプログラムの高速化における問題点
最近では、数値計算用コンピュータ言語の主流が FORTRAN77 から FORTRAN90に移り、スーパーコンピュータ SX シリーズにおいても、コンパイラは f90(またはsxf90) を 使 用 す る こ と が 推 奨 さ れ る よ う に な っ て い る 。 こ の FORTRAN90 はFORTRAN77 と比べて、配列式が使えるという点で非常に便利である。例えば、図 1
PARAMETER N=1000
REAL*8 A(1:N,1:N),B(1:N,1:N),C(1:N,1:N)
DO I=1,N
DO J=1,N
A(I,J)=B(I,J)+C(I,J)
END DOEND DOINTERGER,PARAMETER :: N=1000
REAL(8),DIMENSION(1:N,1:N) :: A,B,C
A(:,:)=B(:,:)+ C(:,:)

FORTRAN77FORTRAN90

図 1 FORTRAN90 における配列式の記述

に示すように、FORTRAN77 では複数個の DO 文を使って 5 行に渡って記述していた 2 次元配列の計算式が、FORTRAN90 では配列の演算式だけで 2 次元配列の演算ができるようになっている。配列式を使用することによってコードが簡潔になり、デバッグやコードの解読が容易になるというメリットがある。
しかしながら、図 1 右のような配列のすべての要素に対して演算を行う配列式を含むコードをコンパイルすると、ベクトル化が可能であるにもかかわらず、自動ベクトル化される場合とされない場合があることが分かった。以下にその例を示す。
1. 自動ベクトル化されない場合
REAL(8),DIMENSION(1:81,‑16:15,0:3,0:7) :: AR,AIREAL(8),DIMENSION(1:81,‑16:15,‑16:15) :: R,I
DO J=0,7AR(:,:,0,J)= R(:,:,‑16+J)+R(:,:,‑8+J)+R(:,:,J)+R(:,:,8+J)
AR(:,:,1,J)=‑R(:,:,‑16+J)‑I(:,:,‑8+J)+R(:,:,J)+I(:,:,8+J)
AR(:,:,2,J)= R(:,:,‑16+J)‑R(:,:,‑8+J)+R(:,:,J)‑R(:,:,8+J)
AR(:,:,3,J)=‑R(:,:,‑16+J)+I(:,:,‑8+J)+R(:,:,J)‑I(:,:,8+J)
END DO

この類のコードを sxf90 コマンドでコンパイルすると自動ベクトル化されない旨のメッセージが表示される。
2. 自動ベクトル化される場合
REAL(8),DIMENSION(1:81,‑16:15,0:3,0:7) :: AR,AIREAL(8),DIMENSION(1:81,‑16:15,‑16:15) :: R,I
DO J=0,3
R(:,:,‑8+J)= AR(:,:,J,0)+AR(:,:,J,4)‑AR(:,:,J,2)‑AR(:,:,J,6)
R(:,:,‑4+J)= AR(:,:,J,0)‑AR(:,:,J,4)‑AI(:,:,J,2)+AI(:,:,J,6)
R(:,:,J)= AR(:,:,J,0)+AR(:,:,J,4)+AR(:,:,J,2)+AR(:,:,J,6)
R(:,:, 4+J)= AR(:,:,J,0)‑AR(:,:,J,4)+AI(:,:,J,2)‑AI(:,:,J,6)
END DO

この類のコードを sxf90 コマンドでコンパイルすると自動ベクトル化される旨のメッセージが表示され、配列の 1 次元から 3 次元目までを一重ループ化※の上、ベクトル化を行う。
どちらのケースでも 1 次元目と 2 次元目はすべての要素に対して演算を行っており、ベクトル化が可能であるはずだが、このように自動的にベクトル化されない場合が存在する。これらの 1 次元目と 2 次元目を自動ベクトル化させるために、ループ一重化を行い、次のように上記のそれぞれのコードを変更した。
1'. 配列式で記述すると自動ベクトル化されなかったコード
REAL(8),DIMENSION(1:81,‑16:15,0:3,0:7) :: AR,AI
REAL(8),DIMENSION(1:81,‑16:15,‑16:15) :: R,I
DO J=0,7
DO K=1,2592
AR(K,‑16,0,J)= R(K,‑16,‑16+J)+R(K,‑16,‑8+J)+R(K,‑16,J)+R(K,‑16,8+J)※ 「一重ループ化」法は、配列の最初の二つ以上の次元を処理するために入れ子になっている複数の DO ループを一体化して取り扱う方法のことである。ベクトル長が長くなり効率が上がるため Fortran における最適化の手法として一般的に使われている。配列変数が記憶領域に割り付けられる際に第一添字、第二添字の順に配置される Fortran の特性を利用している。

AR(K,‑16,1,J)=‑R(K,‑16,‑16+J)‑I(K,‑16,‑8+J)+R(K,‑16,J)+I(K,‑16,8+J)
AR(K,‑16,2,J)= R(K,‑16,‑16+J)‑R(K,‑16,‑8+J)+R(K,‑16,J)‑R(K,‑16,8+J)
AR(K,‑16,3,J)=‑R(K,‑16,‑16+J)+I(K,‑16,‑8+J)+R(K,‑16,J)‑I(K,‑16,8+J)
END DOEND DO

配列式で記述しても自動ベクトル化されたコード
REAL(8),DIMENSION(1:81,‑16:15,0:3,0:7) :: AR,AI
REAL(8),DIMENSION(1:81,‑16:15,‑16:15) :: R,I
DO J=0,3
DO K=1,2592
R(K,‑16,‑8+J)= AR(K,‑16,J,0)+AR(K,‑16,J,4)‑AR(K,‑16,J,2)‑AR(K,‑16,J,6)
R(K,‑16,‑4+J)= AR(K,‑16,J,0)‑AR(K,‑16,J,4)‑AI(K,‑16,J,2)+AI(K,‑16,J,6)
R(K,‑16, J)= AR(K,‑16,J,0)+AR(K,‑16,J,4)+AR(K,‑16,J,2)+AR(K,‑16,J,6)
R(K,‑16, 4+J)= AR(K,‑16,J,0)‑AR(K,‑16,J,4)+AI(K,‑16,J,2)‑AI(K,‑16,J,6)
END DOEND DO

コンパイルを行うと、これらはどちらも最内側の DO ループのみが自動ベクトル化される旨のメッセージが表示された。また、2’のコードは 1 次元目と 2 次元目のみが一重ループ化の対象となり、3 次元目はベクトル化されなかった。これらをまとめると次のようになる。
配列式を採用する配列式を採用しないベクトル化しない場合がある。
ベクトル化した場合、ベクトル長が長くなる。
ほとんどの場合、ベクトル化する。
3 次元目はベクトル化されない。

以上のように、現在の FORTRAN90 コンパイラは、新たに採用された配列式の自動ベクトル化において、必ずしも最善の最適化ができない場合があり注意を要する。
今後のコンパイラの最適化能力の向上に期待したい。
３．ナノ電子デバイスの電気特性の解析や超精密加工プロセスの解析に応用した例
3.1 金属ナノワイヤーの電気伝導特性の解析
近年の電子デバイスにおける小型化の傾向はとどまることを知らず、原子スケールで制御された超小型の電子デバイスまでもが要求されるようになってきている。このような固体中の電子の平均自由行程よりもはるかに小さいサイズのデバイスでは、電子の運動は、原子核に散乱(フォノン散乱)されることのない弾道的な(ballistic)運動となる。この運動には電子の波動性や粒子性などの量子効果が顕著に現れ、従来の古典電磁気学では予想されない奇異な電気伝導特性が生じる(例:トンネル伝導現

図 2 ナノワイヤーにおける電気伝導象・クーロン閉塞等)ことが明らかになってきている。最近では、これらの量子効果を含む電気伝導特性を、デバイスの基本動作原理として利用した新しい電子デバイスの研究・開発も始まっている。

図 3 Na ナノワイヤーモデル

我々は、原子スケールで作成されたナノ構造物の電気伝導特性に対して理論的な面からアプローチするため、量子力学の第一原理に基づく電子状態計算手法を開発し、スーパーコンピュータを用いて実際に解析を行った。電気伝導特性の解析は、図 2 に示すように、電気伝導計算の対象となるナノ構造とその両端に接続した半無限電極からなるモデルに、左右の無限遠から電子が入射すると、ナノ構造を透過または反射する、という一連の過程を計算することによって行った。さまざまな入射エネルギーや入射方

図 4 Na ナノワイヤーのコンダ
クタンスの変化 7)向の電子に関してシミュレートする必要があるため、計算量は膨大になるが、最近の発展著しいスーパーコンピュータの高度な処理能力と並列計算技術、そして計算精度の劣化無しに計算量を節約する新しい計算手法の採用によって、現実的な時間内での計算に成功した 7) 。
ここでは本計算手法によるナトリウム原子鎖(ナノワイヤー)の電気伝導特性の解析について紹介する。実験によるナトリウムナノワイヤーの電気伝導の測定は 1995 年に Krans 等 8)によって初めてなされた。彼等はナトリウムの微小板片を引き伸ばす過程でできるナノワイヤーの電気伝導度(コンダクタンス)が、G0=2e2/h (e: 電荷素量、h: プランク定数)を量子化単位として、6G0, 5G0, 3G0, 1G0 と量子的に減少し、さらに量子化コンダクタンス 4G0,2G0 が現れないという消滅則をも見出した。我々は、この実験事実に対して、図 3 に示すナトリウムナノワイヤーモデルを用いて、ナノワイヤーの長さ(N: ナノワイヤーに含まれる Na 原子数)に対するコンダクタンス G の変化を調べた。その結果、図 4 に示すように、コンダクタンスは N=1 のときのみ 3G0を示し、それよりも長いナノワイヤー(N=2-5)では 1G0 を示した。つまり、量子化コンダクタンス 2G0 についての消滅則が理論的に示された。
次に、N=5 のナノワイヤーを含むナノ構造に注目し、その両端に有限バイアス電圧を印加して、そのときの電流-電圧特性をシミュレートした。図 5 上図は計算結果として得られた電流-電圧特性、図 5 下図はその微分コンダクタンスを示している。印加バイアス電圧が比較的低い場合(0.0~0.5V)には、電流-電圧特性の傾きはほとんど一定であるように見えるが、微分コンダクタンスは急激に減少しており、有限バイアス電圧を印加することでコンダクタンスの量子化が解消されている。これは、本来ナノワイヤーが持つ電子の固有準位に合致しないエネルギーの入射電子は透過率 1 でナノワイヤーを通過できず、有限バイアス電圧下では、その図 5 Na ナノワイヤーの電流-電圧特ような入射電子の電気伝導特性がコンダク性と微分コンダクタンス 7)タンスに反映されるからである。なお、N=5のナノワイヤーにおいてバイアス電圧が 0.0V のときに 1G0 であるのは、ナノワイヤーの固有準位がフェルミレベル付近にあり、電気伝導に寄与するフェルミエネルギーをもった電子の透過率が 1 であるためである。また、バイアス電圧 2.5V 付近では、電圧の上昇にもかかわらず電流値の減少が発生し、負の微分コンダクタンスが得られた。一般的には、両電極間の電位差が大きくなると、一方の電極から他方の電極へと一方的に流れる電子の量は多くなり、コンダクタンスは増加の一途をたどるはずである。しかし、この場合(電圧 2.5V 付近)では、電圧の上昇とともに、負の電圧が印加されている側の Na 表面のポテンシャルが上昇したため、ナノワイヤーの固有準位に合致する入射電子が、ナノワイヤーに到達する前に反射され、その結果、全体的に電流に寄与する電子の数の減少を招いていることが分かった。
以上、量子力学に基づく数値計算によってナノスケール電子デバイスの基本的な特性解析が可能になりつつあることを述べた。今回紹介したナノワイヤーはナノデバイ
スの中でも特に単純なものであるが、今後のスーパーコンピュータのさらなる計算能力の増大により、さらに複雑なデバイスの解析が可能になっていくものと大いに期待している。
3.2 超純水のみによる電気化学的加工プロセスのシミュレーション― Si(001)表面における加工現象の発見―「超純水のみによる電気化学的加工法」は、通常の電解液を利用した電解加工法とは全く異なる加工原理の新しい電気化学的加工法であり、超純水中の OH―イオンの化学的作用を用いて、被加工物表面の原子を除去加工するものである 9)。加工物表面

の物性的損傷や汚染が無く、また化学薬品の廃液を伴わないため、半導体産業を始めとする先端加工分野での応用が期待できる。最近、陰極表面における反応素過程の第一原理分子動力学シミュレーションによる解析が行われ、被加工物を陽極とする従来の方法では加工できなかった Si が、被加工物を陰極とする方法によって加工が可能になることが予測された 10)。ここでは、この結果について紹介する。
超純水中において電解加工を行う場合、陰極表面には常に H+イオンが供給されている。このイオンは、陰極から電子を受け取り、中性の H 原子となる。また、H2Oも陰極表面近傍に存在している。そこで、水素終端化 Si(001)表面に 2 個の H2O 分子と 1 個の H 原子を配置してシミュレーションを行った 10)。図 6(a)(b)に、初期状態と構造緩和後の原子配置と電子密度分布を示す。この結果から、H 原子が H2O 分子と反応し、OH と H2 分子が生成されることがわかる。また、生成した OH は、Si 表面原子と結合し、表面原子の 2 本のバックボンドの 1 本を切断している。以上の結果か
ら、陰極 Si(001)水素終端化表面においては、
H2O + H → OH + H2a(1)

b

H

H 2O

H2H 2O

H2O
バックボンドOH(1)バックボンド

CB

C

Bzzyy
OOxxc

バックボンドd

バックボンド

OH(2)

OH(1)

OH(2)

H

A

HOH(1)AB

CBz

CzyyOOxx

図 6 超純水中の電解加工において、Si(001)を陰極とした場合の表面反応素過程の第一原理分子動力学シミュレーション結果 10)

によって OH が生成され、生成された OH が Si 表面原子と結合し、1 本のバックボンドを切断することがわかった。
そこで次に、Si(001)水素終端化表面に 2 個の OH を作用させてバックボンドを切断し、同時に H も作用させて、切断されたバックボンドの終端化を行った。図 6(c)(d)に、初期状態と構造緩和後の原子配置と電子密度分布を示す。Si 表面原子(A)の 2 本のバックボンド(A−B、A−C)は完全に切断され、Si 原子が上昇した。以上の結果か
ら、電解加工中に陰極表面に供給される H 原子と、H と H2O との反応によって生成された OH が Si 表面原子に作用することよって、表面原子が SiH2(OH)2 分子として除去加工されることがわかった。
現在、超純水中の水分子を解離する触媒材料を用いて OH−イオンや H＋イオンの密度を増加させることにより、電解電流密度を常温・常圧超純水の 100 万倍近くまで増加させることが可能になっている 9)。図 7 に、触媒材料を用いた超純水のみによる電気化学的加工法の概念を示す。このような方法によって、加工物を陰極として加工実験を行った結果を図 8 に示す。シミュレーションで予測されたように、加工物を陰極とすることによって、Si や Al の加工が可能であることが実証された 11)。現在、大阪大学超精密加工研究拠点ウルトラクリーンルームでは、超純水電解加工システムが試作・設置され、実用化を目指した研究・開発が行われている。
陽極
超純水

触媒材料
反応生成物

陰極
加工物原子

酸素原子

水素原子

3.0
陰極
Al
2.5
2.01.51.00.5
Fe0.0Mo-0.5
0.0 0.5 1.0 1.5 2.0
電流密度（A/cm2）

加工速度（μm/min）35
陽極3025Mo201510 CuFe5Al0Si-5
0.0 0.5 1.0 1.5 2.0
電流密度（A/cm2）

加工速度（μm/min）

加工速度（μm/min）

図 7 触媒材料を利用した加工法の概念（被加工材料が陰極の場合）
0.06
0.05

陰極

0.04
0.03
0.02
0.01
0.00
0.0

Si

0.1
0.2
0.3
電流密度（A/cm2）

図 8 各種材料の加工特性 11) （加工速度が負の領域は酸化膜形成を表す）

４．おわりに
本稿では、Fortran90 を用いた高速化で問題となった事例について述べるとともに、本研究で開発されたシミュレーションプログラムを利用したナノワイヤーの電気伝導特性の解析ならびに超純水電気化学加工法における Si の陰極でのエッチング現象の発見について述べた。本年度は東北大学のスーパーコンピュータの機種更新が行われるということで、さらなる大規模計算が実現可能になることを期待している。
本研究は文部省中核的研究拠点(COE)形成プログラムである大阪大学・超精密加工研究拠点における共同研究の一環として行われた。また、計算ソフトウェアの並列化・高速化については東北大学情報シナジーセンター殿に多大な協力をいただいた。

参考文献
1) W.Kohn and L.J.Sham: Self-consistent Equations Including Exchange andCorrelation Effects, Phys. Rev., 140, 4A (1965) 1133.
2) N.Troullier and L.Martins: Efficient Pseudopotentials for Plane-wave
Calculations, Phy.Rev., B43, 3 (1991) 1993.
3) M.C.Pyne, M.P.Teter, D.C.Allan, T.A.Arias and J.D. Joannopoulos: Iterative
Minimization Techniques for ab initio Total-energy Calculations: MolecularDynamics and Conjugate Gradients, Rev. Mod. Phys., 64,4 (1992) 1045.
4) K. Hirose: First-Principles Molecular-Dynamics Simulations Aiding Creationof Perfect Surfaces, Proc. of the 9th ICPE, p. 913 (1999).
5) Tomoya Ono and Kikuji Hirose: Timesaving Double-Grid Method forReal-Space Electronic-Structure Calculations, Phys. Rev. Lett. Vol.82,
No.25 ,p.5016 (1999).
6) K.Hirose and T.Ono: Direct minimization to generate electronic states withproper occupation numbers, Phys. Rev. B, Vol. 64, No.8, 085105.
7) S. Tsukamoto and K. Hirose: Electron-Transport Properties of Na Nanowires
under Applied Bias Voltages, Phys. Rev. B, in print.
8) J. M. Karns, J. M. van Ruitenbeek, V. V. Firsun, I. K. Yanson, and L. J. deJongh: The signature of conductance quantization in metallic point contacts,Nature 375, p.767 (1995).
9) 森 勇藏、後藤英和、広瀬喜久治、小畠厳貴、當間 康、森田健一：超純水のみによる電気化学的加工法の研究―触媒反応を利用した超純水中の OH−イオンの増加方法−、精密工学会誌、67,6 p.932 (2001).
10) 後藤英和、広瀬喜久治、小畠厳貴、當間 康、森 勇藏：超純水のみによる電気化学的加工法の研究―陰極 Si(001)表面における除去加工現象の第一原理分子動
力学シミュレーション−、精密工学会誌、67,10 p.1680 (2001).
11) 森 勇藏、後藤英和、広瀬喜久治、小畠厳貴、當間 康、西村丈人：超純水のみによる電気化学的加工法の研究―陰極表面における加工現象−、精密工学会誌、
68,9(2002).