(57)【要約】【課題】  自己相互作用力が従来の方法より小さい粒子運動シミュレーション方法を提供する。【解決手段】  複数の粒子が存在する空間を計算対象とする粒子運動シミュレーションであって、予め、粒子が存在せず場の勾配が０と仮定した計算対象の空間の中に１個の粒子を単独で置いたと仮定し、該１個の粒子に作用する場をポアソン方程式を解いて計算し、得られた場から粒子に作用する力を計算し、これを自己相互作用力とし、この自己相互作用力を計算対象の空間内の各位置において計算しておき、次に、複数の粒子が存在する空間において、ポアソン方程式を解いて各粒子に作用する場を計算し、得られた場から各粒子に作用する仮の力を計算し、この仮の力から前記自己相互作用力を差し引いて各粒子に作用する力を計算し、得られた力から各粒子の運動を計算する粒子運動シミュレーション方法。【選択図】          図３
2005301651.tif 000002
【特許請求の範囲】
【請求項１】
複数の粒子が存在する空間を計算対象とする粒子運動シミュレーションであって、予め、粒子が存在せず場の勾配が０と仮定した計算対象の空間の中に１個の粒子を単独で置いたと仮定し、該１個の粒子に作用する場をポアソン方程式を解いて計算し、得られた場から該１個の粒子に作用する力を計算し、これを該１個の粒子の位置における自己相互作用力とし、この自己相互作用力を計算対象の空間内の各位置において計算しておき、次に、複数の粒子が存在する空間において、ポアソン方程式を解いて各粒子に作用する場を計算し、得られた場から各粒子に作用する仮の力を計算し、この仮の力から粒子が存在する位置における前記自己相互作用力を差し引いて各粒子に作用する力を計算し、得られた力から各粒子の運動を計算することを特徴とする粒子運動シミュレーション方法。
【請求項２】
粒子運動シミュレーションが、半導体デバイス中の粒子運動のシミュレーションであり、自己相互作用力を、全ての電極間の電位差を０と仮定し、かつ、不純物濃度分布とバンド構造に起因する電場の勾配を０と仮定して計算し、粒子散乱のタンミングと粒子の散乱方向等を、モンテカルロ法を用いて計算する粒子運動シミュレーションである請求項１記載の方法。
【請求項３】
コンピュータによって読み取り可能な記録媒体であって、請求項１または２に記載の方法を実行するためのプログラムを格納した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の粒子が存在する空間を計算対象とし、ポアソン方程式を解くことにより粒子に作用する場を計算して行う粒子運動シミュレーション方法に関する。
【背景技術】
【０００２】
粒子運動のシミュレーション技術は、半導体デバイスの電気特性解析やプラズマ・ディスプレイ・セル内におけるプラズマの挙動解析などに適用可能な、産業上有用な技術である。
【０００３】
この粒子運動シミュレーションを実行する際、粒子同士は相互に力を及ぼし合っているので、すべての粒子同士についてこの相互作用力を算出する方法が粒子運動シミュレーション方法として理想的であるが、計算に時間が掛かり過ぎ、実用的ではないという問題点がある。
【０００４】
この問題を解決して、粒子運動シミュレーションを実用的な時間で実行できるようにするために、空間をメッシュに分割し、各粒子をメッシュに割り付けて粒子密度を計算した後、粒子密度と場の関係を表現する式（１）に示すようなポアソン方程式を立て、そのポアソン方程式を解いて計算した場から各粒子に作用する力を計算し、得られた力から各粒子の運動を計算して行う粒子運動シミュレーション方法が知られている（例えば、非特許文献１参照。）。
【０００５】
Δφ＝−４πρ      （１）
（但し、Δは微分演算子∂2/∂ｘ2＋∂2/∂ｙ2＋∂2/∂ｚ2、φは場を表現するスカラーポテンシャル、πは円周率、ρは粒子密度から計算される湧き出しの強さである。）
【０００６】
このような粒子運動シミュレーション方法は、例えば、半導体デバイス中のキャリアを粒子として取り扱い、この粒子の運動を、電場から受ける力と格子振動や不純物等から受ける散乱を考慮して計算し、電気特性を予測するモンテカルロ法半導体デバイスシミュレーションにおいては、キャリアをメッシュに割り付けてキャリア密度を計算した後、式（２）に示すポアソン方程式を、電極におけるポテンシャル等の境界条件のもとに解いて電場を求めて具体的に実行されている。
【０００７】
∇・（ε∇φ）＝−ｑ（ｐ−ｎ＋ＮD−ＮA）      （２）
（但し、∇は空間微分演算子であり、左辺は式（３）となる。
∂（ε∂φ／∂ｘ）／∂ｘ＋∂（ε∂φ／∂ｙ）／∂ｙ＋∂（ε∂φ／∂ｚ）／∂ｚ
（３）
φは静電ポテンシャル、εは誘電率、ｑは電気素量、ｐはホール密度、ｎは電子密度、
ＮDはドナー型イオン化不純物濃度、ＮAはアクセプター型イオン化不純物濃度である。）
【０００８】
しかし、前記シミュレーション方法によれば、ある粒子に作用する力から、その粒子以外の計算対象の粒子に起因する力と、電界や重力など計算対象空間の外からその粒子に作用する力をすべて差し引いた後にも、自己相互作用力と称される力が計算上は残存する場合があることが知られている。この自己相互作用力はシミュレーション方法に起因する計算誤差であり、自己相互作用力が小さい粒子運動シミュレーション方法が求められていた。
【０００９】
そこで、粒子を、それを囲む近接メッシュ（粒子の周囲の空間で、粒子から各々のメッシュの１辺の半分の距離の空間を含むメッシュを「近接メッシュ」と称する。）の中心に均等な比率で割り付け（図１参照。）、あとは従来と同様にして粒子密度と場の関係を表現するポアソン方程式を解いて場を計算し、粒子を囲む近接メッシュ中心間の場の勾配を平均して得られた値から粒子に作用する力を計算して粒子運動シミュレーションを行うＮＥＣ法（Ｎｅａｒｅｓｔ−Ｅｌｅｍｅｎｔ−Ｃｅｎｔｅｒ法）が提案されており（例えば、非特許文献２参照。）、自己相互作用力が小さくなることが知られている。例えば、図１に示す２次元の場合、粒子は、メッシュＡ、Ｂ、Ｃ、Ｄの中心に囲まれており、メッシュＡ、Ｂ、Ｃ、Ｄに１／４ずつ割り付けて粒子密度を計算する。
【００１０】
しかしながら、自己相互作用力の低減は十分ではなく、さらに自己相互作用力の小さくなる粒子運動シミュレーション方法が求められていた。
【００１１】
【非特許文献１】アール・ダブリュ・ホックニー（Ｒ．Ｗ．Ｈｏｃｋｎｅｙ）、ジェー・ダブリュ・イーストウッド（Ｊ．Ｗ．Ｅａｓｔｗｏｏｄ）著、「コンピューター・シミュレーション・ユージング・パーティクルズ（Ｃｏｍｐｕｔｅｒ  Ｓｉｍｕｌａｔｉｏｎ  Ｕｓｉｎｇ  Ｐａｒｔｉｃｌｅｓ）」、インスティテュート・オブ・フィジックス・パブリッシング（Ｉｎｓｔｉｔｕｔｅ  ｏｆ  Ｐｈｙｓｉｃｓ  Ｐｕｂｌｉｓｈｉｎｇ）、（英国）、１９８８年
【００１２】
【非特許文献２】エス・イー・ラウクス（Ｓ．Ｅ．Ｌａｕｘ）、「アイ・イー・イー・イー・トランザクションズ・オン・コンピュータ−エイディッド・デザイン・オブ・インテグレーティッド・サーキッツ・アンド・システムズ（ＩＥＥＥ  ＴＲＡＳＡＣＴＩＯＮＳ  ＯＮ  ＣＯＭＰＵＴＥＲ−ＡＩＤＥＤ  ＤＥＳＩＧＮ  ＯＦ  ＩＮＴＥＧＲＡＴＥＤ  ＣＩＲＣＵＩＴＳ  ＡＮＤ  ＳＹＳＴＥＭＳ）」、Ｔｈｅ  ＩＥＥＥ  Ｃｉｒｃｕｉｔｓ  ａｎｄ  Ｓｙｓｔｅｍｓ  Ｓｏｃｉｅｔｙ、（米国）、１９９６年、Ｖｏｌ．１５、ｐ．１２６６−１２７７
【発明の開示】
【発明が解決しようとする課題】
【００１３】
そこで、本発明の目的は、自己相互作用力が従来の方法より小さい粒子運動シミュレーション方法を提供することにある。
【課題を解決するための手段】
【００１４】
上記課題を解決するために、本発明者は、複数の粒子が存在する空間を計算対象とする粒子運動シミュレーションにおける力の計算方法について鋭意検討し、空間をメッシュ分割し、各粒子をメッシュに割り付けて粒子密度を計算した後、粒子密度と場の関係を表現するポアソン方程式を解くことにより粒子に作用する場を計算し、得られた場から粒子に作用する力を計算する粒子運動シミュレーションにおいて、計算対象の空間内において、粒子が存在せず場の勾配が０と仮定した空間の中に、ある位置に１個の粒子を単独で置いたと仮定し、該１個の粒子を置いたことにより生じ、該１個の粒子に作用する場をポアソン方程式を解いて計算し、得られた場から該１個の粒子に作用する力を該１個の粒子の位置における自己相互作用力として計算することに思い至り、予め、計算対象の空間内の各位置においてこの自己相互作用力を計算しておき、次に、計算対象の空間において複数粒子の運動を計算するにあたり、ポアソン方程式を解いて各粒子に作用する場を計算し、得られた場から各粒子に作用する仮の力を計算し、その仮の力から粒子が存在する位置における前記自己相互作用力を差し引くことにより各粒子に作用する力を計算し、得られた力から各粒子の運動を計算することにより、従来の方法より自己相互作用力が大幅に小さい粒子運動シミュレーションを行うことができることを見出し、本発明を完成させるに到った。
すなわち本発明は、複数の粒子が存在する空間を計算対象とする粒子運動シミュレーションであって、予め、粒子が存在せず場の勾配が０と仮定した計算対象の空間の中に１個の粒子を単独で置いたと仮定し、該１個の粒子に作用する場をポアソン方程式を解いて計算し、得られた場から該１個の粒子に作用する力を計算し、これを該１個の粒子の位置における自己相互作用力とし、この自己相互作用力を計算対象の空間内の各位置において計算しておき、次に、複数の粒子が存在する空間において、ポアソン方程式を解いて各粒子に作用する場を計算し、得られた場から各粒子に作用する仮の力を計算し、この仮の力から粒子が存在する位置における前記自己相互作用力を差し引いて各粒子に作用する力を計算し、得られた力から各粒子の運動を計算することを特徴とする粒子運動シミュレーション方法を提供する。
【発明の効果】
【００１５】
本発明の粒子運動のシミュレーション方法及びその記録媒体によれば、自己相互作用力の影響が従来の方法より極めて小さい粒子運動シミュレーションを行うことが可能である。この特徴は、計算時間と計算機記憶領域の制約から非等間隔メッシュを用いざるを得ない場合、特に有利となり、本発明の粒子運動のシミュレーション方法は、半導体デバイスの電気特性解析やプラズマ・ディスプレイ・セル内におけるプラズマの挙動解析などに幅広く適用することができるので、工業的に極めて有用である。
【発明を実施するための最良の形態】
【００１６】
以下、本発明に係わる粒子運動のシミュレーション方法、及び、粒子運動のシミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体の実施形態について、図面を参照しながら具体的に説明する。
【００１７】
第１に、本発明のシミュレーション方法は、複数の粒子が存在する空間を計算対象とする粒子運動シミュレーションであって、複数の粒子が存在する空間において、ポアソン方程式を解いて各粒子に作用する場を計算し、得られた場から各粒子に作用する力を仮の力として計算し、この仮の力から後述の自己相互作用力を差し引いて各粒子に作用する力を計算し、得られた力から各粒子の運動を計算することを特徴とする。図２により、本発明の方法の第１の特徴を具体的に説明する。
【００１８】
本発明の粒子運動シミュレーション方法においては、図２に示したように、まず、各粒子の初期状態（位置、運動量ベクトル）を生成する（ステップＳ２０１）。各粒子の初期位置の生成方法は特に限定されず、乱数を用いてランダムに生成させる方法、等間隔に配置させて生成する方法、計算対象の空間の中の一つの面など特定の位置に配置させる方法など、通常用いられ、実行する粒子運動シミュレーションに適した方法を選ぶことができる。そして、計算対象の空間をメッシュに分割する。メッシュ分割方法は特に限定されず、メッシュ同士の境界は等間隔であっても非等間隔であってもよい。
【００１９】
上記の初期設定の後、各粒子をポアソン方程式を適用するために、空間分割メッシュに割り付けて粒子密度を計算し（ステップＳ２０２）、ポアソン方程式を立てて解くことによって場を計算する（ステップ２０３）。粒子のメッシュへの割り付け方法は、粒子が属するメッシュの中心に割り付ける方法、粒子を囲む近接メッシュ中心に定められた比率で割り付ける方法等のいずれの割り付け方法を用いてもよい。
【００２０】
次に、各粒子が場から受ける力を計算し（Ｓ２０４）、その力を用いて、各粒子を所定の時間運動させる（ステップ２０５）。このときに、本発明の方法においては、場から計算される各粒子に作用する仮の力から、後述の自己相互作用力を差し引いて各粒子に作用する力を計算し、得られた力を用いて粒子の運動を計算する。
【００２１】
なお、場から受ける仮の力の大きさは場の勾配（微分値）から計算される。ここで、粒子が散乱を受ける状況を含んだシミュレーションを行う場合は、散乱を受けるタイミング、あるいは散乱方向については、乱数を用いたモンテカルロ法等の通常用いられる方法により設定することができる。また、不純物から受ける散乱については、不純物が一様に分布した状態、あるいは、偏って分布した状態等を仮定して計算することができる。
【００２２】
このサイクルを所定の回数繰り返した後（ステップＳ２０６）、各粒子の運動状態を平均して平均速度、平均密度等の計算値を得る（ステップＳ２０７）ことができる。
【００２３】
第２に、本発明の粒子運動シミュレーション方法は、予め、計算対象の空間であって、粒子が存在せず場の勾配が０と仮定した空間の中に、ある位置に１個の粒子を単独で置いたと仮定し、該１個の粒子を置いたことにより生じ、該１個の粒子に作用する場をポアソン方程式を解いて計算し、得られた場から該１個の粒子に作用する力を計算しておき、これを該１個の粒子の位置における自己相互作用力とすることを特徴とする。
【００２４】
本発明において、自己相互作用力を計算する方法を、フローチャートの例として図３を参照しながら説明する。
まず、計算対象の全空間に１個の粒子（粒子としては、複数の粒子をまとめて１個の粒子として取り扱う超粒子も含む。）のみが存在する条件を仮定し、この粒子を所定の位置に単独で置き、メッシュに割り付けて粒子密度を計算する（ステップＳ３０１）。粒子をメッシュに割り付ける方法は、複数の粒子が存在する空間における場合と同様に、粒子が属するメッシュの中心に割り付ける方法、粒子を囲むメッシュ中心に定められた比率で割り付ける方法等のいずれの割り付け方法を用いてもよい。次に、計算対象の空間内に粒子が存在せず、かつ場の勾配が全く無いと仮定した条件下、例えば、計算対象の空間の外から加わる電界や磁界や重力の影響を無視する条件下、前記１個の粒子が存在することによって生じる場をポアソン方程式を解くことによって計算する（ステップＳ３０２）。この場から前記１個の粒子に作用する力、即ち自己相互作用力を計算する（ステップＳ３０３）。力を算出する方法は、粒子が属するメッシュ中心における場の勾配から計算する方法、粒子を囲むメッシュ中心間の場の勾配を定められた比率で平均した値から計算する方法等のいずれの方法を用いてもよい。
【００２５】
Ｓ３０１からＳ３０３までの作業を計算対象の空間内の各位置について実行した後（ステップＳ３０４）、図２のステップＳ２０１に移る。以上の方法で求めた自己相互作用力を、図２のステップＳ２０４において各粒子に作用する力を計算する際、場から受ける仮の力の計算値から差し引くことにより、計算手法に起因する数値誤差の小さい粒子運動シミュレーションを実行することができる。なお、メッシュへの粒子の割り付け方法、及び、場から粒子が受ける力の計算方法には、自己相互作用力を計算する際に用いる方法と、複数の粒子が存在する空間において作用する仮の力を計算する際に用いる方法は同じ方法を用いることが好ましい。
【００２６】
本発明の粒子運動シミュレーション方法は、電荷キャリアを粒子とし、場を電場とし、半導体デバイスを対象としたシミュレーション、特に、粒子散乱のタイミングや散乱方向等をモンテカルロ法を用いて計算するモンテカルロ法半導体デバイスシミュレーションに好適に適用することができる。このとき、粒子密度はキャリア密度となる。格子振動や不純物等から受ける散乱を数値化して計算に加えて反映させることができる。
また、モンテカルロ法半導体デバイスシミュレーションの場合は、自己相互作用力を計算するために仮定する場の勾配が０の条件には、電極間の電位差や、不純物濃度分布、バンド構造に起因する場の勾配も０とすることを含む。
【００２７】
以上説明した粒子運動のシミュレーション方法は、それを実行するための計算機プログラムを用いて、スーパーコンピュータやエンジニアリングワークステーション（ＥＷＳ）などの計算機システムによって実行される。粒子運動シミュレーションのためのプログラムは、磁気テープやＣＤ−ＲＯＭなどの記録媒体によって、計算機システムに読み込まれる。図４は粒子運動のシミュレーションを実行する計算機システムの構成例を示すブロック図である。
【００２８】
この計算機システムは、中央処理装置（ＣＰＵ）４１と、プログラムやデータを格納するためのハードディスク装置４２と、主メモリー４３と、キーボードやマウスなどの入力装置４４と、ＣＲＴなどの表示装置４５と、磁気テープやＣＤ−ＲＯＭなどの記録媒体４７を読み取る読み取り装置４６とから構成されている。ハードディスク装置４２、主メモリー４３、入力装置４４、表示装置４５及び読み取り装置４６は、いずれも中央処理装置４１に接続している。この計算機システムでは、粒子運動シミュレーションを行うためのプログラムを格納した記録媒体４７を読み取り装置４６に装着し、記録媒体４７からプログラムを読み出してハードディスク装置４２に格納し、ハードディスク装置４２に格納されたプログラムを中央処理装置４１が主メモリー４３上に展開して実行することにより、粒子運動シミュレーションが実行される。
【実施例】
【００２９】
以下、本発明を実施例によりさらに詳しく説明するが、本発明はこれらに限定されるものではない。
【００３０】
実施例１
図５に断面図を示した両端に電極を付けた底辺１×１０-5ｍ、高さ１×１０-7ｍの直方体空箱の空間の中の荷電粒子のシミュレーションを行った。奥行き方向は一様とし、２次元のポアソン方程式を解いた。電極の中心を結ぶ線（図５に示した点線）上に電気素量の２０倍の電荷を有する粒子が１個だけ存在する状態を前提条件とし、電極間の電位差を１Ｖ、側面に垂直な電界を０とする境界条件下、ポアソン方程式を解いて、場から粒子に作用する力を計算した。その際、高さ方向について電極端での初期間隔１×１０-9ｍ、電極端−直方体中心間のメッシュ数１０の等比級数メッシュを設定し、横方向についてメッシュ数２１の等間隔メッシュを設定した。粒子密度を計算する際、粒子のメッシュへの割り付け方法にはＮＥＣ法（非特許文献２参照）を用い、粒子に作用する力はＮＥＣ法から計算される力を仮の力とし、その仮の力から自己相互作用力を差し引いて求めた。この空間については、理論上粒子に作用する力を粒子の電荷と電界の大きさの積として求めることが可能であり、上記の方法で計算した力と、理論上粒子に作用する力［２０×ｑ×１／（１×１０-7）［ｋｇ・ｍ／ｓ2］（ここでｑは電気素量であり１．６０２×１０-19［Ｃ］）］との差の絶対値を後者で割って得られた誤差率の最大値は１．４×１０-14（１．４×１０-12％）であった。
【００３１】
比較例２
実施例１に記載した計算において自己相互作用力を差し引かない場合、実施例１と同様にして計算した誤差率の最大値は３．８×１０-2（３．８％）であった。電極端からの距離に対し、誤差率をプロットすると図６の点線のようになった。自己相互作用力を差し引いた場合（実線、実施例１）よりも誤差率は大きかった。
【図面の簡単な説明】
【００３２】
【図１】ＮＥＣ法における粒子を囲む近接メッシュとその中心を説明するための図である。
【図２】本発明の粒子運動シミュレーション方法のフローチャートの例である。
【図３】本発明の粒子運動シミュレーション方法の内、自己相互作用力を計算するためのフローチャートの例である。
【図４】本発明を実行する計算機システムの構成例を示すブロック図である。
【図５】実施例１において、本発明の適用対象として取り上げた両端に電極を付けた直方体空箱の断面図である。
【図６】実施例１と比較例１の粒子シミュレーション結果を示すグラフである。
【符号の説明】
【００３３】
４１  中央処理装置
４２  ハードディスク装置
４３  主メモリー
４４  入力装置
４５  表示装置
４６  読み取り装置
４７  記録媒体
【図１】
2005301651.tif 000003
【図２】
2005301651.tif 000004
【図３】
2005301651.tif 000005
【図４】
2005301651.tif 000006
【図５】
2005301651.tif 000007
【図６】
2005301651.tif 000008