J. Plasma Fusion Res. Vol.89, No.12 (2013)８
３５‐８４４

解説

分子軌道法を用いたプロセスプラズマへのアプローチ
An Approach of Molecular Orbital Calculations to Process Plasmas
林俊雄

HAYASHI Toshio
名古屋大学プラズマナノ工学研究センター（原稿受付：２
０
１
３年７月１
１日）

分子軌道法の発展とパーソナルコンピュータの発展は目覚ましく，これまで困難と思われていた分子物性の計算が机上で比較的容易にできるようになってきた．一方，エッチング等のプロセス開発はプラズマプロセス中
で起っている多くの現象が解明されないまま進められてきた．今後の開発をさらに発展させるためにはプロセス中で起っている現象解明が必須である．このため，未経験者でも分子軌道法を用いることができるように初歩的な解説を試みる．
Keywords:
molecular orbital calculation, excitation energy, ionization energy, electron attachment energy, potential surface,electron collision, dissociation, Gaussian09, process plasma, chemical reaction in etching process.

１．はじめにの低い誘電体を放電部に用いている場合には２００℃程度に分子軌道法による計算について述べる前にプラズマプロ

なることもあるので表面での熱化学反応が起る可能性はあセスにおける分子の一次解離過程について考えてみる．プる．それでも気相中での衝突頻度が低いのでガス温度は解ラズマ中には熱，UV，イオン，電子，原子，分子が存在す

離するほど高くはならない．振動準位の脱励起速度はる．この内，ガス（分子）の温度は〜450 K 程度［１‐５］であ10−9 s よりも早いと考えられる．脱励起する前に次の衝突るので分子の熱解離は無視できる．気相中におけるイオが起り，より高い準位への励起が起ればやがては解離エネ

ン・分子反応も低圧プラズマでは衝突頻度が低いので無視

ルギーを超えることになる（N2 分子は赤外不活性なので脱できる（粘性流領域では無視できない）．UV 解離の可能性

励起は起り難く赤外活性の分子に比べれば比較的温度は高も考えられるが，電子励起された後の２次過程であり，電くなるが 100 Pa の DC 放電でも 0.26 eV 程度である［６］）．一子衝突断面積に比べれば低いと思われるので無視して差し方，０℃，１気圧における Ar 原子の気相中における衝突頻支えないであろう．２次過程でも，表面吸着物質の UV

度は分子運動論から 1 GHz と計算される．したがって，低励起や化学反応，電子・イオン衝突刺激脱離等は無視でき

圧プラズマでは気相中における分子の温度が高くなることない．総じて，低圧プラズマの気相中における一次解離過はない．以上のことから，気相中における分子の一次解離

程は，主に電子衝突によって引き起こされるので，電子付過程としては電子付着解離，電子励起解離，イオン化解離

着解離，電子励起解離，イオン化解離が中心に起こっていを考えればよいことがわかる．
分子軌道法について詳しく知りたい方はよい解説書［７‐ると考えてよい．
結合エネルギーと電子エネルギーで分子の解離を考える

１１］が出ているのでそれを読んでいただくこととして，こ議論を時折みかけるが，結合エネルギーは基底状態の物理

こでは簡単にその特徴だけを述べるにとどめ，分子軌道法

量で熱解離の対象でしかなく，振動・回転励起の選択則はをプラズマプロセスに応用した例を参考にしてプロセス技±１準位であり，一回の衝突で受けるエネルギーは高くて術者が利用できるようになることをめざす．ただし，３体

も〜0.4 eV（C-H 伸縮振動エネルギー等）である．低圧プラ問題が解析的に解けない以上，多体問題を扱う分子軌道法

ズマでは衝突頻度が低いのですぐに脱励起してしまい検討

は近似計算であることを最初に断っておきたい．したがっする必要はない．もちろん，非常に圧力・温度とも高い条て，計算結果が絶対的に正しいとは思わないでほしい．ま件下であれば検討すべき内容ではある．分子の温度は小さ

た，入力データが正しくなければ間違った計算結果も出てな分子の余剰エネルギーを伴う解離反応，電子衝突によるくることを知っておいてもらいたい．分子軌道法を用いる

振動・回転励起，イオン衝撃によって加熱された壁面からときには実験データを基に分子の物性や構造をよく理解の熱伝達と考えられる．通常壁面は水冷されているか放電

し，反応過程をひとつのルートに絞らずに多面的に検討し部から離れているのでそれほど高くはならない．熱伝導性てから行ってほしい．これは筆者の反省を込めたメッセー

Plasma Nano Technology Research Center, Nagoya University, Nagoya, AICHI 464-8603, Japanauthor’s e-mail: hayashi@plasma.engg.nagoya-u.ac.jp

８３５!2013 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.89, No.12 December 2013

れている．その一つが活性化エネルギーを低く評価してしジでもある．

まうことである．この改良法がいくつか報告されている．

２．分子軌道法について

CAM-B3LYP（Coulomb-attenuating method）法［１７］も そ分子軌道法を大きく分け る と 波 動 関 数 を 基 本 としたのひとつである．遠距離効果を加えたことで B3LYP 法の

Hartree Fock（HF）法に基づく計算方法と電子密度で表わ

活性化エネルギーを低く評価してしまうことを補っていされるポテンシャル汎関数を基本と し た 密 度 汎 関数法

る．B3LYP 法は計算時間が大幅に節約でき，これまでの経（Density Function Theory; DFT）がある．HF 法は１電子験からすると，電子付着した負イオンのエネルギー（電子ハミルトニアンを用い，分子軌道は原子軌道の線型結合

付着エネルギー）やイオン化エネルギーを６‐311G（d, p）基（Linear Combination of Atomic Orbital; LCAO）で表せると底関数系で正確に計算できている．最低三重項励起エネル

して，各原子軌道の係数を未定定数法で，エネルギーと係

ギーは B3LYP/6-311＋G（d, p）
で求めることができる．振動数が自己無撞着になるまで繰り返し解く方法（Self Consis-

エネルギーも少し大きく計算されるが，補正係数を掛ける

tent Field method; SCF 法）である．このようにして解いたことでほぼ正確に予測できる．密度汎関数法では相関汎関

解は近似解であるから真値より高い値を示す．もう一つの数を含むので電子相関を一部取り入れてはいるが HF 法と特徴は原子間距離が離れた場合，電子相関を考慮していな

同様に一重項の反応ポテンシャルを十分に表現できてはいいので各原子のイオン性を強く評価してしまうところにあない．HF 法と同様に解離原子・分子の総エネルギーを個る．特に一重項状態にある分子の解離エネルギーを求める

別に計算して評価することが必要になる．

時に，原子間距離を関数として計算すると解離分子（原子）こ れ ら の 計 算 を 行 う 主 な 汎 用 プ ロ グ ラ ム に は，の総エネルギー値の和よりも高いエネルギー値を示してしGAMESS［１８］，Gaussian［１９］，MOLPRO［２０］な ど が あ

まう．解離エネルギーを求める時にはそれぞれの解離分子

る．GAMESS は無料で頒布されているプログラムで登録（原子）の総エネルギーを求めてその和によって評価する

すればだれでも使用できる．基本的には Linux ベースで使ことが必要になる．ラジカルや三重項状態のエネルギーを用する．Gaussian と MOLPRO は購入する必要がある．一求めるときにはこのようなことはなく，比較的正確に表現

番使いやすいのは Windows 上で使用できるようにもなっできる．このとき基底関数重なり誤差（Basis Set Superpo-ている Gaussian であろう．しかし，Windows 上で使うこと

sition Error; BSSE）が発生するので counterpoise 法によるのできる汎用の３２ビット版では計算上の制約が多く，大き補正［１２］が必要になるが，反応性をみるだけならば和にな分子の計算には向いていない．特に励起状態の計算ではよって表現しても大きな間違いにはならない．電子相関を

メモリー制限のため，c-C4F8 や c-C5F8 の CCSD の計算がで考慮した計算方法として MP2（M!ller-Plesset の２次摂動きなかった．大きな分子の励起状態を計算するのには６４法）や CCSD（Coupled Cluster Single and Double 法），そビ ッ ト 版 の Gaussian09（３２ビ ッ ト 版 の 価 格 の３倍 は すの他の方法がある．

る）やGAMESSかMOLPROが向いている．特にMOLPROHF 法に基づいた方法を用いてより正確に分子の諸量をには大幅に収束条件を速めるアルゴリズムが取り入れられ求めるためにはいろいろな改良を加える必要があり，まており，励起状態を中心に計算したい人には向いていよた，改良された方法を用いて計算すると時間がかかりすぎう．

て実用的ではなくなる．この観点から見直されて普及した波動関数は Slater 型の原子軌道の線形結合として用いらのが密度汎関数（DFT）法である．この方法では波動関数れるのが基本であるが，主量子数によって波動関数の形がの代わりに電子密度で表されるポテンシャル汎関数を用い異なっており，２電子積分項を計算するのに多くの時間をており，多電子系の電子状態計算時間を大幅に短縮してい要する．Slater 型の原子軌道を Gauss 関数で表すと解析的

る．密度汎関数法は古くから議論されているが，１９６５年にに計算できるので，この２電子積分項を計算する時間が大Kohn と Sham［１３］が電子密度の汎関数である電子間相互幅に短縮される．Slater 型の原子軌道は一つの Gauss 関数

作用ポテンシャルについて，変分原理に従って最低エネル

で表すことができないので，複数のGauss関数で表現する．

ギーとそれに対応する分子軌道を求める方法を提唱したこ

３つの Gauss 関数で Slater 型の原子軌道を表現するのをとにより，化学物性や固体物性の定量的な計算が可能にな

STO-3Gと呼ぶ．構造を求めるときによくB3LYP/6-31＋G（d, p）

り見直された．その後，改良が加えられ B3LYP が提案され法が用いられる．この基底関数系は，周期律表の２列目のると一 気 に 化 学物 性 計 算 に 利 用 さ れ る よ う に な っ た．

原子（Li-Ne）に対しては内殻の１s 軌道の原子軌道を６つのB3LYP は１９９３年に Becke によって提案された混成汎関数

primitive Gauss（原始）関数で記述し，価電子の原子軌道に［１４］で，Hartree-Fock の交換積分項，LYP（Lee, Yang and

対しては４つ（３＋１）のprimitive Gauss関数で記述すること

Parr）の相関汎関数［１５］，交換汎関数を含む項に３つのパを示す．＋は defuse 関数（複数の広がった原子軌道で構成

ラ メ ー タ を 加 え て，小 さ な 原 子・分 子 の G2 値
［１６］

される関数で Rydberg orbital を含む）を加えたことを，（MP2-G2 によって求めた値）に合うようにパラメータセッ（d, p）
は結合した時原子軌道が分極するのを表現するため

トを決めたハイブリット密度汎関数法である．HF 法を基に，H 以外の原子に d 軌道を加え，H には p 軌道を加えたにして摂動や電子相関を取り入れた計算方法で分子の構造ことを意味する．後に述べる aug-cc-pVDZ 基底関数系［２１］を求めるとかなりの時間を要するが B3LYP 法では数分程

では，correlation consistent polarized Valence Double Zeta

度で分子構造が求まる．しかし，幾つかの問題点も指摘さに Diffuse 関数を加えたことを意味する augmented の略語
８３６

Commentary

An Approach of Molecular Orbital Calculations to Process Plasmas

T. Hayashi

を先頭に付けている．この関数系は１電子および２電子励起を含めた配置間相互作用で得られるエネルギーを最適化するものとしてよく用いられる．

３．気相化学反応の取り扱い
低圧プラズマでは気相反応を無視できると述べたが，chemical dry etching の圧力領域は 100 Pa（原子・分子の平均自由行程が数 μm）以上の粘性流領域で行われるので気相中の化学反応を無視できない．その上，通常アフターグ
ロー領域からウエハまでの距離が 10 cm 以上と長いのでその間でも化学反応が起こる．
N2 マイクロ波プラズマのアフターグロー領域にNF3 を導入すると Si のエッチングが起ることが偶然に発見された［２２］．本来の実験目的は N2+H2 のアフターグローに NF3を導入して形成される HF+NH3 による SiO2 のダメージレ
スエッチングであった．しかし，SiO2 がエッチングされず
本来エッチングされないはずの Si がエッチングされた．H2の導入を忘れたためである．このように偶然に発見され，

図１

NF3→NF2+F（C3v）に沿ったポテンシャル曲線．点線は frozen core ポテンシャル，実線は非断熱ポテンシャルを示す．

図２

GaussView で Calculation Setup 画面を開いて基底状態の構造を最適化するための optimization を選択した時の画
面．

Si のダメージレスエッチングの可能性が見出された．しかし，Si がエッチングされるので F が生成されることは予想されるがどのような機構を経て F が生成されるのかわからなかった．反応機構がわからないプロセスを使うのはあまり気持ちのよいものではない．プラズマを用いる半導体プ
ロセスの多くはその発生機構がわからないで使用されているのでこの機構がわからなくても実用化には支障はないであろう．しかし，わかればプロセスパラメータの最適化に有利であることは間違いがない．幸いこの機構は気相反応で説明できそうであり分子軌道法を用いて解明することにした．
圧力領域は少し低く（〜100 mTorr（13.3 Pa））なるが N2
アフターグロー中の活性種の測定は堀川等［２３］によって報3! ＋
告されている．それによると，主な活性種である N（A
，
u ）
2

N
（4So），N（2Do）の内，後者の２種の状態にある N 原子は圧力に大きく依存しない．N（4So）は基底状態にある原子なので反応しない限り消滅することはないのでわかるが，励起状態の N（2Do）も下流領域まで大きく減衰しないで残りそうであることが報告されている．N2 ダウンフロープラズマ

で，一番左の F 原子を示す項目を選び，NH3 の H をクリッ

中に NF3 を導入するとこれらの物質と反応するに違いない

クするとFに置き換わる．こうしてNF3 ができたところで，と思われた．図１は NF2-F の距離を関数としたときの基底

Gauss View に戻り Calculate をクリックする．一番上の状態，最低励起三重項状態，イオン化状態のポテンシャル

Gaussian Calculation Setup をクリックして計算の方法を示曲線である［２４］．点線で示したのは基底状態（C３v）の構造

す画面を出す．一番左のJob typeをクリックしてOptimiza-を凍結して N-F の距離のみを延ばして（C３v→CS の対称性

tion を押して，Optimization to Minimum になっていることを保持して）求めた frozen core ポテンシャルである．これを確認し，Calculation Force Constants Never になっているを求めるための計算は，まず Gaussian プログラムの GUIところを Once に変える（図３）．B3LYP や MP2 では構造である Gauss View で file の項目を開き，New→create Mo-

最適化の過程で Force Constant を計算することができるよ

lecular Group で分子構造作成画面を開き，file の下にある

うになっているのでこれを用いると容易に構造最適化がで

6

C をクリックする（図２）．周期表が現れるので N をク

きる．この操作が終わったら，再度 Gauss View に戻り，リックすると N 化合物の構造が幾つか現れる．ここで，Method をクリックして計算方法を選ぶ．最初に基底状態

NH3 の構造をクリックすると予備画面に NH3 が表れるのでの構造を決めなければならないので，Ground State を選

この状態で分子構造作成画面をクリックすると分子構造作

び，その右側の Custom のところをクリックして DFT を選

6

成画面に NH3 が作成される．次に，もう一度， C をクリッ

ぶと Gauss View 5 では B3LYP が表示される．次に Basis

クして周期律表上の F を押すと F の化合物が示されるのSet に移り６‐３１あるいは６‐３１１を選ぶ．右側に新しい□が現
８３７

Journal of Plasma and Fusion Research Vol.89, No.12 December 2013

れるのでそこで Diffuse 関数を選択する．負イオンや正イ計算する場合にはファイル名を残したほうがよい．後にオンのエネルギーやポテンシャルを計算するときにはこの

orbital の形をみることもできる．単に構造を求めるだけな□内を空白にしておく．６‐３１あるいは６‐３１１の Diffuse 関数

らばこの項目で don't save を選択する．％mem は計算に必として用意されているのは＋，＋＋で示される関数形であ

要なメモリー領域を確保する項目である．B3LYP で構造

る．価電子を形成する原子軌道より大きな関数系で構成さ決定するだけならば default 値の 6 MW で十分である．励起れており，Rydberg orbital を含んでいる．Diffuse 関数を選状態を計算するときには 100 MW 以上必要になる場合が多択する□の右側の□をクリックすると分極関数が示され

いが３２ビット，４コア対応の Gaussian09 でも最大 200 MW

る．左側は水素以外の原子に対する分極関数で結合によりしか受け付けない（６４ビット版では 200 GW まで対応する

歪んだ効果をdおよびf関数で補正するものである．高精度ようである）．３２ビット，２コア対応の Gaussian09 では最の計算をしないのであれば，周期表の２列の原子を計算す

大 120 MW 程度である．右の general の画面では計算結果

るのには通常d関数だけで十分である（図４参照）．右側にを用いて次の計算を行う場合には connectivity にチェック

H に対する分極関数が示される．結合することにより Hを入れる．特にその必要がなければチェックを外す．構造の１s 軌道が歪み，p 関数の成分をもつことを補正するため決定だけならば Guess 画面以降は通常必要ないので何もしである．Gaussian Calculation Setup 画面の一番下に charge

なくてよい．以上のチェックや記入を終了したら下のタスと spin の□が表示されているので，charge（中性であればクバーに示された左端の submit を押す．Save を押して０，負イオンであれば−１，正イオンであれば１）を入力

ファイル名を記入すると計算が始まる．これで最も安定なする．この画面の入力が終わったら title 画面に移動し，わ構造が求まる．問題とする結合距離を除く他の原子の座標

かりやすくタイトルを記入する．ここには注釈を入れるのを凍結して（Frozen core の）ポテンシャル曲線を求めるにで何を記入しても計算には関係ない．次に link 画面に移動

は，上記で求まった output file を出して Gauss View 上に表し，％chk＝にファイル名を付けて残す場合には default

示されている座標を示す A をクリックする．ここには Z-を選択する．ここには計算された詳細なデータが残される

matrix で表示された座標が示されている．ここの C-F の距ので多くのハードディスク領域が必要になる．励起状態を離が例えば 1.38 と示されているとして，C-F を 1.18〜2.18の frozen core ポテンシャル曲線を求めるならば，calculation→job type→scan と進み，この画面で rigid を選択する．
Method の画面で一重項基底状態のポテンシャルを計算するのであれば，０と singlet を選択し，次の画面で構造決定の時と同様にファイル名とメモリー容量を記入して Edit
を押す．この時，save 画面が現れるがそのままクリックしないで，左下のCartesian coordinateのところにチェックが入っているので（GaussView０３にはない），それを外し，名前を付けて save する（図５）．これを外さないと Cartesian
coordinate で原子の座標が表示されてしまい，scan で Fの座標を移動しながらポテンシャル曲線を求める計算ができない．Cartesian coordinate を使って一点ずつ計算するの図３であればそれでよい が，C-F の 距 離 を 変 え て 行 う に はGaussView で Calculation Setup 画面を開いて基底状態の構造を最適化するための optimization を選択した時の画面．

図４

Z-matrix 法を使った方が便利である．ファイルの save が終基底関数系の選択画面．

図５

８３８

rigid scan を選択し frozen core ポテンシャルを計算する時に現れる save file 画面．

Commentary

An Approach of Molecular Orbital Calculations to Process Plasmas

T. Hayashi

わると Z-matrix で表示された座標がテキスト画面で現れンシャル曲線が図１である．負イオンと正イオンのポテン

るので変えたい C-F の距離が表示されている数字の右横にシャル曲線は B3LYP/6-311G（d）で，最低三重項励起状態

半角の一字分開けて，scan したい数字を，さらに一字分開のポテンシャル曲線は６‐３１＋G（d）で計算している．同一図

けてどの間隔で計算したいかの数字 を 実 数 で 記 入するに示してあるが，basis function が異なると計算される総エ（図６）．例えば，C-F の距離を 1.18〜2.18 を 0.1Å刻みに１０

ネ ル ギ ー 値 が 異 な る の で，三 重 項 の エ ネ ル ギ ー 値 は点計算したいときには C-F 1.38000000 のところを（テキス

６‐３１＋G（d）を用いて別に計算した基底状態のエネルギー

ト画面では原子の指定順に座標が示され，C1-F2 の距離は値との差を用いて求め重ね合わせてある．Frozen core ポ

B1として表示されるので）B1 1.18000000 10 0.1と記入する．

テンシャルの曲線は基底状態の最低エネルギー付近を除い上書きした後に，この画面で submit を押せば計算が始まて大きな意味をもたない．電子付着，励起，イオン化等がる．このとき，最後の数字が整数ではなく実数でないと er-

起こった後は構造が変わると考えられるからである．このror の表示が出力され動かないので注意する必要がある

構造変化をみるのが非断熱ポテンシャル曲線である．この（例えば，角度を１０°毎に変えてエネルギーを求める時には曲線が frozen core ポテンシャルの下（より低エネルギー）

10.0 とする）．距離を縮めるときには−0.1 と記入すればよに位置し，解離性を示せばその励起により解離すると考え

い．励起や電子付着，イオン化プロセスはフェムト秒の現られる．図１を見ると負イオンの frozen core ポテンシャル

象であるが電子励起等が起こった後は核間距離が安定な方とそれらの非断熱ポテンシャルの間には大きな差があり，向に動くので非断熱ポテンシャルを求める必要がある．こ明らかに解離することを示している．この計算で，Errorの計算を行うには，Gauss View の最初の画面で Job typetermination request processed by link 9999 と表示されるこを開いて scan の項目を選択し，その下の項目として Re-とがある．これは optimization step の default 値が２０回にlaxed（Z-matrix）を選択する．他は前述した scan と同様になっており，それを超えても最適構造が求まらなかったと設定し，最後に一番下に示されている Edit 画面に入る．

きに出力される．このような場合には option 項目欄に２０以Save の画面が出るので scan の時と同様に Cartesian coordi-

上の step 数を指定するか，上記で求まらなかった時の out-nate のチェックを外し，ファイル名を記入して save する．

put file を開いて，再度 optimization を行うと求まる．また，同様に，Z-matrix で表示された座標系がテキスト画面とし

calculation 画 面 の general の 項 目 で ignore symmetry にて現れるので，動かしたい結合距離のところにカーソルをチェックが入っていないと最適構造が求まらない場合があ

移動し，B1 1.3800000 S 10 0.1 のように最初に scan の S をこる．正負のイオン状態や励起状態では基底状態と異なったこに記入する．Job type のところの scan で指定した Re-構造になる場合である．Ignore symmetry にチェックを付

laxed（Z-matrix）はここでは単に Z-matrix で最適構造を求けないと対称性を保持して最適構造を求めることになるのめるという指定でしかない．S を記入することで指定したで，異なった対称性の構造を求めようとしない．NF3 の負結合距離を動かす（scan する）と同時に各点における構造

イオン化状態も最低励起三重項状態もNH3 の励起状態と同最適化を行いながら計算する構造になっている．一点だけ様に傘型構造をした C3v から平面構造のD3h になり，これかの計算をするのであれば，この S の代わりに F（Freeze

ら C2v の対称性に沿ったルートで NF2 と F に解離すると思の意）のみを記入するとその指定した結合距離を凍結してわれる．このような場合，ignore symmetry にチェックを

その他の原子位置の最適座標を探してくれる．F を使って入れないで，D3h 構造で非断熱ポテンシャルを求めると一必要とする構造を探す手法は励起状態のポテンシャルを計

つの F が離れていくポテンシャル曲線が得られる．このよ算するときに必要になる．最後に上書き保存を押して閉じ

うに基底状態と同じ対称性で計算して error termination

ると計算が始まる．

になる場合の多くは，異なった対称性になっている可能性

このようにして計算した NF3 の NF2-F を関数としたポテが高い．しかし，図１で見てわかるように，非断熱ポテン
シャルを計算しなくても，frozen core ポテンシャルを見るだけで解離するかしないかの判断はできる．
図１から NF3 の最低励起三重項状態の励起エネルギーが3! ＋N（Aの エ ネ ル ギ ー（約 7 eV）に 近 い こ と が わ か
u ）
2
3! ＋
が消滅し
る．これから，アフターグロー領域に N（Au ）2ないで流れてくれば共鳴的にエネルギー移動が起こり，
NF3→NF2+F の解離が起こって F が生成されると考えられ
3! ＋
る．しかし，N（Aは励起状態であるから衝突により
u ）2

エネルギーを失って消滅し，生き残って流れてくる確率は低いと考えられる．そこで，N（4So）および N
（2Do）
との反応
を考えることになる．
図７はこの原子種と NF3 と の 反 応 ポ テ ン シ ャ ル で ある．左端は基底状態の N（4So）と基底状態の NF3 の総エネル
ギーの和を基準とした時の N（4So）および N（2Do）と NF3 を
図６

frozen core ポテンシャルを求める時の設定．

別々に計算した時の総エネルギーの差である．系全体を二
８３９

Journal of Plasma and Fusion Research Vol.89, No.12 December 2013

４．電子衝突解離
４．
１ 電子付着状態とイオン化状態の計算
前述したように，低圧プラズマ中の気相中における一次
解離過程は電子付着解離，励起解離，イオン化解離が主である．電子付着解離とイオン化解離は負イオンと正イオンのポテンシャル曲線を求めればわかる．これらの量は B3
LYP/6-311G（d，p）で計算できる．Diffuse function である＋を付加すると実験値を再現しない．多分これは電子付着やイオン化が価電子の virtual orbital への付加と価電子のイオン化であるからと思われる．HF 法を基本とした計算法で求められた軌道エネルギーは Koopmans の定理［２５］（軌
道エネルギーが電子の束縛エネルギーに等しいとする定理
図７で符号を反対にすればイオン化エネルギーになる）をほぼN 原子と NF3 の反応ポテンシャル．

満足するが密度汎関数では分子の総エネルギーと軌道エネ
重項，四重項にして計算すると図に示されるような反応ポ

ルギーが別々に計算されるので，この定理を満足しないもテンシャルが得られる．B3LYP では活性化エネルギーがの と 思 わ れ る．図９は CF4 の 基 底 状 態，三 重 項 励 起 状

4 o低く計算される傾向にあることを考えれば，N（ S ）との反態，電子付着した負イオン化状態，正イオン化状態のポテ

応は考えなくてもよいであろう．N（2Do）と NF3 が反応する

ンシャル曲線である．最低三重項励起状態を除き，B3LYPと不安定な N=NF3 が形成され，N=NF2+F に解離する．これから，F が生成される反応は N（2Do）との反応が主であると推定される．
図１と図７から N2 ダウンフロープラズマ中に NF3 を導入して F が生成される反応経路は→NF2+FNF3+e（6〜7 eV）
3! ＋
1! ＋
3
N（A
u ）＋NF3→N（Xg ）＋NF（
2
2
3 E1）
1! ＋
→N（X
g ）＋NF2+F22oN（ D ）＋NF3→N=NF3→N=NF2+Fの３つの経路であることがわかった．NF3 の導入部にプラズマが伸びていないことは明らかなので電子付着は最初に
3! ＋
除くことができる．N（2Do）か N（Au ）のどちらかがア
2

フターグローまで生き延びている可能性が考えられるので

図８

N2 アフターグロー中にNF3 を導入した時に起きると予想される反応図．

図９

CF4 のイオン化状態，最低三重項状態，負イオン状態，基底状態の C-F を関数とするポテンシャル曲線．

この内どちらかである．分子よりは原子のほうが小さく平均自由行程が大きいので N（2Do）と反応している可能性が高いと考えられる．これらの反応種は圧力が高いと衝突により短時間で消滅してしまう．したがって，エッチング速
度は圧力と流量および NF3 を導入する位置に大きく依存す
ることがわかる．この傾向は実験事実とよく一致する．上記反応が起こった後にさらにどのような反応が起こるかを同様な方法で計算してまとめたのが図８である．上記反応が起こって NF2 や N=NF2 が生成されると，不安定なこれ
らの化合物は容易に基底状態の N（4So）と反応する．これらの反応を計算するとき，時々 angle alpha exceeds 0-180 のerror 表示が出て計算が止まってしまうことがある．このときは Gauss View 画面の A をクリックし Z-matrix で表示されている原子の属性を変更すると計算できるようにな
る．Z-matrix ではどの原子がどの原子に接続しその角度が何度になっているかで表示する．接続している原子を変えると角度が変わり１８０°を超えてしまう角度を１８０°以内にでき再度計算が続行できる．

８４０

Commentary

An Approach of Molecular Orbital Calculations to Process Plasmas

T. Hayashi/6-311G（d）で 計 算 し て い る．最 低 三 重 項 励 起 状 態 はた VUV スペクトルが実測値と非常によい一致を示すからRydberg state なので 6-31＋G（d）で計算した．実験値と非である．

常によく一致していることがわかる．他の化合物でもこの励起状態を計算する場合，先ず前述した方法で基底状態方法で計算した値は実験値をよく再現する．frozen coreの構造を求めるが，求めた構造が本来の対称性からずれてポテンシャルを計算したのち，非断熱ポテンシャルを計算いる場合がある．例えば，CF4 の場合，B3LYP を用いた構することで解離するかしないかがわかる．

造最適化で Td の対称性が得られたとする．この時，一つの４．
２ 励起状態の計算

結合距離が小数点以下７桁以降で他の結合距離と違ってい最低三重項励起状態を除く一重項，三重項励起状態の計たとする．B3LYP では小数点以下６桁で収束条件を満た

算は大変難しく，多くの時間を要する．一般には密度汎関

すと Td と表示して計算を終了する．しかし，この座標系を

数 TD-DFT（Time Dependent Density Function Theory）

使用して EOMCCSD/aug-cc-pVDZ の計算をすると大変が多く用いられるが筆者はあまり使用していない．HF 法（１０倍以上）時間がかかる．場合によっては，計算途中で表を基本とした EOMCCSD/aug-cc-pVDZ を主に使用してい示される軌道対称性の表示が Td の点群に含まれないで C３v

る．後者の方が，時間はかかるが，実測値をよく再現するや他の点群の表示になることがある．この時は Gaussian 09

ことがその主な理由である．しかし，大きな分子には計算の画面で表示されている Process のところをクリックし時間がかかり過ぎるのと前述した制約があるので適用困難

Kill job を実行して終わらせる．間違っても，Gaussian 09

である．この場合には，精度は低いが CIS（Configuration

ウインドウ画面で表示される右上の画面消去ボタン X で止Interaction Single）法 を 用 い て い る．計 算 時 間 は CCSD

めてはならない．X で消したつもりでも，計算の実行を終を用いるよりも TD-DFT の方がはるかに少なくて済む．まわらせたことにならないため，次の計算ができなくなってた，実測された VUV スペクトルをよく再現する場合があしまう．Kill Job で終了させたら，もう一度構造決定のプロ

る．図１
０［２６］は c-C4F8 の TD-DFT で計算された VUV スペ

グラムを起動して，A の座標系を表示させ，A の左の対称操

クトルと実測された VUV スペクトルである．TD-DFT の作のボタンを押す．A の座標系で対称になっている結合距

計算値は一般には実測値よりも低いエネルギー値を示す．

離が小数点以下８桁まで同じになっていることを確認し筆者が TD-DFT を使用しない理由は一重項状態に対称性て，再度構造最適化を行う．結果を View File か Stream

禁制で振動子強度のない低い励起状態が現れることと三重

Output File で求める対称性になっているかを確認すると項状態を計算すると多数の低い励起状態が現れることにあともに A で座標を表示させ小数点以下８桁まで対称になっる．励起状態がどの virtual orbital への遷移で構成されていているかを確認する．一度の再計算で求める対称性の座標

るかは未知なことが多くこれは好みの問題といえる．

値が正確に得られなかったらもう一度同じ操作を行う．通図１１［２７］は CF4 のいくつかの計算方法と基底関数系を用

常，１〜２回で求める対称性をもつ座標系が得られる．特いて，一重項励起状態を計算して得た VUV スペクトルでに一つの結合距離を動かしてポテンシャル曲線を求める場ある．計算する励起状態の数を２０と指定して計算した．励合，前述したように F をテキスト画面で表示される Z-

起状態の計算では分子の対称性が非常に重要になる．そのmatrix の座標表示のところに付けて，計算しようとする結ため，励起状態の計算を始める前に基底状態の構造を精度

合距離を固定し，残りの原子の位置を最適化する．このとよく計算しなければならない．基底状態の構造は B3LYP

き，一度の構造最適化では求める値が得られない場合が多か MP２で通常計算される．筆者は時々 CAM-B3LYP をい．例えば，CF4 で Td→C３v で解離する時の励起状態を計使 っ て 計 算 し て い る．こ の 方 法 で 分 子 構 造 を 求 め，

算する場合，一つの C-F の結合距離を0.1 Åずつ離しながら

EOMCCSD/aug-cc-pVDZ で励起状態を計算すると得られ構造決定をする．F を用いて一つの結合距離を固定し，他の原子の座標系を最適化するが，一度の構造最適化では
C３v の対称性をもつ座標系は得られない．同じ操作を２〜

図１
０ c-C4F8 の実 測 VUV ス ペ ク ト ル と TD-DFT で 求 め た VUVスペクトル．

図１
１ CF4 の VUV スペクトルと各計算法による VUV スペクトル．

８４１

Journal of Plasma and Fusion Research Vol.89, No.12 December 2013のエネルギーなので，基底状態のエネルギーを計算した後３度行って C３v の対称性が得られる．
励起状態の計算の手順は，Gauss View で構造最適化をに，これに計算された励起エネルギーを加える形でポテン行った時の output file を開き，Calculate をクリックしてシャル曲線が得られる．B3LYP で構造最適化した時の全Gaussian Calculation Setup 画 面 を 開 く．Job Type を En-

エネルギー（total energy）を含め，分子軌道計算ではエネ

ergy に し て，Method の 画 面 に 入 り，Method で 例 え ば

ルギー，結合距離等において原子単位（atomic unit）が用

EOMCCSD を選択する．次に，Basis Set で cc-pVDZ を選いられている．励起状態の計算結果は eV 単位で表示され択すると左に空白の□が現れる．ここに Diffuse 関数を加るが，総エネルギーは hartree 単位（1 hartree＝27.211 eV）

える aug‐を表示させる．一重項の計算をするのであればそで示されているので，基底状態の最小値を０とした時の各の 下 の□で default か singlet only を 選 択 す る．Solve for

点でのエネルギー差を hartree 単位から eV 単位に換算し，More States に励起状態をいくつまで計算させたいかの数

基底状態のポテンシャル曲線を求めてから，各点における字を記入する．Default 値は６になっている．筆者は通常イ励起エネルギーを加えて励起状態のポテンシャル曲線を求オン化閾値付近まで計算するようにしている（N＝２０程めることになる．CH2F2 の C-H あるいは C-F の結合距離を度）．この数字はどの励起状態までを計算させるかで決め関数としたポテンシャル曲線を求めるときの例を図１３に示ればよいが，この数が多くなれば高い励起状態に対する計す．C２v→Cs の対称性に沿って解離する場合には，この点群

算値の精度は上がるものの，計算時間はその分長くなる．

に含まれる対称性が，C-H の場合は A1→A'，A2→A"，B1

次にタイトル画面に入り，任意の文字列で計算内容をわか→A'，B2→A"に，C-F の場合は A1→A'，A2→A"，B1→A"，

りやすく記入する．注釈欄なので何を書いても計算には関B2→A'に変わる．前述したように，各点における構造最適

係ない．Link 0 の欄に入り，Memory Limit を３２ビット 4化の時に対称性を十分確認しておかないとポテンシャル曲

core コンピュータ用のプログラムであれば 200 MW まで設線を求めることが大変困難になる．その上に，前述したよ定できる．３２ビット 2 core 用のプログラムでは 120 MW

うに，CCSD 法では求めた基底状態の対称性が低いと１０倍以上にすると計算しないで終了してしまうことがある．

以上の時間を要するので注意が必要になる．解離ポテン

Chkpoint File の欄はできれば save にした方がよい．かなり

シャルはこのモードに沿って結合距離を 0.1 Å変化させな

ハードディスクのメモリー領域を必要とするので計算内容

がら基底状態および励起状態のエネルギーをプロットしてをチェックする必要がなければ Don't save を選択する．こ求めている．

れらの設定が済んだら下の Submit を押すと save 画面が現

図１４は NH3 のポテンシャル曲線である．左図は C３v の傘

れて，file 名を記入し計算をスタートさせることができる．

型構造から Cs の対称性に沿って H2N-H に解離する時の基

図１２は CF4 の励起状態計算に用いた aug-cc-pVDZ 基底関

底状態および一重項励起状態の EOMCCSD/ aug-cc-pVDZ

数系で構成される分子軌道の例である．CF4 の occupied or-を用いて得た曲線，右図は平面構造の D３h における H2N-H

bital の数は２１で，aug-cc-pVDZ を用いた時の virtual（unoccupied）orbital の数は９４である．aug-cc-pVTZ を用いた時の virtual orbital の数はさらに多くなり２０９になる．多くの
virtual orbital を用いればそれだけ計算精度は高くなるが１０倍以上の時間を要する．現時点で最速の CPU（i７‐３７７０）を使用しても，CF4 を EOMCCSD/ aug-cc-pVTZ で計算すると１６時間を超える．この方法で多数点を計算し励起状態のポテンシャル曲線を求めるのは現実的ではない．pVDZと pVTZ を用いた時のエネルギー値は僅かに異なるが低い励起状態のエネルギー値は大きく異ならないので，pVDZ
図１２ aug-cc-pVDZ を用いた時の occupied orbital と Virtual orbital の準位と数（上に行く程エネルギーが低くなる）．

を用いることとした．
ポテンシャル曲線は，基底状態の振動構造解析から，
Td→C3v の対称性で解離するルートと Td→C2v で解離するルートが考えられる．このルートに沿って C-F の距離を0.1 Åずつ変えながら，以前に述べたFによる構造最適化法
を用いて，各結合距離で得られた最適化構造の座標を用いて励起状態の計算を行うことでポテンシャル曲線が得られる．Td→C3v では T1→A2+E，T2→A1+E に，Td→C2v では
T1→A2+B1+B2，T2→A1+B1+B2 に分かれるので，この対称性に従い各点の励起エネルギーをプロットして求める．基底状態のエネルギーは output file の最後に HF，MP2，MP3の値として出力される．基底状態のエネルギーとしては

図１
３ 基底状態の各振動モードおよびそれに沿って解離ポテン
シャルを求める時の対称性の変化．

MP2 か MP3 の値が通常選ばれる．筆者は MP2 の値を用いている．各点における励起エネルギーはこの基底状態から
８４２

Commentary

An Approach of Molecular Orbital Calculations to Process Plasmas

T. Hayashi

影響を与える C3F＋5 イオンが大きなピークとして質量分析
計で測定されている［３１］．このイオンは c-C4F8 のイオン化
閾値付近で F の転移と開裂反応が同時に起ることにより
CF2CFCF+2 +CF3 になることが分子軌道計算で明らかになっている［３２］．また，CF4 プラズマの負イオンを計測す
ると F‐イオンが主となるが，この F‐イオンは電子付着解離
により CF3 ラジカルとともに生成されることが分子軌道法
を用いると容易にわかる（図９）．文中で述べたように N2アフターグロー中に NF3 を導入すると F の生成が起きることも説明できる．このように計測だけではわからないことが分子軌道法を用いることにより明らかになる．分子（原図１
４ NH3 のポテンシャル曲線と最低励起一重項の解離．

子）およびそのイオンの計測でわかることは原料物質のプ
ラズマ分解種が何であるかであり，原料物質と分解種の間に解離する時の基底状態および一重項励起状態のを繋ぐ情報を提供するのが分子軌道法の１つの役割である

EOMCCSD/aug-cc-pVDZ を用いて得た曲線である．基底ともいえる．

状態はC３v なのでこの構造の第１励起状態へ励起が起こる．

この分子軌道法について，簡単に Windows 上で使用で励起されると平面構造の D３h になり，最低一重項励起状態

きる Gaussian09 を例にして，未経験者でも使用できるようは解離性に変わるので NH2+H に解離する．さらに，N-Hに述べてきたつもりであるが，計算方法や基底関数系も多が２Åを超えるとこの励起状態は基底状態と入れ替わる

義にわたっており，用いる際には参考文献に挙げた入門書［２８，２９］．実験では 6.3 eV 付近に吸収極大をもつ振動構造を是非よく読んで間違った使い方をしないように注意してが得られている［２７］が励起されると解離する．計算で求め戴きたい．また，ここで述べた使い方は非常に限られた範たポテンシャル曲線はこの実験結果をよく表現している．

囲のものであり，Gaussian09 や GAMESS，MOLPRO 等に励起状態の計算例をいくつか述べてきたが，励起状態の付随している取扱説明書をよく読んで使用することをお薦

非断熱曲線を求めるのは容易ではない．CF4 の場合，C-F

めする．励起状態の計算では群論の知識が必要になるのでの距離を変えながら各励起状態のエネルギーを

分子軌道法の入門書［７］あるいは群論の専門書を読み基礎

EOMCCSD/aug-cc-pVDZ を用いて求めるが，１点の計算知識を身につけてから計算することをお薦めする．

に最速のCPUを使っても１時間２０分の時間を要する．８点

全くの未経験者には確かに敷居が高いかもしれないが基の計算をするには１０時間以上の時間が必要となり，三重項

本的なことを知れば使えないことはなく，使い易い GUIの計算も行うのであれば２０時間以上になる．入力を間違えが装備されているものもあり，パーソナルコンピュータのると更に時間を要する．また，分光学的にどのような構造性能も高くなって使いやすい環境は整いつつある．エッチになっているかがわかっていれば，小さな分子について

ングプロセスに限らず，CVDやその他のプロセスへの適用

は，比較的容易に計算できるが未知の化合物の励起解離曲も可能である．コンピュータの性能やソフトは改良が続け線を求めるのは現実的に無理といってよい．励起状態にあられ今後さらに使い易くなるものと思われる．現状困難でる分子の構造を求める方法が大変難しい上に，デスクトッあるエッチング等の表面反応へ適用できる日も近いであろ

プ用の３２ビット版 Gaussian では多くの制約があるからでう．プロセスプラズマでどのような気相および表面反応がある．この計算にはよく CASSCF（Complete Active Space

起こっているかの解析を進めるのに分子軌道法は大いに役

Multi Configuration SCF）
［３０］が用いられるが制約のないに立つものと考える．多くの人がプロセス中で起こってい大型コンピュータ用のソフトを使うことが必要となる．デる現象の解明に携わることができればプロセスの開発や効スクトップコンピュータで計算するのならば GAMESS や率向上に役に立つことは間違いがない．是非多くの人がこ

MOLPRO などが利用できると思うが CASSCF を使うにはの分野に参入し議論を戦わすことを期待したい．

分子軌道法に精通していないと使用するのは難しい．まず謝辞

は励起状態の構造がわかっている簡単な分子について計算
し，解析が十分行えるようになった後に挑戦してみてはど

草稿をまとめるにあたり終始議論に加わって下さった名うであろう．

古屋大学の石川健治教授，励起状態の計算や GAMESS について御教授いただいた大阪府大の小関史朗教授に感謝い

５．まとめ

たします．

以上，プラズマプロセスへのアプローチという観点から

参考文献

分子軌道法を用いた計算方法について述べてきた．分子軌
道法は原料物質が気相中でどのような解離や化学反応を起［１］Y. Matsumi et al., J. Appl. Phys. 60, 4102 (1986).
［２］D.B. Oh et al., J. Vac. Sci. Technol. B 13, 954 (1995).
［３］T. Nakano et al., Jpn. J. Appl. Phys. 35, L338 (1996).
［４］J. P. Booth et al., J. Phys. Chem. 100, 47 (1996).

こしエッチング種や堆積性の物質がどのように形成されるかの情報を得る手段として使うことができると考えている．例えば，c-C4F8 プラズマではエッチング反応に大きく
８４３

Journal of Plasma and Fusion Research Vol.89, No.12 December 2013［５］P. Fendel et al., Plasma Sources Sci. Technol. 14, 1 (2005).
［６］K.S. Suraj et al., Surf. Coatings Technol. 202, 301 (2007).
［７］原田義也：量子化学 上下巻（裳華房，２００
７）．
［８］平尾公彦監修，武次徹也編：すぐできる量子化学計算
ビギナーズマニュアル（講談社，２００
６）．
［９］常田貴夫：密度汎関数法の基礎（講談社，２
０
１
２）
．
［１
０］堀 憲次監修：Gaussianプログラムによる量子化学計算
マニュアル（丸善出版，２
０
０
９）
．
［１
１］古賀伸明編：実験科学講座１
２「計算化学」（丸善出版，
２
０
０
４）
．
［１
２］S. Simon et al., J. Chem. Phys. 105, 11024 (1996).
［１３］W. Kohn et al., Phys. Rev. 140, A1133 (1965).
［１４］A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
［１５］C. Lee et al., Phys. Rev. B 37, 785 (1988).
［１６］L.A. Curtiss et al., J. Chem. Phys. 94, 7221 (1991).
［１７］T. Yanai et al., Chem. Phys. Lett. 393, 51 (2004).
［１８］www.msg.ameslab.gov/gamess［１９］www.gaussian.com/g̲prod/g09.htm［２０］www.molpro.net［２１］例えば，M. Kállay et al., J. Chem. Phys., 121, 9257 (2004).
［２２］Y. Tajima and S. Takahashi: Japan Patent JP2010165954,
and private communication (2011).
［２３］Y. Horikawa, et al., Jpn. J. Appl. Phys. 49, 026101 (2010).
［２４］T. Hayashi et al., Jpn. J. Appl. Phys. 51, 026505 (2012).
［２５］T. Koopmans, Physica, 1, 104 (1933).
［２６］E. Va"eková, Ph. D Theses of The Open University, 2006.
［２７］R.Y.L. Chim et al., Chem. Phys. Lett. 367, 697 (2003).
［２８］G. Herzberg, Molecular spectra and molecular structure III,
(Van Nostrand Reinhold Co. New York 1966).
［２９］B. Cheng et al., Astrophys. J. 647, 1535 (2006).
［３０］例 え ば，N. Yamamoto et al., Chem. Phys. Lett. 250, 373(1996).を参照．
［３１］H. Hayashi et al., J. Vac. Sci. Technol. A 17, 2557 (1999).
［３２］T. Hayashi et al., Jpn. J. Appl. Phys. 50, 036203 (2011).

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

はやし

とし

林

俊 雄

お

現在の所属；名古屋大学大学院工学研究科
付属プラズマナノ工学研究センター．
略歴；東北大大学院理学研究科博士課程を
修了し，１
９
７
７年に日本真空技術株式会社（現：!アルバック）に入社．質量分析器，MOCVD 装置の開発に従事した後，１９９０年ころからエッチング装置の開発に従事し，NLD エッチング装置を実用化する．２００
６年秋から名古屋大学・プラズマナノ工学研究センター特任教授に赴任し，現在に至る．現在，エッチング中の化学反応をテーマにプラズ
マ中の解離過程や気相中の反応過程を中心に研究を進めている．趣味はゴルフであるが近隣の先生方は真面目で趣味のゴルフができないのが現在の悩みである．新しい音楽についていけないので主に古い音楽を聴いている．バロックやロマン派の音楽が好みであり，研究が思うように進まないときはバロックを聞いて心を癒している．ドライブも嫌いではなく２
週に一度神奈川県の家に車で帰る．最近はどの程度の低燃費で走れるかを考えており，スピードを出したことはない．現在のところ 29 km/L が最高で 30 km/ L は未だ達成したことがない．

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

８４４