JP2013250662A - 形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラム - Google Patents

Abstract

【課題】物質の表面に直接的、間接的に到達する反応種を考慮した形状シミュレーションを高速化することが可能な形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラムを提供する。
【解決手段】実施形態による形状シミュレーション装置は、物質の表面を複数の計算要素に分割する分割部を備える。さらに、前記装置は、各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部を備える。さらに、前記装置は、前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部を備える。
【選択図】図６

Description

本発明の実施形態は、形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラムに関する。

ＣＶＤ（Chemical Vapor Deposition）やＲＩＥ（Reactive Ion Etching）等による物質の表面の加工において、加工形状のシミュレーションは重要な技術となっている。このシミュレーションでは、物質の表面を複数の計算要素（例えば点、線分、多角形等）に分割して、各計算要素に到達する反応種のフラックスや、物質の局所的な表面成長速度を算出することが一般的である。また、物質の表面の加工では、表面に直接的に到達する反応種だけでなく、一旦表面に接触した後さらに別の表面に間接的に到達する反応種が、加工形状に影響を与えることが知られている。しかしながら、これらの反応種を考慮して、フラックスや表面成長速度をすべての表面にて矛盾なく計算するには、長い計算時間が必要となる。理由は、計算時間が、計算要素の個数の２乗のオーダーで増えていくためである。

特開２００９−４９１１０号公報

G. Kokkoris et al., "Simulation of SiO2 and Si feature etching for microelectronics and microelectromechanical systems fabrication: A combined simulator coupling modules of surface etching, local flux calculation, and profile evolution", J. Vac. Sci. Technol. A 22, 1896 (2004). Arpan P. Mahorowala et al., "Etching of polysilicon in inductively coupled Cl2 and HBr discharges. II. Simulation of profile evolution using cellular representation of feature composition and Monte Carlo computation of flux and surface kinetics", J. Vac. Sci. Technol. B 20, 1064 (2002). Shahram Abdollahi-Alibeik et al., "Analytical modeling of silicon etch process in high density plasma", J. Vac. Sci. Technol. A 17, 2485 (1999).

物質の表面に直接的、間接的に到達する反応種を考慮した形状シミュレーションを高速化することが可能な形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラムを提供する。

一の実施形態による形状シミュレーション装置は、物質の表面を複数の計算要素に分割する分割部を備える。さらに、前記装置は、各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部を備える。さらに、前記装置は、前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部を備える。

第１実施形態の形状シミュレーション方法の手順を示したフローチャート図である。 第１実施形態における物質の初期構造の例を示す斜視図である。 レベルセット関数について説明するための模式図である。 図１のステップＳ３の詳細を示したフローチャート図である。 物質表面を複数の計算要素に分割した様子を示した模式図である。 図４のステップＳ１２の詳細を示したフローチャート図である。 ローカル座標系について説明するための図である。 可視判定値について説明するための模式図である。 可視係数について説明するための模式図である。 入射角度θ_inについて説明するための模式図である。 鏡面境界条件について説明するための模式図である。 周期境界条件について説明するための模式図である。 ２次元における計算要素可視判定値について説明するための模式図である。 ３次元における計算要素可視判定値について説明するための模式図である。 図６のステップＳ２２〜Ｓ２７の手順の変形例を示したフローチャート図である。 グローバル座標系について説明するための図である。 比較例における計算時間の例を示したグラフである。 第１実施形態における計算時間の例を示したグラフである。 第１実施形態と比較例の計算時間を比較したグラフである。 第１実施形態と比較例におけるθ分割数と計算誤差との関係を示したグラフである。 第２実施形態の形状シミュレーション装置の構成を示す外観図である。 図２１の制御部の構成を示すブロック図である。

以下、本発明の実施形態を、図面を参照して説明する。これらの図面では、同一または類似の構成要素に同一の符号を付しており、重複する説明は適宜省略する。

（第１実施形態）
図１は、第１実施形態の形状シミュレーション方法の手順を示したフローチャート図である。本実施形態の形状シミュレーション方法は、パーソナルコンピュータやワークステーション等の情報処理装置を使用して行われる。

本実施形態の形状シミュレーション方法ではまず、物質の初期構造を情報処理装置に入力する（ステップＳ１）。図２は、第１実施形態における物質の初期構造の例を示す斜視図である。図２に示す初期構造は、シリコン基板１と、シリコン基板１上に順に形成されたシリコン窒化膜２およびシリコン酸化膜３と、シリコン窒化膜２とシリコン酸化膜３とを貫通する貫通孔４を含んでいる。初期構造の入力方法の例としては様々なフォーマットが考えられるが、本実施形態では、物質の表面の形状を点列で表現し、情報処理装置がこれを読み取る方法を用いている。

次に、入力された初期構造から、初期レベルセット関数を作成する（ステップＳ２）。図３は、レベルセット関数について説明するための模式図である。レベルセット関数ψは、物質の表面からの距離ｄを用いて定義される関数であり、計算領域内のメッシュごとに値を有する。レベルセット関数ψの値は、物質の表面において０と定義される（ψ＝０）。また、物質の外部（真空中）ではψ＞０であり、物質の内部（物質中）ではψ＜０である。初期レベルセット関数を作成する際には、各メッシュ点から最近接となる表面を探し、その距離ｄを計算し、メッシュ点が真空中であればその符号を正とし、物質中であればその符号を負とする。なお、初期レベルセット関数は、ステップＳ２で作成する代わりに、ステップＳ１で入力してもよい。

次に、物質の局所的な表面成長速度Ｆを計算する（ステップＳ３）。ここで、表面の成長とは、表面への堆積だけでなく、表面のエッチングも含むものとする。なお、表面成長速度Ｆの計算は、タイムステップごとに行う必要はない。また、本実施形態では、後述するように、物質の表面におけるフラックス（総フラックス）から表面成長速度Ｆを計算し、表面成長速度Ｆからレベルセット関数を計算するが、代わりにフラックスからレベルセット関数を計算し、表面成長速度Ｆの計算は省略してもよい。

次に、表面成長速度Ｆを用いて、時間Δｔ経過後のレベルセット関数を計算する（ステップＳ４）。時間ｔにおけるレベルセット関数ψ_tは、以下の式（１）から計算できる。

ただし、∇はベクトル微分演算子を表し、｜∇ψ_t｜は∇ψ_tのノルムを表す。時間Δｔ経過後のレベルセット関数は、式（１）を離散化した式に従い、レベルセット関数を時間発展させることで計算可能である。なお、本実施形態では、表面形状を時間発展させる代わりに、ある表面形状における表面成長速度Ｆやフラックスを計算してもよい。これは、後述するステップＳ５を１ステップ目でＹｅｓと判定する場合に相当する。

次に、予め設定したプロセス時間が経過したか否かを判定する（ステップＳ５）。プロセス時間が終了した場合には、物質の最終形状を出力し（ステップＳ６）、計算終了となる。プロセス時間が終了していない場合には、ステップＳ３に戻る。

なお、本実施形態では、形状表現の手法としてレベルセット法を用いているが、レベルセット法以外のセル法やストリング法などの手法を用いてもよい。

（１）ステップＳ３の詳細
次に、図４を参照し、ステップＳ３の詳細について説明する。

図４は、図１のステップＳ３の詳細を示したフローチャート図である。

まず、レベルセット法で表された物質表面を、複数の計算要素に分割する（ステップＳ１１）。図５は、物質表面を複数の計算要素に分割した様子を示した模式図である。図５の例では、物質表面がメッシュごとに分割されている。その結果、１つのメッシュ内の物質表面が、１つの計算要素となっている。ステップＳ１１の処理を行うブロックは、本開示の分割部の例である。

なお、物質表面の分割方法は、メッシュ単位に限られるものではなく、どのような方法を採用してもよい。また、物質表面の分割は、タイムステップごとに行う必要はなく、例えばステップＳ１の直後に行ってもよい。

また、図５に示す計算領域は、２次元領域となっているが、代わりに３次元領域としてもよい。また、図５に示す各計算要素の形状は、線分となっているが、代わりに点や多角形などとしてもよい。

図５は、第１の計算要素ａと、第２の計算要素Ｂを示している。計算要素Ｂに到達する反応種のフラックスを計算する際には、気層から計算要素Ｂに直接的に到達する反応種のフラックスと、気層から任意の計算要素ａを介して計算要素Ｂに間接的に到達する反応種のフラックスの両方を考慮するのが一般的である。前者のフラックスを直接フラックスと呼び、後者のフラックスを間接フラックスと呼ぶ。また、これらの合計を、総フラックスと呼ぶ。なお、反応種の例としては、堆積種やエッチング種などが挙げられる。

計算要素Ｂにおける総フラックスΓ_Bは、以下の式（２）のように、計算要素Ｂにおける直接フラックスΓ_B,directと、任意の計算要素ａからの間接フラックスΓ_aB,indirectの合計との和で表される。

ここで、間接フラックスΓ_aB,indirectは、例えば以下の式（３）で表すことができる。

Ｓ_a(Γ_a)は、計算要素ａで吸収されるフラックスの割合を示す付着確率を表す。Ｓ_a(Γ_a)の値は、計算要素ａにおける総フラックスΓ_aに依存する。また、ν(ａ,Ｂ)は、計算要素ａと計算要素Ｂが互いに見えるか否かを示す可視係数（面間可視係数）を表す。計算要素ａ、Ｂを結ぶ直線が、計算要素ａ、Ｂ間において物質表面とぶつかる場合には、ν＝０となり、ぶつからない場合にはν＝１となる。また、ｇ(ａ,Ｂ)は、計算要素ａと計算要素Ｂとの位置関係（面関係）を示す形態係数を表す。ｇ(ａ,Ｂ)の値は、計算要素ａ、Ｂの互いの見えやすさの程度を示している。ｇ(ａ,Ｂ)の値は、計算要素ａ、Ｂ間の距離や角度などに依存する。

式（２）に式（３）を代入すると、計算要素Ｂにおける総フラックスΓ_Bは、以下の式（４）で表すことができる。

図４のフローでは次に、任意の計算要素における直接フラックスや、任意の計算要素間における可視係数ν、形態係数ｇを計算する（ステップＳ１２）。

次に、各計算要素ｉの直接フラックスΓ_i,directを仮の総フラックスΓ_iとして使用し、各計算要素ｉにおける付着確率Ｓ_i(Γ_i)を計算する（ステップＳ１３）。この際、このフラックスは、中性分子を含んでいても、志向性を持ったイオンを含んでいてもよいし、その両方を含んでいてもよい。

次に、可視係数ν、形態係数ｇ、直接フラックスΓ_i,direct、付着確率Ｓ_i(Γ_i)を用いて、以下の式（５）から、各計算要素ｉにおける総フラックスΓ_iを計算する（ステップＳ１４）。

次に、ステップＳ１３とステップＳ１４の処理を、付着確率Ｓ_i(Γ_i)の値が収束するまで繰り返す（ステップＳ１５）。なお、２回目以降のステップＳ１３では、前回のステップＳ１４で計算した総フラックスΓ_iを、仮の総フラックスΓ_iとして使用する。また、ステップＳ１５では、Ｓ_i(Γ_i)の値が収束したか否かを、Ｓ_i(Γ_i)の変化が閾値以下になったか否かで判定する。そして、Ｓ_i(Γ_i)の値が収束した際の総フラックスΓ_iを、総フラックスΓ_iの正しい計算結果として取り扱う。

なお、計算要素の個数をＮとする場合、任意の計算要素間における可視係数νや形態係数ｇは、まとめてＮ×Ｎ行列で表すことができる。行列の形で表された可視係数ν、形態係数ｇを、それぞれ可視係数行列、形態係数行列と呼ぶ。また、任意の計算要素におけるフラックスは、まとめてＮ行ベクトルで表すことができる。ベクトルの形で表されたフラックスを、フラックスベクトルと呼ぶ。

この場合、式（５）は、以下の式（６）のように行列方程式で表現することができる。

ただし、Ｉ、Ｊは処理対象の計算要素の個数を表し、例えばＩ＝Ｊ＝Ｎである。行列方程式（６）は、どのような解法で解いてもよい。解法の例としては、反復法（ガウスザイデル法、ＳＯＲ法、ヤコビ法、共役勾配法など）や、直接法（ガウスの消去法、ＬＵ分解法、コレスキー分解法など）が挙げられる。行列方程式（６）を解く際、行列Ａが疎行列である場合には、ＣＲＳなどの保存方法を用いた上で疎行列に適したルーチンを使用することで、計算処理の省メモリ化や高速化を図ってもよい。

図４のフローでは次に、総フラックスΓ_iから、各計算要素ｉにおける局所的な表面成長速度Ｆ_iを計算する（ステップＳ１６）。例えば、Ｋ種類の反応種を使用する場合には、表面成長速度Ｆ_iは、Ｋ個の局所的な総フラックスΓ_1,i〜Γ_K,iに依存した以下の式（１０）のような形でモデリングされる。

ただし、ｋは１≦ｋ≦Ｋを満たす任意の実数である。以上のようにして、ステップＳ３の処理が終了する。

（２）ステップＳ１２の詳細
次に、図６を参照し、ステップＳ１２の詳細について説明する。

図６は、図４のステップＳ１２の詳細を示したフローチャート図である。

図６のフローでは、各計算要素に固有のローカル座標系を使用する。図７は、ローカル座標系について説明するための図である。図７(ａ)は、各計算要素の法線ベクトルを示し、図７(ｂ)は、各計算要素におけるローカル座標系を示す。図７(ｂ)に示すように、ローカル座標系の直交座標（ｘ_local,ｙ_local,ｚ_local）は、＋ｚ_local方向が法線ベクトル方向と一致するように定められる。また、ローカル座標系の極座標（ｒ_local,θ_local,φ_local）は、極角θ_localが動径ｒ_localと＋ｚ_local方向との間の角度となり、偏角φ_localが動径ｒ_localと＋ｘ_local方向との間の角度となるように定められる。

計算要素Ｂにおける直接フラックスΓ_B,directは、以下の式（１１）により計算することができる。

ただし、η(θ_local,φ_local)は、計算要素Ｂからθ_local,φ_localの方向に直線を伸ばした場合の可視判定の結果を示しており、可視判定値と呼ぶことにする。図８は、可視判定値ηについて説明するための模式図である。図８に示すように、上記直線が物質表面とぶつかる場合には、η＝０となり、ぶつからない場合にはη＝１となる。なお、図８に示すように、物質表面の片側の方向にだけ直線を伸ばす場合には、式（１１）におけるθ_localの積分範囲は、０からπではなく、０からπ／２にしてもよい。

可視判定値ηと可視係数νとの違いについては、図９を参照されたい。図９は、可視係数νについて説明するための模式図である。ν(ａ,Ｂ)は、計算要素ａと計算要素Ｂが互いに見えるか否かを示す。計算要素ａ、Ｂを結ぶ直線が、計算要素ａ、Ｂ間において物質表面とぶつかる場合には、ν＝０となり、ぶつからない場合にはν＝１となる。前者の例としては計算要素ｄを、後者の例としては計算要素ｃを参照されたい。

また、ｆ_flatは、平坦面での直接フラックスを示し、入力値として事前に与えられる。さらに、Normは、以下の式（１２）で与えられる規格化定数を示す。さらに、ｆ(θ_local)は、直接フラックスの面積素片の係数を示し、例えば以下の式（１３）で与えられる。

ただし、θ_inは、図１０に示すような入射角度である。図１０は、入射角度θ_inについて説明するための模式図である。入射角度θ_inは、法線ベクトル方向とθ_local,φ_localの方向との間の角度に相当する。よって、ローカル座標系（ｒ_local,θ_local,φ_local）を用いる場合には、θ_in＝θ_localが成り立つ。

以下、図６のフローについて具体的に説明する。

まず、極角θ_localの数列θ_local(ｍ)の値（ｍ＝０、１、…、Ｍ−１）と、偏角φ_localの数列φ_local(ｏ)の値（ｏ＝０、１、…、Ｏ−１）とを計算する（ステップＳ２１）。これは、０からπまでの極角θ_localの範囲をＭ個の領域に分割し、０から２πまでの偏角φ_localの範囲をＯ個の領域に分割することに相当する。後述するように、式（１１）の積分計算は、これらの数列θ_local(ｍ),φ_local(ｏ)を用いて離散化される。

式（１１）の直接フラックスΓ_B,directの計算に、式（１３）に示す面積素片係数を使用する場合には、例えば、以下の式（１４）、式（１５）のような数列θ_local(ｍ),φ_local(ｏ)を用意する。

ここで、数列∂(ｍ)は、以下の式（１６）で与えられる。

式（１４）のθ_local(ｍ)は、ｆ(θ_local)｜sinθ_local｜をθ_local＝０からθ_local＝θ_local(ｍ)まで積分した場合に、積分結果が∂(ｍ)となる角度を表す。この定義から式（１７）の関係が成り立ち、式（１７）から式（１８）が導出され、式（１８）を変形することで式（１４）が得られる。

以上のように、ステップＳ２１では、０からπまでの極角θ_localの範囲を非等間隔に分割し、０から２πまでの偏角φ_localの範囲を等間隔に分割する。なお、本実施形態では、極角θ_localの範囲だけでなく、偏角φ_localの範囲も非等間隔に分割してもよい。また、極角θ_localの積分範囲を０からπ／２とする場合には、０からπまでではなく、０からπ／２までの極角θ_localの範囲をＭ個の領域に分割するようにしてもよい。

次に、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する（ステップＳ２４）。計算要素ａから直線を伸ばす方向は、計算要素ａにおける数列θ_local(ｍ),φ_local(ｏ)により決定する。すなわち、ステップＳ２４では、計算要素ａからθ_local(ｍ),φ_local(ｏ)の方向に直線を伸ばす。よって、各計算要素ａからは、Ｍ×Ｏ本の直線を伸ばすこととなる。ステップＳ２４の処理は、Ｎ個の計算要素ａの各々について行われる。ステップＳ２４の処理を行うブロックは、本開示の判定部の例である。

なお、ステップＳ２４では、鏡面境界条件や周期境界条件を考慮して可視判定を行ってもよい。図１１、図１２はそれぞれ、鏡面境界条件、周期境界条件について説明するための模式図である。このような判定を行っておくと、境界条件を取り込んだフラックス計算をローコストで行うことが可能となる。

以上のように、ステップＳ２４では、複数の計算要素ａからの各直線が、物質表面にぶつかるか否かと、どの計算要素にぶつかるかとを判定する。物質表面にぶつかった直線に関しては、ステップＳ２５の処理を行い、物質表面にぶつからなかった直線に関しては、ステップＳ２６の処理を行う。

ステップＳ２５では、ある計算要素ａからのいずれかの直線が計算要素Ｂにぶつかった場合には、その計算要素ａを、計算要素Ｂの可視計算要素としてカウントする。一方、ある計算要素ａからのいずれの直線も計算要素Ｂにぶつからなかった場合には、その計算要素ａは、計算要素Ｂの可視計算要素としてカウントしない。このような処理をすべての計算要素ａについて行うことで、計算要素Ｂから見ることのできるすべての計算要素ａを特定することができる。なお、この処理は、計算要素Ｂのみに限らず、Ｎ個のすべての計算要素について同様に行う。

一方、ステップＳ２６では、ある計算要素ａからのある直線が物質表面にぶつからなかった場合（すなわち、気層に到達した場合）には、その直線の方向を、計算要素ａの気層可視方向としてカウントする。このような処理をすべての直線について行うことで、気層から各計算要素ａに直接的に反応種が到達する方向をすべて特定することができる。この特定結果は、直接フラックスの計算に利用することができる。例えば、計算要素Ｂにおける直接フラックスの計算には、計算要素Ｂについての気層可視方向のカウント結果が利用される。

図６のフローでは次に、ステップＳ２６のカウント結果を利用して、計算要素Ｂにおける直接フラックスΓ_B,directを計算する（ステップＳ２８）。直接フラックスΓ_B,directは、数列θ_local(ｍ),φ_local(ｏ)により式（１１）を離散化することで、以下の式（１９）のように表される。

ただし、θ_Blocal(ｍ),φ_Blocal(ｏ)は、計算要素Ｂにおける数列θ_local(ｍ),φ_local(ｏ)を表す。式（１９）のη(θ_Blocal,φ_Blocal)は、計算要素Ｂの気層可視方向ではη＝１となり、その他の方向ではη＝０となる。よって、式（１９）は、ステップＳ２６でカウントされた計算要素Ｂの気層可視方向を利用することで計算することができる。

図６のフローでは次に、ステップＳ２５のカウント結果を利用して、計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)と形態係数ｇ(ａ,Ｂ)を計算する（ステップＳ２９）。形態係数ｇ(ａ,Ｂ)は、数列θ_Blocal(ｍ),φ_Blocal(ｏ)により以下の式（２０）のように表すことができる。

ただし、κ(θ_Blocal,φ_Blocal,ａ)は、計算要素Ｂからθ_Blocal,φ_Blocalの方向に計算要素ａが見えるか否かの可視判定の結果を示しており、計算要素可視判定値と呼ぶことにする。計算要素Ｂからθ_Blocal,φ_Blocalの方向に計算要素ａが見える場合にはκ(θ_Blocal,φ_Blocal,ａ)＝１となり、見えない場合にはκ(θ_Blocal,φ_Blocal,ａ)＝０となる。よって、式（２０）は、ステップＳ２５にて計算要素ａが計算要素Ｂの可視計算要素としてカウントされたか否かを参酌することで計算することができる。

計算要素可視判定値κを２次元、３次元で計算する例を、それぞれ図１３、図１４に示す。図１３、図１４はそれぞれ、２次元、３次元における計算要素可視判定値κについて説明するための模式図である。

なお、可視係数ν(ａ,Ｂ)は、式（２０）によるｇ(ａ,Ｂ)の計算結果から算出することができる。具体的には、ｇ(ａ,Ｂ)＝０の場合にはν(ａ,Ｂ)＝０となり、ｇ(ａ,Ｂ)＞０の場合にはν(ａ,Ｂ)＝１となる。

以上のように、ステップＳ２７、Ｓ２８では、ステップＳ２４における判定の結果に基づいて、直接フラックスΓ_B,direct、可視係数ν(ａ,Ｂ)、および形態係数ｇ(ａ,Ｂ)を計算する。ステップＳ２７、Ｓ２８の処理を行うブロックは、本開示の計算部の例である。

図６のフローによるΓ_B,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)の計算結果は、図１や図４のフローにおいて、総フラックスΓ_B、表面成長速度Ｆ_i、レベルセット関数ψ_tなどの計算に用いられる。これらを計算するブロックも、本開示の計算部の例である。

（３）第１実施形態における計算時間
次に、以上の説明を踏まえて、第１実施形態における計算時間について説明する。

従来の方法では、任意の計算要素Ｂの直接フラックスΓ_B,directを計算するのに、計算要素数Ｎに比例する時間がかかる。理由は、計算要素Ｂに関するループ計算をＮ回繰り返し行うからである。また、従来の方法では、任意の計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)や形態係数ｇ(ａ,Ｂ)を計算するのに、Ｎ²に比例する時間がかかる。理由は、計算要素ａに関するループ計算と計算要素Ｂに関するループ計算を、それぞれＮ回繰り返し行うからである。ν(ａ,Ｂ)やｇ(ａ,Ｂ)の計算時間は、鏡面境界条件や周期境界条件を採用するとさらに長くなる。よって、従来の方法における計算時間の多くは、ν(ａ,Ｂ)やｇ(ａ,Ｂ)の計算に費やされる。

一方、本実施形態では、図６に示すように、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、この判定結果に基づいてΓ_B,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を計算する。そのため、ν(ａ,Ｂ)とｇ(ａ,Ｂ)が、Γ_B,directと同様に、計算要素ａに関するループ計算をＮ回繰り返すことで計算される（ステップＳ２２、Ｓ３０を参照）。よって、本実施形態によれば、Γ_B,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)の計算時間を、計算要素数Ｎに比例する時間に抑えることができる。

このような計算時間の短縮の効果は、物質表面に直接的に到達する反応種だけでなく、物質表面に間接的に到達する反応種も考慮する場合に有効である。理由は、式（３）から理解されるように、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)の計算時間の短縮は、間接フラックスΓ_aB,indirectの計算時間の短縮をもたらすからである。よって、本実施形態によれば、物質表面に直接的、間接的に到達する反応種を考慮した形状シミュレーションを高速化することが可能となる。従来の方法に対する本実施形態の計算時間短縮の効果は、鏡面境界条件や周期境界条件を採用する場合にさらに顕著となる。

また、本実施形態では、式（１９）と式（２０）によりΓ_B,directとｇ(ａ,Ｂ)を計算するため、式（１１）や式（１３）のようなsinやcosの計算を深いループの中で行う必要がない。よって、本実施形態によれば、sinやconのように長い計算時間を要する工程を深いループの中からカットして、計算時間をさらに短縮することが可能となる。

また、本実施形態によるｇ(ａ,Ｂ)の計算によれば、Ｎ²回のループ計算でｇ(ａ,Ｂ)を計算する従来の方法に比べて、ｇ(ａ,Ｂ)行列における０要素の個数が多くなる傾向にある。よって、本実施形態によれば、式（６）のような行列方程式を繰り返し解きながら化学反応計算を行う場合に、これらの０要素に着目した計算アルゴリズムを採用することで、計算時間をさらに短縮することが可能となる。さらには、ＣＲＳなどの疎行列保持アルゴリズムを採用することで、０要素が多くなるほど省メモリ化を図ることも可能となる。

また、本実施形態では、Γ_B,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を計算するために、ステップＳ２２〜Ｓ２７において、極角θ_localと偏角θ_localについてのループ計算を行う。そして、ステップＳ２８、Ｓ２９では、このループ計算の結果からΓ_B,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を計算する。このように、本実施形態では、Γ_B,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を、ステップＳ２２〜Ｓ２７という同一のループ計算により並列的に計算するため、計算時間をさらに短縮することができる。

なお、本実施形態では、面積素片係数ｆ(θ)を、式（１３）以外の式で与えてもよい。この場合には、この式で与えられるｆ(θ)に応じた数列θ_local(ｍ)を、ステップＳ２１において決定する。

（４）第１実施形態の変形例
次に、第１実施形態の変形例について説明する。

（４．１）ループ計算の省略
図１５は、図６のステップＳ２２〜Ｓ２７の手順の変形例を示したフローチャート図である。

図１５におけるステップＳ３３〜Ｓ３５は、図６におけるステップＳ２４〜Ｓ２６に対応している。図１５では、ステップＳ３３〜Ｓ３５の処理を、極角θ_localの大きい方向から極角θ_localの小さい方向へと順に行う（ステップＳ３１、Ｓ３２、Ｓ３６、Ｓ３８）。例えば、θ_local＝θ₁となる方向が（θ₁,φ₁）、（θ₁,φ₂）、（θ₁,φ₃）の３つあり、θ_local＝θ₂となる方向が（θ₂,φ₄）、（θ₂,φ₅）の２つある場合には（θ₁＞θ₂）、まず前者の３方向についてステップＳ３３〜Ｓ３５の処理を行い、次に後者の２方向についてステップＳ３３〜Ｓ３５の処理を行う。

このループ計算では、極角θ_localが同一のすべての偏角φ_localにおいてステップＳ３３〜Ｓ３５の処理が終わるごとに、これらの偏角φ_localにおける可視判定値η(θ_local,φ_local)がすべて１であるか否かを確認する（ステップＳ３７）。

そして、これらの可視判定値ηがすべて１である場合（すなわち、これらの偏角φ_localの方向の直線がすべて物質表面にぶつからなかった場合）には、それ以降のループ計算をすべて省略し、図１５のフローを終了する。例えば、極角θ_localがθ₀（θ₀は定数）となる方向において、任意の偏角φ_localで可視判定値ηが１である場合には、極角θ_localがθ₀よりも小さい方向についてのループ計算をすべて省略する。

このような処理を行う理由は、極角θ_localがθ₀となる全方向の直線が気層に到達する場合には、極角θ_localが０〜θ₀の範囲内に物質表面は存在しないことが多く、極角θ_localがθ₀よりも小さい全方向の直線も気層に到達することが多いからである。よって、本変形例では、これらの方向のループ計算を省略し、これらの方向をすべて気層可視方向としてカウントする。これにより、無駄なループ計算を減らし、計算時間をさらに短縮することが可能となる。

（４．２）グローバル座標
また、本実施形態では、図６や図１５のフローにおいて、各計算要素に固有のローカル座標系を使用したが、代わりに、すべての計算要素に共通のグローバル座標系を使用してもよい。

図１６は、グローバル座標系について説明するための図である。図１６(ａ)は、各計算要素の法線ベクトルを示し、図１６(ｂ)は、グローバル座標系の直交座標（ｘ,ｙ,ｚ）と極座標（ｒ,θ,φ）を示す。

グローバル座標系を使用する場合、直接フラックスΓ_B,directと形態係数ｇ(ａ,Ｂ)は、グローバル座標系の数列θ_B(ｍ),φ_B(ｏ),θ_Bin(ｍ)により、それぞれ以下の式（２１）、式（２２）のように表すことができる。

ただし、θ_B(ｍ)、φ_B(ｏ)、θ_Bin(ｍ)はそれぞれ、計算要素Ｂにおける極角θ、偏角φ、入射角度θ_inの数列を表す。

なお、ローカル座標系を使用することには、計算がシンプルになる、誤差が少なくなるなどの利点がある。例えば、図８に示すように物質表面の片側の方向にだけ直線を伸ばす場合、ローカル座標系を使用すれば、極角θ_localの範囲を０〜πから０〜π／２に変更するだけでこれに対処することができる。よって、この場合にはローカル座標系を使用することで計算がシンプルになり、その結果、誤差も少なくなる。一方、グローバル座標系を使用することには、計算要素ごとの座標系の違いを考慮する必要がなくなるという利点がある。

（４．３）偏角の分割数
また、本実施形態では、偏角φ_localの分割数Ｏは極角θ_localによらず一定としたが、代わりに、偏角φ_localの分割数Ｏは極角θ_localに応じて変化させてもよい。すなわち、偏角φ_localの分割数Ｏは、極角θ_localの刻みｍに依存する変数としてもよい。

よって、極角θ_local(ｍ)における偏角φ_localの刻みｏ、分割数Ｏをそれぞれｏ_m、Ｏ_mで表すと、式（１５）、式（１９）、式（２０）はそれぞれ以下の式（２３）、式（２４）、式（２５）のように書き換えられる。

前述または後述のような第１実施形態の効果は、このような偏角φ_localの分割方法を用いる場合にも享受することができる。なお、このような分割方法は、ローカル座標系だけでなく、グローバル座標系にも適用可能である。

（５）第１実施形態の効果
最後に、第１実施形態の効果について説明する。

以上のように、本実施形態では、各計算要素から複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、この判定結果に基づいて直接フラックスと形態係数を計算する。さらには、この判定結果に基づいて可視係数を計算する。

よって、本実施形態によれば、直接フラックスや形態係数の計算時間を、計算要素数に比例する時間に抑えることができる。よって、本実施形態によれば、間接フラックスの計算時間に影響を与える形態係数の計算時間を短縮することで、物質表面に直接的、間接的に到達する反応種を考慮した形状シミュレーションを高速化することが可能となる。

図１７、図１８はそれぞれ、比較例、第１実施形態における計算時間の例を示したグラフである。比較例では、従来の方法を用いて直接フラックス、可視係数、形態係数を計算した。図１７、図１８は、直接フラックスの計算時間、可視計算（可視係数と形態係数の計算）の計算時間、化学反応収束計算の計算時間、全計算時間の合計を示している。

図１７、図１８に示すように、第１実施形態によれば、全計算時間を比較例よりも短縮することが可能となる。これらの比較結果を図１９に示す。図１９は、計算要素数が４万の場合の第１実施形態と比較例の計算時間を比較したグラフである。

図２０は、第１実施形態と比較例におけるθ分割数と計算誤差との関係を示したグラフである。ただし、図２０における計算では、ローカル座標系を使用した。また、ｎ＝１とｎ＝１０００のグラフは、比較例による計算結果を示す。図２０に示すように、第１実施形態によれば、計算誤差を抑制できることが分かる。

なお、第１実施形態の形状シミュレーション方法は、どのような情報処理装置を用いて実行してもよい。第２実施形態では、このような情報処理装置の一例である形状シミュレーション装置について説明する。

（第２実施形態）
図２１は、第２実施形態の形状シミュレーション装置の構成を示す外観図である。

図２１の形状シミュレーション装置は、制御部１１と、表示部１２と、入力部１３とを備えている。

制御部１１は、形状シミュレーション装置の動作を制御するモジュールである。制御部１１は、例えば、第１実施形態の形状シミュレーション方法を実行する。制御部１１の詳細については後述する。

表示部１２は、液晶モニタなどの表示デバイスを有している。表示部１２は、例えば、形状シミュレーション用の設定情報の入力画面や、形状シミュレーションの計算結果などを表示する。

入力部１３は、キーボード１３ａやマウス１３ｂなどの入力デバイスを有している。入力部１３は、例えば、形状シミュレーション用の設定情報の入力用に使用される。設定情報の例としては、計算式に関する情報、実験値や予測値に関する情報、物質の構造に関する情報、フラックスに関する情報、形状シミュレーションの条件や手順に関する指示情報などが挙げられる。

図２２は、図２１の制御部１１の構成を示すブロック図である。

制御部１１は、ＣＰＵ（Central Processing Unit）２１と、ＲＯＭ（Read Only Memory）２２と、ＲＡＭ（Random Access Memory）２３と、ＨＤＤ（Hard Disk Drive）２４と、ＣＤ（Compact Disc）ドライブなどのメモリドライブ２５と、メモリポートやメモリスロットなどのメモリＩ／Ｆ（interface）２６とを備えている。

本実施形態では、第１実施形態の形状シミュレーション方法用のプログラムである形状シミュレーションプログラムが、ＲＯＭ２２またはＨＤＤ２４内に格納されている。入力部１３から所定の指示情報が入力されると、ＣＰＵ２１は、ＲＯＭ２２またはＨＤＤ２４からプログラムを読み出し、読み出したプログラムをＲＡＭ２３に展開し、このプログラムにより形状シミュレーションを実行する。この処理の際に生じる各種データは、ＲＡＭ２３内に保持される。

なお、本実施形態では、形状シミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体を用意し、この記録媒体からＲＯＭ２２やＨＤＤ２４内に形状シミュレーションプログラムをインストールしてもよい。このような記録媒体の例としては、ＣＤ−ＲＯＭやＤＶＤ（Digital Versatile Disk）−ＲＯＭなどが挙げられる。

また、本実施形態では、形状シミュレーションプログラムを、インターネットなどのネットワーク経由でダウンロードすることで、ＲＯＭ２２やＨＤＤ２４内にインストールしてもよい。

以上のように、本実施形態によれば、第１実施形態の形状シミュレーション方法を実行するための形状シミュレーション装置や形状シミュレーションプログラムを提供することが可能となる。

なお、第１及び第２実施形態では、形状シミュレーションの適用対象の例として半導体デバイスを取り上げたが、この形状シミュレーションは半導体デバイス以外のデバイスにも適用可能である。このようなデバイスの例としては、ＭＥＭＳ（Micro Electro Mechanical Systems）デバイスやディスプレイデバイスなどが挙げられる。

以上、第１及び第２実施形態について説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することを意図したものではない。これらの実施形態は、その他の様々な形態で実施することができる。また、これらの実施形態に対し、発明の要旨を逸脱しない範囲内で、種々の省略、置換、変更を行うことにより、様々な変形例を得ることもできる。これらの形態や変形例は、発明の範囲や要旨に含まれており、特許請求の範囲及びこれに均等な範囲には、これらの形態や変形例が含まれる。

１：シリコン基板、２：シリコン窒化膜、３：シリコン酸化膜、４：貫通孔、
１１：制御部、１２：表示部、１３：入力部、
２１：ＣＰＵ、２２：ＲＯＭ、２３：ＲＡＭ、２４：ＨＤＤ、
２５：メモリドライブ、２６：メモリＩ／Ｆ

Claims (6)

1. 物質の表面を複数の計算要素に分割する分割部と、
   各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
   前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部とを備え、
   前記計算部は、前記直接フラックスと前記形態係数とを使用して、各計算要素に直接的または間接的に到達する反応種のフラックスである総フラックスと、前記物質の局所的な表面成長速度の少なくともいずれかを計算し、
   前記計算部は、前記直接フラックスと前記形態係数を、前記直線の方向を示す極角θと偏角φについてのループ計算により計算し、
   前記計算部は、前記ループ計算を前記極角θの大きい方向から前記極角θの小さい方向へと順に行い、前記極角θがθ₀となる方向において任意の前記偏角φで各直線が前記物質の表面にぶつからないと判定された場合には、前記極角θがθ₀よりも小さい方向についての前記ループ計算を省略する、
   形状シミュレーション装置。
2. 物質の表面を複数の計算要素に分割する分割部と、
   各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
   前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部と、
   を備える形状シミュレーション装置。
3. 前記計算部は、前記直接フラックスと前記形態係数とを使用して、各計算要素に直接的または間接的に到達する反応種のフラックスである総フラックスと、前記物質の局所的な表面成長速度の少なくともいずれかを計算する、
   請求項２に記載の形状シミュレーション装置。
4. 前記計算部は、前記直接フラックスと前記形態係数を、前記直線の方向を示す極角θと偏角φについてのループ計算により計算し、
   前記計算部は、前記ループ計算を前記極角θの大きい方向から前記極角θの小さい方向へと順に行い、前記極角θがθ₀となる方向において任意の前記偏角φで各直線が前記物質の表面にぶつからないと判定された場合には、前記極角θがθ₀よりも小さい方向についての前記ループ計算を省略する、
   請求項２または３に記載の形状シミュレーション装置。
5. 物質の表面を複数の計算要素に分割し、
   各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、
   前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する、
   ことを含む形状シミュレーション方法。
6. 物質の表面を複数の計算要素に分割し、
   各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、
   前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する、
   ことを含む形状シミュレーション方法をコンピュータに実行させる形状シミュレーションプログラム。

Images