SX-4を使ったGussian94による
水素終端Si（１１１）表面初期酸化過程の検討
金島岳、弓達新治、奥山雅則
大阪大学基礎工学研究科
平成12年3月21日
１．はじめに
近年の計算機の非常な高速化により、材料探索や、反応解析のためにad-initio量子化学計算によるシュミレーションが広く用いられるようになっている。また計算方法の改良も合わせ、反応座標や遷移状態の大規模な計算も行なうことが可能となりつつある。しかし、大規模計算であり、多量の計算機資源を必要とすることに変わりはない。そこで実用的な計算を行ない、どの程度の計算機資源が必要であるかを調べた。ここでは、Si表面初期酸化過程の計算を行なった。Ｓｉ酸化膜はＬＳＩのゲート絶縁膜として使われており、その信頼性はＬＳｉのデバイスとしての信頼性に直結する。また、高集積化のために２mm程度の非常に薄いSi酸化膜が求められている。そのような極薄シリコン酸化膜を作製するためには、Si表面レベルの平坦化、およびその表面酸化過程の解明が必要である。Si表面の初期酸化過程の計算はいくつか報告があるが3,4,7-9、より現実的な系―水素終端Ｓｉ表面と酸素分子については殆ど報告がない。そこで、本研究では、水素終端Ｓｉ（111）と酸素分子との反応クラスター計算を行ない、初期化の様子を調べた。
２．計算方法、
計算は、すべて大阪大学大型計算機センターで行なった。使用した分子軌道計算ソフトウェアはＧａｕｓｓｉａｎ９４で、ＳＸ-４上で動作させた。プリ・ポストプロセッサはCerius２を使用した。本文中で使用しているクラスターの図は全てCerius２でPostscriptファイルに出力したものである。ｐ４クラスのbatchで実行し、以下のオプションにより実行情報を収集した。- F_PRINDIFF=NO
export F_PRINDIFF
F_PROGINF＝DENTAL
export F_PROGINF
また、半経験的計算により反応の概略を調べ、その後、密度汎関数法（DFT）により精度を上げた計算を行なった。基底関数は3-21Gを使用した。
図１に使用したクラスターを示す。これは、水素終端Ｓｉは全て自由に動けるようにし、周りの水素を固定することで結晶構造を維持するようにした。Si-Si距離の初期値は実測値2.35Aとし、Si-H距離は、1.48Aに固定した。（Si-t中心の Si-中心のＨ）面に対して対称とし、酸素および表面水素の２面角は固定した。
３．結果
３．１ Ｓｉ－Ｏ－Ｓｉ結合の計算
Ｓｉの表面初期酸化においては、ｍバックボンドから酸化されるといわれている。そこで、まずバックボンドが１つ酸化された構造をab-inition密度汎関数法により求めた。Becke'ｓ three-parameter fomulation（B3LYP)を使用して計算を行なった。（計算時間は付録A参照。）
図２に安定構造を示す。Si-O距離は1.6Aと、バルクシリカガラス（アモルファルSiO2）において約144°といわれているが、Si/SiO2界面ではストレスやサブオキサイドの影響などによりこの値より小さくなるとの報告がある。本計算は135°となっており、実験結果と良い一致を示している。よって、B3LYP/3-21Gで水素終端Si表面酸化状態（安定構造）は再現可能であると考えられる。
３．２ 半経験的計算によるO2分子反応の計算
Siバックボンドが酸化された状態は実験と良く一致することが示されたので、次に酸素分子による表面初期酸化過程についての計算を行なった。非経験的計算は、計算時間がかかるため、まずおよその構造を計算時間の短い半経験的計算で探した。計算は全てGaussian９４で行ない、パラメータはAM1とした。（計算時間は付録B参照。）図３に計算された安定状態を、図４に遷移状態を示す。安定状態では、O原子がSi-Siの間に入り込んでおりO-O距離はほぼ１.２Aと実験値と近い値になっている。また、酸素が３配位になっていることが計算された。
また、反応前の構造を求めた。SiクラスターとO２の距離を十分離して、両者のクラスターの構造安定化をさせた結果を図５に示す。Si-Si距離は２.４Aと、Si結晶における結合長の実験値と良い一致を示しており、この状態でO2からの影響を殆ど受けていないと考えられる。
３．３ 密度汎関数法によるO2反応計算
活性化エネルギーを求めるために密度汎関数法を用いて水素終端Si表面クラスターとO2分子との反応を調べた。初期構造はAMIで求めたものを利用した。Gauusian９４ではスピン状態を指定する必要があるが、酸素分子で基底状態のスピン状態がtripletであることより、tripletを指定し、非制限計算を行なった。
まず、反応する前のエネルギーを求めた。AM1の結果（図５）を初期構造としてB３LYP/3-21Gにより構造の最適化を行なった。なお、SiクラスタとO２の距離は固定した。図６は計算結果を示す。Si-Si結合距離は2.33Aとほぼ実験値と一致しており、O２からの影響を殆ど受けていないと考えられる。次に、AM１の結果（図３）を初期構造として、構造安定化を行なった。（計算時間は付録C参照。）結果を図７に示す。酸素分子が解離し、Si-O-Si結合が形成されて安定化することが示された。
最後に遷移状態の計算および、振動解析を行なった。図８に遷移状態の計算結果を示す。酸素およびSi間の距離が実験値よりわずかに大きくなっている。振動解析を行なったところ、酸素が乖離する方向に虚の振動を持つことが示され、遷移状態である可能性が示唆されたが、IRC計算により反応経路の確認を行なう必要があると思われる。（計算時間は付録D参照。）
活性化エネルギーの変化をまとめると、図９のようになる。酸素分子がSiバックボンドを１つ酸化するためには1.11eV必要であり、Si-O-Si＋Ｏで安定することが示された。（計算時間は付録Ｅ参照。）水素終端されておらず、ダングリングボンドを持つSiアドアトムへO２分子が反応するときの活性化エネルギーは既に報告されている。これによると、O２分子が表面に接近し、物理吸着するときの活性化エネルギーが2.5eV、Si-Si結合をきるためのエネルギー1.2eVである。さらに、Si-Hを切る反応は起こりにくいことにより、水素終端Si(111)表面においては、今回計算した水素終端Siのバックボンドを直接攻撃する反応が比較的起こり易いことが示唆された。
A Si-O-Si結合の計算時間
Si-O-Si結合計算に使用したG94のコマンド、ファイルサイズおよび計算時間を示す。構造最適化した変数は結合距離５、結合角５、2面各1の合計11個である。#B3LYP/3-21G* FChk=OptInt Opt=(Z-Matrix,，
MaxCycle=50) SCF=MaxCycle=256
File lengths (MBytes) ： RWF= 42 Int＝ 0
D2E= 0 Chk= 2 SRC= １
****** Program Information ******
| User Time (sec) | ： | 3434．969448 |
| Sys Time (esc) | ： | 213．581399 |
| Vector Time (esc) | ： | 1597．167053 |
| Inst．Count | ： | 20221451806 |
| V．Inst．Count | ： | 19402393258 |
| V．Element Count | ： | 2177811139301 |
| FLOP Count | ： | 828725030854 |
| MOPS | ： | 687．232546 |
| MFLOPS | ： | 241．261252 |
| VLEN | ： | 112．244459 |
| V．Op．Ratio (%) | ： | 92．255787 |
| Memory Size (MB) | ： | 958．031250 |
| MIPS | ： | 58．869243 |
| I-Cache (sec) | ： | 110．064745 |
| O-Cache (sec) | ： | 87．958163 |
| Bank (sec) | ： | 4．736364 |
B 半経験的計算によるO２分子反応の計算時間
半経験的計算によるO２分子反応の計算の流れを述べる。構造最適化した変数は結合距離６、結合角６、２面各１の合計１１個である。全てGaussian９４のAM1で計算を行なった。１）前節で求めたSi-O-SiのOの位置にO２分子をおき構造最適化する。２）十分離れた点から構造最適化した位置までをGaussian９４ののSCANコマンドにより反応座標を計算する。３）遷移状態をGaussian９４のSCANコマンドにより計算する。４）遷移状態が正しい方向を向いてるかどうか、Gaussian９４のFREQコマンドで振動解析を行なう。以下の例は、構造最適化を行なったときの、実行情報を示す。この系においては典型的な出力である。
#AM１ FCHK=OPTINT OPT=(Z-MATRIX，
MAXCYCLE=99) SCF=(MAXCYCLE=256)
File lengths (MBytes) ： RWF＝ 5 Ｉｎｔ＝ 0
Ｄ２Ｅ＝ 0 Chk＝ 2 Scr＝ １
****** Program Information ******
| User Time (sec) | ： | 28.151421 |
| Sys Time (esc) | ： | 53．451633 |
| Vector Time (esc) | ： | 6．995417 |
| Inst．Count | ： | 1646169700 |
| V．Inst．Count | ： | 99896968 |
| V．Element Count | ： | 3787800383 |
| FLOP Count | ： | 331474353 |
| MOPS | ： | 189．477933 |
| MFLOPS | ： | 11．774693 |
| VLEN | ： | 37．917071 |
| V．Op．Ratio (%) | ： | 71．011407 |
| Memory Size (MB) | ： | 958．031250 |
| MIPS | ： | 15．475545 |
| I-Cache (sec) | ： | 1．760194 |
| O-Cache (sec) | ： | 0．854065 |
| Bank (sec) | ： | 0．000097 |
計算時間は、26-140秒でVector TimeはUser Timeの約20%であった。これらの計算に使用した時間の総計は1,500秒であった。Vector化は悪いものの、密度汎関数法に比較して約2桁高速であり、試行錯誤には都合よいと思えわれる。
Ｃ 密度汎関数法によるO２反応計算時間
密度汎関数法で安定構造を求めたときの実行情報を示す。
#B3LYP/3-21G* FCHK=OPTINT OPT=(Z-MATRIX，
MAXCYCLE=99) SCF=(MAXCYCLE=256)
File lengths (MBytes) ： RWF＝ 89 Ｉｎｔ＝ 0
Ｄ２Ｅ＝ 0 Chk＝ 2 Scr＝ １
****** Program Information ******
| User Time (sec) | ： | 16978．81742 |
| Sys Time (esc) | ： | 588．353947 |
| Vector Time (esc) | ： | 9539．876277 |
| Inst．Count | ： | 917894199878 |
| V．Inst．Count | ： | 121869909505 |
| V．Element Count | ： | 15977102296322 |
| FLOP Count | ： | 6287842959684 |
| MOPS | ： | 987．826546 |
| MFLOPS | ： | 370．334564 |
| VLEN | ： | 131．091443 |
| V．Op．Ratio (%) | ： | 95.253886 |
| Memory Size (MB) | ： | 958.031250 |
| MIPS | ： | 54．061138 |
| I-Cache (sec) | ： | 537．328854 |
| O-Cache (sec) | ： | 369．029995 |
| Bank (sec) | ： | 20．929611 |
D 振動解析計算時間
振動解析を行なったときの実行情報を示す。構造最適化は行なわず､FREQコマンドにより振動解析のみを行なった。
#B3LYP/3-21G* FCHK=OPTINT EREQ
SCF=(MAXCYCLE=256)
File lengths (MBytes) ： RWF＝ 70 Ｉｎｔ＝ 0
Ｄ２Ｅ＝ 0 Chk＝ 2 Scr＝ １
****** Program Information ******
| User Time (sec) | ： | 5655．551585 |
| Sys Time (esc) | ： | 23．946464 |
| Vector Time (esc) | ： | 3805．245699 |
| Inst．Count | ： | 293472691827 |
| V．Inst．Count | ： | 10241190763 |
| V．Element Count | ： | 7456240395731 |
| FLOP Count | ： | 3899623849783 |
| MOPS | ： | 1352．176011 |
| MFLOPS | ： | 689．521401 |
| VLEN | ： | 72．806381 |
| V．Op．Ratio (%) | ： | 97．501592 |
| Memory Size (MB) | ： | 958.031250 |
| MIPS | ： | 51．891082 |
| I-Cache (sec) | ： | 62．999507 |
| O-Cache (sec) | ： | 420．639962 |
| Bank (sec) | ： | 2．320988 |
E 密度汎関数法計算時間のまとめ
DFTによる活性化エネルギーの計算を行なうのに投入したjobの回数および合計時間を示す。遷移状態は収束しないことがあり、初期構造を多少変化させながら試行を行なった。
- 安定構造計算（1回）17,000秒
- 遷移状態計算（5回）45,000秒
- 離れた状態計算（3回）29,000秒
一般に密度汎関数法は計算時間が短いといわれているが、この程度の系においてはB3LYP/3-21G*とMP2/3-21G*の間に大きな差は見られなかった。Siバックボンドを１つ酸化したときの安定状態の計算を示す。
制限密度汎関数法（B3LYP/3-21G*）
#B3LYP/3-21G* FCHK=OPTINT OPT=（Z-MATRIX，
MAXCYCLE=99） SCF=(MAXCYCLE=256)
File lengths (MBytes) ： RWF＝ 49 Ｉｎｔ＝ 0
Ｄ２Ｅ＝ 0 Chk＝ 2 Scr＝ １
****** Program Information ******
| User Time (sec) | ： | 6461．447791 |
| Sys Time (esc) | ： | 318．158986 |
| Vector Time (esc) | ： | 3067．704229 |
| Inst．Count | ： | 379340356253 |
| V．Inst．Count | ： | 37080793106 |
| V．Element Count | ： | 4325944324259 |
| FLOP Count | ： | 1654690316529 |
| MOPS | ： | 722．470263 |
| MFLOPS | ： | 256．086619 |
| VLEN | ： | 116．662670 |
| V．Op．Ratio (%) | ： | 92．668282 |
| Memory Size (MB) | ： | 958．031250 |
| MIPS | ： | 58．708260 |
| I-Cache (sec) | ： | 191．902709 |
| O-Cache (sec) | ： | 175．077964 |
| Bank (sec) | ： | 14．084807 |
非制限密度汎関数法（UB3LYP/3-21G*）
#UB3LYP/3-21G* FCHK=OPTINT OPT=（Z-MATRIX，
MAXCYCLE=99） SCF=(MAXCYCLE=256)
File lengths (MBytes) ： RWF＝ 89 Ｉｎｔ＝ 0
Ｄ２Ｅ＝ 0 Chk＝ 2 Scr＝ １
****** Program Information ******
| User Time (sec) | ： | 7968．087843 |
| Sys Time (esc) | ： | 336．108290 |
| Vector Time (esc) | ： | 4416．166406 |
| Inst．Count | ： | 432174056421 |
| V．Inst．Count | ： | 58230301930 |
| V．Element Count | ： | 7482076983119 |
| FLOP Count | ： | 2953939444955 |
| MOPS | ： | 985．935508 |
| MFLOPS | ： | 370．721245 |
| VLEN | ： | 128．491125 |
| V．Op．Ratio (%) | ： | 95．240036 |
| Memory Size (MB) | ： | 958．031250 |
| MIPS | ： | 54．238114 |
| I-Cache (sec) | ： | 247．366595 |
| O-Cache (sec) | ： | 187．891118 |
| Bank (sec) | ： | 14．815780 |
非制限摂動計算（MP2/3-21G*）
#UMP2/3-21G* FCHK=OPTINT OPT=（Z-MATRIX，
MAXCYCLE=99） SCF=(MAXCYCLE=256)
File lengths (MBytes) ： RWF＝ 89 Ｉｎｔ＝ 0
Ｄ２Ｅ＝ 0 Chk＝ 2 Scr＝ １
****** Program Information ******
| User Time (sec) | ： | 5261．573834 |
| Sys Time (esc) | ： | 188．718329 |
| Vector Time (esc) | ： | 2612．641821 |
| Inst．Count | ： | 305342073567 |
| V．Inst．Count | ： | 31872348429 |
| V．Element Count | ： | 375470309679 |
| FLOP Count | ： | 1399103395853 |
| MOPS | ： | 765．583255 |
| MFLOPS | ： | 265．909676 |
| VLEN | ： | 117．804407 |
| V．Op．Ratio (%) | ： | 93．211073 |
| Memory Size (MB) | ： | 958．031250 |
| MIPS | ： | 58．032460 |
| I-Cache (sec) | ： | 48．869679 |
| O-Cache (sec) | ： | 65．142536 |
| Bank (sec) | ： | 17．710712 |
参考文献
- Ｊ．Ｂ．Foresuman and AE．Frisch： Exploring Chemistry with Electronic Structure Methoods，Gaussian Inc．，Pittsburgh，second edition（1993）．
- M．J．Frisch，G．W．Trucks，H．B．Schlegel，P．M．W．Gill，B．G．Johnson，M．A．Robb，J．R．Cheeseman，T．Keith，G．A．Petersson，J．A．Montgomery，K．Raghavachari，M．A．Al-Laham，V．G．Zakrzewski，J．V．Ortiz，J．B．Foresman，J．Cioslowski，B．B．Stefanov，A．Nanayakkara，M．Chllcombe，C．Y．Peng，P．Y．Ayala，W．Chen，M．W．Wong，J．L．Andres，E．S．Replogle，R．gomperts，R．L．Martin，D．J．Fox，J．S．Binkley，D．J．Defrees，J．Baker，J．P．Stewart，M．Head-Gordon，C．gonzalez，and J．A．P．and：Gaussian 94，Gaussian，Inc．，Pittsburgh PA（1995）
- T．Hoshino：Phsy．Rev．B 59（1999）2332-2340．
- T．Hoshino：Phsy．Rev．B 59（1999）2332-40．
- H．Ikeda，K．Hotta，T．Yamada，S．Zaima，H．Iwano and Y．Yasuda：J．Appl．Phys．77（1995）5125-5129．
- H．Ikeda，K．Hotta，T．Yamada，S．Zaima，Y．Yasuda：Jpn．J．Appl．Phsy．34（1995）2191-2195．
- K．Sakata，A．Tachibana，S．Zaima and Y．Yasuda：Jpn．J．Appl．Phsy．37（1998）4962-4973．
- B．B．Stefanov and K．Raghavachari：Appl．Phys．Lett．73（1998）824-826．
- S．Yudate，N．Okada，T．Kanashima，M．Okuyama and K．Nihihara：Precision Science and Technology for Perfect Surfaces （eds．Y．Furukawa，Y．Mori and T．Kataoka），Osaka，Japan（1999），The Japan Society for Precision Engineering（JSPE），9th Int．Natl．Conf．on Production Engineering．
- 国立天文台（編）：理科年表，63，丸善，東京都（1989）．
- 義人，真司：固体物理２３（1988）888-899．