JP4516030B2 - 半導体装置の製造プロセスの設計方法、半導体装置の製造方法、及び半導体装置の製造システム - Google Patents

Description

本発明は、半導体装置の製造方法、製造装置、シミュレーション方法、及びシミュレータに係り、特に、枚葉式又はクラスタ化した半導体装置製造装置における限られた量又は質的内容の“その場”測定値を最大限に拡張し、テストピースを用いることなく、半導体装置製造プロセスを、所望の工程通りに又は修正しながら進行することを可能とする半導体装置の製造方法、製造装置、シミュレーション方法、及びシミュレータに関する。

従来、半導体装置の製造装置は、次のような変遷を経て現在に至っている。これを各世代の流れのなかで見てみると、ｄＲＡＭの世代とその世代毎のチップコストは、６４Ｋから２５６Ｍまで容量が増加するとともに、図１２８に示すように変化する。なお、図１２８の縦軸は、１Ｍにおけるチップコストを１として表現している。また、図１２８において、プロットされた点はチップコストを示し、横軸はｄＲＡＭの世代を示している。

図１２８に示すように、６４Ｋから２５６Ｋまでは、比較的工程数も安定しており、チップコストも安定している。それが、４Ｍから世代を経る毎に、工程数が増えており、それに伴ってチップコストも上昇している。この工程数の増大によるコストの増加を少しでも低減させるため、図１２８にも示したように、例えば、６４Ｍでは、（１）ムダの排除、（２）作業の標準化、（３）チップサイズの縮小、（４）微細化推進、（５）生産性向上、（６）ＦＡ化推進、（７）大口径化等の手段を講じて、コスト低減をはかっている。

２５６Ｍになると、新材料の適用や、高誘電体の採用、新構造の適用、立体化等を図ることにより、バイルールに乗せてゆくのである。しかし、１Ｇクラスになると、状況は一転し、バイルールに従うかどうか不明となる。これまでの２５６Ｍクラスまでの製造装置個々の高性能化を極限まで追求して行く従来の方法では、もはや事業としての解を求めることは困難となる。製造装置個々の高性能化を極限まで追求して行くことは、即ちコスト高を招いてしまう。

また、図１２９に示すように、各工程毎に、高性能化された装置が、例えば、工程ＡからＢ→Ｃ→Ｄ→Ｅ→Ｆの順序で使用されて行くものとする。この場合、設備領域の増大もコスト高を招くことは明らかである。

また、１Ｇクラス以上で、もう一つ忘れてはならないことは、各工程において物理的必然性による精密制御が必須になってくることである。例えば、汚染の遮断、自然酸化膜の排除、原子レベルでの平坦化の追求等である。現在では、これらのことを考慮して、研究段階ではあるが新しい概念として、図１３０に示すようなクラスタ化した製造装置が提案されている。このクラスタ化した製造装置（クラスタツール）は、バッチ式もあるが、主として枚葉式であり、また、図１３０からも分かるように、例えば中央に搬送系があり、各工程は、それぞれＡからＢ→Ｃ→Ｄの如く、矢印に示すように順次進行して行くものである。このとき、各工程及び各工程内では、多くの場合高真空に保たれ、自然酸化膜の被着も減少でき、また、汚染の遮断、原子レベルでの表面制御ができるようになってきた。

一方、半導体上に形成した電界効果ＭＩＳ素子について、昨今の高速化及び高集積化の技術開発の中で、勢い、強誘電体は薄膜化に向かっている。強誘電体膜を薄膜化すると、勿論閾値電圧は浅くなり、この分だけ主動作速度は早くなり、これにより特にＡＣ特性は格段に改良させる。また、ＥＥＰＲＯＭ等を考えてみると、微細化されてくると、非常に過酷な使用状況になってくる。このような場合には従来の通常の製法の強誘電体膜では、もはや、充分な信頼性がえられない。

しかし、ＭＩＳ素子が微細化されても、強誘電体膜自体の特性の改善策は余り省みられず、電源電圧がさほど下がらない状態では、特に主動作時にゲート強誘電体膜に高電界が掛かる。また、チャネル領域からインパクトイオン化等で発生した電子正孔は、それぞれゲート電極の極性及びドレイン電圧等の境界条件により、強誘電体膜中に注入される。そしてこれらのキャリアは、強誘電体膜中にトラップされ長期信頼性を低下させるだけでなく、耐圧低下等を引き起こす。

もう少し原子レベルでみてみる。熱強誘電体膜においては、例えばシリコン強誘電体膜に高電界を印加すると、ネットワークを構成するＳｉ−Ｏ結合が、外部から印加された高電界と相互作用する。そして、結合が切断されて、電子や正孔が捕獲される捕獲中心が形成され、続いて通過する電子や正孔は、この捕獲中心に捕獲されて、膜厚方向の電界強度分布が局部的に平均電界より高くなり、やがて絶縁破壊に至る。

これらの状況を改善しようとして、最近では、強誘電体膜を単結晶化するアイデアが学会レベルではあるが提唱され出した。例えば、シリコン（１１１）面上には、酸化セリウムＣeＯ₂が単結晶成長することが、J.Appl.Phys.,Aol.69(12),p8313(1991)に報告されている。あるいは、弗化カルシウムＣaＦ２がシリコン単結晶上に単結晶成長することが、Japan.J.Appl.Phys.Suppl.,Aol.21-1,p187(1982)に報告されている。これらは、一部であり、まだアイデアの状態で有るばかりではなく、その指針のもとになる計算手法にも多々疑問が残っている。それで、一つには、例えば強誘電体膜を例に上げても、それを用いた強誘電体膜の単結晶の姿自体が充分認識されていない状況にある。この単結晶強誘電体膜の構造に付いて報告したものを見てみよう。一つには、A. Miyamot, K.Takeichi, T.Hattori, M.Kubo, and T.Inui, "Mechanism of Layer-Layer Homoepitaxial Growth of SrTiO3(100) as Investigated by Molecular Dynamics and Computer Graphics". Jpn. J. Appl. Phys. Vol 31, (1992) 4463-4464.がある。ここでは、ペロブスカイトを単結晶強誘電体膜として、それを下地Ｓｉ基板と接着させるとき、安定構造を計算している。しかし、Ｓｉ−Ｏ−Ｔｉ叉は、Ｔｉ−Ｓｉ−Ｓｒ間の角度や距離の初期配置を誤って用いている。また、これらは、単結晶強誘電体膜の構造を単純化しており、必ずしもその指針は正しくはない。

これらのことから分かる様に単結晶強誘電体膜をどの様に設計するかに至っては殆ど正しい議論できていないのが現状である。また他方、強誘電体膜そのものについては、その特性の向上に関しては、形成プロセス上のリファインに未だ終始しているのが現状である。例えば、出来るだけ清浄な面を予め用意するとか、或いは、膜生成の為のスパッタ温度やスパッタ雰囲気の議論に終始している。これらは、単結晶強誘電体膜や強誘電体膜に関する構造を含めての認識がまだ極めて低いと言わざるを得ない。また強誘電体の自発分極発現の理由や、その結晶構造との関連が分かっていないので、現象論的追求にしか過ぎない。

しかし、以上のような新しいクラスタツールにおいて、現在直面している最も大きな困難は、その状況観察手法である。即ち、単一行程を行う装置の場合では、まずテストピースと称するウエハを、各工程に導入していた。このテストピースを用いて、たとえば、膜厚や、屈折率等を測定し、所定の条件通り工程が進行しているかどうかをテストしていた。

かかるテストピースの使用による各クラスタツール内の条件のモニタリングには、量的、質的理解上において問題がある。例えば、酸化膜の成膜工程又は成長工程においては、ＸＰＳやＦＴ−ＩＲの装置を用いるのが通例である。これをｉｎ−ｓｉｔｕで測定していたとしても、そのデ−タの意味するところは、現状では十分に検討出来ていない。せいぜいそのスペクトルのピーク波数位置から判断して、所定の物質が生成しているかどうかという判定程度である。このような状況では、研究者の判断に頼っているところが多い。また、定量的な対応が出来ず、更に、迅速な対応も出来ない。また、研究者の経験に頼ることのできない新しい未知の物質があるときは、全く判断の仕様がない。

又、一方で分子動力学シミュレータの利用も提案されている。しかし、一般に、分子動力学シミュレータは、プロセス変動算出プログラムに比べ、（１）非常に計算が遅く、また（２）対象とする物質の実行領域が、非常に小さいと言われている。それで、分子動力学シミュレータと、プロセス変動算出プログラムを直接つなげると言うことは、殆ど考えられないことであった。即ち、分子動力学部分は、極めて計算が遅いので、合体化させても、結局は、分子動力学部分の遅いスピードに全体が引っ張れ、実用に向かないと思われていた。また、分子動力学部分では、計算領域が極めて小さく、せいぜい数十オングストローム程度であり、一方、プロセス変動算出プログラムでは、数ミクロン領域を計算しているため、計算領域のずれの点からもこれらを直接つなげることは難しい。

たとえば、本発明者等は、単結晶強誘電体膜の設置条件として、界面を含む系の全体の自由エネルギに着目し、印加条件下でこのエネルギが、最小になるように、位置を決定し、これにより、高信頼性のＭＩＳ素子を作るとともに、界面に準位を形成しないようにした。特に、系の全体のエネルギを求めるにあたり、単結晶強誘電体膜の構造を確実に再現することと、これらの構成している各原子がどの様に運動するのかを厳密にしかも正確に掴む必要があったが、本発明者等は、ab-initio分子動力学を完成させ、さらに、分子動力学においては独自のポテンシャル積算法を確立し、初めてこれらの計算を可能にした。しかも、これを用いて、電気的特性等の向上点を十分保証する意味から、単結晶強誘電体膜の品質の度合いまでを明らかにしようと言うものである。

例えば、トンネル絶縁膜に、高電界を印加すると、トンネル絶縁膜の基本構造を構成している結合（例えば、シリコン強誘電体膜の場合にはＳｉ−Ｏ結合）と相互作用し、結合が切断されて、切断された結合に、続いてトンネル絶縁膜を通過する電子や正孔が捕獲されて、膜中の電荷分布に応じて、局所的に電界が平均電界より大きくなり、劣化が更に進行する。劣化の進行を抑えるには、本発明者によれば、この電界による絶縁膜の結合（例えばＳｉ−Ｏ結合）の切断を抑えることである。

電界と絶縁膜の結合との相互作用は、電界の向きと絶縁膜の結合の向きに依存する。したがって、相互作用を弱くするには、絶縁膜の結合の向きが相互作用が強い向きのものをより少なくなるようにすることである。非晶質のように、全くランダムとした場合、絶縁膜の結合の一部は必ず電界との相互作用が強い向きとなる。

本発明の目的は、強誘電体絶縁膜の結合の向きを制御し、電界との相互作用を抑え、電界による誘電性の劣化の抑制を図ることを目的としている。

これは、高信頼誘電体膜作成技術に関するもので、本発明者等が、新しいａｂｉｎｉｔｉｏ分子動力学システムによる計算結果と独自の実験データを基に生み出したものである。単結晶強誘電体膜は、これを強誘電体膜等に適用した場合、材料そのものの選択の他、基板や基板電極との位置関係等を最適化しないと、かえって非晶質誘電体膜よりも悪くなり得る。この時の材料及び結晶方位のクライテリアを示そうと言うものである。本発明者による新しい指針としては、あえて一言で記せば「主印加電界下で系の自由エネルギが最小になる様に結晶配置を行うこと等」である。また実際上の問題として「単結晶強誘電体膜に存在し得る結晶欠陥濃度の許容範囲」等をも記している。さらに、「主印加時の系の自由エネルギが最低になる配置のほかに、劣化のトリガを抑止するため、局所的にみた場合の自由エネルギ分布の凸凹の許容範囲等」をも記している。

単結晶強誘電体膜一つ取り上げても、４回対称のペロブスカイトなど非常に複雑な構造を持っていることや、基板との接合面の組み合わせによって大きく特性がことなる事など、大変複雑であるが、本発明者等によれば、結晶学的には「損」をしても主動作電圧下で「得」をしようとするである。

又、本発明は、クラスタ化におけるモニタ機能の問題点を解決し、半導体装置製造装置、特に枚葉式あるいはクラスタ化製造装置における限られた量又は質的内容の“その場”測定値を最大限に拡張し、半導体装置製造プロセスを、テストピースなしに、所望の工程通り又は修正しながら進行することを可能とする半導体装置の製造方法、製造装置、シミュレーション方法、及びシミュレータを提供することを目的とする。

本発明者らは、ａｂｉｎｉｔｉｏ分子動力学を用い、かつ新たな変分計算法を開発し、（ｉ）各プロセス要素の揺らぎや、各要素のずれの重畳現象も診断し、（ii）所望のシーケンス通り実行されているかのチェックや、（iii)変動因子の発見も可能とするシミュレーションシステムを開発した。かかるシステムのさらなる特徴は、（vi）プロセスの変動を数学的・統計的に、しかも効率よく取り扱えるようにした点と、（ｖ）新しい推論エンジンを付加した点と、（ｖi)新しい逆方向システムを付加した点にある。

即ち、本発明は、複数の工程からなる半導体装置の製造方法において、前記複数の工程の少なくとも１つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも１つに関する複数の特性の少なくとも１つを測定して実観測データを得る工程と、ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、前記複数の工程の少なくとも１つにおける、複数の製造工程因子の少なくとも１つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも１つの予測データを得る工程と、前記予測データと実観測データとを逐次、実時間で比較検定する工程と、前記比較検定により、前記複数の製造工程因子の少なくとも１つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも１つとの間に有意差が認められた場合、前記製造工程因子を逐次実時間で修正処理する工程とを具備することを特徴とする半導体装置の製造方法を提供する。

本発明は、上述の半導体装置の製造方法（請求項１）において、前記ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータは、前記複数の工程において形成又は処理される半導体材料、絶縁材料、又は導電性材料を構成する原子の運動が平衡に達したところで、時刻を０に設定し、その後、時刻ｔまでをサンプリング時間とし、各原子の分極率、電荷、及び位置ベクトルとの時間相関を求め、所定時間における光学的スペクトルを求める機能を有し、前記ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータから導かれる予測データは、個々の原子の時々刻々の運動から導かれる振動挙動、その固有振動数、及び固有振動数の分布情報を含むことを特徴とする半導体装置の製造方法を提供する。

本発明は、上述の半導体装置の製造方法（請求項１）において、前記ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータは、前記複数の工程において形成又は処理される半導体材料、絶縁材料、又は導電性材料を構成する原子の運動から、該材料の材料強度定数を逐次算出し、これを用いて該材料に誘起される応力を求める機能を有し、該応力が規定値を越える場合には、前記製造工程因子の少なくとも１つが逐次修正処理されることを特徴とする半導体装置の製造方法を提供する。

本発明は、上述の半導体装置の製造方法において、前記半導体材料、絶縁材料、又は導電性材料を構成する原子は、シリコン原子及び酸素原子である半導体装置の製造方法を提供する。

本発明は、複数の工程を実施する半導体装置の製造装置において、前記複数の工程の少なくとも１つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも１つに関する複数の特性の少なくとも１つを測定して実観測データを得る手段と、ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、前記複数の工程の少なくとも１つにおける、複数の製造工程因子の少なくとも１つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも１つの予測データを得る手段と、前記予測データと実観測データとを逐次、実時間で比較検定する手段と、前記比較検定により、前記複数の製造工程因子の少なくとも１つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも１つとの間に有意差が認められた場合、前記製造工程因子を逐次実時間で修正処理する手段とを具備することを特徴とする半導体装置の製造装置を提供する。

本発明は、ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、半導体装置の製造プロセスの複数の工程の少なくとも１つにおける、複数の製造工程因子の少なくとも１つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の工程の少なくとも１つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも１つに関する複数の特性の少なくとも１つの予測データを得る工程と、この予測データと前記複数の工程の少なくとも１つにおける前記複数の特性の少なくとも１つを測定して得られた実観測データとを逐次、実時間で比較検定して、前記製造工程因子の診断及び判定を行う工程とを具備するシミュレーション方法を提供する。

本発明は、上述のシミュレーション方法（請求項６）において、前記比較検定により前記複数の製造工程因子の少なくとも１つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも１つの間に有意差が認められた場合、前記製造工程因子の修正値を実時間で求める工程を具備することを特徴とするシミュレーション方法を提供する。

本発明は、半導体装置の製造プロセスの複数の工程の少なくとも１つにおける製造装置内、半導体基板表面、及び／又は半導体基板上に形成された膜に関する複数の特性の少なくとも１つを測定して実観測データを得る手段と、ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学ミュレータにより、前記複数の工程の少なくとも１つにおける、複数の製造工程因子の少なくとも１つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも１つの予測データを得る手段と、前記予測データと実観測データとを逐次、実時間で比較検定する手段とを具備する、半導体装置の製造工程因子を診断・判定するシミュレータを提供する。

本発明は、ａｂｉｎｉｔｉｏ分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、半導体装置の製造プロセスの複数の工程の少なくとも１つにおける、複数の製造工程因子の少なくとも１つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の工程の少なくとも１つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも１つに関する複数の特性の少なくとも１つの予測データを得る手段と、この予測データと前記複数の工程の少なくとも１つにおける前記複数の特性の少なくとも１つを測定して得られた実観測データとを逐次、実時間で比較検定して、前記製造工程因子の診断及び判定を行う手段とを具備する、半導体装置の製造工程因子を診断・判定するシミュレータを提供する。

一方、本発明者等は分子動力学シミュレータについて、前記したような非常に計算が遅いという問題、及び対象とする物質の実行領域が、非常に小さいという問題に対してそれぞれ新しい考えを導入して克服し、分子動力学シミュレータとプロセス変動算出プログラムとを新しい方法で合体させ、しかも、この合体により、今までにない新しい機能を引き出すことが出来た。以下に順に説明する。

まず（１）の計算速度の問題に関して少し詳しく触れてみる。分子動力学では、原子一個一個について、いわばニュートン力学に従って、運動方程式を立てて、これを、着目している１００個または１０００個の原子全部に渡って連立させて解くものである。

この運動方程式の基になる力の計算が一番複雑である。これが、従来の計算手法では遅かったゆえんである。各原子は、その廻りの原子から力を受けている。

近くの原子からも、遠くの原子からも、力を貰っている。これを考慮する必要がある。一般にはテレソフの式に従って導出されるポテンシャルエネルギＶ_ijをもとめ、これを基に、力を算出して行く。しかし、従来では、概ね、ポテンシャルを、距離や方向等をもとに、予め各原子に加わる力について方向毎のテーブルを作っておく手法を取っていた。この方法では、上記の１００個または１０００個の原子について、毎回、テーブルを参照し、ポテンシャルの大きさを見て行くものである。これは、非常に時間の掛かるものである。またポテンシャルそのものからは、力を算出することは出来なくて、ポテンシャルの大きさのそれぞれの方向に対する傾きが、それぞれの方向の力に成るわけである（図２参照）。これゆえ、非常に計算機の時間がかかることが分かる。

なぜこの様に、数値テーブルを作るかと言うと、ポテンシャルの形状を表した式が非常に複雑であるからである。上にも述べた様に、ポテンシャル（Ｖ）は、廻りの原子の位置はお互いの角度（θ）、さらには、距離（ｒ）等が複雑にからみ合っている。これ故に、これを、数学的に、微分して関数を作るのは非常に難しい。例えばθはｒの関数であるのに、また廻り回って再びｒの関数に成ったりする。

本発明者らは、この様な一般に非常に複雑で、上記の様にポテンシャルの表現式の引数が、巡回する様な場合について、従来、偏微分が試みられなかった部分において、新しい高階巡回偏微分式を作成した。これにより、式自体の変形は一見非常に複雑になるが、計算は非常にはやく出来るようになった。従来のように、数表を参照しなくて良いようになった。

では、この高階巡回偏微分式の一番の概念を以下に示す。これは、一般に３個の原子を想定して、考えている。すなわち、原子３個に関して一般的に記述すれば、あとは、この手法を広げて行けばよい。３個の原子があれば、各原子との距離や、その着目している原子の間の角度等も算出できる。これらにより、ポテンシャルを一義的に求められ、これをもとに、何れの方向にも、ポテンシャルの「傾き」を容易に算出することが出来る。本発明者等は、巡回して用いている変数について全て偏微分式を作っておき、これらを、連立させて解くことで、正しく、微分型が作られていることを確認した。

以下本発明によるシステムの概要を具体的に説明する。最初に、図１を参照して、本発明の応用の実際を説明する。即ち、プロセス装置Ｖは、実際に稼働している半導体製造システムの一部である。ここでは、スペクトルメータＶＩによって光の吸収スペクトルをみている。

スペクトルメータＶＩによって得られた光の吸収スペクトルは、一旦ファイル４として格納され、スケジューラＶＩＩに送られる。スケジューラＶＩＩには、後述する本発明による予測値も入力され、これらは対比診断検定部ＩＶで比較される。

対比診断検定部ＩＶでの比較結果が、許容値を超えていれば、スペックを満たさないものとしてラインから外される。また、満たしていれば、その変動分を推論エンジンとしてのプロセス変動算出ユニットＩに送られる。

プロセス変動算出ユニット１では、基本的には原子レベルまで考慮しないマクロ的な性質を計算し、ファイル２に格納する。ファイル２の内容は、分子動力学シミュレータＩＩで原子間位置ベクトルと結合角が算出される。

分子動力学シミュレータＩＩで算出された原子間位置ベクトルと結合角は、膜特性算出ユニットＩＩＩに入力され、ここで基本膜特性が算出されてファイル３に戻される。膜特性算出ユニットＩＩＩで算出された基本膜特性は、プロセス変動算出ユニットＩにも戻され、パラメータの一部とされる。

このように、本発明によって、半導体装置製造プロセスを、テストピースなしに、所望の工程どおり、または修正しながら進行する事を可能とする半導体装置の製造方法が実現される。

実際には、分子動力学シミュレータＩＩの内容が特に重要である。具体的には、図７のブロック図に示したように、２つの部分からなっている。第１の部分ＩＩａは、本件発明者によって導入された高階巡回偏微分法を利用するものであり、第２の部分ＩＩｂは、三次以上の高次項を考慮した解析を行うものである。なお、以後、高階巡回偏微分法や三次以上の高次項を考慮した解析の詳細が順次説明される。

即ち、本発明者らは、純粋数学理論にたち帰り、新たな変分計算法を確立するとともに、新しい計算化学理論により、各種スペクトルをも予測する技術を開発した。以下、本発明における新しいａｂｉｎｉｔｉｏ分子運動学について説明する。

本発明者等が開発した新しいシミュレーションの全体構成を図３に示す。分子動力学部分（ＭＤ）と密度汎関数（ＤＦＴ）とを合成した新しい手法を提案した。まず本発明者らによる分子動力学法と密度汎関数法との使い分けを図３を用いながら以下に示す。

図３に示すように、分子動力学（ＭＤ）部分はフローチャートの下半分を占め、密度汎関数部分（ＤＦＴ）は上半分を占めている。分子動力学法は零Ｋでない所謂有限温度でのやや大きな系を扱い、原子（イオン）間のポテンシャルには、予め第一原理から求めたポテンシャルを用いる。他方、密度汎関数（ＤＦＴ）は、小さな系の基底状態を求める時に使うことにした。そして、当然ＤＦＴの部分は温度は零Ｋとした。ポテンシャル部分は、電子系の基底状態から求めた。具体的には波動関数（ＫＳ方程式）を解き、局所密度汎関数を併用した。

ＤＦＴ部分で、縮重のない電子系の基底状態の全エネルギーを求めて置く。全エネルギーは（波動関数までさかのぼらずに）電子密度ｎ（ｒ）の汎関数Ｅ［ｎ（ｒ）］である。

即ち、Ｅ［ｎ（ｒ）］は、下記式（１２０１）により表される。

右辺第１項のＶz （ｒ）は核のクーロン場などの外場、第２項は電子間のクーロン相互作用エネルギー、第３項は、電子の１体近似を想定して、電子間相互作用のない時の運動エネルギーＴｓ[ｎ（ｒ）]であり、下記式（１２０２）により表される。

第４項は交換相関エネルギーＥｘｃ[ｎ（ｒ）]で、その形は局所密度近似による。電子密度ｎ（ｒ）は下記式（１２０３）により表される。

基底状態は、式（１２０１）のＥ［ｎ（ｒ）］の変分を０とおいて得る、下記式（２１），（２２）に示すＫＳ方程式（Ｋｏｈｎ−ｓｈａｍ ｅｑｕａｔｉｏｎ）を解いて、固有値Ｅの小さい方からＮ個採用する。式（１２１１）にｎ（ｒ）が入っているから、ｎ（ｒ）、ψ（ｒ）、有効ポテンシャルＶeff （ｒ）を自己無撞着（ｓｅｌｆ−ｃｏｎｓｉｓｔａｎｔ）に解かなければならない。

〔数４〕
［−∇²＋Ｖ_eff（ｒ）］Ψ（ｒ）＝ＥΨ（ｒ） …（1211）
Ｖ_eff（ｒ）＝Ｖ_Z（ｒ）＋∫２ｎ（ｒ′）／（｜ｒ−ｒ′｜）・
ｄｒ′＋δＥ_xc［ｎ（ｒ）］／δｎ（ｒ）
…（1212）
交換相関エネルギーＥxc［ｎ（ｒ）］は、局所密度近似（ＬＤＡ）で１電子当たりのエネルギーεxc（ｎ）に電子密度ｎを掛けて積分することにより、下記式（１２１３）を得た。

〔数５〕
Ｅ_xc［ｎ（ｒ） ］＝∫ε_xc（ｎ（r））ｎ（ｒ）ｄｒ …（1213）
εxc（ｎ）の形は、Ｗｉｇｎｅｒ他の式など、いろいろ提案されている。更に、内殻電子は考えずに価電子だけ扱う時には、式（１２１１）のＶz （ｒ）として核のクーロン場をそのまま使うのでなく、核の位置での特異性を消したノルム保存擬ポテンシャルでおきかえる。このようにならすことによって、フーリエ展開での成分の個数をおさえることができた。

全エネルギーとして、本発明者等は厳密を期するため、原子（イオン）間のクーロン・ポテンシャルも含むものを使った。即ち、下記式（１２２１）である。

ここで｛Ｒ_I｝はイオンの核の座標、｛αν｝は外部の制約条件（体積Ω、ひずみεＬνなど）。式（１２２１）の右辺第１項は、ＤＦＴの式（１２０１）から右辺のＴs ［ｎ（ｒ）］であり、式（１２０１）右辺の残りが式（１２２１）のＵの１部である。その部分は、ＤＦＴの場合と同様に１体近似の式（１２１１）からＶ_eff（ｒ）で表され、やはりＶz （ｒ）を擬ポテンシャルで置き換えて使う。

ここで新しい考え方を入れた。即ち、Ｅをポテンシャル・エネルギーと見なして、他に、下記式（１２２２）に示す仮想的な運動エネルギーを導入する。

従ってラグランジアンは、下記式（１２２３）により表される。

Ｌ＝Ｋ−Ｅ …（1223）
式（１２２２）中のＭ_Iは、本当の原子（イオン）の質量だがμ、μνは仮想的な質量であり、後述のように目的に応じて変えた。束縛条件として、下記式（１２２４）に示す正規直交条件
〔数８〕
∫ψ_i ^*（ｒ，ｔ）ψ_k（ｒ，ｔ）ｄｒ＝δ_ik …（1224）
を課するから、オイラー方程式を作る時に未定乗数Λ_ikを導入して、Ｌに下記式（１２３１）を加えたものの変分をとる。その結果、下記式（１２３２），（１２３３），（１２３４）が得られる。

式（１２２２）の運動エネルギーの温度Ｔが対応する。古典統計力学のエネルギー等分配則によって各自由度が温度を持っているのだから、Σ1/2 Ｍ_I Ｒ_I 2 との関係でいえば、仮想的でなく現実の温度である。またμ、μνの大小によってψi 、ανを抑制したり増強したりすることが出来る。

例えば、μ<<Ｍ_Iとして高温からＴ→０とすれば、Ｒ_I をとめたままψi が変化して電子系の基底状態を得る。この時、式（１２３１）の左辺が０になり、右辺は式（１２２１）とその下の注意により、下記式（３０）となるから、ψi の線形結合を適当に作る｛Λ_ik｝が対角線化されてＫＳ方程式（１２３５）を得る。

〔数１０〕
［∇²−Ｖ_eff（ｒ）］Ψ_j（ｒ）＋ΣΨ_k（ｒ） …（1235）
つまりＤＦＴでＫＳ方程式を解いたのと同じ結果をｓｉｍｕｌａｔｅｄ ａｎｎｅａｌｉｎｇで得たことになる。ただし、温度Ｔの状態を作るのに、モンテカルロ法によらず、仮想的な力学系の運動を追っているから、ｄｙｎａｍｉｃａｌ ｓｉｍｕｌａｔｅｄ ａｎｎｅａｌｉｎｇ（ＤＳＡ）と呼ばれる。

式（１２３１）の積分をする時、必ずしもψi をそのまま変数として扱うのでなく、それを平面波で展開した展開係数（つまりフーリエ係数）を変数として扱うことができる。この展開係数の個数をむやみにふやさないためには、ＤＦＴの場合と同様に、擬ポテンシャルで核の位置での特異性を消しておくことが必要である。特に、結晶の周期性を利用できるときには、上記の｛Λ_ik｝を対角線化するψi の線形結合をφnkと記す。ｋはブリユアン帯域内の波数ベクトル、ｎは帯域インデクスである。

φnkは、下記式（１２４１）を満たすはずである。

〔数１１〕
μＳ_nk(ｒ，ｔ)＝［∇²−Ｖ_eff(ｒ)］φ_nk(ｒ)＋λ_nkφ_nk(ｒ，ｔ)
…（1241）
λnkは行列｛Λ_ik｝の固有値であり、エネルギー準位である。φnk（ｒ，ｔ）を展開して、下記式（１２４２）が得られる。

擬ポテンシャルが局所的であるならば（この条件は一般に満たされない）、有効ポテンシャルＶ_effも展開して、下記式（１２４３）が得られる。

式（１２４２）と式（１２４３）を式（１２４１）に代入して、下記式（１２５１）が得られる。

右辺の最後の頃でＫ＝Ｋ′＋Ｋ″とおくと、Σexp ［ｉ（Ｋ＋Ｇ）ｒ］が全部の項に共通となるから外して、下記式（１２５２）を得る。

〔数１５〕
μＴ_n,K+G(ｔ)＝［−｜Ｋ＋Ｇ｜²＋λ_nk］Ｃ_n,K+G(ｔ)
−ΣＶ_G-G′(ｔ) Ｃ_n,K+G′(ｔ)
＝−［｜Ｋ＋Ｇ｜²＋Ｖ₀(ｔ)−λ_nx］Ｃ_n,K+G(ｔ)
−ΣＶ_G-G′(ｔ) Ｃ_n,K+G′(ｔ)
…（1252）
ここで周波数を、下記式（１２５３）に示すように定義し、
〔数１６〕
ω＝((｜Ｋ＋Ｇ｜²＋Ｖ₀(ｔ)−λ_nx)／μ)^1/2
…（1253）
式（１２５１）の左辺を差分になおすと、下記式（１２６１）が得られる。

〔数１７〕
Ｃ_n,K+G(t＋Δt)＝２cos(ωΔｔ）Ｃ_n,K+G（t)−Ｃ_n,K+G(t−Δｔ）
−２［１−cos(ωΔｔ）］｛Σ_G´Ｖ_G-G´(t)・Ｃ_n,K+G´(t)｝
…（1261）
この式は、振動部分の積分を解析的にすましているので、純粋に数値的な方法よりもΔｔを大きくすることにより、計算量を減少させることが出来る。

クーロンポテンシャルは、以下のようになる。即ち、本発明者は、クーロンポテンシャルを分解してゆくと４項に分かれる事を見い出した。特に、従来は３項しかなかったが、新しく第４項があることを確認した。これらの方法を順次説明する。

本発明者らは、まず厳密な式の出発として、誘電率を加味した基本式から解き始めた。まず基本式である下記式（１２６２）から解き始めた。

導体ε＝∞と真空ε＝１の場合とでクーロンポテンシャルが異なり、Ｌは（立方体の）単位結晶の一稜、Σは単位結晶内でとり、Ｚi 、ｒi は、それぞれ第ｉ粒子の荷電と位置である。これは球内の荷電によって球の内面に分極が生じることによる。導体でない球の内面に双極子の層が出来るが、上記式の最終項がちょうどその効果を打ち消す働きをする。Ｅｗａｌｄの方法は左辺、ε＝∞のものを与えるから、真空内の値を求めるには上記式の最終項を加えなければいけない。

結果だけ掲げる。

クーロン・ポテンシャルを下記式（１２７１）により表わすものとする

式（１２７１）において、Ｎは単位結晶内の原子数、その中の第１原子の位置と荷電がｒi ，Ｚi であり、ｎは単位結晶とそれを周期的にずらしたものを指定するベクトルであって、下記式（１２７２）に示すように設定した。

ｎ＝ｎ_xＺ_x＋ｎ_yＺ_y＋ｎ_zＺ_z …（1272）
と表される。ここに、Ｚx ，Ｚy ，Ｚz がそれぞれｘ，ｙ，ｚ方向の稜のベクトルで、ｎx ，ｎy ，ｎz はそれぞれ（バルク結晶の場合）−∞から＋∞にわたる整数である。Σ′の′はｎ＝０の時のｊ＝ｉを除くことを意味するものとする。

ここで本発明者らは下記式（１２７３）に示す新しいＦ関数を導入する。

上記式（１２７３）において、Ｇ（ｒ，ｔ）はｔの周期関数である。これはフーリエ級数で表現できることを見い出した。

Ｇ（ｒ，ｔ）は、更に下記（１２８１）に示すように変形することが出来る。

上記式（１２８１）において、ｍは逆格子ベクトルであり、下記式（１２８２）により表わされる。

〔数２２〕
ｍ＝ｍ_x・(１/Ｌ_x，０，０)＋ｍ_y・(０，１/Ｌ_y，０)
＋ｍ_z・(０，０，１/Ｌ_z)
…（1282）
指数ｍx ，ｍy ，ｍz はそれぞれ−∞から∞にわたる整数である。ｍ＝０の時は、exp ［…］＝１だから、ΣＺi ＝０であり、単位結晶内の総荷電が０であれば、ｍ＝０の項は消える。あらかじめｍ＝０の項を除くと、下記式（１２８３）となる。

本発明者等は積分範囲をｋで分け、Ｇ（ｒ，ｔ）の２つの形を使い分けて、下記式（１２９１）を得る。

ここで本発明者等は、公式として下記式（１２９２）を用いた。

これにより、始めのクーロンポテンシャルは、下記式（１２９３）により表わされる。

Σ′の右上の′印は、ｎ＝０の時、ｊ＝ｉを除くことを意味する。

更に、下記式（１３０１）が成立する。

この式（１３０１）はｒを含まないから、力には関係ない。この式（１３０１）は、更に下記式（１３０２）に示すように変形することが出来る。

以上の結果をまとめると、クーロンポテンシャルは、結局、下記式（１３１１）〜（１３１５）により表わされる。

ここに、ｎは、ｎ＝ｎx ・（Ｌx ，０，０）＋ｎy ・（０，Ｌy ，０）＋ｎz ・（０，０，Ｌz ）であり、、ｍは、ｍ＝ｍx ・（１／Ｌx ，０，０）＋ｍy ・（０，１／Ｌy ，０）＋ｍz ・（０，０，１／Ｌz ）である。上記の表現の具合のよい点は、もとのΣの項が逆数のオーダーでしか減衰しないのに対して、Φ(1) ではΣの項がerfcの因子によって、Φ(2) ではΣの項がexp の因子によって、急速に減衰することである。Φ(3) での減衰との遅速に逆向きに効くから、両者のバランスがとれるような適当なκを選ぶ必要がある。これらは、距離の近いところから順にクーロン力の寄与を計算した結果であり、しかも周囲が導体である場合の結果である。

周囲が真空である場合によればもう１項が加わる。これが特に従来省みられなかった部分である。

一般に分子動力学では、原子間ポテンシャルが最も大切な量である。これには、現状では、単純には、２体間でのみ記述するものも結構多い。これは、分子動力学計算が非常に、時間のかかるもので、その時間節減のために、ポテンシャルを簡便なものを使っているわけである。この様に、ポテンシャルを簡便にしてしまうと、計算自体は速くなるが、実際を反映したものではありえない。

上に記したように重要な原子間ポテンシャルは、一般には、３体で記述できる。これの意味合いは以下の通りである。即ち、第一の原子（ｉ）と、第二の原子（ｊ）との作用に、さらに第三の原子（ｋ）から第二（ｊ）を通して入ってくる効果がある。これを図示したものが図４である。図に有るように、ｒを粒子間の距離とし、θを粒子角度とする。これらの諸量の主従の関係は、図５に示す通りである。特に図５で分かるように、非常に複雑であることが分かる。一般に、ポテンシャルを微分したものが、力になるわけであるが、これで分かるように計算は、非常に複雑になる。ここでもう一度、従来の計算手法と本発明の位置つけを図６を用いて示した。

本発明者らは、図６にある高階巡回偏微分法を新しく開発して、これにより、高精度に且つ高速に演算できる様にした。本発明者等の方法がどれくらい速いかを以下に実際の比較例を入れて記す。まず、図６にある簡便な２体間ポテンシャル計算に関しては、ここでは、触れない。なぜなら、半導体ＬＳＩの物理現象を記述するには殆ど使用に耐える精度が得られないので、ここでは、議論しないことにした。一方、３体問題を入れたポテンシャルで、従来の方法として、テーブル参照方法のプログラムを、本発明者等は作成した。これによると、ｒ_ijと,θ_ijkと、ｒ_jkに関して、表を作成した。実際には、block dataとして、配列にして格納した。そして、この参照項目からずれていた時のことを考え、補完プログラムをも作成した。補完には、cubic spline法を用いた。ある値、即ち、ｒ_ijと, θ_ijkと、ｒ_jkとに関して、テーブルを引用出来るようなプログラムにした。し かし、所望の条件（ｒ_ijと,θ_ijkと、ｒ_jk）が、配列内のどの行列に入るかを探る必要がある。本発明者等は、出来るだけ高速に成るように気を配りながらプログラムを作成したが、やはり、数字の大きさによる場合分けをする部分が必要で、そこにはif then else文を用いた。この様にして、所望の条件（ｒ_ijと,θ_ijkと、ｒ_jk）がテーブルのどの枠にはまるかを、見いだすとともに、その中での補完をおも行った。これで、まず、所望の条件（ｒ_ijと,θ_ijkと、ｒ_jk）でのポテンシャルが求まったわけであるが、実は欲しいのは、この点における力である。

本発明者等は、微小変化させた部分でもう一点についても、ポテンシャルを同様の手続きを用いて求め、両者の差から、平均変化率の考えを用いて、微分、即ち力を求めた。

これで分かるように、非常に多くのif then else手続きを用いた。電子計算機では、if then else手続は比較的苦手で、遅い。事実、本発明者等が作成した、参照型プログラムを用いた場合、力を計算する部分のみについて、詳細について調べた。その結果、１６００個のＳｉ原子からなる単結晶基板で、絶対温度３００Ｋで、外圧１気圧の場合であるが、参照型プログラムでは、約３５０倍の時間を要した。これで分かるように、本発明者等が提案する、高階巡回偏微分法のほうが、有利であることがわかる。また参照型プログラムでは、原子が異なる度に、その都度データ表を作成する必要がある。

次に分子動力学から材料特性値等をどの様に抽出して、それを汎用シミュレータにどの様に入れているのか」に関しても簡単に説明する。

通常、ＬＳＩプロセスに於いては、種々様々な評価技術を使っている。赤外吸収による膜の物性値評価法、さらには、ラマン分光による膜中の応力評価、Ｘ線による膜の構造解析等々がある。またこれら光学的なものに限らず、膜の誘電率等電気的な評価技術もある。

本発明者等は、それぞれの評価手法の原理を、一旦、鋭意、理論的に個々原子の運動にまで分解し詳細に追跡した。その上で、本発明者等の開発した分子動力学を用いて、個々原子の動きを追跡し、それぞれの評価量を初めて理論的に算出し、たとえばスペクトルの様な形に直接表示することが出来た。これにより、例えば、実測値との定量的な比較が出来る等、従来にない大きな技術の進捗を得たわけである。

本発明者等が、初めて克服した理論的手段による原子運動の追跡と、具体的な測定評価量の導出までを以下に記す。幾つかのものがあるが、基本的には、本発明者が導出した基本概念は以下の通りである。即ち、着目する原子のモーメントを計算する。そして、これを、着目する時間の窓の中で積算し、しかる後に、フーリエ変換し、固有スペクトルを得るものである。これが大きな骨子である。

次に、Ｓｉを例に、経験的なポテンシャルに関して、本発明者がどの様にして、そのポテンシャルを力に変換したかを示す。これは、一例であるが、Ｓｉに限らず、イオン結合による効果を考慮しなければ、例えば酸化膜のポテンシャルも、このＳｉの問題と同様に、３体問題として、ここに述べる高階巡回偏微分を用いることが出来る。本発明者等は、Ｓｉに関してどんなポテンシャルが良いかも吟味した。

ここで用いるポテンシャルはTersoffによるものである。（Phys.Rev.Lett.56.632 (1986), Phys.Rev.B37,6991(1988))。

今、Tersoffによれば、i番目のＳｉに関する全ポテンシャルは、

で記述できる。これは３体で記述しているので、本発明者等が注意を喚起したいのは、上記式（０００１）に於いて、Ｖ_ij≠Ｖjiである。着目するＳｉ原子の位置番号をｉとし、その周辺の他の粒子番号をｊとすると、上記Ｖ_ijは、
〔数３１〕
Ｖ_ij＝ｆ_c(ｒ_ij)［ａ_ijｆ_R(ｒ_ij)＋ｂ_ijｆ_A(ｒ_ij)］ …（0002）
である。ここにｒは粒子間の距離である。また、ｆc（ｒ_ij）は、カットオフ関 数で、ｆR（ｒ_ij）は斥力を示し、ｆA（ｒ_ij）は引力を示す。ａ_ijは配位数を考慮したカットオフ係数、ｂ_ijも配位数を考慮したカットオフ係数である。配位数を表現するパラメータを用いて３体的表現を行っていることになる。これは、他のものでも基本形は同じであるので、ここでは、本発明者等が初めて開発した巡回手法の思想を交えながら以下に示す。

ｆR（ｒ_ij）とｆA（ｒ_ij）は、Morse型のポテンシャルを変形したもので、ｆR（ｒ_ij）＝Ａexp(-λ1ｒ)、ｆA（ｒ_ij）＝−Ｂexp(-λ2ｒ)である。

この内、λ1とλ2は、定数であり、その大きさは、原子間距離程度の値の逆数である。

ところで、カットオフ関数ｆc（ｒ_ij）とは、
〔数３２〕
Ｖ_ij＝ｆ_c(ｒ_ij)［ａ_ijＡexp(−λ₁ｒ_ij)−ｂ_ijＢexp(−λ₂ｒ_ij)］
ｆ_c(ｒ)＝１ （ｒ≦Ｒ−Ｄ）
ｆ_c(ｒ)＝1/2−1/2sin［π/2(ｒ−Ｒ)/Ｄ］(Ｒ−Ｄ＜ｒ＜Ｒ＋Ｄ）
ｆ_c（ｒ）＝０ （ｒ≧Ｒ＋Ｄ）
…（0003）
であり、ここに、Ｒは通常対象とする構造の第一隣接ゾーンだけを含む様にその寸法を選ぶ。その値は、本発明者の詳細な吟味によると大体２〜３Åである。次に、ｂ_ijは実配位数を示すことになるが、ここでも上記カットオフ関数を使う。

その定義は、Tersoffによれば、 ｂ_ij＝（１＋βⁿζ_ij ⁿ）^-1/2n …（0004）
である。ここに、
〔数３３〕
ζ_ij＝Σｆ_c(ｒ_ik)ｇ(θ_ijk)ｅｘｐ［λ₃ ³(ｒ_ij−ｒ_ik)³］ …（0005）
である。Σ記号はk≠i,jで回す。ここで分かるように、ζ_ijの意味は第３の原子kが入ることによる環境因子であるので、iから見た場合と、jから見た場合とで 、互いに大きさは異なる。即ち
ζ_ij≠ζjjであり、さらに、上記式（０００１）で述べた様に、Ｖ_ij≠Ｖjiである。従って、ｂ_ij≠ｂjiである。

またｇ（θ）はボンド角因子であり、
〔数３４〕
ｇ(θ)＝１−(ｃ²/ｄ²)−ｃ²/(ｄ²＋cosθ²) …（0006）
である。ここで、θは図４の様に取るものとする。θを求めるにあたり、実際の直交座標を用いて表現してみる。即ち、
〔数３５〕
ｒ_ij＝√｛(ｘ_j−ｘ_i)²＋(ｙ_j−ｙ_i)²＋(ｚ_j−ｚ_i)² ｝ …（0007）
であり、ｒ_ikも同様の手続きで求められる。そうすると、内積をＰ_ijとすると、
〔数３６〕
Ｐ_ijk＝(ｘ_j−ｘ_i)(ｘ_k−ｘ_i)＋(ｙ_j−ｙ_i)(ｙ_k−ｙ_i)
＋(ｚ_j−ｚ_i)(ｚ_k−ｚ_i) …（0008）
である。これらを用いて、
cosθ_ijk＝Ｐ_ijk/(ｒ_ij・ｒ_ik) …（0009）
となる。ここで、本発明者等が用いた上記各式における定数を示す。即ち、
〔数３７〕
Ｒ＝3.0Å，Ｄ＝0.2Å，Ａ＝3264.7ev，Ｂ＝95.373ev，Ｃ＝4.8381，
λ₁＝3.2394Å，λ₂＝1.3258Å，λ₃＝λ₂ ，
β＝0.33675，ｎ＝22.956，ｄ＝2.0417 …（0010）
である。

これらの準備を経ていよいよ本発明者の巡回部分を詳しく説明する。ここからいよいよ力を計算して行くわけである。ポテンシャルの（０００２）式を位置の座標で微分すると力になる。即ち、

がそれぞれ、粒子i、jに働く力のベクトルのｘ成分である。しかし、実際にそれを求めるのは、そう簡単ではない。本発明者等の数学理論に立ち戻った工夫と発明が以下に有るわけである。それが、高階巡回偏微分と言う分けである。

（００１１）式及び（００１２）式で述べた式を、具体的に偏微分すれば良い。

しかし、そんなに単純では無いことが直ぐ分かる。

即ち、

また、j及びｋに関する偏微分方程式の変形は、同様に以下の様になる。また 、上記は主にｘについて記述したが、ｙやｚについても行う必要がある。ここでは、特に、上記との対応が分かるように、空白部分は、空白のままにして置いた。このような複雑な変形は、本発明者等が、調べた限りでは見あたらなかった。

この様な複雑な変形を今まで克服しなっかたため、逆に、単純なテーブル参照方式に逃げたりしていたわけである。この様な発明を克服したため、本発明者等は、分子動力学と汎用流体シミュレータとを合体出来たわけである。

これらの式の各項は次のように変形できる。

〔数４１〕
∂Ｖ_ij/∂ｒ_ij
＝∂ｆ_c(ｒ_ij)/∂ｒ_ij・［Aexp(−λ₁/ｒ_ij)−ｂ_ij Bexp(−λ₂/ｒ_ij)］＋
ｆ_c(ｒ_ij)［−λ₁ Aexp(−λ₁/ｒ_ij)＋λ₂ｂ_ijBexp(-λ₂/ｒ_ij)］
＝Aexp(−λ₁/ｒ_ij)・［∂ｆ_c（ｒ_ij）/∂ｒ_ij−λ₂ｆ_c(ｒ_ij)］−
Bexp(−λ₂/ｒ_ij)・［∂ｆ_c（ｒ_ij）/∂ｒ_ij−λ₂ｆ_c(ｒ_ij)］ｂ_ij
…（0016）

∂ｆ_c(ｒ_ij)/∂ｒ_ij＝−π/４Ｄ・cos［π/２・(ｒ−Ｒ)/Ｄ］
（Ｒ−Ｄ＜ｒ＜Ｒ＋Ｄ）
∂ｆ_c(ｒ_ij)/∂ｒ_ij＝０ （ｒ≧Ｒ＋Ｄ）
…（0017）

∂Ｖ_ij/∂ζ_ij＝(∂Ｖ_ij/∂ｂ_ij)・(∂ｂ_ij/∂ζ_ij)
＝−Ｂfc(ｒ_ij)exp-(λ２ｒ_ij)(−1/2ｎ)(１＋βⁿζ_ij ⁿ)−^1/2n-1βⁿζ_ij ^n-1
＝Ｂfc(ｒ_ij)exp-(λ２ｒ_ij)ｂ_ij(βζ_ij)ⁿ/［２｛１＋(βζ_ij)ⁿ｝ζ_ij］
…（0018）

∂ζ_ij／∂ｒ_ij
＝３λ₃ ³Σｆ_c(ｒ_ik)ｇ(θ_ijk)exp(λ₃ ³(ｒ_ij−ｒ_ik)³)(ｒ_ij−ｒ_ik)²
…（0019）
〔数４２〕
∂ζ_ij／∂ｒ_ik＝ｄｆ_c(ｒ_ij)/ｄｒ_ij・ｇ(θ_ijk)exp(λ₃ ³(ｒ_ij−ｒ_ik)³)
−３λ₃ ³ｆ_c(ｒ_ik)ｆ_c(ｒ_ij)(ｒ_ij−ｒ_ik)²
…（0020）

∂ζ_ij/∂cosθ_ijk
＝ｆ_c(ｒ_ik)exp(λ₃ ³(ｒ_ij−ｒ_ik)³)ｄｇ(θ_ijk/ｄcosθ_ijk
ｆ_c(ｒ_ik)exp(λ₃ ³(ｒ_ij−ｒ_ik)³)[2c²cosθ_ijk/{ｄ²cos²θ_ijk}²]
…（0021）

∂ｒ_ij／∂ｘ_i＝（ｘ_i−ｘ_j）／ｒ_ij＝∂ｒ_ij／∂ｘ_i
…（0022）
∂ｒ_ik／∂ｘ_i＝（ｘ_i−ｘ_k）／ｒ_ik＝∂ｒ_ij／∂ｘ_k
…（0023）
〔数４３〕
∂cosθ_ijk/∂ｘ_i＝１/(ｒ_ijｒ_ik)・∂Ｐ_ijk/∂ｘ_i ⁺
Ｐ_ijk{1/(ｒ_ik)∂/∂ｘ_i[１/(ｒ_ik)]⁺
１/(ｒ_ij)∂/∂ｘ_i[１/(ｒ_ik)]｝
＝１/(ｒ_ijｒ_ik)・(ｘ_i−ｘ_k＋ｘ_i−ｘ_j)−
Ｐ_ijk{(ｘ_iｘ_j)/(ｒ_ikｒ_ij ³)＋(ｘ_i−ｘ_k)/(ｒ_ikｒ_ij ³)
＝１/(ｒ_ik)・{(ｘ_i−ｘ_j)/(ｒ_ij)−
(ｘ_i−ｘ_k)/(ｒ_ik)cosθ_ij}
…（0024）

∂cosθ_ijk/∂ｘ_i＝−１/(ｒ_ij)・{(ｘ_i−ｘ_k)/(ｒ_ik)
−(ｘ_i−ｘ_j)/(ｒ_ik)cosθ_ij}
…（0025）

∂cosθ_ijk/∂ｘ_k＝−１/(ｒ_ik)・{(ｘ_i−ｘ_j)/(ｒ_ij)
−(ｘ_i−ｘ_k)/(ｒ_ik)cosθ_ij}
…（0026）

Ｌ(ｒ_i∂ｒ_i/∂ｔ，Ｖ，∂Ｖ/∂ｔ)
＝1/2Σｍ(∂ｒ_i/∂ｔ)²−Ｕ(ｒ_i)＋1/2Ｍ(∂Ｖ/∂ｔ)²−Ｐ_EＶ
(∂Ｌ(ｑ_j,ｑ′_j)/∂ｑ_j)−ｄ(∂Ｌ/∂ｑ′_j)/ｄｔ＝０
…（0027）
図７及び図８に、分子動力学シミュレータＩＩと、上記プロセス変動算出ユニットＩとの連携を説明する。分子動力学シミュレータＩＩには、第一の入力ａとしては、図８に示したファイル（２２）、（３２）、（４２）、（５２）、（６２）、（７２）等から転送されるものがある。内容は、境界条件や、領域の大きさ、応力、歪み、温度などである。また必要に応じて第二のものとして補助的な入力ｂがあり、分子動力学の計算時間やアンサンブルの指定等がある。

他方、第一の出力としては、位置ベクトルｒ（ｔ）が出てくる。そして、図８に示したファイル（２３）、（３３）、（４３）、（５３）、（６３）、（７３）、（８３）等へ転送される。叉、膜特性算出プログラムＩＩＩにより、光学的・電気的・結晶学的特性等の膜特性を算出し、図８のファイル（３）へ転送される。

また、図９に、汎用２・３次元プロセスシミュレータ（図８の（Ｉａ））の入力シーケンスの一例を示す。また、図１０に、同じく汎用２・３次元プロセスシミュレータ（図８の（１０））のもうひとつの入力シーケンスを示す。図９に示した内容は、例えば、それぞれ用いるクラスタツールをも指定するとともに、例えばその中にある酸化工程の詳細には、ガス量や、昇温シーケンス、さらには、バルブシーケンス等がある。

他方、図１０の入力シーケンスは、図８で示した未知材料を含む場合である。

この図１０の場合は、例えば、強誘電体膜であり、膜をスパッタする工程を含んでいる。このとき、強誘電体膜の詳細データがないので、分子動力学シミュレータＩＩに、一旦、飛び、そこで計算を行われる。

又、本発明者らは、経験的ポテンシャルを与えた分子動力学シミュレーションを行い、(1)各プロセス要素の揺らぎや各要素のずれの重畳現象も診断し、(2)所望のシーケンス通り実行されているかのチェックや、(3)変動因子の発見も可能 とするシミュレーションシステムを開発した。

具体的には、単結晶強誘電体膜とＳｉ基板との位置関係や、金属電極との位置関連を提案するとともに、さらに、単結晶強誘電体膜として、実際上どの程度の品質のものを用意すれば良いかについても、初めて指針を示すことが出来た。かかる単結晶強誘電体膜構造を厳密に独自の計算手法により、計算機を用いて再現するとともに、特性評価関数を付与し、これに基ずき、印加条件を考えて基板単結晶との位置関係をどの様にすれば良いかと同時にその欠陥密度の許容範囲の指針を示すものである。

すなわち、本発明は、複数の工程からなる半導体装置の製造方法において、前記複数の工程の少なくとも一つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも一つに関する複数の特性の少なくとも一つを測定して実観測データを得る工程と、経験的ポテンシャルを与えた分子動力学シミュレータにより、前記複数の工程の少なくとも一つにおける、複数の製造工程因子の少なくとも一つの設定値、及びこの設定ちの揺らぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも一つの予測データを得る工程と、前記予測データと実観測データとを、逐次、実時間で比較検定する工程と、前記比較検定により、前記複数の製造工程因子の少なくとも一つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも一つとの間に有為差が認められた場合、前記製造工程因子を逐次実時間で修正処理する工程とを具備することを特徴とする半導体装置の製造方法を提供する。

本発明は、上述の半導体装置の製造方法において使用する経験的ポテンシャルを与えた分子動力学シミュレータは、前記複数の工程において形成または処理される半導体材料、絶縁材料、または導電性材料を構成する原子の運動が平行に達した所で、時刻を０に設定し、その後、時刻ｔまでをサンプリング時間とし、各原子の電荷及び位置ベクトルとの時間相関を求め、所定時間における光学的スペクトルを求める機能を有し、前記経験的ポテンシャルを与えた分子動力学シミュレータから導かれる予測データは、個々の原子の時時刻々の運動から導かれる振動挙動、その固有振動数、及び、固有振動数の分布情報を含む半導体装置の製造方法を提供する。

本発明は、上述の経験的ポテンシャルを与えた分子動力学シミュレータにおいて、原子の運動を計算する際に、原子の位置ベクトルと各元素の電気陰性度、もしくは各元素のイオン化ポテンシャルと電子親和力とを入力し、これらのデータから、個々原子の振動及び変位による電荷の変化を計算し、さらにこの電荷の変化を加味して前記ポテンシャルを計算する半導体装置の製造方法を提供する。

本発明者らは新しい計算化学理論により、種々元素を含む系における様々な現象を高速に且つ高精度に計算予測する方法を発明した。これによって、半導体製造プロセスに使用される様々な材料の特性が高精度に計算できるようになった。

従来、種々元素が含まれている結晶や分子の構造及び特性を予測する手段としては、ａｂｉｎｉｔｉｏ分子動力学法や第一原理分子軌道法が使用されてきた。

ａｂｉｎｉｔｉｏ分子動力学法は主に結晶の構造や電子状態の調査に使用され、分子軌道法は通常の分子や、結晶の一部を切り出したクラスタを対象として、その構造や電子状態の計算を行うのに使用されてきた。ところが、これらの方法では、系の各電子に対してシュレテ゛ィンガー方程式を解くため、膨大な計算時間が 必要であった。そのため、経験的ポテンシャルを与えた分子動力学法のように、数10万ステップもの計算回数を必要とする手法の中で、これら第一原理手法を取り込む事は現実的に不可能であった。一方、分子動力学の分野では、1991年に、A.K.RappeとW.A.Goddardが"Charge Equillibration for Molecular Dynamics Simulations"という論文の中で、多原子分子に対する新しい電荷計算手法を 提案している。本発明者らは、A.R.Rappeらによって提案されている多原子分子 における電荷計算手法をさらに大きく拡張し、種々の元素や結晶系に適用できるように大幅に変更し、半導体装置の材料として一般に用いられるSiやSiO₂はもちろん、PbTiO3やPZTといった強誘電体や、高誘電体にも適用できる材料設計シミ ュレーション手法を構築した。まずはじめに、Rappeらの手法について以下に述 べる。

Rappeらは、まず、孤立原子のエネルギーEをイオン化ポテンシャルと電子親和力の関数として考えた。すなわち、ある原子Aにおいて、電気的に中性な状態をA0で示す。テイラー展開の2次の項まで考えると、

原子が電気的に中性な状態のエネルギーは、(３０３１)式のQAに0を代入して得 られる。

Ｅ_A（０）＝Ｅ_A0 …（3041）
原子が+1価イオン化した時のエネルギーは、(３０３１)のQAに1を代入して得ら れ、-1価にイオン化した時のエネルギーは、-1を代入して得られる。

(３０４２)と(３０４１)の差はイオン化ポテンシャルIPになり、(３０４３)と( ３０４１)の差は電子親和力EAになる。

(３０４５)式のEAに−が付いているのは、電子を取るのを正の仕事と考えると、電子を加える事は負の仕事をした事になるためである。(３０４４)と(３０４５)の差をとると、電気陰性度χA0となる。

また、(３０４４)と(３０４５)の和をとると、次式が成り立つ。

ここで、(３０５２)式の物理的意味を考える。水素原子の場合のイオン化時のクーロン相互作用の模式図を図１２乃至図１４に示す。図から分かる様に、IPとEAの差は電子間クーロン相互作用bになっている。従って、φA軌道中の2電子間 のクーロン斥力をJAA0と表すと、
ＩＰ−ＥＡ≡Ｊ_AA ⁰ …（3053）
JAA0をidempotentialと呼ぶことにする。電子を加える事により、軌道形状が 変化するので、Hartree-Fock等の第一原理計算で計算されたJAAと(３０５３)式 の値は多少異なってくるが、ここでは(３０５３)式で近似する。

Idempotential JAA0と原子直径RA0の関係を考えてみる。図１５に示した簡単 な原子構造モデルを考えると、電子間のクーロン相互作用JAA0は以下の様になる。

ここで、 JAA0はeV単位、RはÅ単位である。(３０５３)式の妥当性を確認するため、J値とRとの関係を調べ、表1に纏めた。表1に示した様に、Jは原子サイズに ほぼ逆比例し、(３０５４)式から求められた原子の直径は等極結合距離とおおよそ一致していることが分かった。

(３０５１)(３０５３)式を用いると、孤立原子のエネルギー(３０３１)は、次式で書ける。

多原子分子の場合には、原子間の静電エネルギーも考慮する必要がある。静電エネルギーは次式で表せる。

ここで、JABは原子A,Bの中心に1電子電荷がある時のクーロン相互作用である。JABはAB間の距離に依存する。(３０５５)式を用いると、全静電エネルギーは(３ ０５７)式で表せる。

原子スケールの化学ポテンシャルχI(Q1,・・・・,QN)は次式で表せる。

電荷平衡の条件(３０５９)から、(N-1)ヶの連立方程式が得られる。

χ₁＝χ₂＝………＝χ_N …（3059）
すなわち、

となる。これに、全電荷の条件(３０６２)式を含めて、N個の連立方程式を解く 。

(３０６１)(３０６２)をマトリックスで表示すると、以下の式で表せる。

(３０６３)式を解くに当たって、各イオンが取り得る電荷には、上限下限があるので、電荷がその範囲外ならば、その電荷を境界値に固定する。(３０６１)式を境界値に固定した電荷と、そうでない電荷とに分離し、(３０６３)式を修正してマトリックスを解き、収束させる。すると、電荷のばらつきが計算できる。

この方法の中で問題なのはJABの決定である。原子AとBとの距離RABが大きいときには

でよいが、電子雲が重なる距離においては斥力を考える必要がある。ここでは、1つのスレーター軌道の項で電子密度を近似してしまう。すなわち、外側の原子 価軌道がns,np,あるいはndである原子において、何れの場合にも次の型の規格化されたnsスレーター軌道からなるとする。

(３０６５)式を用いて原子の平均サイズを求めると、

となる。下記の関係によって、原子価軌道指数ζAを選ぶ。

λは調整用のパラメタで、(３０６６)式で与えられた平均的な原子サイズと結晶の共有結合半径RAとの差を計算するために含まれている。

スケーリングパラメタλは、図１６に示したようなアルカリ金属ハライド分子を用いて設定する。アルカリ金属ハライド分子をＭＸとし、Ｍをアルカリ金属( Ｎａ、Ｋ、Ｒｂ、Ｃｓ)、Ｘをハライド(Ｃｌ、Ｂｒ、Ｉ)とする。それぞれの電 荷をQM、QXと表すと、(３０６２)式から、
Ｑ_total＝Ｑ_M＋Ｑ_X＝０ ∴Ｑ_X＝−Ｑ_M …（3068）
(３０６１)式から

(３０６９)式に(３０６８)式を代入すると、
〔数６２〕
（Ｊ_MM−Ｊ_MX−Ｊ_XM＋Ｊ_XX）Ｑ_M＝χ⁰ _x−χ⁰ _M …（3071）ところで明らかに
Ｊ_MX＝Ｊ_XM …（3072）
であるから、

(３０７３)式の変数はλだけである。計算されたQから(３０７４)式を用いて双 極子モーメントμMXを求め、実験値と比較してλを決定する。

μ_MX＝4.80324Ｑ_MＲ_MX …（3074）
RMXは結合距離の実験値である。(３０７４)式はQがe単位、RがÅ、μがデバイの単位になっている。双極子モーメントの実験値との比較から、λの最適値として0.4913が得られている。このλを用いると(３０６７)式より原子価軌道指数ζが図１７の表の値に求まる。

水素に関しては、Rappeらは、他の元素と異なる取り扱いをしている。

水素の電気陰性度をMullikenとPauling他の実験上の値と比較すると一致して いない。Paulingスケールでは、水素は炭素より電気陽性で、ホ゛ロンよりわずかに電気 陰性であるのに対し、Mullikenスケールでは、水素は炭素や窒素より電気陰性である。この電気陰性度の問題は結合している水素の軌道は自由な水素イオンのように広がる事ができないので、実行的な電子親和力EAが原子の時の値よりも小さくなるためである。

従って、EAHを変数にして、水素のχ0HとJ0HHとを再定義する。

実験との比較から、最小二乗法によって、
χ⁰ _H＝4.528eV
Ｊ⁰ _HH(0)＝13.8904eV …（3076）
実験でなく、HF波動関数の静電ポテンシャルから計算された電荷との比較から、 χ⁰ _H＝4.7174eV
Ｊ⁰ _HH(0)＝13.4725eV …（3077）
を使用しても良いとしている。

また、電荷に関する境界値を設定する。

χ_A ⁰，Ｊ_AA ⁰は、明らかに、原子価殻が空っぽ、あるいは満杯に対応する範囲 外では妥当でない。例えば、Ｌｉ，Ｃ，Ｏの場合には、下図に示したような範囲に電荷を制限する。この範囲外ではＥ_A（Ｑ）＝∞をとる。

マトリックス(３０６３)を解いて求めた電荷を、図１８乃至図２０に示した不等式の範囲内にあるか否か判定し、その範囲外である場合には、その電荷を境界値に固定する。

以上がRappeらの手法のあらましである。Rappeらは、前述の手法を有機系の分子に適用している。本発明者らは、Rappeらの提案した原子スケールの化学ポテンシャ ルという概念だけを採用し、結晶系における個々原子の電荷計算手法を新たに開発した。以下に本発明者らによる電荷計算手法を述べる。

(３０５７)式右辺第二項は、すべての結合の静電エネルギーの和を意味する。本発明者らは、分子から結晶に拡張し、周期境界条件を取り入れようと試みた。すると、静電エネルギーを計算する際、計算に用いている数百〜数千の原子だけでなく、その外側にある無数の仮想原子からうける静電エネルギーも考慮する必要があることを発見した。すなわち、静電エネルギーは減衰が小さく、遠方の原子からの影響も蒸し出来ないことに気がついた。そこで、計算に用いているセルの外側に仮想的に無数のセルが存在すると仮定して、仮想セルからの静電エネルギーも考慮するようにした。

この場合、(３０５７)式仮想セル中の原子から受ける静電エネルギーJABは 、計算セルの大きさを考慮すると、RABが通常10オングストローム以上の大きな値 となるので、クーロンエネルギーだけと考がえて良い。従って、JAB=1/RABと一 律に扱う。仮想セルをνで表すと、(３０５１)式は次の様に変更される。

ここで、ν=0は計算セル自体を意味し、ν≠0は仮想セルを意味する。(３０７８)式 第3項にエワルト゛法を適用する。エワルト゛法を用いると、第3項は次の様に求めることが出来る。

Gは逆格子ベクトル、Ωは計算セルの体積、κはエワルド計算の収束を調節する パラメタである。更に、(３０８２)式右辺第一項はν=0且つA=Bの項は含まず、 第二項はG=0の項は含まない。これらの式から、原子スケールの化学ポテンシャルχIを求めると、

となる。電荷平衡の条件から、N-1個の方程式が得られる。

χ₁＝χ₂＝χ₃＝……＝χ_N …（3092）
(３０９１)式を代入して、

(３０９３)式で表せるN個の連立方程式を解くことになる。

(３０９３)式からCQ=-Dのマトリックスを作ると、次式が得られる。

本発明者らは、SiO₂の様な結晶や非晶質を扱う場合には、周期境界条件を考慮しながら(３０９４)式を解けば良い事を発見した。(３０９４)式は電荷Qに対し て非線形となっているため、計算プログラムにはQを収束させるためのルーフ゜を持 ってくる。

さらに、本発明者らは、実計算にあたってのクーロン積分のモデル化にも注意を払った。JABは、ζAとζBに対するスレーター関数の2原子クーロン積分として求め た。クーロン積分JAB。はRoothernの方法を用いて計算した。Roothernらは、ク ーロン積分JABには、 という標識を使用する。 (３１０１)式で定義されたハ゜ラメタを用いて、クーロン積分は(３１０２)式で表せる。

ここで
〔数７１〕
（α，α′，β，β′）
＝（１＋τ_a）α（１−τ_a）α′（１＋τ_b）β（１−τ_b）β′
…（3103）である。さらに、

ここで、上式の単位系は原子単位系であり、距離はbohr、エネルギーはHartreeを使用している。

Roothernらは1s,2s,2pに対するクーロン積分の解析式を与えているが、本発明者らは、3s以降の原子軌道に対するクーロン積分まで計算し、半導体装置で一般的に用いられるSi及びSiO₂が扱えるように、最適な軌道の選択を行った。

JABの値は電荷計算に大きく影響する。そのためJを算出する際に用いるスレーター軌道や原子価軌道指数ζの設定は重要である。そこで本発明者らは、Si-OやSi-Siを扱う場合のJの値について調べた。

Rappeらは「分極電荷がns軌道にあると仮定」しているが、種々元素に対して どの軌道を考慮するのが最適であるかという点は不明である。 Siの主量子数nが3であることから、Siの電荷が3sにあるとモデル化した場合と、2sにあると仮定 した場合の計算を行い、比較した。このとき、n=3におけるJの解析式はRoothernの方法に従い、本発明者らが、別途用意した。例えば、Si-Siに対しては、次式 を用いた。

図２１はSi-Si間のクーロン相互作用Jの計算結果である。この場合には、2sと3sの何れを使用してもJの値に殆ど変化が無いことが分かる。ところがSi-O間の クーロン相互作用Jは2sと3sを使用した場合で大きく異なっている。図２２にOの電荷が2sにあり、Siが2sあるいは3sにあるとモデル化した場合のJABの計算結果 を示した。Siの電荷を3sだけにおくと、3s軌道が非常に広がっている影響を受けて、Jの値がr=1.6Å近辺で極小値を持ってしまう。Si-Oの結合距離は丁度この当たりにあるので、Jが非常に小さく見積もられ、結果、SiとOの何れの電荷も極めて過小に評価されてしまう事が分かった。Si-Hの場合も同様であった。従って、Siの電荷のモデルとしては、敢えて3s軌道で仮定すると、Jの見積もりが全く精 度を欠いてしまう。このことから、スレーター軌道を用いるのは、電荷の広がりを考慮するためだけとし、モデルとしては、1sあるいは2sのスレーター軌道を用いるのが最適である事が分かった。これは、Si以降の原子番号の元素に対しても同様であり、3s以降の広がった軌道をモデルとすると、原子間距離が近い場合のクーロン力の計算精度が非常に悪くなることが分かった。

水素に対する補正にも、本発明者らは、プログラム上の注意を払った。 (３０３１)式から、

水素補正した時の原子スケールの化学ポテンシャルχHlを求める。まず、(３０ ５７)式に対応する多原子分子のエネルギーの式は、全原子数N個の内、水素M個 とすると、

ここで、(３１１４)式の右辺第一項、第二項、第五項のΣ′の記号は、それぞれ の和を求めるに当たってHは除いている事を意味している。原子スケールの化学 ポテンシャルχHlは、

(３１２１)式右辺第四項は次の様に解きほぐす事が出来る。

(３１２２)式より、(27)式の第四項と第五項との和は、Hl以外の全原子との間の原子間クーロン相互作用を意味する事が分かる。従って、化学ポテンシャルは、

(３１２３)式を用い、水素補正時のマトリックス(３０６３)に変更を加える事になる。今、χ1が水素以外の原子であるとすると、(３０５９)式に(３１２３)式 を代入し、

マトリックス(３０６３)と比較すると、水素のJii0を

に変更する事になる。

次に原子間クーロンポテンシャルJAHの水素補正を考える。水素の場合のスレーター軌道関数は、(３０７５)式より、

水素が含まれる場合には、電荷を求めた後に(３１２６)式を用いてJAHを修正 し、再度電荷を求める。そして電荷の値が収束するまで、これを反復する。この際、I=1の行が水素の場合、水素でない行と入れ替えてマトリックスを作成する ように、本発明者らはアルゴリズムを作成した。

電荷境界条件の扱いにも、本発明者らはマトリックスの具体的な変更方法について注意を払った。

(３０６１)式を、境界値に固定した電荷QBと、そうでない電荷QCに分けて考える。

(３０６２)式を(３１３１)式に変更すると、マトリックスの要素は、のとき、

上式の様にマトリックスを変更して、電荷を境界条件の範囲内に収める。この時、I=1の元素が境界値に設定される場合、マトリックスの行を入れ替え、境界 値設定されていない元素に対する行がI=1になるように行の選択(ピボット選択) をするようにした。

本発明者らは、結晶系及び非晶質系に新たに開発した上記電荷計算法を用い、これを分子動力学法に追加して、個々原子の動きを逐次計算し、双極子モーメントを時々刻々計算して記録し、これをフーリエ変換してIRスペクトルの算出を行った。図２３は、以上の手続きをまとめたものである。

一方、１６Ｇ以降のＤＲＡＭに使用されるゲート酸化膜やＮＡＮＤ ＥＥＰＲＯＭのトンネル酸化膜においては、高信頼化が従来にも増して重要な課題となっている。特に、ＥＥＰＲＯＭにおいては、高電界下でトンネル電流を流しながら絶縁性を維持させるという過酷な条件下で酸化膜を使用するため、絶縁材料としての限界を原子・電子レベルで探索することが重要となっている。

このような状況から、Ｓｉ／ＳｉＯ₂界面構造やＳｉＯ₂ネットワーク構造の原子レベルでの探索実験や、ＳｉＯ₂のトラップ機構の実験、及び理論計算による 予測等が種々の研究機関で進められている。しかし、本発明のように、ネットワーク構造とＩＲフルスペクトルとの関連を原子レベルで関連づけ、利用する手法は、従来にはなかった。僅かに、ＩＲスペクトルのＳｉ−Ｏの非対称伸縮ピークと酸化界面の応力とを関連づける実験結果が報告さているだけであった。

また、実験結果の解釈としても、「中心力モデル」を用い、現実のネットワークを非常に簡略化したＳｉ−Ｏ−Ｓｉ分子構造の振動解析からの議論しかなされておらず、非晶質特有のＳｉ−Ｏ結合長の分布、Ｓｉ−Ｏ−Ｓｉ、Ｏ−Ｓｉ−Ｏ結合角度の分布等を全て考慮した議論は、従来全く行われていなかった。本発明においては、非晶質特有のこれら様々な分布を考慮することにより、中心力モデルとは異なった視点で、電気的特性と構造との相関が得られたのである。

ここで、本発明者らが克服した点について付言する。本発明者らは、古典的分子動力学法を使用するにあたり、そのポテンシャルを吟味した上、さらに、改良を加えた。本発明者らは、最初、常行氏らによって開発されたＴＴＡＭポテンシャルの使用を試みたが、Ｓｉ−Ｏ非対称伸縮振動が非常に低波数側に表されるという欠点を見い出した。そこで、ｖａｎ Ｂｅｅｓｔ等によって開発されたＢＳＫポテンシャルを使用した。彼らのポテンシャルは、ＴＴＡＭポテンシャルよりも高波数側にピークが表れたが、尚、実測値よりも低い波数に伸縮ピークが表された。

本発明者らは、このピーク値を実測に近づけるため、これら既存ポテンシャルで使用されている、電荷一定という条件を取り除き、Ｓｉ−Ｏ−Ｓｉ角等が大変位した場合の電荷の分布を、古典的分子動力学に取り込んだ。この電荷分布の計算は、理想的には分子軌道法を用いて算出できるが、現実問題としては、数百原子もの構造を一度に取り扱って電荷分布計算を行うことは現状では不可能である。スーパーコンピュータを用いても現実的には不可能である。

そこで、本発明者らは、Ａ．Ｒ．ＲａｐｐｅとＷ．Ａ．Ｇｏｄｄａｒｄ III によって提案されている手法を用いて電荷分布計算を行った。引用論文は“Ｃｈａｒｇｅ Ｅｑｕｉｌｌｉｂｒａｔｉｏｎ ｆｏｒ Ｍｏｌｅｃｕｌａｒ Ｄｙｎａｍｉｃｓ Ｓｉｍｕｌａｔｉｏｎ”，Ｊ．Ｐｈｙｓ．Ｃｈｅｍ．，９５，１９９１，ｐ．３３５８−３３６３である。本発明者らは、彼らの論文を参考に数式を展開し、論文中の誤りを訂正の上、必要なパラメータ算出を独自に行って、彼らの手法をＳｉＯ₂系に適用した。Ｒａｐｐｅらの手法の概略について以下に 述べる。

原子ｉの個別の電荷によるエネルギーＥi （Ｑ）を、テーラー展開の２次まで考えて、下記式（１６５１）で近似する。

ここで、Ｘi 0 は電気陰性度であり、下記式（１６５２）で表される。

Ｘ_i ⁰＝１／２（ＩＰ＋ＥＡ） …（1652）
Ｊiiは、原子Ｉの軌道上の電子間のクーロン斥力であり、下記式（１６６１）で定義される。

Ｊii＝ＩＰ−ＥＡ …（1661）
ここでＩＰはイオン化ポテンシャル、ＥＡは電気親和力を示す。Ｎ個の原子で構成される分子の全エネルギーＥ（Ｑ1 ，…，ＱN ）は、個別の原子のエネルギーと原子間の相互作用静電エネルギーＪ_ijＱi Ｑj の和である。ただし、Ｊ_ijは原子_ijの中心にある単位電荷のクーロン相互作用である。従って、下記式（１６６２）が成立する。

原子スケールの化学ポテンシャルは、下記式（１６６３）で与えられる。

この式（１６６３）を、電荷平衡の条件、Ｘ1 ＝…＝ＸN に代入すると、下記式（１６７１）が得られる。

すなわち、下記式（１６７２）が得られる。

さらに、下記式（１６７３）に示す全電荷の条件を考慮する。

Ｊ_ijは既知であって、およそ下記式（１６７４）である。

ここでＲはＩ＝ｊのときは原子の大きさ、Ｉ≠ｊのときは原子間距離である。従って、行列Ｃとベクトルｄを導入すると、Ｑに関する下記式（１６８１）に示す連立方程式が成立する。すなわち、

と定義すると、
ＣＱ＝−ｄ …（1682）
となる。論文ではＣ1i＝Ｑi 、Ｃd ＝ｄとあったが、本発明者らは、これらを上述の様に修正して使用した。さらに、Ｓｉ−Ｏ−Ｈ等の原子の大きさ等、必要なパラメータを独自に算出し、連立方程式を解いて電荷分布を求めた。これによって、ピーク波数が格段に実測と一致するようになった。

分子動力学シミュレータＩＩで算出された「戸籍簿」、つまり単位時間毎の結合長（位置ベクトルｒ（ｔ））と結合角を示す膨大なデータの一例を、図２４に示す。このデータは、膜特性算出ユニットＩＩＩに入力され、ここで基本膜特性（光学的、電気的、結晶学的特性）が算出される。

又、プロセス変動算出ユニットＩについては、その２次元汎用システムの一例が、特開平３−１６４９０４号に記載されている。

３次元のプロセスシミュレータは、今回本件発明者が構築に成功した。以下これを、図８をもちいながら説明する。まず技術とシミュレータの構造の概観をする。本願は、原子レベルに主体を置いているが、他のプロセス変動算出ユニットＩとの新しい合体方法をも提案しており、より効果的かつ正確で、従来にない新しい機能をも提供できることをも示している。合体する相手のプロセス変動算出ユニット（１０）は、例えば３次元プロセスシミュレータＩｃであり、いわゆる流体モデルに属するものである。とくに、本発明者らが新しく構築した３次元プロセスシミュレータＩｃの内容を以下に記す。要点は２つあり，まず、本発明者らが構築した「３次元プロセスシミュレータ」は、原子レベルシミュレータとの接続が可能であること。次に、この「３次元プロセスシミュレータ」の部分には、本発明者等が新しく克服した化学モデルや、数学式を駆使している。これらにより、初めて、原子レベルとの合体が可能になったのである。原子レベルシミュレータから出力されるものとして例えば、粘弾性係数、応力定数、拡散定数、Ｃ_ij、等があり、これらが、流体的汎用３次元プロセスシミュレータに転送されるものである。図８では、（２３）、（３３）、（４３）、（５３）、（６３）、（７３）、（８３）である。例えば、酸化剤等の拡散定数は、（２３）に転送されるわけである。転送手続きの方法等は後述する。

汎用３次元プロセスシミュレータ（１０）側で、入力データに従って、数ミクロンレベルの広領域の工程誘起応力の計算を、（２）や（３）で、行い、この結果を、ファイル（２２）やファイル（３２）を介して、一旦分子動力学シミュレータＩＩに転送し、あるいは、再び汎用シミュレータ（１０）に戻し、これを繰り返すことにより、より正確でより多くの機能を利用することもできる。計算手続きとしては、例えばある時刻で、プロセス変動算出ユニットで求めたある程度広領域の工程誘起応力計算の結果を、一旦分子動力学シミュレータＩＩに転送し、ここで、物性定数の修正量を計算しながら、改めて、原子の動きをみ、再び、修正された物性定数を汎用シミュレータに送り込み、汎用シミュレータで広領域の工程誘起応力等の結果を、進める。これを繰り返すことにより、もし、十分大きな応力が局所的に生成されれば、分子動力学シミュレータＩＩ側で原子の動きをみ、塑性変形の可否を見ることもできる。

勿論上記の汎用プロセスシミュレータ部分での不純物再分布計算や、応力計算の時には、図８に記した、点欠陥挙動計算部分（９）が呼ばれている。点欠陥としてはＳｉ基板中の空格子と格子間Ｓｉを扱っている。不純物の拡散にはこれら点欠陥の挙動を正しく取り入れないと、不純物の再分布は正確には計算されない。イオン注入計算部分（４）では、もちろん、注入損傷による多量の空格子と格子間Ｓｉが生成される。これらは、熱工程中に再分布するが、その時、不純物とも相互作用を持ち、不純物分布の計算に反映させなければならない。またシリサイド・形状・応力計算部分（７）でも、勿論、Ｓｉ基板中の空格子と格子間Ｓｉの挙動に十分配慮しなければならない。

ところで、図８のその他（８）についても、少し触れる。ＬＳＩ素子試作工程で、酸化・拡散・イオン注入等は、既に多くのデータを基に、このプロセス変動算出ユニット（１０）には、各必要物性定数が格納されている。例えば、Ｓｉの酸化に関しては、後に掲げる式（１１３０）や（１１３１）のＰＢ１やＰＢ２等がそれである。また、ＳｉＯ２膜中の酸化剤の拡散定数などは、後に掲げる式（１５３１）のＤ’などがそれである。一般には、この様に、概ねデータがあるものである。

しかし、ＬＳＩ開発の中では、必ずしもデータ等が揃っているとは限らない。

たとえば、今、強誘電体膜を用いたメモリーセルの試作開発を考える。ＣＤＶ或いはスパッタ工程を用いて、ＰｂＴｉＯ₃を形成する場合を考える。現状の汎用 プロセスシミュレータでは、多くはまだＰｂＴｉＯ₃膜中での種々の不純物の拡 散係数も、さほど分かっていないし、また、その応力定数もさほど分かっていない。また、界面での不純物の偏析定数もさほど分かっていない。この様な場合、まず（８）の工程に入り、すぐさま、（８４）に示す様に、分子動力学シミュレータＩＩに一旦行き、そこで、必要な元素のポテンシャル等を算出し、それを基に、分子動力学部分が動き、そして、物性定数を算出する。このあとは、既に述べてきた部分と非常に良く似ているが、（８３）からこれらの物性定数がファイルに格納される。そして、プロセス変動算出ユニットと同じようにして、計算部分が進められる。

ここで、図９を参照して、具体的な処理に即した説明を行う。ここでは、酸化工程を取り上げる。即ち、ウエハーの投入を行い、これを基板処理工程に移し、しかる後に、酸化工程に入る。

この酸化工程では、例えば、高純度酸素ガスの量を例えば、Poxとし、また、 これの分散をσtoxとして構成する。叉、昇温シーケンスはＴempとし、またこれの分散をσtemp、バルプシーケンスをtox、σtox等の様に構成する。特に、バルプシーケンスは、ここでは、例えばパルプを開けている時問をさしているわけであるが、詳細は実施例で述べたようなものである。

たとえば、図８では、（２）の酸化・形状・応力計算部分をアクセスし、ここに分圧Poxや、温度Ｔemp等が入力される。そして、酸化成長計算や、同時に応力計算も進められる。この時、（２１）を用いて、分散・変動計算も行う。ここではσpoxやσte即が入力して、中心値からの変分を効率よく計算するものである 。

これらの計算結果は、例えば、酸化膜の成長量や、酸化膜とＳｉの界面の形状や、さらには、不純物量も計算される。さらに、本発明者によれば、応力も求められるわけである。また分散・変動からも同じように Pox＝Pox＋σpoxや、Ｔemp＝Ｔemp十σtemp等の部分についても、酸化膜の成長量や、酸化膜とＳｉの界 面の形状や、さらには、不純物量が計算される。そして、これらが、一旦メモリまたはファイル（２２）にたくわえられ、そして、分子動力学部分に送られる。

分子動力学部分では信号ａとしてこれらが入力され、位置ベクトルｒ（ｔ）が算出される。

図１０は、図９と同様の行程図を示しているが、ここでは強誘電体膜の形成に関するものである。

ところで、ここで本発明での診断部分に関する説明も行っておく。図１１に、本実施例による２進木の一部を示す。ここでは酸化工程に着目した２進木を示している。図中にも記したように、ここではまず、in-situモ二夕から送付される ＳｉＯ₂膜厚と、計算から予測されるＳｉＯ₂膜厚を比較する。そして、膜厚が許容範囲に収まっている時を取り上げ議論している。そうすると、次の手続きとして、inlituモニタから送付される１Ｒスベクトルカープと、計算から予測される１Ｒカーブとを比較する。そして、図に示した様に、まず、スペクトルを比較して、少しずれていた場合で、しかも、スペクトルが高波数側にずれた場合を取り上げる。そうすると、次の手続きとして、ＳｉＯ₂中に残留応力が大きく誘起し ていることが考えられる。それで、この場合、ゆっくり降温して、残留応力を軽減させることが出来る。それで、図９に記載したアニール部分を修正して、図８の（２）に入力される。他方、スペクトルを比較し、許容範囲に入っていた場合は、当初のプロセスシーケンスのまま続行される。

図８に関して、概ねデータの流れがこれで分かったと思う。では、もう少し、図８の出力に関して言及したい。即ち、実行結果は、一旦、ファイル（３）に蓄えられる。ここでは、各工程で求めた、工程に対応するスペクトルや、その他の光学的・電気的特性結果が格納される。これらは、あとで、詳細を述べるが、実測データとしてモニタをしている量と対比できる様になっている。

一方、本発明者らの１人は、以前、プロセス変動の変分計算を用いたプロセスシミュレータシステムを提案している。即ち、特開昭６３−１７４３３１号および特開平３−１６４９０４号である。これらの提案では、流体的な取扱いにより、酸化／拡散、イオン注入、化学的堆積（ＣＶＤ）、エッチング、リソグラフィ等の各工程を取り扱うことができる。例えば、流体モデルによる拡散工程では、不純物総量と電気的活性化不純物量とを区別しながら、電気効果を取り入れつつ、不純物の再分布を求めている。

特開昭６３−１７４３３１号および特開平３−１６４９０４号に示すシミュレーションシステムの概略を、図２５〜図４９を参照して説明する。これらの提案では、プロセス変動については次のように扱っている。例えば、酸化、イオン注入、リソグラフィの工程を例にとると、平均の温度（Ｔ）とその許容バラつき（δＴ）、平均時間（ｔ）とその許容バラつき（δｔ）、平均圧力（Ｐ）とその許容バラつき（δＰ）、更に平均マスク寸法（Ｌ）とその許容バラつき（δＬ）、平均ドーズ量（Ｄ）とその許容バラつき（δＤ）などが入力値である。

プロセスシミュレータには、図２５に示すように、イオン注入、酸化／拡散、リソグラフィ等の各サブルーチンに分散変動（ｄｅｖｉａｔｉｏｎ）を計算する部分が付設されている。計算手順としては、これらの中心値について計算を進めながらその変分を計算して行く方法を採用している。

この変分計算について若干説明する。まず拡散工程における温度変分について説明すると、拡散工程では、拡散係数をＤ、不純物分布による内部電界をε、電気的活性化不純物数をＮとして、Ｆｉｃｋの第２法則に基づく方程式を解いて、不純物濃度分布Ｃを求めている。図２６を用いて変分計算の手順を示す。ここで、図中の左側を正常プロセス（Ｎｏｒｍａｌ ｐｒｏｃｅｓｓ）とする。

非線形拡散方程式は、下記式（００７１）により表わすことが出来る。

ここでｋ、ｋ＋１は時刻であり、時刻ｋ＋１のときの不純物濃度Ｃk+1 が未知数である。右辺のＢは、下記式（１０７２）の第２項に由来するものである。

∂Ｃ／∂ｔ→（Ｃ_k+1−Ｃ_k）／ｔ …（1072）
ここでは陰解法を用いているので、上式に現れる係数行列はすべてＤや平均温度Ｔ等を用いて表現されている。正常プロセスの場合は、所定の拡散時間ｔn+1 番目までを適当に分割し、ｔ1 ，ｔ2 ，…，ｔn+1 まで順次解を求めて行く。ここではそれぞれの時刻、ｔ1 ，ｔ2 ，…，ｔn+1 においてＮｅｗｔｏｎ法を用いてこれらを線形化して計算をすすめている。

次に、変動プロセス（Ｐｒｏｃｅｓｓ ｄｅｖｉａｔｅｄ）について説明する。ここでは拡散温度に変動δＴが発生した場合を示している。上式の変分をとり、時刻ｋにおける濃度の変化量δＣK がわかっていると、時刻ｋ＋１の時の変分δＣk+1 は、下記式（１０８１）により表される。

すなわち時刻ｋ＋１において正常プロセスの収束時の値Ｄ，ψ，Ｎ，Ｃ等を保存しておき、これらを用いてＦ′を作成し、またその温度Ｔに対する変化量を用意しておき、δＣk+1 を求める。ここでＦ′を作成するための行列要素の一部を図２６に示す。拡散係数の導出には、基本的には空孔モデルを用いている。この場合、温度変化によってフェルミレベルが変化し、これによって荷電空孔の濃度が変化し、拡散係数が変化する。従って、まずその基になる真性キャリア濃度ｎi に対してはその変分は、下記式（１０８２）に示すようになる。

〔数９４〕
δｎ_i＝ｎ_i（１．５／Ｔ＋０．６０５／ｋＴ²）δＴ …（1082）
また、不純物の内部電界の温度による変化も、これらを用いて表現できる。更に、酸化性雰囲気における拡散係数の増大現象（Ｏｘｉｄａｔｉｏｎ Ｅｎｈａｎｃｅｄ Ｄｉｆｆｕｓｉｏｎ Ｅｆｆｅｃｔ）の変分に対しても考慮した。酸化工程では線形酸化速度定数（Ｌｉｎｅａｒ ｒａｔｅ ｃｏｎｓｔａｎｔ）及び放物線酸化速度定数（ｐａｒａｂｏｒｉｃ ｒａｔｅ ｃｏｎｓｔａｎｔ）も変化するのでこれも考慮した。これらを用いて先のＦ′を作成する。

ここで、ある典型的な拡散工程を取り上げ、その１つの拡散時間における収束の様子を図２７に示す。図２７の横軸は反復回数を示し、縦軸は残差を示している。図２７では、正常プロセスにおける収束状況が示されている。正常プロセスでは、３回のＮｅｗｔｏｎ ｌｏｏｐで収束していることがわかる。この収束時のＤ，ψ，Ｎ等を用いて先に示した変分行列Ｆ′を作成し、温度変分を計算する。図２７に示すように、約５０回で収束し、全体の１／１０の計算量で済んでいることがわかる。

また、図２８、図２９、図３０は、変分の精度を示す。図２８は、正常プロセスとして１０００℃の時の不純物分布を２次元格子にわたって表現したものである。図２９は、変動プロセスであって、δＴ＝１℃の場合で図２８をもとに計算したものである。他方、図３０は、１００１℃を正常プロセスとして打ち出したものである。

これらの図から、図３０の値＝図２８の値＋図２９の値なる関係があることがわかる。また、改めて図３０の値をプロットしてみると、図３１に示すようになり、高温にズレると高濃度側では不純物濃度が下がる傾向にあり、低濃度側では増大する傾向にあることがわかる。

半導体素子を作成する場合、多くのプロセスがあり、各プロセスにおいては複数個のパラメータがある。例えば、ＣＭＯＳ工程では主なものだけでも、（ｉ）Ｎ−ｗｅｌｌ領域へのイオン注入、（ii）Ｐ−ｗｅｌｌ領域へのイオン注入、（iii)ＬＯＣＯＳ工程のためのＳｉ₃Ｎ₄のパターニング、（iv）ＬＯＣＯＳ工程、（ｖ）ｎ+ イオン注入、（iv）ｐ+ イオン注入、（vii)アニール工程などからなる。プロセスの揺らぎを考える場合、これらの数多くのパラメータの変動とその度数（或は頻度）を取り扱う必要がある。

特開昭６３−１７４３３１号および特開平３−１６４９０４号では、変分の度数と実行の組み合わせについて考えている。今、簡単のためにα，β，γの３つの工程を考えてみる。それぞれにおいてプロセスパラメータの中心値をＣ（Ｃｅｎｔｅｒ）、正に変位した値をＵ（Ｕｐｐｅｒ）、そして負に変位した値をＬ（Ｌｏｗｅｒ）と名づけ、それらの度数を、特開昭６３−１７４３３１号および特開平３−１６４９０４号から、図３２に示す様に定義してみた。

工程αとβとを考えた場合、組み合わせは、ＵＵやＬＬなどの計算実行を無視すれば、ＣＵ，ＵＣ，ＣＣ，ＬＣ，ＣＬの５個の組み合わせが考えられる。さらにγのプロセスを引き続き経るときは、図中に示した７つの組み合わせが考えられる。しかし、特開昭６３−１７４３３１号および特開平３−１６４９０４号では、これらを全て取り扱うのではなく、ＵとＬの対称性を考慮し、Ｌの場合はＵから計算するようにしている。

このシステムの中心部である推論エンジンの探索の木の一部を、図３３に示す。ここでは、酸化・拡散工程について検討している。図中の“予”は、統計解析システムによる予測値を示し“実”は、実測値を示す。

このとき用いる許容値は、酸化・拡散・エッチング等全工程を経た値である。

このときε値は、２次元プロセス統計解析システムを実行して得られるものである。“実”が“予”±εの外に出た場合は、図３３に示すように、酸化工程に異常があった旨記述する。この場合、温度又は時間がどれだけズレていたかを必要に応じて、統計解析システムのロードモジュールを実行させる。上記検定は、酸化・拡散工程を中心に説明しているが、リソグラフィ工程についても行う。リソグラフィ工程でのゆらぎはＳｉＯ₂膜の形状の横ズレ、接合形状の横ズレなどに 相当する。

推論エンジンには、Ｃｏｍｍｏｎ Ｌｉｓｐを用いているが、これは数値処理のみならず文章（Ｌｉｓｔ）の処理、形状の処理に適していること、さらにシステムの保守が極めて簡単であるためである。本システムは、推論エンジンから実測値のデータを読むだけではなく、統計解析システムのロードモジュールをも実行している。

特開昭６３−１７４３３１号および特開平３−１６４９０４号に示すシミュレーションシステムによると、各工程における正常反復計算の収束時の行列の係数を保存しておき、且つこれを用いて一次変分行列を作成し、前記個々のプロセスパラメータの変動、及びこれらの重畳効果をも考慮することにより、従来の約１／１００の計算量で変動部分が処理できるようになった。また、たとえばＣｏｍｍｏｎ Ｌｉｓｐによる推論エンジンと上記システムを連動させることにより、プロセスマージンさらにはプロセス異常の診断さらに最適化まで出来ることがわかった。

特開昭６３−１７４３３１号および特開平３−１６４９０４号では、基本拡散方程式は、各不純物に対して下記式（１１１１）に示すように立てている。

ここで、ｘ，ｚは２次元空間の水平方向と鉛直（下向き）方向の座標であり、ｔは時間を示し、Ｃは不純物濃度、Ｎは電気的活性化不純物量である。またＤは濃度依存性と増速拡散効果を含む拡散係数であり、下記式（１１１２）により表わすことが出来る。

式中、Ｄi は濃度に依存しない基本拡散係数であり、その温度依存性は、下記式（１１１３）で表現される。

Ｄ_i＝ＢＢ・_exp(−ＢＡ/ｋＴ) …（1113）
ＢＢとＢＡは不純物の種類による定数、ｋはボルツマン定数である。ＦＥはｎi を真性キャリア濃度として、下記式（１１２１）で与えられる。

またｎは１／２[(Ｃ_D −Ｃ_A )＋（Ｃ_D −Ｃ_A ）² ＋４ｎ_i ² ]である。さらにｎi は、
〔数９８〕
ｎ_i＝３．８７×１０¹⁶×Ｔ^1.5×_exp（−０．６０５／ｋＴ）
…（1122）
とした。

酸化膜成長速度は、放物性酸化速度定数Ｒ１と線形酸化速度定数Ｒ２とを用いてそれぞれ下記式（１１２３）、（１１２４）に示すように表現した。

Ｒ１＝ＰＢ３・_exp（ＰＢ１／ｋＴ）
…（1123）
Ｒ２＝ＰＢ４・_exp（ＰＢ２／ｋＴ）
…（1124）
上記放物性酸化速度定数をあらためてＢ、Ｂ／Ａと記すと、水平方向の座標ｘの点で、酸化膜厚をＺo 、酸化時間をｔ、τを定数として、下記式（１１２５）に示す関係がある。

Ｚ_o＝Ａ／２・［１＋４Ｂ（ｔ＋τ）／Ａ２−１］
…（1125）
これらを用意しておき上記各係数の温度δＴに対する変分を求めてみる。

δｎi ＝ｎi ・（１．５／Ｔ＋０．６０５／ｋＴ2 ）δＴとしてｎの変分は、下記式（１１３１）に示すようになる。

また、下記式（１１３２）が成立する。

成長係数の変分は、下記式（１１３３），（１１３４）により表わされる。

Ｂ＝Ｒ１，Ａ＝Ｒ１／Ｒ２に書き改めると、下記式（１１４１），（１１４２）に示すようになる。

そして、ｔ＋τの変分は、下記式（１１４３）のようになる。

以上説明したシミュレーションシステムは、当時としては画期的なものであった。しかし、最近のクラスタ化ツールにおける工程制御やモニタ技術には、もはや使えなくなっている。

そこで本発明者等は、更に検討を重ね、本件発明を完成するに至った。ここでは、ペロブスカイト単結晶を用いた電界効果型ＭＩＳ素子の強誘電体膜の最適設計を行い、これに従い、実際に素子を試作し評価し発明完成を確認した。

即ち、１つは、半導体単結晶基板上に電界効果型ＭＩＳ素子を形成するにあたり、該ＭＩＳ素子の強誘電体膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、対象特性を表現するａｂｉｎｉｔｉｏ分子動力学理論に基ずいた評価関数を導入し、これにより所望の主動作条件下での系の最適解に従った単結晶強誘電体膜および基板の設計配置とすることを特徴とする電界効果型ＭＩＳ半導体装置の製造方法である。

叉、他には、半導体単結晶基板上に、電界効果型ＭＩＳ素子を形成するにあたり、該ＭＩＳ素子の強誘電体膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するａｂｉｎｉｔｉｏ分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、主動作条件下で自由エネルギが最小になるように前記単結晶強誘電体膜を配置することを特徴とする電界効果型ＭＩＳ半導体装置の製造方法である。

更に、他には、Ｓｉ半導体単結晶基板上に、電界効果型ＭＩＳ素子を形成するにあたり、該ＭＩＳ素子のゲート絶縁酸化膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するａｂｉｎｉｔｉｏ分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、この自由エネルギが、主動作条件下で、最小になるように前記単結晶強誘電体膜を配置することを特徴とする電界効果型ＳｉＭＩＳ半導体装置の製造方法である。

更に、他には、Ｓｉ半導体単結晶基板上に、電界効果型ＭＩＳ素子を形成するにあたり、該ＭＩＳ素子のゲート絶縁酸化膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するａｂｉｎｉｔｉｏ分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、この自由エネルギが、主動作条件下で最小になるように前記単結晶強誘電体膜を配置するとともに、その単結晶強誘電体膜の酸素欠陥濃度を０．０１％以下に抑えたことを特徴とする電界効果型ＳｉＭＩＳ半導体装置の製造方法である。

更に、他には、Ｓｉ半導体単結晶基板上に、電界効果型ＭＩＳ素子を形成するにあたり、該ＭＩＳ素子のゲート絶縁酸化膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するａｂｉｎｉｔｉｏ分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、この系の自由エネルギが、主動作条件下で最小になるように前記単結晶強誘電体膜を配置するとともに、特に局所的なエネルギの凸凹が偏差値３σ以内に納め、その単結晶強誘電体膜の酸素欠陥濃度を０．０１％以下に抑えたことを特徴とする電界効果型ＳｉＭＩＳ半導体装置の製造方法である。

以下にこれらの最適化手順ならびに、試作結果を記す。

ここでは、単結晶強誘電体膜として、例えばＰＺＴ系ペロブスカイトを例に取り、またＳｉ表面を例えば（１００）面とした。この基板面は例えば金属電極としてＦＣＣ構造のＣｕでも良い。一番重要な問題は、単結晶強誘電体膜ＰＺＴ系ペロブスカイトとＳｉ（１００）面とをどのような結晶位置関係に配置すれば良いかという問題に置き換えられる。ここでは、本発明の骨子である評価関数として系の全エネルギを取り上げた。しかも、この系としては、Ｓｉ／ＰＺＴ界面を取り上げた。本発明者等は厳密な評価関数としてのエネルギ表現式を予め作成しておき、これを用いることにより、印加条件下でのＰＺＴとしてのペロブスカイトと、Ｓｉ（１００）の位置関係は以下の様にするのが最適であることを見い出した。即ち、ペロブスカイトのＣ4ｖ構造表現のａ１軸とａ２軸のなす角の中線 を、Ｓｉ（１００）面上の［１１０］軸方向に合わせ、且つｃ軸を傾け、第１Ｓｉと第３Ｓｉ位置をＳｉ（１００）面上のＳｉに一致させることが最適値であることを見い出した。またこの時、ＰＺＴ中には０．０１％以内の酸素欠損が許容範囲であることも同時に示した。この様に作成した本発明の提案する概念図を図３４に示した。

図３４に付いて詳細に説明する。図３４aの上半分は、Ｓｉ側の単一格子の平 面図を示しており、下半分はＰＺＴ面のスケッチである。また、図３４cは、［ １１０］Ｓｉ方向に対して、ペロブスカイトを結ぶ線を並行に配置させるのが最も良いと指摘している。この時、後で詳細は説明するが、シリコンと結ぶ直線が［１１０］Ｓｉに平行になる為には、ペロブスカイトのＣ軸が傾く必要がある。

こうすることによって、本発明者らの計算によれば、単に，ペロブスカイトとＳｉを結ぶ線を並行に配置させて、Ｃ軸を傾けない場合にくらべ、全エネルギがさらに１０％低下することをも見い出した。この上で更に、実際素子への利用を考え、ＰＺＴ中に０．０１％以内までで有れば、充分単結晶ペロブスカイトの本来の特徴を示す事を見い出した。

ここで指摘する界面エネルギの問題は以下に説明するように、素子特性には多々影響を及ぼすものである。強誘電体膜と基板とが形成する界面に発生する不整合によるエネルギは、不対原子を形成しやすくし、これにより、界面準位の形成を助長することになる。この界面準位は、伝導現象の妨げや、信頼性の低下につながる。また膜中においても、異常ホ゛ント゛等の欠陥ができ、準位になる可能性を はらんでいる。ここでは、これらの欠陥の許容範囲を初めて示すものである。

本発明によって作成した、ＭＩＳ特性と従来例との比較をも行った。ここでは、特に界面準位を効果的に評価できるものとして、ＣＶ法を用いた。本発明による特性には、データのばらつきも正確に示した。本発明によった場合、その界面準位に値は従来例に比べ、約１/１０に低下していることがわかる。また本発明 によるＭＩＳ素子の移動度も同一電界状況下で評価すると約１２％の向上が見られた。また信頼性も向上した。

本発明者等は今までにない新しい軸関係を発明するに至った手順をさらに詳しく説明する。図３５が本発明者等が作成した、全く新しい厳密ａｂｉｎｉｔｉｏ/分子動力学シミュレータの構造図である。図３６の上半分は、ａｂｉｎｉｔｉ ｏの部分で、下半分は分子動力学シミュレータの部分である。このシミュレータを用いて、初めて、発明者らは、例えば ペロブスカイトの構造が如何なるものかを示した。図３７が、本発明者らが、初めて厳密な構造を示すことができた、ペロブスカイトの透視図である。図中の大きい玉は酸素を示し、小さい玉は、Ｐｂを示している。ここでＺ軸はペロブスカイトのＣ軸を示しており、図でわかる様に、ペロブスカイトのどの面とＳｉのどの面が巧く接合するかは、この様に複雑な面であるので、極めて推測は困難であることが良く分かる。

また、上記の分子動力学シミュレータの精度チェックとして、酸化膜に着目し、分光データの対比作業を進めた。即ち、時々刻々のＳｉ原子や酸素原子の位置等から、O-Si-OやSi-O-Siの、ロッキンモードやストレッチングモード等を読み とり、これらの数値から固有周波数を求めた。これと実測値とを比べたところ、それぞれ、440cm−１と1100cm−１に対し、490cm−１及び1111cm−１であり、妥当な範囲である。また圧力を印加した場合、O-Si-Oモードの強度が増大することも実測値と良く一致していることが分かった。

これらの検証作業をもすすめながら、本発明者らは、ＰＺＴに非常に大きな圧縮応力を印加した場合についても計算してみた。界面としては、この様な構造をとったほうが、全エネルギの利得があることもわかった。

計算に使った領域等について若干詳しく論じてみる。計算ではＳｉ基板側に、深さ方向に数十原子層をとり、また、表面での奥行き方向に数十原子層をとり、また、左右方向にも数十原子層をとり、計算領域をとった。またペロブスカイトにおいても、縦横高さをそれぞれ数十原子層ずつとった。

今までにない手法で、新しい、位置関係を見いだせたのは、この方法で摘めている。本発明者らは、強誘電体膜及びＳｉ単結晶を含む領域に対して、厳密なイオン結合ポテンシャルの厳密積算式を開発した。従来の手法では、厳密性に欠け、特に系のエネルギを正確に求めることが出来なかった。本発明者はこの明細書のなかで従来の欠点や不正確な場所に付いて順次明らかにして行く。

本発明者によると、以下のように説明できる。本発明者らは、従来にない新しいシミュレータを構築し、一例として、強誘電体膜の単結晶、具体的にはペロブスカイトについて計算を行った。

以下に、本発明者が作成した、シミュレータについても従来との比較をしながら、若干説明を加える。たとえば、まず、酸化膜について述べてみる。Ｓｉ−Ｏ間、Ｏ−Ｏ間、Ｓｉ−Ｓｉ間の３種類のポテンシャルを厳密に表現しなければならない。Ｓｉ−Ｏ間とＯ−Ｏ間、Ｓｉ−Ｓｉ間のtotalのポテンシャルは実際は すこぶる難しく、３種とも距離ｒの項により積算量は、無限大に発散することになる。また計算自体もクーロンポテンシャルは遠くまで裾を引くので打ち切ることができず、少し厳密な計算では従来の手法ではEwaldの方法を使っていた。し かし、後述するが、計算に厳密性を欠いている。

実際に配置するに当たり本発明者らが用いた手法を以下に示す。ここでは、本発明者らが初めて作成したシミュレータを以下に示す。

本発明者らは、鋭意検討し、計算物理学に則りながらも、新しい独自のａｂｉｎｉｔｉｏ分子動力学システムを構築した。この新しいシステムがあったからこそ、新しい試みが可能になった訳である。以下に本発明者等が作成したものを詳解する。

ここではＳｉＯ₂を例に採っているが、他のものについては、単に元素の番号 等を入れ換えるだけで使用できることも確認した。

本発明者等が提案した新しいシミュレータの全体構成を図３８に記す。分子動力学部分（ＭＤ）と密度汎関数（ＤＦＴ）とを合成した新しい手法を提案した。

まず本発明者らによる分子動力学法と密度汎関数法との使い分けを図３９を用いながら以下に示す。

図３６に示す様に、分子動力学（ＭＤ）部分はフローチャートの下半分を占め、密度汎関数部分（ＤＦＴ）は上半分を占めている。分子動力学法は零Ｋでない所謂有限温度でのやや大きな系を扱い、原子（イオン）間のポテンシャルには、予め第一原理から求めたポテンシャルを用いる。他方、密度汎関数（ＤＦＴ）は、小さな系の基底状態を求める時に使うことにした。そして、当然ＤＦＴの部分は温度は零Ｋとした。ポテンシャル部分は、電子系の基底状態から求めた。具体的には波動関数（ＫＳ方程式）を解き、局所密度汎関数を併用した。

本発明者等によるシミュレータの計算手続きを実際にＳｉとＯを例にとり説明してみる。ＤＦＴ部分で、縮重のない電子系の基底状態の全エネルギを求めて置
〔外１〕

]である。

即ち

〔外２〕

ン相互作用エネルギー、第３項は、電子の１体近似を想定して、電子間相互作用のない時の運動エネルギーであり、

と書かれる。第４項は交換相関エネルギーで、その形は局所密度近似による。電
〔外３〕

〔外４〕

方程式(Kohn-Sham equation)

〔外５〕

を自己無撞着(self-consistant) に解かなければならない。

〔外６〕

ルギーε_xc(n)に電子密度nを掛けて積分して、

とした。ε_xc(ｎ)の形はWigner他の式など、いろいろ提案されている。さらに、
〔外７〕

のクーロン場をそのまま使うのでなく、核の位置での特異性を消したノルム保存擬ポテンシャルでおきかえる。このようにならすことによって、フーリエ展開での成分の個数をおさえることができた。

全エネルギーとして、本発明者等は厳密を期するため、原子（イオン）間のクーロン・ポテンシャルも含むものを使う。即ち、

〔外８〕

ε など）。式（２０７１）から右辺第１項はＤＦＴの式（２０６５）から右
〔外９〕

である。その部分はＤＦＴの場合と同様に１体近似の式（２０６５）からＶeff
〔外１０〕

ここで新しい考え方を入れた。Ｅをポテンシャル・エネルギーと見なして、他に仮想的な運動エネルギー

を導入する。従ってラグランジアンは
Ｌ＝Ｋ−Ｅ …（2073）
式（２０７３）中のＭ_Iは本当の原子（イオン）の質量だがμ、μ は仮想的な 質量であり、後述のように目的に応じて変えた。束縛条件として正規直交条件

を課するから、オイラー方程式を作る時に未定乗数Λ_ikを導入して、Ｌに

を加えたものの変分をとる。その結果

式（２０７３）の運動エネルギーの温度Ｔが対応する。古典統計力学のエネルギー等分配則によって各自由度が温度を持っているのだから、Σ1/2Ｍ_IＲ_I ²との関係でいえば仮想的でなく現実の温度である。またμ、μ の大小によってν_iν を抑制したり増強したりできる。

例えば、μ Ｍ_Iとして高温からＴ→０とすれば、Ｒ_Iをとめたままψ_iが変化 して電子系の基底状態を得る。この時式（２０７６）の左辺が０になり、右辺は式（２０７１）とその下の注意により

となるからψ_iの線形結合を適当に作る｛Λ_ik｝が対角線化されてＫＳ方程式（ ２０６４）を得る。つまりＤＦＴでＫＳ方程式を解いたのと同じ結果をsimulated annealing で得たことになる。ただし、温度Ｔの状態を作るのに、モンテカルロ法によらず、仮想的な力学系の運動を追っているから、dynamical simulated annealing(DSA)と呼ばれる。

式（２０７６）の積分をする時、必ずしもψ_iをそのまま変数として扱うので なく、それを平面波で展開した展開係数（つまりフーリエ係数）を変数として扱うことができる。この展開係数の個数をむやみにふやさないためには、ＤＦＴの場合と同様に擬ポテンシャルで核の位置での特異性を消しておくことが必要である。特に結晶の周期性を利用できるときには、上記の｛Λ_ik｝を対角線化するψ_iの線形結合をφ_nkと記す。ｋはブリユアン帯域内の波数ベクトル、ｎは帯域イ ンデクス。

φ_nkは

を満たすはずである。λ_nkは行列｛Λ_ik｝の固有値であり、エネルギー準位である。φ_nk(ｎ、t)を展開して

〔外１１〕

満たされない）有効ポテンシャルＶ_effも展開して

式（２０７５）と式（２０７６）を式（２０７４）に代入して

〔外１２〕

共通になるから外して

を得る。この式は振動部分の積分を解析的にすましているから、純粋に数値的な方法よりもΔｔを大きくして、計算量を減少することが出来る。

クーロンポテンシャルは、以下の様になる。本発明者は分解してゆくと４項に分かれる事を見い出した。特に従来は３項しかなかったが、新しく第４項があることを確認した。これらの方法を順次しるす。これらは以下に示す様に、本発明者らは、まず厳密な式の出発として、誘電率を加味した基本式から解き始めた。

まず基本式である。

導体ε＝∞と真空ε＝１の場合とでクーロンポテンシャルが異なり、Ｌは（立方体の）単位結晶の一稜、Σは単位結晶内でとり、Ｚｉ、酸素ｉ第ｉ粒子の荷電と位置である。これは球内の荷電によって球の内面に分極が生じることによる。

導体でない球の内面に双極子の層が出来るが、上記式の最終項がちょうどその効果を打ち消す働きをする。Ewaldの方法は左辺、ε＝∞のものを与えるから、真 空内の値を求めるには上記式の最終項を加えなければいけない。結果だけ掲げる。この部分でも従来の導出はしばしば誤りがある。

クーロン・ポテンシャルを

〔外１３〕

、以下のように設定した。

と表される。ここに、ｎ_x，ｎ_y，ｎ_zがそれぞれｘ，ｙ，ｚ方向の稜のベクトル で、ｎ_x，ｎ_y，ｎ_zはそれぞれ（バルク結晶の場合）−∞から＋∞にわたる整数 である。Σ′の´はｎ=oの時のj=iを除くことを意味するものとする。

ここで新しいＦ関数を導入して、

〔外１４〕

を見い出した。

これを本発明者等はさらに下の様に変形して、

指数ｍ_x，ｍ_y，ｍ_zはそれぞれ−∞から∞にわたる整数である。 = の時は、
〔外１５〕

〔外１６〕

を使った。

これにより、始めのクーロンポテンシャルは、

〔外１７〕

これは を含まないから力には関係ない。さらに次ぎのように変形できる。

今までの結果をまとめて、クーロン・ポテンシャルは、結局、

〔外１８〕

/L_x,0,0)+m_y（0,1/L_y,0)+m_z（0,0,1/L_z)である。上記の表現の具合のいい点は、もとのΣの項が逆数のオーダーでしか減衰しないのに対して、Φ（１）ではΣの項がerfcの因子によって、Φ（２）ではΣの項がexpの因子によって、急速に減 衰することである。Φ（３）での減衰との遅速に逆向きに効くから、両者のバランスがとれるような適当なκを選ぶ。距離の近いところから順にクーロン力の寄与を計算した結果であり、しかも周囲が導体である場合の結果である。周囲が真空である場合によればもう１項が加わる。これが特に従来省みられなかった部分である。

本発明者等が作成した回転点群を図４０を用いて以下に示す。Ｓｉ原子を半径rの球の中心(0,0,0)に置き、球面上正四面体の頂点の位置に4個の酸素原子を置 く。Ｃ原子の位置は、オイラー角(θ、φ)で指定できる。それを(x,y,z)座標で表すには、北極に向かうベクトル(0、0、γ)に対して、まずy軸まわりθの回転、つ いでz軸まわりφの回転をしてやればよい。

その回転行列は、

ここで、4個の酸素原子のうち1つを北極(0、*)に置いた時、しかもφは不定で*としておくと、他の３つの酸素原子の位置は、未定の角ψによって、それぞれ、(109°28´,ψ)、(109°28´,ψ+120°)，(109°28´,ψ+240°)と表される。従って、それらの(x,y,z)座標は、式（２０６１）で、cosθ=-1/3、sinθ=2√(2)/3、そしてφをそれぞれψ、ψ+120°、ψ+240°としたものを、ｔ(0,0,r)に作用させて、

となる。北極を指していた第１の酸素原子を(θ、φ)の向きにするには、式（２ ０６１）のＲ(φ、θ)をそのまま作用させればよい。その結果、本発明者等は以 下のように求めた。4個の酸素原子の(x,y,z)座標は、

となる。単位格子内の4個のＳｉ原子の(x,y,z)座標をパラメタa,b,cで表す。図 ４０と照合すると、２ａが縦方向、横方向に共通の周期である。

I ( 0, 0, 0)
II ( a, b, c)
III(a-b,a+b, 2c)
IV ( -b, a, 3c)
V ( 0, 0, 4c)
すなわち、第１原子を原点におき、ｚ座標が一定値ｃずつ上り、(x,y)面に正 射影すると、I 、II、III、IVで正方形をなしている。第5原子は次の単位格子の第1原子でｚ座標が4cふえて、第1原子の真上にある。各Ｓｉ原子に属する4個の 酸素原子の配置は、Ｓｉ原子に相対的な配置がｚ軸まわりに90°ずつ回転して行く。それに上述のようなＳｉ原子の移動が加算される。第１Ｓｉ原子周囲の4個 の酸素原子の(x,y,z)座標は式（２０６１）をそのまま使えばよく、第２、第３ 、第４Ｓｉ原子の周囲の酸素原子の(x,y,z)座標は、第１原子周囲の酸素原子に 対して、それぞれ
Ｒ(0, 90°)の回転 と( a, b, c)の平行移動
Ｒ(0,180°)の回転 と(a-b,a+b, 2c)の平行移動
Ｒ(0,270°)の回転 と( -b, a, 3c)の平行移動を作用させれば得られる。

その他に、I II1とIII2、II III1とIV2、III IV1とV2の組がそれぞれ同一の酸素原子であるが、これらは（当然のことだが）上と同じ関係式を与える。酸素原子Ｉ1とII4は、横方向に１周期分2aだけずれている。

この様に調べてきた上で、さらに実際プロセスを考えてどの程度の欠陥密度まで許容されるかについて調べた。即ち、今、例えば３２００個のＳｉと３２００個のＴｉを用意し、これに０個から１００個の酸素欠陥を作った。この場合の時々刻々のＳｉ原子や酸素原子の位置等から、次の事が明らかに成って来た。まず、酸素原子及びＴｉ原子の動的挙動とＰＺＴ構造因子との関連を調べてみた。単結晶ＰＺＴ中に酸素欠陥が存在すると、数原子先の遠方のＴｉ及び酸素までも、無欠陥の場合と比較すると明らかに擾乱を受けており、しかも、点欠陥の移動度はすこぶる大きいことである。これらの動的挙動を、単結晶Ｓｉの場合と比べると、ＰＺＴの場合の方が極めて影響も範囲も大であることが分かった。この様な活発な動きは、ＰＺＴにはイオン結合成分がかなり有り、分極しやすい性格を持っており、従って、局所的な欠損が引き金に成ってポテンシャル変化に連動するものとも考えられる。本発明者らは、このような手続きで、詳細に調べたところ、０．０１％以下であれば充分単結晶強誘電体膜の利点を引き出すことが出来ることを見い出した。

本発明者による強誘電体材料の原子レベル計算について以下にします。本発明者等は、ＰＴＯやＰＺＯ、ＰＺＴの原子レベルでの構造把握や、構造と膜特性との把握計算を進めた。計算の結果を報告してみる。原子レベル計算で用いている対象構造としては、(1)ＰＴＯの完全結晶で、酸素、Ｔｉがそれぞれ高対称位置 にあるもの、(2)ＰＴＯの完全結晶で、酸素が、0.3Aずれたもの、(3)前記(2)に 酸素欠損を、ｃ軸にそった部分に１個入れたもの（図３７）、(4)前記(2)に酸素欠損を、ａ軸にそった部分に１個入れたもの（図４３）、(5)ＰＴＯのＴをＺｒ に置き換えたもの等である。特に、(5)においては、Ｚｒを、人工的に層状に置 換したものも実行した。酸素やＰｂの働きは、ＰＴＯやＰＺＯ等で異なっていることが初めて分かった。これらから、強誘電体膜の「自発分極」のきっかけが掴めた。また、強誘電体膜の原子配置や歪み、欠損構造等と、基本膜特性との関連もつかめた。

最も基本的なＰｂＴｉＯ₃膜の原子レベルでの計算についてまず示す。引っ張 り歪みを受けた場合や、試みにＰｂ位置を平衡点からずらせた場合（図４４、図４５）や、酸素欠損を有した場合等について、膜特性を求める上で最も基本になる電荷分布が上手く求まるかを調べた。図４４はＰｂＴｉＯ₃の単位胞の絵で、 さらにＰｂを中心位置から、少しずらせた場合である。この時の電荷分布の等濃度図を求めたものである。また図４５はやはり単位胞であるが、Ti-酸素結合の 部分をひずませた場合である。この結果から以下の点が分かってきた。(i)全体 に非常にイオン結合成分が多いこと。また、(ii)電気陰性度の強い酸素に殆ど大部分の電荷が引き寄せられ、次にTiの廻りに電荷があることが分かった。また(iii)このままでは電荷分布が対称的で自発分極は起こり得ないが、酸素「欠損」 や、Tiの位置変動を考えると、局所的電気的中性が破られることが、容易にわかる。さらに引っ張りにより、誘電率が低下することも確認出来た。

さらに例えば、２５℃で熱擾乱がないとして見たときの計算結果では、以下のことが分かった。即ち、点群Ｃ4vからなる理想的なＰｂＴｉＯ３膜の安定構造を求めた。その結果、酸素を極わずかずらしても、エネルギの安定な元の点群Ｃ4vに戻ることが分かった。さらに大きく、酸素位置を0.3Åずらしたところにも、 通常言われている様に、もっと安定な平衡点があるか。これを調べた。出発点0.3Åより極わずか離れた位置に向かってさらに安定になりつつある。この辺りが 通常言われている位置で、このずれ故に自発分極が生じるとも言われている。安定点が少なくとも二つあり、１つは単純立方晶の位置にあるもの、もう１つは若干ずれた位置であることが初めて分かった。

分極の起源は、局所的に電荷量の「中性」を、くずす様な原子位置が、自由エネルギ的に安定点であることとも置き換えられる。また、以上の計算では、酸素の存在は極めて大きいことが分かった。それ故、酸素欠損は、ひょっとしたら、結晶的にも、電子論的にも大きなトリガーになりうることがわかる。即ち、
(1)各３種類の原子にそれぞれ電荷分布は局在している様である。

(2)周期律表にある順序に従って、Ｏ及びＴｉ及びＰｂは、それぞれ（1s2,2s2,2p4)，(1s2,2s2,2p6,3s2,3p6,3d2),(1s2,2s2,2p6,3s2,3p6,3d10,4s2,4p6,4d10,4f14,5s2,5p6,5d10,6s2,6p2)である。これを反映して、特にＰｂの電荷分布領域が大きいのが分かる。

(3)良く見ると、とりわけ、Ｔｉの電荷分布は直ぐ隣のＯによって変調されて いるのが見える。

これを抑止するため、Ｆ，Ｃｌ、Ｂｒ、Ｉ等のハロゲン元素を上手く置き換えてやる手も有ることが分かった。

例２ 上の例では、Ｓｉ・ＰＺＴについて記したが、単結晶強誘電体膜として、スピネルＭｇＡｌ₂Ｏ₄膜、酸化セリウムＣｅＯ₂ 、チタン酸ストロンチウムＳｒＴｉＯ₃ を用いることもできる。また、酸化アルミニウムＡｌ₂Ｏ₃、酸化ジルコニウムＺｒＯ₂ 、酸化イットリウムＹ₂Ｏ₃、イットリウム安定化ジルコニウムＹＳＺ、ＰｒＯ₂ 、弗化カルシウムＣａＦ２ など、および、その積層膜をも用 いることができる。また、該単結晶強誘電体膜と下地シリコンウエハあるいは下地電極との界面にシリコンＰＺＴを挟む変形は本例以外の構造においても応用することは可能である。

また、電圧下での自由エネルギを小さくすることが大切であり、歪のエネルギだけに着目すれば、そのときは必ずしも最低になっていなくても良い。我々がペロブスカイトについて示したのも、思想に基づいている。

本発明による実際の測定系についても説明する。図４７は、ソーヤータウワの回路である。この方法の回路について説明する。Ｔは、高圧トランス、Ｃoは既 知固定容量、ＣＲＯは陰極線オシロスコープである。Ｃoを、強誘電体試料の容 量より、充分大きく選ぶと、Ｔの２次電圧は殆ど試料にかかるが、それは、抵抗Ｒｄによって、適当に分割されて、ＣＲＯの横軸に加えられるから、ＣＲＯの横軸の振れは、試料内の電場Ｅに比例する。一方、ＣＲＯの縦軸の振れは、Ｃoの 両端子間の電圧を表すが、これは、ＤＳ／Ｃoに等しい。即ち、ＣＲＯの図形は 、Ｄ−Ｅ曲線を表し、縦軸の振れからＤ，Ｅの絶対値も容易に得られる。自発変位Ｄsは、飽和部分をＥ＝０へ外捜したときのＤとして得られ、自発分極Ｐsは、Ｐs＝Ｄsできまる。損失のない常誘電体をこの回路で観測すると、原点を通る直線になるが、損失のある場合や、強誘電体はヒステレシスをもつ。図４８に本発明になるＤ−Ｅ曲線を実線で表し、比較のために従来のＰＺＴの結果を破線で示した。これから分かる様に、本発明になるものは、非常に誘電率が大きいこともわかる。

ここで指摘するところは、強誘電体膜と基板とが形成する界面に発生する不整合によるエネルギを印加時で最低にしようとするものである。しかも、この条件で実際利用上のことを考えその品質の許容範囲も示した。この界面エネルギを縮小することは、伝導現象の妨げになる準位の形成抑止や、信頼性の向上につながる。この構造を本発明によれば、新しいａｂｉｎｉｔｉｏ/分子動力学シミュレ ータによって予測するものである。所望の評価関数を搭載したａｂｉｎｉｔｉｏ/分子動力学シミュレータは、経験的なモデル等を一切用いていないので、この システムに指摘するところは、全て自然現象を完全に予測できるものである。ここでは、系の印加時の全エネルギに着目した例を示しているが、例えば誘電率を所望の値にすべく設計についても、さらには、劣化抑止を所望の問題であれば、それに適した設計が可能になる。

以上２次元の物性を計算する方法を示した。３次元の物性の計算方法の一例については、図８と共に既に説明した。

叉、実際の対比や診断がどの様にして行われているかをしめそう。即ち、シミュレーション結果は、スケジューラ（１３）に転送される。この時、他方、実際の実験系からは、モニターデータが、時々刻々送られ来る。そして、これらは、ファイル（４）に格納されるととものに、ファイル（４、３）からのデータ転送は、スケジューラによって同期させておく。そして、対比診断検定部ＩＶにある対比・検定・診断ブロックに転送される。

図４９は、本発明に使用されるクラスタ化ツールにおいて、ｉｎ−ｓｉｔｕモニタ測定装置がどのようにして付加されているかを示す概念図である。図４９に示すように、モニタ量は、例えば、入力ガス構成分析やＸＰＳシグナル、ＦＴ−ＩＲシグナル等であり、これらは、所定の時間空間内における平均値であり、ほんの一部分の情報である。本発明者らは、これら限られた測定量と、他方計算から得られた十分な情報とを対比させる方法を提供するものである。これらは、図５０の様に、モニター上のグラフィックな画像として表示される。

図５１は、本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムへの入力データの一例としての温度シーケンスを示すものである。また図５２は、本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムによる、原子レベル実行予測の概念を示すものである。このように、例えば原子レベルで、基板上での反応や、分子の解離などが記されている。

図５３は、本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムにより得られた原子レベル実行予測の可視化図である。これは、非晶質Ｓｉが堆積されたあと、６００℃近傍の高温でアニールされた状態を原子レベルでみたものである。

また、本発明のａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムは、Ｓｉのみならず、酸化膜に対しても適用することができる。即ち、図５４は、本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムのＳｉＯ₂につ いての実行例を示す。ここでは、ＳｉＯ₂におけるＳｉ−Ｏ−Ｓｉ角が時々刻々 求められており、これから、フーリエ変換し、固有振動数をもとめ、これを重畳することにより、ＦＴ−ＩＲのスペクトルが計算のみで求められる。この予測結果を用いて、実測値の時々刻々データと比較検定すれば良い。

図５５は、さきに示した図５１の入力データを用いた場合の、本発明のシミュレーションシステムを用いた原子レベルでの実行例を示す。図５５から、予想どおり実行されていることがわかる。

図５６から図６９は、図５１の入力データを用いた場合のＳｉ膜の成膜から、アニール、さらに固相成長までにおけるＳｉ−Ｓｉ間距離、結合角を、時間と温度を変化させてもとめたものである。

図１に示したシステムでは、基本膜特性の算出を膜特性算出ユニットＩＩＩで行った。しかし、図７０に示したシステムのように、特別な膜特性算出ユニットＩＩＩを省略し、プロセス変動算出ユニットＩ内部で必要な特性の算出を行ってもよい。それらは設計の要請に応じて、適宜変形可能である。

以上のように、本発明によって、ハードウエアや計算量を従来に比較して大幅に軽減することが可能となった。例えば、計算量について言えば、プロセスシミュレータＴＯＰＡＺでは、図９８に示すように、メモリは約２４Ｍバイトで、実行時間は数分間となる。これに対して、有名なプロセスシミュレータＳＵＰＲＥＭ４では、粘弾性等を考慮しているので、メモリは２００Ｍバイト以上、実行時間はクレイで１日にも達する。さらに分子動力学やａｂ−ｉｎｉｔｉｏ計算になると、メモリが数倍、実行時間も数倍にも達している。

このような状況では、分子動力学やａｂ−ｉｎｉｔｉｏツールをフルに用いて、原子の動的挙動の予測や、光学的諸特性の予測には、計算機のコストが大幅にかかかってしまう。本発明者らは、これら不利益を、新しい数学的な技巧を駆使し、克服する方法を見出すとともに、この技巧により、シミュレータの新たな利用技術をも見出したものである。

従来でも、上記に若干の手続きを加えたものもある。それが図９９に示すフローシートであり、従来の分子動力学シミュレータにおける密度汎関数の流れの主要部分を示す。即ち、この手法では、仮想質量を与え、全系のラグランジアンを用いる。これは、本発明者等が先に記した式（１２２３）に該当する。これはすでに知られているＣａｒ−Ｐａｒｒｉｎｅｌｌｏ法に相当するものである。この方法は一見高速に計算できるとされている。しかし、それでも例えば、クレイ社レベルのスーパーコンピュータでも、１つの場合の計算でも１昼夜を要する。従って、種々のパラメータの変動を調べるため、１回１回計算をしていたのでは、到底、実時間での診断ツールとはなり得ない。

以上説明したように、本発明によると、クラスタ化におけるモニタ機能の問題点を解決するだけでなく、ＭＯＳＬＳＩ等の半導体装置を製造する装置、とりわけ枚葉式あるいはクラスタ化製造装置において、各装置における限られた量の“その場”測定値や、限られた質的内容の“その場”測定値を最大限に拡張し、テストピース等が無くても、所望の工程通り、あるいは修正しながら進行するものである。更に、本発明によると、（ｉ）各プロセス要素の揺らぎや、各要素のずれの重畳現象を診断し、（ii）所望のシーケンス通り実行されているかのチェックや、（iii)変動因子の発見も可能である。更にまた、このシステムの特徴は、（iv）プロセスの変動を数学的・統計的に、しかも効率よく取り扱えるようにした点、（ｖ）新しい推論エンジンを付加させた点、及び（vi）新しい逆方向システムを付加した点にある。

以下、本発明の種々の実施例について説明する。

実施例１
この実施例は、本発明の大きな特徴の１つである光学的スペクトルの直接予測機能に係るものである。ここでもう一度、本発明者らが考え出した手法から詳細に記述してみる。

まず、高信頼化酸化膜設計に関する技術およびその評価関連技術等について、従来の技術的レベルから少し記述してみる。公知刊行物としては、例えば「The Physics and Technology of Amorphous ＳｉＯ₂ 」Plenum Publishing Corporation,ISBN 0-306-42929-2がある。この刊行物には、種々の酸化膜の光学的測定結果として、ＩＲの結果やラマンの結果が記載されている。この刊行物は、特に、原子レベルの議論としては優れたものと言われている。

この刊行物のｐ．６３には、ポール・マクミリアン（Paul McMillan ）によって非晶質酸化膜のラマン信号が、工程によってどの様に変化するかが論じられている。図７１及び図７２は、ポール・マクミリアンによって測定された結果である。例えば、図７１の曲線ａは、通常の酸化膜のラマン測定結果を示し、曲線ｂは、デンシファイした酸化膜のラマン測定結果を示す。

従来では、この両者のスペクトルを比較し、例えば、デンシファイにより水分が脱離し、これによってスペクトルが変化したものとされていた。また図７２は、単結晶水晶に圧力印加する前後でのスペクトルをみたものである。特に６１０と５００ｃｍ-1のところに変化があるとしている。しかし、これでわかるように、まず一見して、スペクトルの変化量はきわめて小さく、どこに着目すべきか、またそれを有為差とみるべきかについては大変技術を要するものである。またさらに、そのスペクトル変化量から物質の変化量を同定するにはさらに複雑な経験を要するものである。

この様な状況なので、今後さらに必要になってくる新しい材料の開発においては、スペクトルの同定テーブルがないものが多々あり、スペクトルの解釈が大きな障害になってくる。

また、上記刊行物のｐ．７１からは、小林等が高応力下での酸化物ファイバーのラマンスペクトルの変化をみている。その結果の一部を図７３に示す。ここではＳｉＯ₂ファイバに応力を無印加の場合と３．５％の引っ張りひずみを印加し た場合について示している。図７３から分かるように、３．５％の引っ張り歪みを加えた場合、４５０ｃｍ^-1近傍のところに変化がある。しかし、ここでも分かる様に、スペクトルの変化量はきわめて小さく、やはりどこに着目すべきか、またそれを有為差とみるべきかについては、大変技術を要する。また、そのスペクトルの変化量から物質の変化量を同定するにはさらに経験を要する。

上記刊行物において小林は、図７４を用いてそのスペクトルの解釈を試みている。図７４から分かる様に、引っ張りひずみを受けるとＳｉＯ₂構造の１０９， ５度の４面体部分が歪むとしている。しかし、彼は、たった１組の４面体構造を持ち出しているにすぎず、となりの４面体部分とはどのようなつながりになっているか全く議論していない。実際には、ある有限の温度であるので、必ず角度の分布を伴うはずである。

本発明者がこの明細書で何度も記述してきたように、定量的な分布を把握しないと、およそスペクトルの解析には程遠いものとなる。また図３３も同じ議論が小林によって展開されている。図３３（ａ）は、Ｏ−Ｓｉ−Ｏ結合角の関数としてのＳｉＯ4 四面体の通常モード振動の波長を示す特性図である。また、図３３（ｂ）、（ｃ）は、ストレスにより生じたＶ４モードのデカップリングによる、強度の減少を示す特性図である。

図７６は、上記刊行物のｐ．１０７に記載されたW.B.Fowlerによる酸化膜中の欠陥モデル図である。図７６の（ａ）は、完全α型クォーツの一部の結合モデルを示し、図中Ｌは長結合、Ｓは短結合である。（ｂ）はＥ1 ´中心のモデル、（ｃ）はＥ4 ´モデル、（ｄ）は、Ｅ2 ´モデルをそれぞれ示す。（ｂ）〜（ｄ）における不対電子は、ｅ- で表されている。

図７６からも分かるように、上記刊行物では、１０９．５度の４面体構造を２個または１個用意して欠損構造を議論している。しかし、上述したように、実際には、ある有限の温度で必ず角度の分布を伴うはずである。従って、ある種の多数個の計算を行わなければ、およそスペクトルの解析には程遠い。

図７７は、本発明者等が作成した新しい分子動力学シミュレータによる実行例である。Ｓｉは３２個、酸素は６４個を用いた。また逆格子点は５２０点であった。縦軸は体積の変化率、横軸は時間をそれぞれ示す。図７７から、１０００ｆｔ後でほぼ外圧と平衡になることが確められた。例えば、このときのＯ−Ｓｉ−Ｏ角、Ｓｉ−Ｏ−Ｓｉ角、及びＯ−Ｓｉ距離を図２４に示す。

図２４から分かる様に、各原子の「戸籍簿」に相当する膨大データが必要なのである。本発明者等は、計算手法にも幾つかの工夫を用意している。例えば、赤外吸収では、時々刻々の原子の振動運動の時系列データから、赤外吸収の選択律を満たす振動モードを抽出し、その頻度を求めればよい。いずれにしてもこの様な試みはまだ発表がなく、これはプロセス開発の新しい「兵器」と言うことが出来る。

本発明者等は、分子動力学システムを用いて、十分、平衡に達したところで、時刻を０に設定し、その後、時刻ｔまでを今サンプリング時間とする。時間相関は以下の通り求めた。即ち、
Ｆ（ω）＝∫＜Ｍ（ｔ）・Ｍ（０）＞ｅ^-1ω^wtｄｔ
…（1361）
である。ここでＭ（ｔ）は分極率であり、下記式で表される。なお、上記式の積分値は−∞から∞まで積分した値である。

Ｍ（ｔ）＝Σｅ_i・ｒ（ｔ） …（1362）
また＜Ｍ（ｔ）・Ｍ（０）＞は、Ｍの自己相関関数を表わす。

ここで、ｅi はその点における電荷数であり、ｉは粒子の番号である。またｒはベクトルであり、その測定起点は特に制約はないが、計算セルの角にとるのがよい。

時間相関を求めるには、それぞれの時刻におけるｓ（ｔ1 ）・ｓ（ｔ2 ）・ｓ（ｔ3 ）・ｓ（ｔ4 ）・ｓ（ｔ5 ）・ｓ（ｔ6 ）・ｓ（ｔ7 ）・ｓ（ｔ8 ）を求めればよい。

例えば、サンプリング時間全体に渡っての時々刻々の座標位置の計算値やその運動のエネルギ、さらには運動量等が、全粒子に渡って存在している。しかも、さらに互いの粒子間の距離や、それらの角度位置関係等も計算されている。

本発明者等は、幾つかの結晶性酸化膜や、また非晶質酸化膜についても計算した。その結果の一部を図７８に示す。即ち、図７８は、１つの例として単結晶クリストバライトのＳｉ−Ｏ距離の計算結果を示す。また、図７９は、欠損がある場合のＳｉ−Ｏ−Ｓｉ角の計算結果を示す。図７９から、欠損の回りでは運動が大きく乱れていることが分かる。これらの欠損の濃度を種々変えて同様の挙動を求めることも行った。

図８０は本発明者等が作成した新しい分子動力学の実行例であり、原子の動きをシミュレートしたものである。図中、大きい玉は酸素を示し、小さい玉はＳｉを示している。図８０の（ａ）は、ＸＹＺ空間における原子を示し、（ｂ）は、そのＸＹ面への投影を示す。ここではその軌跡まで示せなかったが、この単結晶例では原子は主にＸＹ軸の中線方向に原子が優先的に運動していることも分かった。このことは上記式に示した
Ｍ（ｔ）＝Σｅ_i・ｒ（ｔ） …（1371）
の部分に大きな偏りを示し、結果として、単結晶特有のスペクトルを示すことになる。

本発明者らは、測定温度を変えて、個々原子の挙動や、これを基にした固有振動の様子、さらには、歪みを加えた場合についても、シミュレーションを実行した。図８１及び図８２はそれらの例を示す。即ち、図８１は本発明による新しい分子動力学シミュレータによるＸ線スペクトル予測図、図８２は、同様の応力予測結果を示す図である。

実施例２ 本発明者らは、酸素欠損がある場合の酸化膜のスペクトルの予測も行えるようにし、実際の酸化膜との対比を可能とした。実施例２として、酸素欠損酸化膜の光学的性質まで計算可能とした例を示す。

ａｂ−ｉｎｉｔｉｏ分子動力学法を使用した場合には問題にならないが、ＴＴＡＴポテンシャルやＢＫＳポテンシャルを使用した分子動力学法の場合には、酸素欠損を扱う場合に、各原子のチャージをどのように設定するのが最良であるかが問題となる。これは、非常に歪みが大きくなると、結合が伸びたり、Ｓｉ−Ｓｉ結合が生じ、一定電荷を仮定することができなくなるからである。例えば、Ｓｉ−Ｏ結合の伸びにより、実際には電荷の分配が変化し、分極が大きくなると、伸縮振動の静電エネルギーは大きくなり、ＩＲのピーク振動数は大きくなると予想される。

そこで本発明者らは、Ｒａｐｐｅ−Ｇｏｄｄａｒｄにより、アルカリや水素を対象に試みられたことのあるイオン化ポテンシャルと電気親和力を用いた平衡電荷近似（Ｃｈａｒｇｅ ｅｑｕｉｌｉｂｕｒａｔｉｏｎ ａｐｐｒｏａｃｈ）を、初めてＳｉＯ₂に適用し、電荷の補正を取り込んで静電ポテンシャルを計算し 、個々の原子の運動とＩＲを求める手法を行った。

なお、平衡電荷近似とは、系の静電エネルギーを各原子のイオン化エネルギーと各原子対の電気親和力と各原子の電荷の関数と仮定して扱い、平衡状態では各原子の原子化学ポテンシャルが等しいとの条件から、原子数に相当する数の連立方程式を解き、各原子の電荷を求める手法である。

この方法をβ−石英に用いた場合のＩＲラインスペクトルと従来例を図８３に示した。本方法を用いることにより、ピーク周波数が約一割、低周波数側にシフトした。この値は実測とよりよい精度で一致する。ａ−ＳｉＯ₂でも本方法を用 いることにより、より一層、構造とＩＲスペクトルとが高精度で関連づけることが可能となり、実測との対比が効果的に行えるようになった。

図８４に、本発明者等が編み出した方法を分かりやすく概念図で示した。また探索手法を図８５に示した。

なお、本発明は、上記実施例に限定されず。液相デバイスや超電導デバイス等、他のデバイスを特にクラスタ化製造装置により製造することに対して有効に適用可能である。その他、本発明の要旨を逸脱しない範囲で種々変形して実施可能である。

実施例３
ここでは、第３の実施例として、光学的スペクトルの直接予測機能についての詳細を示す。はじめにａ−ＳｉＯ₂（アモルファスＳｉＯ₂）の計算機上で作成し、その光学スペクトルを予測する方法について示す。ａ−ＳｉＯ₂は、単結晶Ｓ ｉＯ₂を計算機上で溶融した後、急冷することにより作成した。図８６は、ａ− ＳｉＯ₂の作成の際の温度設定を示す。溶融はゆっくりと行ってもよいが、ここ では温度を急上昇させて急激に溶融させた。この後、同程度の降温速度で降温した。図８７はその際の圧力の変化であり、図８８は計算領域の体積変化である。

本発明者らは、この過程において、経験的分子動力学法を用いた場合でも、ａｂ−ｉｎｉｔｉｏ分子動力学法を用いた場合でも、昇降温により温度を大きく変動させる場合には、外圧と体積の扱いが重要であることを見いだした。例えば、本実施例では、体積をほぼ一定に保って相変態が起こるように外圧を温度の関数として扱い、温度が上昇すると圧力が大きくなるように設定している。その様子が図８７，５５に示されている。

しかし、体積に全く注意を払わず、外圧を一定に保った場合には、図８９に示したような体積の大膨張が生じ、現実の高温溶融ＳｉＯ₂とは全くことなり、著 しく低密度のＳｉＯ₂となる。更に、冷却過程でも、固相に入っても有りうべか らざる低密度のＳｉＯ₂が形成される。もちろん、ポテンシャルエネルギーを調 べると非常に高く、不安定な状態であった。こうなると実測との対比は困難となってしまう。従って、外圧に十分注意を払ってａ−ＳｉＯ₂を作成することが望 ましいことをここで見出し、実行した。

図９０には、溶融時のポテンシャルエネルギーの変動を示した。外圧に注意を払っているので、ポテンシャルは振動しながらも緩やかに上昇している。急冷の際にも同様に緩やかにポテンシャルエネルギーは減少してゆく。図９１にはこのときの計算領域のセルｃ軸に垂直な面の２辺の時間変化を記した。セルは、相変態に応じて形状が変えられるように、Ｐａｒｉｎｅｌｌｏ−Ｒａｈｍａｎの手法を採用しているので、溶融によって変化している。

最初は石英が菱形であるため、２辺のなす角度は６０度となっている。ところが溶融すると、結晶構造は壊れ、各辺に均等に内圧がかかるので、２辺のなす角度は９０度に近づいている。急冷しても石英には戻らず、約９０度で直方体を保つようであった。

更に、本発明者らは、ａ−ＳｉＯ₂を作成する際、もう一点注意すべき点があ ることを見いだした。それは、ＳｉとＯのポテンシャルに関するものである。ａｂ−ｉｎｉｔｉｏでも経験的ポテンシャル（ＡＴＴＡポテンシャルもしくはＢＫＳポテンシャル）の分子動力学法でも、ＳｉＯ₂系を扱う場合には、原子の配置 の変動によるポテンシャルエネルギーの変動が非常に大きく、そのために、原子の運動を扱うには、通常のＶｅｌｒｅｔのアルゴリズムよりも高精度に原子の運動を求める必要があるということである。

図９２は、通常のＶｅｒｌｅｔのアルゴリズムで原子の運動を求めた場合のポテンシャルエネルギーの変動と運動エネルギーの変動、及び両エネルギーの和を示す。外界とのエネルギーのやりとりが無い場合、総エネルギーは一定になるべきである。しかし、通常のＶｅｒｌｅｔのアルゴリズムでは、ラグランジャンから導かれる運動方程式の積分が十分な精度で行われないため、このようにエネルギーの保存が行われない。積分精度を向上させるために、積分間隔である時間間隔を小さくとっても、状況はほぼ変わらなかった。これは調べてみると、酸素のポテンシャルが非常に急峻であることに起因する材料に特有のものであることが分かった。

そこで、酸素を正確に扱うために、Ｖｅｒｌｅｔのアルゴリズムを改良し、時時刻々の個々原子の運動の加速度の高次のまで求めて、これを積分して速度や位置を求めて見た。すると、図９３に示したように、ポテンシャルエネルギーと逆位相で運動エネルギーが変動し、総エネルギーは一定に保たれるようになった。

このようなエネルギーの正確な扱いがなされない場合、本発明者らが実行した結果によると、高温で原子の運動が激しい場合や、液体状のＳｉＯ₂などで、計算 精度の悪さが原因で、構造が無秩序に壊れ、ａ−ＳｉＯ₂の作成ができなかった 。また、ＩＲスペクトルも、全く実測と異なるものになってしまった。

以上のように、本発明者らは、ａ−ＳｉＯ₂の作成において独自に工夫をこら すことにより、初めてａ−ＳｉＯ₂を計算機上で作成することを可能とした。こ のようにして作成したａ−ＳｉＯ₂は、石英に比較して、Ｓｉ−Ｏ結合長が若干 長く、しかもＳｉ−Ｏ長は一定の範囲内に収まっていた。溶融ＳｉＯ₂では、Ｓ ｉ−Ｏ長は石英でのＳｉ−Ｏ長よりも小さな結合長のＳｉ−Ｏが多く存在していたが、冷却するとこれが無くなり、石英よりも長い結合ばかりになった。配位数を調べると、第一近接を２Ａとすると、Ｓｉは４配位が殆どで、３と５が少し混じり、それ以外は生じなかった。酸素は２配位が殆どであった。

次に、ａ−ＳｉＯ₂に静水圧もしくは一軸性の応力をかけ、種々の応力場に応 じたａ−ＳｉＯ₂の構造変化を予測計算するため、ポテンシャルエネルギーが平 衡に達するまで、温度を９００℃に保って、２ナノ秒の時間、ＳｉＯ₂を構造緩 和させた。できあがった種々のａ−ＳｉＯ₂の構造を、個々、原子位置、速度、 加速度、その一次微分値、セル基本ベクトル、セル基本ベクトル変位速度、等の数値に変換し、それぞれ保存した。

ここまで行って、種々のａ−ＳｉＯ₂構造を作成した後、それぞれに対してＩ Ｒスペクトルを求めた。ＩＲスペクトルは、ダイポールモーメントから求めた。

ダイポールモーメントＭは、下記式
Ｍ＝Σｑi ・ｒi …（1461）
から求めることができる。ここで、ｑｉは原子ｉの電荷、ｒｉは位置である。この位置を求めるとき、系が中性の場合（Σｑi ＝０）には、原点ｒ0 を何処にとってもダイポールモーメントの値は同じになる。それは、
〔数１３４〕
Ｍ＝Σｑi （ｒi −ｒ0 ）＝Σｑi ・ｒi −Σｑi ・ｒ0
＝Σｑi ・ｒi −ｒ0 Σｑｉ＝Σｑi ・ｒi
…（1462）
となるためである。また、系が非中性の場合には、ｒ0 の掛かる項が残るので、原点の取り方によってダイポールモーメントの絶対値は変化するが、フーリエ変換した場合には周波数がゼロの成分に現れるだけで、他には影響しないことが分かる。このダイポールモーメントを求めるに当たり、原子位置ｒの取り方には一点注意が必要なことに本発明者らは気づいた。それは、分子動力学法を用いる際に周期境界条件を採用しているため、セルの−ｘ方向に出ていった原子が＋ｘ方向から入ってくることに起因するダイポールモーメントの不連続な変動を無くすことである。

実際にＩＲを測定する温度では、ａ−ＳｉＯ₂は熱振動をしているのみであり 、さほど大きく拡散をしないことを考慮して、ダイポールモーメント算出の際に用いるｒは、周期境界条件を考えず、−ｘ方向にセルから出ていっても、セル外のその位置に電荷があると扱って、モーメントを計算した。つまり、個々原子運動計算の際には周期境界条件を考慮した原子位置と、周期境界条件を考慮しない原子位置の２つの原子位置を両方取り扱うことにした。このようにして求めたダイポールモーメントの時間変化の例を図９４に示した。３次元計算であるので、ｘ，ｙ，ｚ各方向のダイポールモーメントが得られた。

ＩＲスペクトルを求めるには、第一段階として、このダイポールモーメントをフーリエ変換し、ＩＲラインスペクトルであるパワースペクトルを求める。その一例を図９５に示した。この段階でＩＲスペクトルのピーク位置が予測できる。

すなわち、約１０００ｃｍ^-1に大きなピークが出現すること、約５００ｃｍ^-1に大きなピークが出現することである。この結果は、ａ−ＳｉＯ₂のＩＲ実測と一 致した。

ＩＲスペクトルを求めるには、このラインスペクトルＩ（ω）にωと温度Ｔの関数を掛けて次の式で求める。

〔数１３５〕
α（ω）＝(4π 2／3hcn）・ω・(1−exp(−ｈω／ｋＴ))・Ｉ（ω）
…（1463）
種々の応力場でのａ−ＳｉＯ₂のＩＲを求めることにより、応力場とａ−ＳｉＯ₂の構造とＩＲスペクトルの詳細（ピークシフト、ピークショウルダー等）が関連づけられた。これにより、実測との対比が可能になった。

ここで、図９６に、本発明者らが新しく開発した分子動力学シミュレータをα−ＳｉＯ₂に適用し、計算されたＩＲスペクトルを示し、図６５に本発明者らが 新しく開発した分子動力学シミュレータをａ−ＳｉＯ₂に適用し、外圧が加わっ た場合の予測されたＩＲスペクトルを示す。また、図１００及び図１０１に本発明者らによる計算手順を示す。なお、図１００におけるブロックＢ内の、個々の原子の運動の式におけるＡ、Ｂは、下記式に示される項である。

実施例４
本発明に係るシステムを用いることにより、工程中の誘起応力の大きさやその成分の予測、さらには、その応力と材料物性との詳細な関連の把握により、基板や層間膜の塑性変形等の現象を予知することができる。この実施例では、特に、応力問題を例にあげる。応力誘起をいかに予測し、プロセス変動幅を考慮しながら、その低減をはかるかを検討した。

昨今の微細素子プロセス設計には、プロセス進行時の応力を考慮することが不可欠である。応力については、例えば以下の現象が生ずる。

（１）埋込み素子分離法では、Ｓｉ基板に溝を掘り、ここにＳｉＯ₂を埋め込 む。この場合、Ｓｉ基板とＳｉＯ₂との熱膨張係数の違いにより、歪や応力が残 留し、熱履歴によっては、結晶欠陥を誘起する場合もある。従って、できるだけ残留応力を低減できる素子構造や製造プロセスを立案する必要がある。

（２）高品質の酸化膜を作成する目的から、格子間酸素濃度の低いＳｉ基板を用いることが頻繁になってきた。このような状況では、Ｓｉ基板強度が減少し、塑性変形を誘起しやすい。

（３）浅い接合層の形成維持やウエル領域の維持を目的に、アニールや酸化工程等の熱工程が、ＲＴＰの採用により、一層基板に負担が多くなってきた。

このような状況のなかで、如何に欠陥の発生を抑止するかとともに、工程途中で、モニタしながら、応力が規定値やその変動予測値よりも大きい値が観測されたなら、その原因を即座に明確にするとともに、今後の工程において、どのように修正を加えるべきかを示すことが重要である。このことは、クラスタ化ツールでは一層重要である。

さらにまた、応力の種類や同定は極めて重要である。最大応力の主応力と主面を追跡することも大切であるが、分解剪断応力に着目することも重要である。例えば、Ｓｉ（１００）基板を用いている場合、その最優先滑り面は図１０２に示すように（１１１）面である。この（１１１）面上で［１００］方向が滑り方向であるが、この方向に剪断力が働いたときが最も危険であると考えられる。

本発明者らによるシミュレーションシステムは、これらの同定の出力も可能である。このように、塑性変形に注意をはらいながら、ここではまずトレンチ部分の酸化を取り上げ、その酸化形状と応力分布、さらには所定の形状を得るための設計指針について記述してみる。例えば、トレンチ角の酸化の「尖り」や「丸め」現象は、応力が寄与していると言われているが、まだその仕組みが分かっていないのが現状である。

これらの現象をつぶさに理解するとともに、上記「尖り」を抑え「丸め」を進行させる的確なプロセス設計を進めることがまず必要である。また、誘起応力を的確に把握しなければ酸化膜形状も正確には把握できない。応力の問題は、単に塑性変形上の問題だだけでなく、形状把握、さらには、応力場下での不純物分布把握にも欠かせない項目でもある。また、今後の１Ｇ以降の素子で展開されるであろうＳＯＩ基板のＬＯＣＯＳ工程後の酸化膜厚の予測等は、通常のＳｉ基板と異なり、酸化進行時にとりわけ応力の影響を受け易いが、その定量的予測はもっと複雑な問題である。

このような背景において、従来の酸化拡散モデルに代わり、応力下での酸化剤の挙動把握や、応力下での酸化現象の理解、応力下での不純物再拡散にまで作業を広げなければならない。しかも短時間で、高速に計算されなければならない。

本発明者等は、応力の発生とその再分布の予測、さらには応力下での点欠陥２次元拡散部分まで計算可能な新しい高速汎用プロセスシミュレータ部分を作成し、これを、原子レベルシミュレータに付加した。即ち、独自開発による高速２次元酸化／拡散／応力／変形プロセスシミュレータ部を、原子レベル材料設計システムの母体に付加した。

応力のその場観察の実測は、顕微ラマン及びＦＴ−ＩＲを用いた。これらの計測器のキャリブレーションは、あらかじめ行った。ＳｉやＳｉＯ₂に関する物性 諸量は、工程表が入力された時点で、あらかじめ原子レベル材料設計システムにより即座に実行され、求められる。これが、高速２次元酸化／拡散／応力／変形プロセスシミュレータ部に送られる。

本発明者らは、まず、この新しい高速２次元酸化／拡散／応力／変形プロセスシミュレータ部を用いて突き止めた「尖り」と「丸め」酸化の現象分析を記述し、更には、適用例として、埋め込み素子分離工程の基礎検討結果等を示す。Ｓｉが酸化されて、ＳｉＯ₂になると、体積膨張を起こす。このとき、ＳｉＯ₂には粘弾性があり、その応力の一部は緩和されるものの、依然ＳｉとＳｉＯ₂の両方に は応力が残留する。

トレンチコーナーの酸化を例にとり、このような状況を調べ、図示したのが図１０３である。まず、濃度Ｃの酸化剤がＳｉＯ₂に到達し、ＳｉＯ₂膜中に侵入して行くが、この時、酸化剤の拡散Ｄ´はＳｉＯ₂膜中の応力σの影響を受ける、 さらに、Ｓｉ／ＳｉＯ₂界面で、酸化剤がＳｉに反応する時にも基板内の残留蓄 積応力の影響を受ける。

具体的な効果としては、例えば、ＳｉＯ₂中に圧縮応力σが存在すると、ネッ トワークの原子間隔が狭まり、酸化剤の拡散を抑制することになる。また、Ｓｉ基板の面上に圧縮応力σが蓄積すると、ここでも、Ｓｉ原子間隔が狭まり、酸化剤の侵入を抑え、結果的に酸化速度を抑えることになる。

今、主応力Ｓのｘ方向の成分をσxxとし、また、ｙ方向の成分をσyyとすると、２次元面上の応力に関しては、下記式（１５３１）〜（１５３４）の関係式が成立する。ここにμは粘性係数、νはポアソン比であり、ｕは変位である。

ところで、一般に粘性率μの温度依存性は上記式（１５３４）で表現できる。

ここで、ＶＩＳＣ．Ｏ、ＶＩＳＣ．Ｅは、それぞれ比例定数及び活性化エネルギーを示す。したがって、上記角応力式にμの温度依存性を加味すると、温度上昇とともに、各主応力の成分が減少することがわかる。また、Ｓｉ／ＳｉＯ₂ 界面にまで到達した酸化剤がＳｉ表面で反応する時の反応定数κs ′は、同様に上記式（５６）のように書き表せる。ここにκs は、応力が存在しない場合の値である。ここで、σn は、σn ＝−（σxxｎx 2 ＋σyyｎy 2 ＋２σxyｎx ｎy ）であり、στ＝−（σxxｎy 2 ＋σyyｎx 2 −２σxyｎx ｎy ）である。ＶＲはパラメータである。この式から、圧縮応力が増大するとともに、酸化反応定数が減少することがわかる。上記のＳｉ／ＳｉＯ₂界面で、酸化剤がＳｉに反応する時 の基板内の残留蓄積応力の影響がこれである。

さらに、酸化膜中の酸化剤の拡散係数Ｄ′は、上記式（１５３１）で表現できる。ここに、ｐは次式で表現される。即ちｐ＝−１／２（σxx＋σyy）である。またＤは応力が存在しない場合の拡散係数である。この式の意味合いは、以下の様に解釈できる。すなわち、応力σxxやσyyは、引っ張りを正、圧縮を負に取っている。従って、膜内にもし圧縮応力が存在し、物質が密な状態であれば、拡散係数は小さくなることを示している。上記式は、このことを表現したものである。

本発明者らが見いだした、この部分のみの計算手順を図１０４に示す。まず、Ｓｉ／ＳｉＯ₂界面での各節点での酸化剤の濃度を求める。しかる後に、この濃 度を基に、１回目の各節点での仮の酸化速度を算出する。その上で、各節点でのＳｉＯ₂膜の成長量ｕを求める。このｕを用いて、上記式（５７）から応力Ｓを 求める。この応力を用いて、式（５６）、式（１５３１）、さらに式（１５３４）を用いて酸化速度を再度計算し直す。そして再び、成長量を求め、応力を算出する。そして、十分応力が収束すると、改めて、新しい酸化時間の進行を始める。即ちニュートン法を用いているわけである。

この図１０４からも分かるように、酸化時の応力計算は非常に長時間の演算時間を要する。典型的なトレンチ構造で、初期節点１３００点で、９００℃の酸化を試みた。酸化雰囲気はドライ酸素１００％とし、また酸化時間は３００分とした。ＥＷＳ Ｓｕｎ４を用いて、計算時間を調べてみた。従来技術であれば、概ね１５９００分の時間を要する。これはほぼ１０日間にも達し、このままでは、実用には極めて不向きである。

本発明者等は、ニュートン法部分の計算に関しては、偏分原理を導入し、前回の収束状況を保存しておき、これの１次微分を利用して収束計算を進めた。この方法により、９００℃〜１１００℃を例にとり、従来手法の１０００分の１以下に短くすることが出来た。これにより、本シミュレータは、材料物性値の適正化と、演算部の高速化を進め、実時間用ツールとして利用できるところまで到達した。本発明者らはさらに、シミュレータの基本精度として、酸化膜厚の時間依存性値の精度や、トレンチ角部の形状再現性等を調べてみた。

図１０５〜７７は、上記シミュレータを用いて、成長膜厚の温度、時間依存性を調べた結果をグラフにまとめたものである。ここには表示しなかったが、本発明者らは、酸素分圧が１０％やさらに１％の時の減圧挙動についても調べた。その結果、酸化速度が遅いときには、酸化膜の膨張に伴う応力が十分緩和しながら進行するので、基板には応力が蓄積せず、従って、凸部分にも均一な酸化膜が形成される。本発明者らの慎重な検討の結果、特に酸化初期に、酸素分圧をさげて、酸化速度を遅くして、十分緩和させながら徐々に酸素分圧を挙げて行く方法が有効であることがわかった。

本発明者らは、本システムを用い、トレンチ角の「尖り」と「丸め」がなぜ起きるのかを調べた。既に述べたように、トレンチ角の酸化ではいくつかの要素を考慮する必要がある。以下に、対応する５個の重要要素を記す。

（１）酸化膜粘弾性モデル
（２）シリコン弾性モデル
（３）酸化剤拡散の応力依存性モデル
（４）酸化速度の面方位依存性モデル
（５）酸化剤／Ｓｉ反応のＳｉ表面応力依存性モデル
本発明者らはこれら５個のモデルの内、もっとも支配的なものは、意外にも第５のモデル即ち、Ｓｉ／ＳｉＯ₂界面でのＳｉ基板の応力であることを突き止め た。ここで、９００℃でのドライ酸化を取り上げることにする。この時の酸化膜の形状及びＳｉ／ＳｉＯ₂界面の形状をよくみると、殆どまだ「尖り」形状は認 められない。９００℃ドライ酸化で３００分経過後になると、９００℃近傍の酸化温度では、トレンチ角が尖る現象を次のように分析した。

即ち、１００分程度までであれば、さほど応力は発生しない。しかし、酸化とともに、ＳｉＯ₂体積膨張により、圧縮応力が角部のＳｉＯ₂膜内とＳｉ先端部側に集中する。とりわけ、９００℃近傍では、粘性流動効果が小さく、酸化膜の形状緩和効果が小さい。そして、Ｓｉ表面での圧縮応力により、酸化反応が抑止され、結果として尖るわけである。しかし、１１００℃になると、粘性流動効果がおおきく、酸化膜の形状緩和効果がおおきい。そして、Ｓｉ表面での圧縮応力の蓄積がなく、結果としてまるくなる。この状況を示したのが図１０８である。これらの応力の予測値は、ラマン測定やＩＲ測定の結果、±１０％以内の精度で測定確認できた。

以上の知見を準備しておき、さらに本発明者らは、上記シミュレーションシステムを、素子分離領域の面積低減化を目的に、ＳＴＩ（Ｓｈａｌｌｏｗ Ｔｒｅｎｃｈ Ｉｓｏｌａｒｉｏｎ）の基礎検討の一部に適用した。ここでは特に、トレンチの上部角の丸め工程の基礎検討を進めた。酸化温度は１１００℃、酸素分圧は３％で、１００分経過後について見た。その結果、以下のことがわかる。

（１）はじめは、ｐｏｌｙ Ｓｉの酸化と、バッファ酸化膜の後退部分の隙間が酸化される。この時、角部は尖り気味である。

（２）引き続き酸化が進み、隙間が酸化され、埋め尽くされると、酸化剤の供給が急激に減少し、側面からの酸化が律速になり、角部は結果的には丸みを帯びてくる。また主応力は隙間近傍が大きい。

これに対して、１０００℃酸素分圧５０％での酸化時の応力分布について検討した。時間は３０分経過後である。その結果、先の１１００℃酸化の応力分布に比べ、若干応力が大きくなっていることがわかる。以上の考察から、トレンチ角部の尖りを抑止するには、Ｓｉ表面に残留応力を蓄積させないことが第一の条件であることがわかる。このための基本的な施策としては、低温でも十分酸化速度を抑え、酸化膜の変形緩和を十分起こさせること等が考えられる。

本発明者らは、さらに塑性変形シミュレータとの合体化を進めた。その結果、酸化のみならず、ＴＥＯＳ膜堆積後の変形等も調べられることがわかった。

以上の結果をもとに、本発明者らは、特にすべり面である（１１１）面上での分解剪断応力に着目し、２次元シミュレーションシステムでは（１１０）面上に投影した値をみた。また熱履歴として、３００度／分であれば良いが、ウエハ内で熱不均一が発生し、部分的にもし、６００度／分にまで達すると、部位によっては、応力により、転位が発生することも突き止めた。

以上の予備知見を別途用意しておいて、実時間でのバイポーラＣＭＯＳの製造工程について本発明を適用した。すなわち、シミュレーションによる逐次予測と、各工程での必要最小限度のその場測定量との対比、さらには、各工程の因子のばらつき量と許容量について検討した。ここではとくに、犠牲酸化の酸化時間ばらつきと、形状等の許容量とについて実施した例を以下に示す。

基本工程は、図１０９に示す通りである。さらに引き続き、素子分離用ＲＩＥエッチング→酸化→多結晶Ｓｉ埋め込み→ＣＭＰ→多結晶Ｓｉ部分の犠牲酸化→の各工程がある。それぞれの工程の因子に変動幅がある。この実施例では、多結晶Ｓｉ部分の犠牲酸化の部分の特に酸化量に着目した。酸化量の変動幅として、酸化量零の場合と所定量の２つについて実時間計算結果を示す。酸化量と、その途上の酸化形状、応力分布さらには、この後に続く、ＬＯＣＯＳ工程での酸化膜形状及びその応力の成分特性の変化を示す。

図１１０は、１０００℃、ＬＯＣＯＳ工程終了時のσxxの分布図、図１１１は、同じく１０００℃、ＬＯＣＯＳ工程終了時のσyyの分布図、また図１１２は、同じく１０００℃、ＬＯＣＯＳ工程終了時のσxyの分布図、図１１３はＬＯＣＯＳ工程終了後、２５℃に温度を下げた時のσxxの分布図、図１１４は、同じくＬＯＣＯＳ工程終了後、２５℃に温度を下げた時のσyyの分布図、また図１１５は、同じくＬＯＣＯＳ工程終了後、２５℃に温度を下げた時のσxyの分布図をそれぞれ示す。

これらの図で解るように、ＣＭＰ後、多結晶Ｓｉ膜に対する犠牲酸化が全くない場合は、多結晶Ｓｉ膜の表面と窒化膜の表面とが揃っている。この状態でＬＯＣＯＳ工程を行うと、酸化によって膨張した酸化膜が窒化膜の横から押すことになる。一方、犠牲酸化を若干なりとも行った場合は、上手く窒化膜が曲がって行くので、応力が緩和して行く。その結果を図１１６に示した。

また、本発明者らは、本発明者らが提唱してきた上記、経験的原子間ポテンシャルを用いた高速分子動力学シミュレータと汎用プロセスシミュレータとを合体させたシミュレーションシステムを用いて、酸化工程や、層間絶縁膜形成工程、さらには、配線工程を経た後に、ＭＯＳ素子のＳｉ基板や酸化膜に残留した応力が、酸化膜の電気的特性に特に大きな影響を及ぼすことを突き止めた。

本発明者らは、上記シミュレータを用いて、逐次応力生成工程と生成部位をリアルタイムで予測するとともに、例えば、ＥＥＰＲＯＭ素子部分のトンネル酸化膜が最終的にどのような保持特性等を示すかも出力させた。さらに、このシステムを用いることにより、最終的に得られる素子の保持特性の許容幅にたいして、現在進行中の素子がどのような特性を与えるかも予測できる。これが、本発明になるシステムの効力である。さらに、工程進捗中に、上手くこのシステムを利用し、最適な特性を得ることができる様に、引き続く工程を修正あるいは回避したりすることができる。この実施例を以下に示す。

実施例５
この実施例は、本発明をシリコン熱酸化膜中の応力問題に適用した例を示す。

応力問題は、シリコン基板中での結晶欠陥や転位の発生を抑制するために重要であるが、それだけではなく、様々な材料において、膜質に影響を与える要因としても重要である。

シリコン熱酸化膜には、シリコンが酸化される際の体積膨張に起因した酸化応力と、昇降温によって生じる熱応力とが生じており、これらが酸化膜質に影響を与えている。すなわち、原子レベルで言うと、応力はシリコン酸化膜のネットワーク構造に影響を与え、ＳｉＯ₂系の単位構造であるテトラヘドロン構造にまで 影響する可能性がある。引いては、電気的な絶縁性にも影響をもたらす。そこで、本発明者らは、シリコン酸化膜におけるＩＲスペクトルや原子レベルでの構造、及び電気的特性への知見を準備しつつ、酸化膜質を制御しながら酸化工程を進めた。

最初に、本発明者らは、前段階として、シリコン酸化膜のＩＲスペクトルと原子レベルでの構造及び電気的特性を調べた。図１１７に示すように、シリコン（１００）基板を用意し、素子分離工程を行って、素子分離領域により分離された素子領域を作成した。次に、シリコン基板を酸化装置に導入し、図１１８に示すようにゲート熱酸化膜を形成した。この時、酸化温度を８００℃から１１００℃まで種々変化させた。

所定の膜厚のシリコン酸化膜を形成した後、基板温度を降温させた。この時、基板温度の降温は、図１１９に示すように２段階で行った。第１段階は基板温度を９００℃まで降温する工程であり、基板温度を急冷させた。第２段階は９００℃から室温まで降温させる工程であり、降温速度には特に留意はしなかった。酸化温度を９００℃以下に設定した場合には、第２段階だけで降温した。

この降温条件の設定は、別途、酸化膜の粘性緩和挙動の解析から決定した。一般的に、酸化膜は高温では柔らかくなり、粘性流動によって応力が緩和しやすくなると言われている。そこで、本発明者らは、粘性による応力緩和挙動を詳細に調べた。降温速度による応力緩和の特性の変化を調べた例を、図１２０及び図１２１に示した。本例では、酸化温度は１０００℃とした。

図１２０は典型的な酸化工程を行った後、それぞれ１０００℃で保持した場合、６００℃／ｍｉｎで急冷した場合、６℃／ｍｉｎで徐冷した場合における、表面に沿った方向の主応力の最大値を調べた結果である。図１２０に示す結果から、１０００℃保持の場合と徐冷の場合では約１分で応力が緩和し、急冷では粘性緩和は僅かで、熱応力が重畳されるのみであることが分かった。

図１２１は、同じ降温過程での応力挙動を、横軸を温度にしてプロットした図である。図１２１から、９５０℃以上の温度領域では応力緩和が顕著であり、粘性体もしくは粘弾性体の挙動を示し、９５０℃以下では応力緩和が殆ど起こらず、酸化膜は弾性体の挙動を示すことが分かった。

このような計算例から推察されるように、本実施例で行った２段階降温を用いれば、粘性緩和が生じる温度領域を可能な限り短時間で通過させることにより、粘性緩和を抑えることが可能となり、しかも酸化温度を種々変化させることにより、酸化膜とシリコン基板の熱膨張係数の差から生じる熱応力を種々設定することが出来る。

このようにして作成した試料を用い、シリコン酸化膜直下での基板表面方向の引っ張り応力を顕微ラマン分光を用いて断面から測定したところ、シリコン基板中に種々の熱応力が誘起していることが確認された。また、本発明者らは、これら基板中の応力誘起測定値と図１２２に例示した計算結果とを比較して、酸化膜中に誘起されている応力値を推察した。尚、本発明者らは、これら応力計算を行うに当たって、本発明者らが構築した分子レベルから汎用レベルまで含めた材料、工程予測シュミレーションシステムを使用した。その結果、酸化温度と酸化膜中の圧縮応力との間に、図１２３に示す関係が得られた。図１２３から、酸化温度の上昇とともに、酸化膜中の圧縮応力の最大値が直線的に増加することがわかる。

次の準備として、以上のように形成された圧縮応力を生じた酸化膜に対して、ＩＲスペクトルを測定した。また一方、この圧縮応力が付加された場合の酸化膜構造とＩＲぺクトルの計算機予測を、古典的分子動力学法を用いて行った。図１２４は、圧縮応力付加の有無による計算ＩＲぺクトルの変化を示している。

図１２５は、ＩＲスペクトルの測定値からＳｉ−Ｏ−Ｓｉ非対称伸縮ピークの応力による波数シフトを測定した結果である。実測値と計算値とを比較すると、圧縮応力でＳｉ−Ｏ−Ｓｉの非対称伸縮ピークの出現波数が、応力付加無しの場合よりも低波数側にシフトする様子が計算ＩＲスペクトルから得られており、実験と計算は良く一致する事が分かった。そこで、応力を様々に付加した場合の非晶質構造とＩＲフルスペクトルを計算し、準備した。

さらに次の準備として、応力と電気特性との関連を調べた。前述の種々の圧縮応力を生じている酸化膜をゲート絶縁膜として使用し、その上にゲート電極を堆積して、Ｑｂｄの測定を行った。ゲ一卜電極の堆積等は、酸化膜の粘性緩和が生じないように、９００℃以下の温度で全ての工程を行つた。そして、前述の分子動力学計算から予測された酸化膜構造とＱｂｄとの相関を調べた。

図１２６は、種々の応力場での酸化膜のＳｉ−Ｏ原子の動径分布関数を算出した後、第一近接のＳｉ−Ｏ原子間隔が１．８オングストローム以上となっている結合の出現頻度とＱｂｄとの相関を求めた結果である。図１２６から、Ｓｉ−Ｏ原子間隔が１．８Ａ以上となっている出現頻度が５％を越えると、Ｑｂｄが顕著に減少しており、Ｓｉ−０原子間隔と酸化膜絶縁破壊とが関連していることが分かった。また、ＩＲスペクトルシフト量は約１５ｃｍ^-1となつていた。

酸化膜に圧縮応力が加わった場合の非晶質ネフトワーク構造の変化は、（１）酸化膜の密度が増加することから、Ｓｉ−０結合長は減少するという予想と、（２）薄膜熱酸化膜のＩＲ測定と中心カモデルから、Ｓｉ−Ｏ−Ｓｉ結合角度が減少するという予想の２通りの予想がなされてきた。しかしながら、本発明者らは、独自に開発したシミュレータにより、図１２７に示したように、酸化膜の密度増加はＳｉ−Ｏ結合長の単純減少にはつながらず、むしろ、結合長分布が広がり、伸びた結合の割合が増加することを新たに見い出した。加えて、Ｑｂｄとの相関まで見い出した。

Ｑｂｄと伸びたＳｉ−Ｏ結合の相関がある物理的理由については、本発明者らはまだ詳細を解明していないが、その理由のーつを最近報告された論文から以下のように推察する。引用論文は、「薄膜・表面物理分科会研究報告、極薄シリコン酸化膜の形成・評価・信頼性」で金田氏によって報告された「酸化膜中の局所的変形および不純物による正孔捕獲」ｐ．１０１−１０４である。金田氏が分子軌道法プログラムＧＡＵＳＳＩＡＮを使用して計算したところによると、「伸びたＳｉ−Ｏ−Ｓｉがあれば、その中心の酸素への正孔の捕獲が誘発される。捕獲された正孔は酸素に強く局在するので、電子に対しては捕獲中心として働く。この結果、酸素の多くは電子を捕獲して元の状態に戻る。しかし、一部は電子と正孔の再結合で放出された数ｅＶのエネルギーを得て、伸びて弱くなっているＳｉ−Ｏ結合を切って格子間へ飛び出し、例えば酸素欠損のような２次的な欠陥を生成しうる。」と報告されている。

本発明者らによると、この報告のように、応力付加により伸びた結合が増加し、酸素欠損欠陥が形成されることが予想される。さらに、酸素欠損欠陥には様々な構造が考えられるが、酸素欠損両端のＳｉ原子が大きく引き離された構造では、電子状態も大きく変化して、バンドギャップ内に準位が形成されたり、電気伝導を生じる電子分布に変化すると推測される。そして酸化膜の絶縁性が劣化するものと推察される。

以上のような準備を行うことにより、応力場中での酸化膜のＩＲスペクトルが予測でき、さらに、原子レベルでの構造と電気的特性との相関が得られたことから、本発明者らは、「製造工程の修正」を試みた。

所定工程まで進んだＳｉウエハを、酸化膜形成装置に導入し、１０５０℃、５％塩酸を含むｄｒｙ Ｏ₂雰囲気で、所望の膜厚の酸化膜を形成した。酸化膜形 成後に室温までウエハを急冷し、形成した酸化膜のＩＲスペクトルを測定した。測定値より、ＩＲスペクトルのＳｉ−Ｏ非対称伸縮振動ピーク波数を読みとり、応力が生じていない場合の同波数からのシフト量を計算したところ、１７ｃｍ^-1のシフトが測定された。

このシフト量から、酸化膜に約０．２ＧＰａの圧縮応力が酸化膜中に誘起されていることが予測された。このシフト量をシフト臨界値１５ｃｍ^-1と比較した。

この熱工程では臨界値を越えていることが検知されたため、プロセスを停止した。そして、製造工程の修正として、高温アニールによる応力緩和を行った。

また、本実施例の効果を比較するため、一部のウエハには高温アニール工程を付加せず、そのまま工程を通過させた。応力修正用高温アニールのパラメータであるアニール時間とアニール温度、降温速度等は、予め用意された粘性緩和による緩和時間の解析、及び、基板中に既に形成されている不純物拡散層の再拡散挙動を考慮して算出した。この場合は１０００℃、３０秒のアニールを行った。

アニール終了後、再度ＩＲスペクトルを測定し、シフト量を算出したところ、臨界値を下回ったので、そのまま、プロセスを継続した。プロセスを修正した場合と、修正しなかった場合とで、Ｑbdの測定を行った。その結果、修正した場合では、４５Ｃ／ｃｍ² 、修正しなかった場合で５Ｃ／ｃｍ² となり、本発明によるプロセス修正の効果が明確に得られた。

また、別の構造を有するウエハにおいて、９５０℃のｄｒｙ Ｏ₂雰囲気で所 望の膜厚のゲート酸化膜を形成した。酸化膜形成後に室温までウエハを急冷し、形成した酸化膜のＩＲスペクトルを測定したところ、ピークシフト量が１６ｃｍ^-1であった。これは臨界値を越えていたため、応力緩和のためのアニールの追加が必要であったが、不純物拡散層の再拡散を考え合わせると、適切なアニール条件の解が存在しなかった。そこで、このウエハは工程継続が不可能と判断し、ウエハを廃棄した。

本発明によるシステムの全体を示す図。 ポテンシャルの大きさのそれぞれの方向に対する傾きが、それぞれの方向の力となることを示す図。 本発明に用いた新しいａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムを示す図。 第一の原子（ｉ）と、第二の原子（ｊ）との作用に、さらに第三の原子（ｋ）から第二（ｊ）を通して入ってくる効果を示す図。 ｒを粒子間の距離とし、θを粒子角度とし、これらの諸量の主従の関係を示す図。 従来の計算手法と本発明の高階巡回偏微分法との位置づけを示す図。 分子動力学シミュレータＩＩと、プロセス変動算出ユニットＩとの連携を説明する図。 分子動力学シミュレータＩＩと、プロセス変動算出ユニットＩとの連携を説明する図。 汎用２・３次元プロセスシミュレータの入力シーケンスの一例を示すを示す図。 汎用２・３次元プロセスシミュレータの入力シーケンスの一例を示すを示す図。 本件発明の実施例による２進木の一部を示す図。 イオン化時のクーロン相互作用の模式図を示す。 イオン化時のクーロン相互作用の模式図を示す。 イオン化時のクーロン相互作用の模式図を示す。 簡単な原子構造モデルをを示す。 アルカリ金属ハライド分子を示す。 原子価軌道指数ζの表を示す。 マトリックス(３０６３)を解いて求めた電荷を示す。 マトリックス(３０６３)を解いて求めた電荷を示す。 マトリックス(３０６３)を解いて求めた電荷を示す。 Si-Si間のクーロン相互作用Jの計算結果を示す。 Oの電荷が2sにあり、Siが2sあるいは3sにある状態をモデル化した場合のJABの計算結果を示す。 結晶系及び非晶質系に新たに開発した上記電荷計算法の手続きをまとめたものを示す。 本発明による分子動力学シミュレータの実行結果であるＯ−Ｓｉ−Ｏ角、Ｓｉ−Ｏ−Ｓｉ角、及びＯ−Ｓｉ距離を示すデータ図。 従来のシミュレーションシステムにおける変分計算を含む各サブルーチンを示す図。 従来のシミュレーションシステムにおける変分計算の手順を示す図。 従来のプロセスと変動プロセスの収束状況を示す特性図。 変分の精度を示すデータ図。 変分の精度を示すデータ図。 変分の精度を示すデータ図。 従来のシミュレーションシステムによる不純物分布とその変動例を示す図。 従来のシミュレーションシステムによるプロセスの変動とその実行組み合わせを示す図。 従来のシミュレーションシステムによる推論エンジンの探索の木の一部を示す図。 本発明者らが開発した新しいabinitioＭＤシミュレータによるペロブスカイトとＳｉ系結晶の最適化の図。 本発明者等による独自開発によるシミュレータの構成概念図。 本発明者等による独自開発によるシミュレータの流れ図。 本発明者等が開発したシミュレータによるＰＴＯペロブスカイトのＣ４ｖ対象例の図。 ペロブスカイトＢａＴｉＯ₃とＰｂＴｉＯ₃の構成図。 本発明者等が新しく作成したシミュレータによるペロブスカイトＰＴＯのポテンシャルの出力図。 本発明者等が開発したシミュレータにおけるSiO₂の作成例。 本発明者等が新しく作成したシミュレータによるペロブスカイトＰＴＯのチャージの出力図。 本発明者等が新しく作成したabinitio/ＭＤシミュレータによるペロブスカイ トＰＴＯの歪みの例。 本発明者等が新しく作成したabinitio/ＭＤシミュレータによるペロブスカイ トＰＴＯの酸素欠損の例。 本発明者等が新しく作成したabinitio/ＭＤシミュレータによるペロブスカイ トＰＴＯのＰｂの移動の例。 本発明者等が新しく作成したabinitio/ＭＤシミュレータによるペロブスカイ トＰＴＯのＰｂの移動の例。 本発明者らが作成したシミュレータのデータの一部。 ソーヤータウワの図。 本発明の実施例の図。 本発明に使用されるクラスタ化ツールにおいて、ｉｎ−ｓｉｔｕモニタ測定装置がどのようにして付加されているかを示す概念図。 限られた測定量と、他方計算から得られた十分な情報とを対比させる方法を、モニター上のグラフィックな画像として表示される様子を示す図。 本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムへの入力データの一例としての温度シーケンスを示す図。 本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムによる原子レベルの実行予測を示す概念図。 本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムによる原子レベル実行予測を可視化した概念図。 本発明によるａｂｉｎｉｔｉｏ分子動力学シミュレーションシステムによるＳｉＯ₂についての実行例であるＳｉ−Ｏ−Ｓｉ角の経時変化を示す図。 図５１の入力シーケンスに対する出力例である温度の測定値変化の予測図。 本発明によるシミュレーションシステムによる原子レベル実行例である成膜Ｓｉ中の配位角の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ配位角の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例。Ｓｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ配位角の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ配位角の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ配位角の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によるシミュレーションシステムによる原子レベル実行例であるＳｉ−Ｓｉ間距離の分布図。 本発明によって構成された別のシステムの全体を示す図。 従来のラマン測定結果を示す特性図。 従来のラマン測定結果を示す特性図。 従来のラマン測定結果を示す特性図。 従来のラマン測定結果を説明するためのモデルを示す図。 Ｏ−Ｓｉ−Ｏ結合角の関数としてのＳｉＯ₄四面体の通常モード振動の波長を 示す特性図及びストレスにより生じた強度の減少を示す特性図。 ＳｉＯ₂中の欠陥を示すモデル図。 本発明による分子動力学シミュレータの実行例を示す特性図。 本発明による分子動力学シミュレータの実行例を示す特性図。 本発明による分子動力学シミュレータの実行例を示す特性図。 本発明による分子動力学シミュレータの実行例である原子の動きを示す図。 本発明による分子動力学シミュレータの実行例であるＸ線スペクトル予測結果を示す特性図。 本発明による分子動力学シミュレータの実行例である応力予測結果を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータによる算出結果を示す特性図。 本発明による変分手法を用いた分子動力学シミュレータの変動を含む一貫計算を示す図。 本発明による分子動力学シミュレータによる探索手法を示す図。 本発明者らによる非晶質ＳｉＯ₂作成時の温度設定の様子を示す特性図。 本発明者らによる計算機中での非晶質ＳｉＯ₂作成時の圧力変化を示す特性図 。 本発明者らが作成した分子動力学シミュレータによる非晶質ＳｉＯ₂作成時の 体積変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータによる非晶質ＳｉＯ₂作 成時の体積変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータによる非晶質ＳｉＯ₂作 成時のポテンシャルエネルギー変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータによる非晶質ＳｉＯ₂作 成時の計算セル中の２辺の時間変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータに通常のVerletのアルゴリズムを用い、計算したエネルギーの変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータに改良型Verletのアルゴリズムを用い、計算したエネルギーの変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータを用い計算したダイポールモーメントの時間変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータを用い、ダイポールモーメントから更に求めたパワースペクトルを示す特性図。 本発明者らが新しく開発した分子動力学シミュレータをα−ＳｉＯ2 に適用し、計算されたＩＲスペクトルを示す図。 本発明者らが新しく開発した分子動力学シミュレータをα−ＳｉＯ2 に適用し、外圧が加わった場合の予測されたＩＲスペクトルを示す図。 本発明者らが新しく開発した分子動力学シミュレータによる非晶質ＳｉＯ、作成時の体積変化を示す特性図。 本発明者らが新しく開発した分子動力学シミュレータに通常のＶｅｒｌｅｔのアルゴリズムを用い、計算したエネルギーの変化を示す特性図。 本発明者らによる計算手順を示す図。 本発明者らによる計算手順を示す図。 Ｓｉ（１００）における滑り面、分解剪断応力等を示す図。 ＳｉＯ₂及びＳｉ／ＳｉＯ₂界面近傍での反応を示す模式図。 粘弾性計算の逐次計算を模式的に示す図。 酸化膜の膜厚の温度、時間依存性を示す特性図。 酸化膜の膜厚の温度、時間依存性を示す特性図。 酸化膜の膜厚の温度、時間依存性を示す特性図。 合わせ込まれた係数を用いた場合のトレンチ酸化における主応力分布及び最大剪断応力を示す図。 典型的な熱工程における温度と蓄積応力の変化を示す特性図。 トレンチ角におけるσxxの分布を示す計算予測図。 トレンチ角におけるσxxの分布を示す計算予測図。 トレンチ角におけるσxxの分布を示す計算予測図。 トレンチ角におけるσxxの分布を示す計算予測図。 トレンチ角におけるσxxの分布を示す計算予測図。 トレンチ角におけるσxxの分布を示す計算予測図。 最適化された酸化予測の形状を示す図。 実施例７におけるゲート熱酸化膜の形成工程を示す断面図。 実施例７におけるゲート熱酸化膜の形成工程を示す断面図。 実施例７におけるゲート熱酸化膜形成後の昇温工程を示す特性図。 昇温工程の相違による応力の変化を示す特性図 昇温工程の相違による応力の変化を示す特性図。 応力分布の計算例を示す図。 実施例７における酸化温度と膜中の圧縮応力との関係を示す図。 応力付加によるＩＲスペクトルの変化の計算結果を示す図。 酸化膜中の圧縮応力とＳｉ−Ｏ−Ｓｉ非対称伸縮ピークの波数シフトとの関係を示す図。 第１近傍のＳｉ−Ｏ原子間隔が１．８オングストローム以上となっている結合の出現頻度とＱbdとの関係を示す図。 応力付加の有無によるＳｉ−Ｏ原子間隔分布の変化を示す図。 ｄＲＡＭの世代とチップコストとの関係を示す図。 従来例による半導体装置の製造工程を示す模式図。 従来例によるクラスタ化ツールを示す図。

Claims (30)

1. 半導体装置の製造プロセスの製造装置内物理現象のシミュレーションを用いた半導体装置の製造プロセスの設計方法において、前記製造装置内に形成された電子材料の特性に関する製造工程因子を入力値とし、この入力値を基に３次元のシミュレーションを行ない製造装置内に形成された電子材料の特性の予測データを得る方法であって、
   前記シミュレーションは、分子動力学シミュレーションと流体モデルに基づいて行われるシミュレーションとからなり、前記分子動力学シミュレーションは、分子動力学法と密度汎関数法との組み合わせで行われ、当該分子動力学法は、原子間ポテンシャルを３体問題として取り扱い、一点と微小変化させたもう一つの点のポテンシャルの差から平均変化率により力を求める手法を含む高階巡回偏微分法を用い、前記分子動力学シミュレーションの結果として出力される３次元物理量と前記流体モデルに基づいて行われるシミュレーション結果として出力される３次元物理量の少なくとも一方は、他方に転送されそこでのシミュレーションに利用されることを特徴とする半導体装置の製造プロセスの設計方法。
2. 前記分子動力学シミュレーションから前記流体モデルのプロセスシミュレーションへ、粘弾性係数、ポアソン比、拡散定数、Ｃijの内少なくとも１つ以上の物理量が転送されることを特徴とする請求項１に記載の半導体装置製造プロセスの設計方法。
3. 更に、設定値のゆらぎを入力する段階を備え、前記分子動力学シミュレーションは、前記入力値とそのゆらぎを基に変分計算を行うことを特徴とする請求項１に記載の半導体装置製造プロセスの設計方法。
4. 前記分子動力学シミュレーションは、経験的ポテンシャルを用いて行われることを特徴とする請求項１に記載の半導体装置製造プロセスの設計方法。
5. 前記経験的ポテンシャルはＡＴＴＡポテンシャルであることを特徴とする請求項４に記載の半導体装置製造プロセスの設計方法。
6. 前記経験的ポテンシャルはＢＫＳポテンシャルであることを特徴とする請求項４に記載の半導体装置製造プロセスの設計方法。
7. 半導体装置の製造工程において、製造装置内に形成された電子材料の特性を測定して実観測データを得る工程と、前記製造装置内に形成された電子材料の特性に関する製造工程因子を入力値とし、この入力値を基に３次元のシミュレーションを行ない前記製造装置内に形成された電子材料の特性の予測データを得る工程と、前記予測データと実観測データとを比較検定する工程と、前記比較検定により、前記製造工程因子の設定値と、前記実観測データから推測される前記製造工程因子との間に有意差が認められた場合、前記製造工程因子を修正処理する工程とを備え、
   前記シミュレーションは、分子動力学シミュレーションと流体モデルに基づいて行われるシミュレーションとからなり、分子動力学シミュレーションは、分子動力学法と密度汎関数法との組み合わせで行われ、当該分子動力学法は、原子間ポテンシャルを３体問題として取り扱い、一点と微小変化させたもう一つの点のポテンシャルの差から平均変化率により力を求める手法を含む高階巡回偏微分法を用い、前記分子動力学シミュレーションの結果として出力される３次元物理量の少なくとも一方は、他方に転送されそこでのシミュレーションに利用されることを特徴とする半導体装置の製造方法。
8. 前記実観測データから推測される前記製造工程因子との間の有意差が許容範囲を越えていて、修正処理が不能である場合には、不良品として処理することを特徴とする請求項７に記載の半導体装置の製造方法。
9. 前記製造装置はクラスタ化していることを特徴とする請求項７に記載の半導体装置の製造方法。
10. 前記クラスタ化した製造装置は枚葉式であることを特徴とする請求項９に記載の半導体装置の製造方法。
11. 前記分子動力学シミュレーションから前記流体モデルのプロセスシミュレーションへ、粘弾性係数、ポアソン比、拡散定数、Ｃijの内少なくとも１つ以上の物理量が転送されることを特徴とする請求項７に記載の半導体装置の製造方法。
12. 前記特性は、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも１つに関する特性であることを特徴とする請求項７に記載の半導体装置の製造方法。
13. 前記予測データと実観測データとの比較検定は逐次実時間で行われることを特徴とする請求項７に記載の半導体装置の製造方法。
14. 前記製造工程因子を修正処理する工程は逐次実時間で行われることを特徴とする請求項７に記載の半導体装置の製造方法。
15. 前記分子動力学シミュレーションは、経験的ポテンシャルを用いて行われることを特徴とする請求項７に記載の半導体装置の製造方法。
16. 前記経験的ポテンシャルはＡＴＴＡポテンシャルであることを特徴とする請求項１５に記載の半導体装置の製造方法。
17. 前記経験的ポテンシャルはＢＫＳポテンシャルであることを特徴とする請求項１５に記載の半導体装置の製造方法。
18. 前記分子動力学シミュレーションで得られた前記特性の少なくとも１つの予測データを用いて、更に流体モデルのプロセスシミュレーションが行われることを特徴とする請求項７に記載の半導体装置の製造方法。
19. 前記分子動力学シミュレーションと流体モデルのプロセスシミュレーションは、共に３次元シミュレーションであることを特徴とする請求項１８に記載の半導体装置の製造方法。
20. 半導体装置の製造装置と、前記製造装置内に形成された前記半導体装置の電子材料の特性を測定する測定装置と、前記半導体装置の製造工程因子の設定値を入力値とし、この入力値を基に３次元のシミュレーションを行ない前記特性の予測データを得るコンピュータシステムと、前記シミュレーションは、分子動力学シミュレーションと流体モデルに基づいて行われるシミュレーションとからなり、分子動力学シミュレーションは、分子動力学法と密度汎関数法との組み合わせで行われ、当該分子動力学法は、原子間ポテンシャルを３体問題として取り扱い、一点と微小変化させたもう一つの点のポテンシャルの差から平均変化率により力を求める手法を含む高階巡回偏微分法を用い、前記分子動力学シミュレーションの結果として出力される３次元物理量の少なくとも一方は、他方に転送されそこでのシミュレーションに利用されることを特徴とし、
    前記測定装置及び前記製造装置と前記コンピュータシステムを接続する手段とを備え、前記コンピュータシステムは、前記予測データと実観測データとを比較検定を行ない、前記比較検定により、前記製造工程因子の設定値と、前記実観測データから推測される前記製造工程因子との間に有意差が認められた場合、前記製造装置の制御手段を介して製造工程因子を修正処理することを特徴とする半導体装置の製造システム。
21. 前記製造装置はクラスタ化していることを特徴とする請求項２０に記載の半導体装置の製造システム。
22. 前記クラスタ化した製造装置は枚葉式であることを特徴とする請求項２１に記載の半導体装置の製造システム。
23. 前記特性は、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも１つに関する特性であることを特徴とする請求項２０に記載の半導体装置の製造システム。
24. 前記予測データと実観測データとの比較検定は逐次実時間で行われることを特徴とする請求項２０に記載の半導体装置の製造システム。
25. 前記製造工程因子を修正処理する工程は逐次実時間で行われることを特徴とする請求項２０に記載の半導体装置の製造システム。
26. 前記入力値として前記設定値のゆらぎが入力され、前記分子動力学シミュレーションは、前記入力値とそのゆらぎを基に変分計算を行うことを特徴とする請求項２０に記載の半導体装置の製造システム。
27. 前記分子動力学シミュレーションは、経験的ポテンシャルを用いた方法を併用することを特徴とする請求項２０に記載の半導体装置の製造システム。
28. 前記経験的ポテンシャルはＡＴＴＡポテンシャルであることを特徴とする請求項２７に記載の半導体装置の製造システム。
29. 前記経験的ポテンシャルはＢＫＳポテンシャルであることを特徴とする請求項２７に記載の半導体装置の製造システム。
30. 前記分子動力学シミュレーションで得られた前記特性の少なくとも１つの予測データを用いて、更に流体モデルのプロセスシミュレーションが行われることを特徴とする請求項２０に記載の半導体装置の製造システム。

Images