(57)【特許請求の範囲】
【請求項１】
物質の表面を複数の計算要素に分割する分割部と、
各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部とを備え、
前記計算部は、前記直接フラックスと前記形態係数とを使用して、各計算要素に直接的または間接的に到達する反応種のフラックスである総フラックスと、前記物質の局所的な表面成長速度の少なくともいずれかを計算し、
前記計算部は、前記直接フラックスと前記形態係数を、前記直線の方向を示す極角θと偏角φについてのループ計算により計算し、
前記計算部は、前記ループ計算を前記極角θの大きい方向から前記極角θの小さい方向へと順に行い、前記極角θがθ0となる方向において任意の前記偏角φで各直線が前記物質の表面にぶつからないと判定された場合には、前記極角θがθ0よりも小さい方向についての前記ループ計算を省略する、
形状シミュレーション装置。
【請求項２】
物質の表面を複数の計算要素に分割する分割部と、
各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部とを備え、
前記計算部は、前記複数の計算要素のうちの第１計算要素における前記直接フラックスを、前記第１計算要素から伸ばして前記物質の表面にぶつからなかった直線を用いて計算し、
前記計算部は、前記複数の計算要素のうちの前記第１計算要素および第２計算要素に関する前記形態係数を、前記第１計算要素から伸ばして前記第２計算要素にぶつかった直線を用いて計算し、
前記計算部は、前記直接フラックスと前記形態係数とを使用して、各計算要素に直接的または間接的に到達する反応種のフラックスである総フラックスと、前記物質の局所的な表面成長速度の少なくともいずれかを計算する、
形状シミュレーション装置。
【請求項３】
前記判定部は、前記第１計算要素から複数の方向に直線を伸ばす際に、各直線を伸ばす方向を、前記第２計算要素の位置に基づかずに決定する、
請求項２に記載の形状シミュレーション装置。
【請求項４】
前記計算部は、前記直接フラックスと前記形態係数を、前記直線の方向を示す極角θと偏角φについてのループ計算により計算し、
前記計算部は、前記ループ計算を前記極角θの大きい方向から前記極角θの小さい方向へと順に行い、前記極角θがθ0となる方向において任意の前記偏角φで各直線が前記物質の表面にぶつからないと判定された場合には、前記極角θがθ0よりも小さい方向についての前記ループ計算を省略する、
請求項２または３に記載の形状シミュレーション装置。
【請求項５】
形状シミュレーション装置が、物質の表面を複数の計算要素に分割し、
前記形状シミュレーション装置が、各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、
前記形状シミュレーション装置が、前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算し、
前記形状シミュレーション装置が、前記直接フラックスと前記形態係数とを使用して、各計算要素に直接的または間接的に到達する反応種のフラックスである総フラックスと、前記物質の局所的な表面成長速度の少なくともいずれかを計算する、
ことを含み、
前記形状シミュレーション装置は、前記複数の計算要素のうちの第１計算要素における前記直接フラックスを、前記第１計算要素から伸ばして前記物質の表面にぶつからなかった直線を用いて計算し、
前記形状シミュレーション装置は、前記複数の計算要素のうちの前記第１計算要素および第２計算要素に関する前記形態係数を、前記第１計算要素から伸ばして前記第２計算要素にぶつかった直線を用いて計算する、
形状シミュレーション方法。
【請求項６】
物質の表面を複数の計算要素に分割し、
各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、
前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算し、
前記直接フラックスと前記形態係数とを使用して、各計算要素に直接的または間接的に到達する反応種のフラックスである総フラックスと、前記物質の局所的な表面成長速度の少なくともいずれかを計算する、
ことを含む形状シミュレーション方法をコンピュータに実行させる形状シミュレーションプログラムであって、
前記複数の計算要素のうちの第１計算要素における前記直接フラックスは、前記第１計算要素から伸ばして前記物質の表面にぶつからなかった直線を用いて計算され、
前記複数の計算要素のうちの前記第１計算要素および第２計算要素に関する前記形態係数は、前記第１計算要素から伸ばして前記第２計算要素にぶつかった直線を用いて計算される、
形状シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラムに関する。
【背景技術】
【０００２】
ＣＶＤ（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は、以下の式（１）から計算できる。
【数１】
0005918630.tif 000002
ただし、∇はベクトル微分演算子を表し、｜∇ψt｜は∇ψtのノルムを表す。時間Δｔ経過後のレベルセット関数は、式（１）を離散化した式に従い、レベルセット関数を時間発展させることで計算可能である。なお、本実施形態では、表面形状を時間発展させる代わりに、ある表面形状における表面成長速度Ｆやフラックスを計算してもよい。これは、後述するステップＳ５を１ステップ目でＹｅｓと判定する場合に相当する。
【００１４】
次に、予め設定したプロセス時間が経過したか否かを判定する（ステップＳ５）。プロセス時間が終了した場合には、物質の最終形状を出力し（ステップＳ６）、計算終了となる。プロセス時間が終了していない場合には、ステップＳ３に戻る。
【００１５】
なお、本実施形態では、形状表現の手法としてレベルセット法を用いているが、レベルセット法以外のセル法やストリング法などの手法を用いてもよい。
【００１６】
（１）ステップＳ３の詳細
次に、図４を参照し、ステップＳ３の詳細について説明する。
【００１７】
図４は、図１のステップＳ３の詳細を示したフローチャート図である。
【００１８】
まず、レベルセット法で表された物質表面を、複数の計算要素に分割する（ステップＳ１１）。図５は、物質表面を複数の計算要素に分割した様子を示した模式図である。図５の例では、物質表面がメッシュごとに分割されている。その結果、１つのメッシュ内の物質表面が、１つの計算要素となっている。ステップＳ１１の処理を行うブロックは、本開示の分割部の例である。
【００１９】
なお、物質表面の分割方法は、メッシュ単位に限られるものではなく、どのような方法を採用してもよい。また、物質表面の分割は、タイムステップごとに行う必要はなく、例えばステップＳ１の直後に行ってもよい。
【００２０】
また、図５に示す計算領域は、２次元領域となっているが、代わりに３次元領域としてもよい。また、図５に示す各計算要素の形状は、線分となっているが、代わりに点や多角形などとしてもよい。
【００２１】
図５は、第１の計算要素ａと、第２の計算要素Ｂを示している。計算要素Ｂに到達する反応種のフラックスを計算する際には、気層から計算要素Ｂに直接的に到達する反応種のフラックスと、気層から任意の計算要素ａを介して計算要素Ｂに間接的に到達する反応種のフラックスの両方を考慮するのが一般的である。前者のフラックスを直接フラックスと呼び、後者のフラックスを間接フラックスと呼ぶ。また、これらの合計を、総フラックスと呼ぶ。なお、反応種の例としては、堆積種やエッチング種などが挙げられる。
【００２２】
計算要素Ｂにおける総フラックスΓBは、以下の式（２）のように、計算要素Ｂにおける直接フラックスΓB,directと、任意の計算要素ａからの間接フラックスΓaB,indirectの合計との和で表される。
【数２】
0005918630.tif 000003
【００２３】
ここで、間接フラックスΓaB,indirectは、例えば以下の式（３）で表すことができる。
【数３】
0005918630.tif 000004
Ｓa(Γa)は、計算要素ａで吸収されるフラックスの割合を示す付着確率を表す。Ｓa(Γa)の値は、計算要素ａにおける総フラックスΓaに依存する。また、ν(ａ,Ｂ)は、計算要素ａと計算要素Ｂが互いに見えるか否かを示す可視係数（面間可視係数）を表す。計算要素ａ、Ｂを結ぶ直線が、計算要素ａ、Ｂ間において物質表面とぶつかる場合には、ν＝０となり、ぶつからない場合にはν＝１となる。また、ｇ(ａ,Ｂ)は、計算要素ａと計算要素Ｂとの位置関係（面関係）を示す形態係数を表す。ｇ(ａ,Ｂ)の値は、計算要素ａ、Ｂの互いの見えやすさの程度を示している。ｇ(ａ,Ｂ)の値は、計算要素ａ、Ｂ間の距離や角度などに依存する。
【００２４】
式（２）に式（３）を代入すると、計算要素Ｂにおける総フラックスΓBは、以下の式（４）で表すことができる。
【数４】
0005918630.tif 000005
【００２５】
図４のフローでは次に、任意の計算要素における直接フラックスや、任意の計算要素間における可視係数ν、形態係数ｇを計算する（ステップＳ１２）。
【００２６】
次に、各計算要素ｉの直接フラックスΓi,directを仮の総フラックスΓiとして使用し、各計算要素ｉにおける付着確率Ｓi(Γi)を計算する（ステップＳ１３）。この際、このフラックスは、中性分子を含んでいても、志向性を持ったイオンを含んでいてもよいし、その両方を含んでいてもよい。
【００２７】
次に、可視係数ν、形態係数ｇ、直接フラックスΓi,direct、付着確率Ｓi(Γi)を用いて、以下の式（５）から、各計算要素ｉにおける総フラックスΓiを計算する（ステップＳ１４）。
【数５】
0005918630.tif 000006
【００２８】
次に、ステップＳ１３とステップＳ１４の処理を、付着確率Ｓi(Γi)の値が収束するまで繰り返す（ステップＳ１５）。なお、２回目以降のステップＳ１３では、前回のステップＳ１４で計算した総フラックスΓiを、仮の総フラックスΓiとして使用する。また、ステップＳ１５では、Ｓi(Γi)の値が収束したか否かを、Ｓi(Γi)の変化が閾値以下になったか否かで判定する。そして、Ｓi(Γi)の値が収束した際の総フラックスΓiを、総フラックスΓiの正しい計算結果として取り扱う。
【００２９】
なお、計算要素の個数をＮとする場合、任意の計算要素間における可視係数νや形態係数ｇは、まとめてＮ×Ｎ行列で表すことができる。行列の形で表された可視係数ν、形態係数ｇを、それぞれ可視係数行列、形態係数行列と呼ぶ。また、任意の計算要素におけるフラックスは、まとめてＮ行ベクトルで表すことができる。ベクトルの形で表されたフラックスを、フラックスベクトルと呼ぶ。
【００３０】
この場合、式（５）は、以下の式（６）のように行列方程式で表現することができる。
【数６】
0005918630.tif 000007
【数７】
0005918630.tif 000008
【数８】
0005918630.tif 000009
【数９】
0005918630.tif 000010
ただし、Ｉ、Ｊは処理対象の計算要素の個数を表し、例えばＩ＝Ｊ＝Ｎである。行列方程式（６）は、どのような解法で解いてもよい。解法の例としては、反復法（ガウスザイデル法、ＳＯＲ法、ヤコビ法、共役勾配法など）や、直接法（ガウスの消去法、ＬＵ分解法、コレスキー分解法など）が挙げられる。行列方程式（６）を解く際、行列Ａが疎行列である場合には、ＣＲＳなどの保存方法を用いた上で疎行列に適したルーチンを使用することで、計算処理の省メモリ化や高速化を図ってもよい。
【００３１】
図４のフローでは次に、総フラックスΓiから、各計算要素ｉにおける局所的な表面成長速度Ｆiを計算する（ステップＳ１６）。例えば、Ｋ種類の反応種を使用する場合には、表面成長速度Ｆiは、Ｋ個の局所的な総フラックスΓ1,i〜ΓK,iに依存した以下の式（１０）のような形でモデリングされる。
【数１０】
0005918630.tif 000011
ただし、ｋは１≦ｋ≦Ｋを満たす任意の実数である。以上のようにして、ステップＳ３の処理が終了する。
【００３２】
（２）ステップＳ１２の詳細
次に、図６を参照し、ステップＳ１２の詳細について説明する。
【００３３】
図６は、図４のステップＳ１２の詳細を示したフローチャート図である。
【００３４】
図６のフローでは、各計算要素に固有のローカル座標系を使用する。図７は、ローカル座標系について説明するための図である。図７(ａ)は、各計算要素の法線ベクトルを示し、図７(ｂ)は、各計算要素におけるローカル座標系を示す。図７(ｂ)に示すように、ローカル座標系の直交座標（ｘlocal,ｙlocal,ｚlocal）は、＋ｚlocal方向が法線ベクトル方向と一致するように定められる。また、ローカル座標系の極座標（ｒlocal,θlocal,φlocal）は、極角θlocalが動径ｒlocalと＋ｚlocal方向との間の角度となり、偏角φlocalが動径ｒlocalと＋ｘlocal方向との間の角度となるように定められる。
【００３５】
計算要素Ｂにおける直接フラックスΓB,directは、以下の式（１１）により計算することができる。
【数１１】
0005918630.tif 000012
ただし、η(θlocal,φlocal)は、計算要素Ｂからθlocal,φlocalの方向に直線を伸ばした場合の可視判定の結果を示しており、可視判定値と呼ぶことにする。図８は、可視判定値ηについて説明するための模式図である。図８に示すように、上記直線が物質表面とぶつかる場合には、η＝０となり、ぶつからない場合にはη＝１となる。なお、図８に示すように、物質表面の片側の方向にだけ直線を伸ばす場合には、式（１１）におけるθlocalの積分範囲は、０からπではなく、０からπ／２にしてもよい。
【００３６】
可視判定値ηと可視係数νとの違いについては、図９を参照されたい。図９は、可視係数νについて説明するための模式図である。ν(ａ,Ｂ)は、計算要素ａと計算要素Ｂが互いに見えるか否かを示す。計算要素ａ、Ｂを結ぶ直線が、計算要素ａ、Ｂ間において物質表面とぶつかる場合には、ν＝０となり、ぶつからない場合にはν＝１となる。前者の例としては計算要素ｄを、後者の例としては計算要素ｃを参照されたい。
【００３７】
また、ｆflatは、平坦面での直接フラックスを示し、入力値として事前に与えられる。さらに、Normは、以下の式（１２）で与えられる規格化定数を示す。さらに、ｆ(θlocal)は、直接フラックスの面積素片の係数を示し、例えば以下の式（１３）で与えられる。
【数１２】
0005918630.tif 000013
【数１３】
0005918630.tif 000014
ただし、θinは、図１０に示すような入射角度である。図１０は、入射角度θinについて説明するための模式図である。入射角度θinは、法線ベクトル方向とθlocal,φlocalの方向との間の角度に相当する。よって、ローカル座標系（ｒlocal,θlocal,φlocal）を用いる場合には、θin＝θlocalが成り立つ。
【００３８】
以下、図６のフローについて具体的に説明する。
【００３９】
まず、極角θlocalの数列θlocal(ｍ)の値（ｍ＝０、１、…、Ｍ−１）と、偏角φlocalの数列φlocal(ｏ)の値（ｏ＝０、１、…、Ｏ−１）とを計算する（ステップＳ２１）。これは、０からπまでの極角θlocalの範囲をＭ個の領域に分割し、０から２πまでの偏角φlocalの範囲をＯ個の領域に分割することに相当する。後述するように、式（１１）の積分計算は、これらの数列θlocal(ｍ),φlocal(ｏ)を用いて離散化される。
【００４０】
式（１１）の直接フラックスΓB,directの計算に、式（１３）に示す面積素片係数を使用する場合には、例えば、以下の式（１４）、式（１５）のような数列θlocal(ｍ),φlocal(ｏ)を用意する。
【数１４】
0005918630.tif 000015
【数１５】
0005918630.tif 000016
ここで、数列∂(ｍ)は、以下の式（１６）で与えられる。
【数１６】
0005918630.tif 000017
【００４１】
式（１４）のθlocal(ｍ)は、ｆ(θlocal)｜sinθlocal｜をθlocal＝０からθlocal＝θlocal(ｍ)まで積分した場合に、積分結果が∂(ｍ)となる角度を表す。この定義から式（１７）の関係が成り立ち、式（１７）から式（１８）が導出され、式（１８）を変形することで式（１４）が得られる。
【数１７】
0005918630.tif 000018
【数１８】
0005918630.tif 000019
【００４２】
以上のように、ステップＳ２１では、０からπまでの極角θlocalの範囲を非等間隔に分割し、０から２πまでの偏角φlocalの範囲を等間隔に分割する。なお、本実施形態では、極角θlocalの範囲だけでなく、偏角φlocalの範囲も非等間隔に分割してもよい。また、極角θlocalの積分範囲を０からπ／２とする場合には、０からπまでではなく、０からπ／２までの極角θlocalの範囲をＭ個の領域に分割するようにしてもよい。
【００４３】
次に、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する（ステップＳ２４）。計算要素ａから直線を伸ばす方向は、計算要素ａにおける数列θlocal(ｍ),φlocal(ｏ)により決定する。すなわち、ステップＳ２４では、計算要素ａからθlocal(ｍ),φlocal(ｏ)の方向に直線を伸ばす。よって、各計算要素ａからは、Ｍ×Ｏ本の直線を伸ばすこととなる。ステップＳ２４の処理は、Ｎ個の計算要素ａの各々について行われる。ステップＳ２４の処理を行うブロックは、本開示の判定部の例である。
【００４４】
なお、ステップＳ２４では、鏡面境界条件や周期境界条件を考慮して可視判定を行ってもよい。図１１、図１２はそれぞれ、鏡面境界条件、周期境界条件について説明するための模式図である。このような判定を行っておくと、境界条件を取り込んだフラックス計算をローコストで行うことが可能となる。
【００４５】
以上のように、ステップＳ２４では、複数の計算要素ａからの各直線が、物質表面にぶつかるか否かと、どの計算要素にぶつかるかとを判定する。物質表面にぶつかった直線に関しては、ステップＳ２５の処理を行い、物質表面にぶつからなかった直線に関しては、ステップＳ２６の処理を行う。
【００４６】
ステップＳ２５では、ある計算要素ａからのいずれかの直線が計算要素Ｂにぶつかった場合には、その計算要素ａを、計算要素Ｂの可視計算要素としてカウントする。一方、ある計算要素ａからのいずれの直線も計算要素Ｂにぶつからなかった場合には、その計算要素ａは、計算要素Ｂの可視計算要素としてカウントしない。このような処理をすべての計算要素ａについて行うことで、計算要素Ｂから見ることのできるすべての計算要素ａを特定することができる。なお、この処理は、計算要素Ｂのみに限らず、Ｎ個のすべての計算要素について同様に行う。
【００４７】
一方、ステップＳ２６では、ある計算要素ａからのある直線が物質表面にぶつからなかった場合（すなわち、気層に到達した場合）には、その直線の方向を、計算要素ａの気層可視方向としてカウントする。このような処理をすべての直線について行うことで、気層から各計算要素ａに直接的に反応種が到達する方向をすべて特定することができる。この特定結果は、直接フラックスの計算に利用することができる。例えば、計算要素Ｂにおける直接フラックスの計算には、計算要素Ｂについての気層可視方向のカウント結果が利用される。
【００４８】
図６のフローでは次に、ステップＳ２６のカウント結果を利用して、計算要素Ｂにおける直接フラックスΓB,directを計算する（ステップＳ２８）。直接フラックスΓB,directは、数列θlocal(ｍ),φlocal(ｏ)により式（１１）を離散化することで、以下の式（１９）のように表される。
【数１９】
0005918630.tif 000020
ただし、θBlocal(ｍ),φBlocal(ｏ)は、計算要素Ｂにおける数列θlocal(ｍ),φlocal(ｏ)を表す。式（１９）のη(θBlocal,φBlocal)は、計算要素Ｂの気層可視方向ではη＝１となり、その他の方向ではη＝０となる。よって、式（１９）は、ステップＳ２６でカウントされた計算要素Ｂの気層可視方向を利用することで計算することができる。
【００４９】
図６のフローでは次に、ステップＳ２５のカウント結果を利用して、計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)と形態係数ｇ(ａ,Ｂ)を計算する（ステップＳ２９）。形態係数ｇ(ａ,Ｂ)は、数列θBlocal(ｍ),φBlocal(ｏ)により以下の式（２０）のように表すことができる。
【数２０】
0005918630.tif 000021
ただし、κ(θBlocal,φBlocal,ａ)は、計算要素ＢからθBlocal,φBlocalの方向に計算要素ａが見えるか否かの可視判定の結果を示しており、計算要素可視判定値と呼ぶことにする。計算要素ＢからθBlocal,φBlocalの方向に計算要素ａが見える場合にはκ(θBlocal,φBlocal,ａ)＝１となり、見えない場合にはκ(θBlocal,φBlocal,ａ)＝０となる。よって、式（２０）は、ステップＳ２５にて計算要素ａが計算要素Ｂの可視計算要素としてカウントされたか否かを参酌することで計算することができる。
【００５０】
計算要素可視判定値κを２次元、３次元で計算する例を、それぞれ図１３、図１４に示す。図１３、図１４はそれぞれ、２次元、３次元における計算要素可視判定値κについて説明するための模式図である。
【００５１】
なお、可視係数ν(ａ,Ｂ)は、式（２０）によるｇ(ａ,Ｂ)の計算結果から算出することができる。具体的には、ｇ(ａ,Ｂ)＝０の場合にはν(ａ,Ｂ)＝０となり、ｇ(ａ,Ｂ)＞０の場合にはν(ａ,Ｂ)＝１となる。
【００５２】
以上のように、ステップＳ２７、Ｓ２８では、ステップＳ２４における判定の結果に基づいて、直接フラックスΓB,direct、可視係数ν(ａ,Ｂ)、および形態係数ｇ(ａ,Ｂ)を計算する。ステップＳ２７、Ｓ２８の処理を行うブロックは、本開示の計算部の例である。
【００５３】
図６のフローによるΓB,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)の計算結果は、図１や図４のフローにおいて、総フラックスΓB、表面成長速度Ｆi、レベルセット関数ψtなどの計算に用いられる。これらを計算するブロックも、本開示の計算部の例である。
【００５４】
（３）第１実施形態における計算時間
次に、以上の説明を踏まえて、第１実施形態における計算時間について説明する。
【００５５】
従来の方法では、任意の計算要素Ｂの直接フラックスΓB,directを計算するのに、計算要素数Ｎに比例する時間がかかる。理由は、計算要素Ｂに関するループ計算をＮ回繰り返し行うからである。また、従来の方法では、任意の計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)や形態係数ｇ(ａ,Ｂ)を計算するのに、Ｎ2に比例する時間がかかる。理由は、計算要素ａに関するループ計算と計算要素Ｂに関するループ計算を、それぞれＮ回繰り返し行うからである。ν(ａ,Ｂ)やｇ(ａ,Ｂ)の計算時間は、鏡面境界条件や周期境界条件を採用するとさらに長くなる。よって、従来の方法における計算時間の多くは、ν(ａ,Ｂ)やｇ(ａ,Ｂ)の計算に費やされる。
【００５６】
一方、本実施形態では、図６に示すように、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、この判定結果に基づいてΓB,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を計算する。そのため、ν(ａ,Ｂ)とｇ(ａ,Ｂ)が、ΓB,directと同様に、計算要素ａに関するループ計算をＮ回繰り返すことで計算される（ステップＳ２２、Ｓ３０を参照）。よって、本実施形態によれば、ΓB,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)の計算時間を、計算要素数Ｎに比例する時間に抑えることができる。
【００５７】
このような計算時間の短縮の効果は、物質表面に直接的に到達する反応種だけでなく、物質表面に間接的に到達する反応種も考慮する場合に有効である。理由は、式（３）から理解されるように、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)の計算時間の短縮は、間接フラックスΓaB,indirectの計算時間の短縮をもたらすからである。よって、本実施形態によれば、物質表面に直接的、間接的に到達する反応種を考慮した形状シミュレーションを高速化することが可能となる。従来の方法に対する本実施形態の計算時間短縮の効果は、鏡面境界条件や周期境界条件を採用する場合にさらに顕著となる。
【００５８】
また、本実施形態では、式（１９）と式（２０）によりΓB,directとｇ(ａ,Ｂ)を計算するため、式（１１）や式（１３）のようなsinやcosの計算を深いループの中で行う必要がない。よって、本実施形態によれば、sinやconのように長い計算時間を要する工程を深いループの中からカットして、計算時間をさらに短縮することが可能となる。
【００５９】
また、本実施形態によるｇ(ａ,Ｂ)の計算によれば、Ｎ2回のループ計算でｇ(ａ,Ｂ)を計算する従来の方法に比べて、ｇ(ａ,Ｂ)行列における０要素の個数が多くなる傾向にある。よって、本実施形態によれば、式（６）のような行列方程式を繰り返し解きながら化学反応計算を行う場合に、これらの０要素に着目した計算アルゴリズムを採用することで、計算時間をさらに短縮することが可能となる。さらには、ＣＲＳなどの疎行列保持アルゴリズムを採用することで、０要素が多くなるほど省メモリ化を図ることも可能となる。
【００６０】
また、本実施形態では、ΓB,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を計算するために、ステップＳ２２〜Ｓ２７において、極角θlocalと偏角θlocalについてのループ計算を行う。そして、ステップＳ２８、Ｓ２９では、このループ計算の結果からΓB,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を計算する。このように、本実施形態では、ΓB,direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)を、ステップＳ２２〜Ｓ２７という同一のループ計算により並列的に計算するため、計算時間をさらに短縮することができる。
【００６１】
なお、本実施形態では、面積素片係数ｆ(θ)を、式（１３）以外の式で与えてもよい。この場合には、この式で与えられるｆ(θ)に応じた数列θlocal(ｍ)を、ステップＳ２１において決定する。
【００６２】
（４）第１実施形態の変形例
次に、第１実施形態の変形例について説明する。
【００６３】
（４．１）ループ計算の省略
図１５は、図６のステップＳ２２〜Ｓ２７の手順の変形例を示したフローチャート図である。
【００６４】
図１５におけるステップＳ３３〜Ｓ３５は、図６におけるステップＳ２４〜Ｓ２６に対応している。図１５では、ステップＳ３３〜Ｓ３５の処理を、極角θlocalの大きい方向から極角θlocalの小さい方向へと順に行う（ステップＳ３１、Ｓ３２、Ｓ３６、Ｓ３８）。例えば、θlocal＝θ1となる方向が（θ1,φ1）、（θ1,φ2）、（θ1,φ3）の３つあり、θlocal＝θ2となる方向が（θ2,φ4）、（θ2,φ5）の２つある場合には（θ1＞θ2）、まず前者の３方向についてステップＳ３３〜Ｓ３５の処理を行い、次に後者の２方向についてステップＳ３３〜Ｓ３５の処理を行う。
【００６５】
このループ計算では、極角θlocalが同一のすべての偏角φlocalにおいてステップＳ３３〜Ｓ３５の処理が終わるごとに、これらの偏角φlocalにおける可視判定値η(θlocal,φlocal)がすべて１であるか否かを確認する（ステップＳ３７）。
【００６６】
そして、これらの可視判定値ηがすべて１である場合（すなわち、これらの偏角φlocalの方向の直線がすべて物質表面にぶつからなかった場合）には、それ以降のループ計算をすべて省略し、図１５のフローを終了する。例えば、極角θlocalがθ0（θ0は定数）となる方向において、任意の偏角φlocalで可視判定値ηが１である場合には、極角θlocalがθ0よりも小さい方向についてのループ計算をすべて省略する。
【００６７】
このような処理を行う理由は、極角θlocalがθ0となる全方向の直線が気層に到達する場合には、極角θlocalが０〜θ0の範囲内に物質表面は存在しないことが多く、極角θlocalがθ0よりも小さい全方向の直線も気層に到達することが多いからである。よって、本変形例では、これらの方向のループ計算を省略し、これらの方向をすべて気層可視方向としてカウントする。これにより、無駄なループ計算を減らし、計算時間をさらに短縮することが可能となる。
【００６８】
（４．２）グローバル座標
また、本実施形態では、図６や図１５のフローにおいて、各計算要素に固有のローカル座標系を使用したが、代わりに、すべての計算要素に共通のグローバル座標系を使用してもよい。
【００６９】
図１６は、グローバル座標系について説明するための図である。図１６(ａ)は、各計算要素の法線ベクトルを示し、図１６(ｂ)は、グローバル座標系の直交座標（ｘ,ｙ,ｚ）と極座標（ｒ,θ,φ）を示す。
【００７０】
グローバル座標系を使用する場合、直接フラックスΓB,directと形態係数ｇ(ａ,Ｂ)は、グローバル座標系の数列θB(ｍ),φB(ｏ),θBin(ｍ)により、それぞれ以下の式（２１）、式（２２）のように表すことができる。
【数２１】
0005918630.tif 000022
【数２２】
0005918630.tif 000023
ただし、θB(ｍ)、φB(ｏ)、θBin(ｍ)はそれぞれ、計算要素Ｂにおける極角θ、偏角φ、入射角度θinの数列を表す。
【００７１】
なお、ローカル座標系を使用することには、計算がシンプルになる、誤差が少なくなるなどの利点がある。例えば、図８に示すように物質表面の片側の方向にだけ直線を伸ばす場合、ローカル座標系を使用すれば、極角θlocalの範囲を０〜πから０〜π／２に変更するだけでこれに対処することができる。よって、この場合にはローカル座標系を使用することで計算がシンプルになり、その結果、誤差も少なくなる。一方、グローバル座標系を使用することには、計算要素ごとの座標系の違いを考慮する必要がなくなるという利点がある。
【００７２】
（４．３）偏角の分割数
また、本実施形態では、偏角φlocalの分割数Ｏは極角θlocalによらず一定としたが、代わりに、偏角φlocalの分割数Ｏは極角θlocalに応じて変化させてもよい。すなわち、偏角φlocalの分割数Ｏは、極角θlocalの刻みｍに依存する変数としてもよい。
【００７３】
よって、極角θlocal(ｍ)における偏角φlocalの刻みｏ、分割数Ｏをそれぞれｏm、Ｏmで表すと、式（１５）、式（１９）、式（２０）はそれぞれ以下の式（２３）、式（２４）、式（２５）のように書き換えられる。
【数２３】
0005918630.tif 000024
【数２４】
0005918630.tif 000025
【数２５】
0005918630.tif 000026
【００７４】
前述または後述のような第１実施形態の効果は、このような偏角φlocalの分割方法を用いる場合にも享受することができる。なお、このような分割方法は、ローカル座標系だけでなく、グローバル座標系にも適用可能である。
【００７５】
（５）第１実施形態の効果
最後に、第１実施形態の効果について説明する。
【００７６】
以上のように、本実施形態では、各計算要素から複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、この判定結果に基づいて直接フラックスと形態係数を計算する。さらには、この判定結果に基づいて可視係数を計算する。
【００７７】
よって、本実施形態によれば、直接フラックスや形態係数の計算時間を、計算要素数に比例する時間に抑えることができる。よって、本実施形態によれば、間接フラックスの計算時間に影響を与える形態係数の計算時間を短縮することで、物質表面に直接的、間接的に到達する反応種を考慮した形状シミュレーションを高速化することが可能となる。
【００７８】
図１７、図１８はそれぞれ、比較例、第１実施形態における計算時間の例を示したグラフである。比較例では、従来の方法を用いて直接フラックス、可視係数、形態係数を計算した。図１７、図１８は、直接フラックスの計算時間、可視計算（可視係数と形態係数の計算）の計算時間、化学反応収束計算の計算時間、全計算時間の合計を示している。
【００７９】
図１７、図１８に示すように、第１実施形態によれば、全計算時間を比較例よりも短縮することが可能となる。これらの比較結果を図１９に示す。図１９は、計算要素数が４万の場合の第１実施形態と比較例の計算時間を比較したグラフである。
【００８０】
図２０は、第１実施形態と比較例におけるθ分割数と計算誤差との関係を示したグラフである。ただし、図２０における計算では、ローカル座標系を使用した。また、ｎ＝１とｎ＝１０００のグラフは、比較例による計算結果を示す。図２０に示すように、第１実施形態によれば、計算誤差を抑制できることが分かる。
【００８１】
なお、第１実施形態の形状シミュレーション方法は、どのような情報処理装置を用いて実行してもよい。第２実施形態では、このような情報処理装置の一例である形状シミュレーション装置について説明する。
【００８２】
（第２実施形態）
図２１は、第２実施形態の形状シミュレーション装置の構成を示す外観図である。
【００８３】
図２１の形状シミュレーション装置は、制御部１１と、表示部１２と、入力部１３とを備えている。
【００８４】
制御部１１は、形状シミュレーション装置の動作を制御するモジュールである。制御部１１は、例えば、第１実施形態の形状シミュレーション方法を実行する。制御部１１の詳細については後述する。
【００８５】
表示部１２は、液晶モニタなどの表示デバイスを有している。表示部１２は、例えば、形状シミュレーション用の設定情報の入力画面や、形状シミュレーションの計算結果などを表示する。
【００８６】
入力部１３は、キーボード１３ａやマウス１３ｂなどの入力デバイスを有している。入力部１３は、例えば、形状シミュレーション用の設定情報の入力用に使用される。設定情報の例としては、計算式に関する情報、実験値や予測値に関する情報、物質の構造に関する情報、フラックスに関する情報、形状シミュレーションの条件や手順に関する指示情報などが挙げられる。
【００８７】
図２２は、図２１の制御部１１の構成を示すブロック図である。
【００８８】
制御部１１は、ＣＰＵ（Central Processing Unit）２１と、ＲＯＭ（Read Only Memory）２２と、ＲＡＭ（Random Access Memory）２３と、ＨＤＤ（Hard Disk Drive）２４と、ＣＤ（Compact Disc）ドライブなどのメモリドライブ２５と、メモリポートやメモリスロットなどのメモリＩ／Ｆ（interface）２６とを備えている。
【００８９】
本実施形態では、第１実施形態の形状シミュレーション方法用のプログラムである形状シミュレーションプログラムが、ＲＯＭ２２またはＨＤＤ２４内に格納されている。入力部１３から所定の指示情報が入力されると、ＣＰＵ２１は、ＲＯＭ２２またはＨＤＤ２４からプログラムを読み出し、読み出したプログラムをＲＡＭ２３に展開し、このプログラムにより形状シミュレーションを実行する。この処理の際に生じる各種データは、ＲＡＭ２３内に保持される。
【００９０】
なお、本実施形態では、形状シミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体を用意し、この記録媒体からＲＯＭ２２やＨＤＤ２４内に形状シミュレーションプログラムをインストールしてもよい。このような記録媒体の例としては、ＣＤ−ＲＯＭやＤＶＤ（Digital Versatile Disk）−ＲＯＭなどが挙げられる。
【００９１】
また、本実施形態では、形状シミュレーションプログラムを、インターネットなどのネットワーク経由でダウンロードすることで、ＲＯＭ２２やＨＤＤ２４内にインストールしてもよい。
【００９２】
以上のように、本実施形態によれば、第１実施形態の形状シミュレーション方法を実行するための形状シミュレーション装置や形状シミュレーションプログラムを提供することが可能となる。
【００９３】
なお、第１及び第２実施形態では、形状シミュレーションの適用対象の例として半導体デバイスを取り上げたが、この形状シミュレーションは半導体デバイス以外のデバイスにも適用可能である。このようなデバイスの例としては、ＭＥＭＳ（Micro Electro Mechanical Systems）デバイスやディスプレイデバイスなどが挙げられる。
【００９４】
以上、第１及び第２実施形態について説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することを意図したものではない。これらの実施形態は、その他の様々な形態で実施することができる。また、これらの実施形態に対し、発明の要旨を逸脱しない範囲内で、種々の省略、置換、変更を行うことにより、様々な変形例を得ることもできる。これらの形態や変形例は、発明の範囲や要旨に含まれており、特許請求の範囲及びこれに均等な範囲には、これらの形態や変形例が含まれる。
【符号の説明】
【００９５】
１：シリコン基板、２：シリコン窒化膜、３：シリコン酸化膜、４：貫通孔、
１１：制御部、１２：表示部、１３：入力部、
２１：ＣＰＵ、２２：ＲＯＭ、２３：ＲＡＭ、２４：ＨＤＤ、
２５：メモリドライブ、２６：メモリＩ／Ｆ
【図１】
0005918630.tif 000027
【図２】
0005918630.tif 000028
【図３】
0005918630.tif 000029
【図４】
0005918630.tif 000030
【図５】
0005918630.tif 000031
【図６】
0005918630.tif 000032
【図７】
0005918630.tif 000033
【図８】
0005918630.tif 000034
【図９】
0005918630.tif 000035
【図１０】
0005918630.tif 000036
【図１１】
0005918630.tif 000037
【図１２】
0005918630.tif 000038
【図１３】
0005918630.tif 000039
【図１４】
0005918630.tif 000040
【図１５】
0005918630.tif 000041
【図１６】
0005918630.tif 000042
【図１７】
0005918630.tif 000043
【図１８】
0005918630.tif 000044
【図１９】
0005918630.tif 000045
【図２０】
0005918630.tif 000046
【図２１】
0005918630.tif 000047
【図２２】
0005918630.tif 000048