学習目標

このページを読み終えると、以下のことができるようになります：

OghmaNanoが半導体デバイスをモデル化する際に解く方程式として、ドリフト拡散方程式、ポアソン方程式、およびSRHモデルに基づくキャリアの捕捉・放出方程式が挙げられる。
ニュートン法がこれらの方程式をどのように反復的な根探索問題へと変換するかを説明してください。
収束状況を示すプロットを読み取り、初期の平衡状態の解および各電圧ステップにおける解を識別してください。
ソルバーの設定、すなわち反復計算回数の上限値、目標誤差、およびクランピング値を適切に調整することで、物理的精度を損なうことなく実行時間を短縮できます。

デバイスシミュレーションの高速化

多数のシミュレーションを実行する場合、各シミュレーションが完了するまでの速度が重要となります。本ページでは、OghmaNanoの電気的ソルバーの実行速度を向上させる方法について説明します。適切な設定を有効活用するためには、まずソルバーが実際にどのような処理を行っているかを理解する必要があるため、ページの前半ではモデルの基礎、すなわち解かれる方程式およびその解法について述べ、後半では実行時間を短縮するためのソルバーの調整方法を解説します。これらの設定によって物理法則自体は変わらず、数値解が得られるまでの速度のみが変化します。

This page has a little more maths than most. The key idea to carry through is that solving the device equations is really a root-finding problem: OghmaNano repeatedly improves a guess until an error function f(x) falls to zero. Everything in the tuning section is about reaching that point in fewer iterations.

1. モデルの基礎事項

シミュレーションの高速化方法を説明する前に、OghmaNanoが実際に何を計算しているのかを理解しておくことが重要です。デバイス内の電荷輸送を記述するため、OghmaNanoは基本的なドリフト拡散方程式およびポアソン方程式（ガウスの法則）を解きます。これらの方程式はデバイス全体にわたる有限差分メッシュ上で解かれるため、一方の電極から他方の電極へ向かう電荷の移動が伝導帯（LUMO）および価電子帯（HOMO）の両方において追跡可能となります。

1.1 デバイス内における電荷輸送

ポアソン方程式は静電ポテンシャルと自由電荷の関係を表し、2つの電流方程式は電子および正孔の駆動項（ドリフト拡散による寄与）を記述する。

\[ \nabla \cdot \varepsilon_0 \varepsilon_r \nabla \phi = q\,(n - p) \] \[ J_n = q\,\mu_e\,n\,\nabla E_c + q\,D_n\,\nabla n \] \[ J_p = q\,\mu_h\,p\,\nabla E_v - q\,D_p\,\nabla p \]

これらの電流密度は、電子および正孔の連続の式に代入される。この式では、電流の流れと再結合、キャリアのトラップ現象、ならびに自由キャリア密度の時間変化率が釣り合っている。

\[ \nabla \cdot \mathbf{J}_n = q\left( R_n + T_n + \frac{\partial n_{\text{free}}}{\partial t} \right) \] \[ \nabla \cdot \mathbf{J}_p = -q\left( R_p + T_p + \frac{\partial p_{\text{free}}}{\partial t} \right) \]

ここで\( \phi \)は静電ポテンシャル、\( n \)および\( p \)は自由電子密度と正孔密度、\( \mu_e \)および\( \mu_h \)は移動度、\( D_n \)および\( D_p \)は拡散係数、\( E_c \)および\( E_v \)はバンドエッジ、\( R \)は再結合率、\( T \)はキャリアがトラップに捕獲される正味の速度を表す。詳細な導出過程および符号規則についてはMacKenzieらによるJ. Phys. Chem. C（doi:10.1021/jp200234m）を参照されたい。

1.2 キャリアトラップと再結合

Density-of-states diagram plotting electron energy against log(DOS). A green parabolic band of free electrons sits at the top at the electron mobility edge mu_e^0, a magenta parabolic band of free holes at the bottom at the hole mobility edge mu_h^0. Red exponential trap states descend from the electron trap edge E_e^d and blue exponential trap states rise from the hole trap edge E_h^d. Curved arrows labelled r1 and r2 connect the free electrons to the traps (capture and emission), and arrows r3 and r4 connect the traps to the free holes. — OghmaNanoにおいて解かれるSRHモデルに基づくキャリアの捕捉・放出過程 *すべて* メッシュ点上における自由キャリアなし_無料, p_無料これらは移動度端付近において放物線型のバンド構造を形成している \( \mu_e^0 \) / \( \mu_h^0 \)赤色および青色の三角形は、電子および正孔のトラップ状態が指数関数的に分布していることを示している。これらの速度値は \( r_1 \)–\( r_4 \) トラップへのキャリアの捕獲およびトラップからの放出である。

電荷トラップおよび再結合は有機デバイスやその他多くの材料系において重要な現象であるため、これらを明示的にモデル化する必要がある。通常は、解きやすい形としてSRH再結合を単一の定常状態式に簡略化するが、OghmaNanoではそうはしない。時間領域での方程式求解を目的としているため、有限差分格子上の各メッシュ点においてSRHモデルに基づく捕獲・放出率の全方程式を解いているのである。

各メッシュ点上において、??に示す状況が解かれる。上部および下部にある放物線は、電子および正孔の移動度端における自由キャリアバンドである。電子トラップ端から下方へは電子トラップの指数分布（赤色）が広がり、正孔トラップ端から上方へは正孔トラップの指数分布（青色）が広がっている。\( E_e^d \)および\( E_h^d \)で示されるこれらの分布内の各離散準位は個々のトラップであり、それぞれに対して完全な捕捉・放出率方程式が解かれる。ある準位におけるトラップされた電子数の変化は、以下の4つの率のバランスによって決まる：

\[ \frac{\partial n_t}{\partial t} = r_1 - r_2 - r_3 + r_4 \]

ここで\( r_1 \)は電子の捕獲過程、\( r_2 \)は電子の放出過程、\( r_3 \)は正孔の捕獲過程、\( r_4 \)は正孔の放出過程を表し、それぞれトラップ密度、キャリア密度、捕獲断面積、熱速度によって決定される標準的な捕獲・放出式に従う（上記の参考文献を参照）。すべてのエネルギー準位において詳細釣り合いが保たれている。トラップ分布全体にわたって純捕獲量を合計することで、自由キャリアの再結合率が得られる。

\[ R_n = \sum_{0}^{n_{\text{band}}} \left( r_1^{\,e} - r_2^{\,e} \right) \]

速度向上にとってなぜこれが重要なのか。位置空間上の各点においてエネルギー空間内のトラップ分布を解くことで、エネルギーおよび位置の両方の関数としての電荷密度の完全な分布が得られる。これによりOghmaNanoは通常の定常状態近似を用いずに過渡解析を実行できるのだが、同時に同時に解くべき連立方程式の数も増加するため、用いる数値計算手法やその設定が実行時間に大きく影響するのである。

2. 方程式の解法：ニュートン法

上記のすべての方程式は、アイザック・ニュートンが1671年に『流数法』において発表したニュートン法を用いて解かれます。この手法は、この種の連立方程式系を解くための最も効率的かつ堅牢な手法の一つです。手計算では解けないような方程式であっても、ニュートン法ならほぼ確実に解を求めることが可能です。

2.1 根と誤差関数

まずは簡単な問題から始めましょう。どのようなxの値のときに、次の方程式はゼロになりますか？

\[ 0 = (x - 3)(x - 40) \]

見てすぐに答えは \( x = 3 \) および \( x = 40 \) であることがわかる。これらがその方程式の根となる。この問題は関数の形でも表現可能である。

\[ f(x) = (x - 3)(x - 40), \]

\( f(x) \)を\( x \)に対してプロットすると、根とは単に曲線が軸と交わる点、すなわち\( f(x) = 0 \)となる点である。この方程式をゼロに向かって減少していく誤差関数と見なすことが重要な考え方だ。デバイスモデルにおいては、\( x \)の代わりに電位\( \phi \)や電子・正孔の準フェルミ準位\( E_{Fn} \)、\( E_{Fp} \)といった物理的な未知変数が用いられ、これらの値が調整されることで最終的に\( f(x) = 0 \)となるのである。

根を見分けるには、単純な方程式であれば目視でも可能だ。しかしより複雑な場合は別である。

\[ f(x) = \frac{(x^2 + 10x - 3)\,\cos^2(4x)}{\tan(x)}. \]

この関数の根を求めるためにグラフを描くのは非常に手間がかかり、さらに半導体デバイス方程式ははるかに複雑である。まさにこのような状況でニュートン法が真価を発揮するのである。

2.2 手順の詳細

具体的な例として\( f(x) = x^2 + 10x - 3 \)を挙げます。ニュートン法の手順は以下の通りです：

解くべき変数について方程式を微分する：\( f'(x) = 2x + 10 \)
根の初期値を適当に設定しましょう。どの値でも構いません。ここでは\( x = 10 \)を採用します。この値は適切でない値です（\( f(10) = 197 \)、実際にはゼロからほど遠い値です）が、問題ありません。
その推測値において関数とその微分値を評価する：\( f(10) = 197 \)および\( f'(10) = 30 \)。
ニュートンの更新式を用いて、より良い近似値を計算する：

\[ x_{+1} = x - \frac{f(x)}{f'(x)} \]

値を代入することでより良い推定値が得られ、\( x_{+1} = 10 - 197/30 = 3.433 \)となります。その後はこの手順を単純に繰り返すだけです。つまり、改善された推定値を再び代入して再度評価し（\( f(3.433) = 43.115 \)、\( f'(3.433) = 16.866 \)）、さらに良い\( 0.877 \)の推定値を得る、という流れです。各反復計算を重ねるごとに誤差関数はゼロに近づき、推定値も真の根に近づいていきます。

2.3 実践例

この手法全体はわずか数行のコードで実装可能だ。以下のループが\( f(x) = x^2 + 10x - 3 \)に対してこれを実行し、各反復計算における現在の近似値と誤差を出力する。

x = 10;                       % our guess for the answer
for i = 1:10                  % loop ten times
    y  = x*x + 10*x - 3;      % f(x)
    dy = 2*x + 10;            % f'(x)
    new_x = x - (y/dy);       % calculate a better x
    x = new_x;                % update x
    sprintf("%f %f", x, y)    % print the estimate and the error
end                           % repeat

A code editor on the left showing the Newton iteration loop for f(x) = x*x + 10x - 3, and a terminal on the right showing the printed output. The two printed columns are the running estimate of the root and the error f(x): 3.433333 and 197.0, then 0.876746 and 43.12, then 0.320644 and 6.54, then 0.291582 and 0.309, then 0.291503 and 0.000845, then 0.291503 and 0.000000 for the remaining iterations. — 反復計算の実行状況。左列が根の現在の推定値、右列が誤差である \( f(x) \)推定値は安定して一定の値に収束する \( 0.291503 \) 数回の反復計算のみで誤差がゼロに収束する。

出力結果は非常に参考になります。推定値が\( 10 \rightarrow 3.433 \rightarrow 0.877 \rightarrow 0.321 \rightarrow 0.292 \rightarrow 0.291503 \)へと移動し、誤差も\( 197 \rightarrow 43 \rightarrow 6.5 \rightarrow 0.31 \rightarrow 0.00085 \rightarrow 0 \)まで減少します。各反復計算ごとに正しい小数点以下の桁数がほぼ倍増している点に注目してください。これがいわゆる超線形収束であり、ニュートン法が方程式の根を見つけるための極めて効率的な手法となっている理由です。

2.4 同様の考え方を行列形式で表現した場合

このデバイスは単一の方程式ではなく、多数の連立方程式から成る。すなわち、各メッシュ点上における電位、電子および正孔の準フェルミ準位、そしてトラップ占有状態がそれに該当する。行列形式で表現すると、ニュートン法は次のようになる：

\[ [J]\,[\Delta x] = -[f(x)] \]

\( [J] \)はヤコビアン行列であり、すべてのデバイス方程式（ポアソン方程式、ドリフト拡散方程式、SRHモデル）の各未知数に対する偏微分値から成る行列です。\( [f(x)] \)は現在の近似値における方程式の残差誤差を表し、\( [\Delta x] \)はそれに基づいて算出される更新量、すなわちデバイス全体における電位分布、準フェルミ準位、およびトラップ占有状態の改善分となります。この更新量が現在の解に加算されるのです。

\[ [x] = [x] + [\Delta x], \]

この処理を反復して行い、\( f(x) \approx 0 \)に達した時点で半導体デバイス方程式が解かれる。ここで重要な点が2つある。まず、\( f(x) \)はゼロに近づけるべき量であり、収束度を測るための自然な指標となる。また、初期値の精度が高ければ高いほど、ソルバーはより正解に近い状態から計算を開始できるため、必要な反復回数も少なくて済む。

3. 実際のシミュレーションの高速化方法

基礎が確立された上で、実際の設定内容も理解しやすくなる。ここで用いる例は二次元の電界効果トランジスタ（すべてのトラップを含む二次元デバイス）であるが、同様の考え方はあらゆるシミュレーションに適用可能だ。

3.1 ディスクへの書き込みを無効化する

最も簡単に行える改善策は、ディスクへのデータ書き込み量を減らすことです。ファイルの書き込み処理は時間がかかり、大量の出力データを生成するシミュレーションでは実行時間の大部分を占める可能性があります。OghmaNanoにはワンクリックでディスクへの書き込み量を最小限に抑える機能が備わっており、出力設定画面では個々の出力ファイルオプション（例えば時間領域データのダンプ機能など）も個別に有効化・無効化可能なため、必要なファイルのみを保存できます。最速のシミュレーション実行を実現するためには、ディスクへの書き込みをできるだけ少なくすることが重要です。

まずはこれをお試しください。以下の数値設定に触れる前に、不要な出力機能を無効にして再実行してください。ディスク書き込み量が多いシミュレーションにおいては、この操作だけで最も大きな速度向上が得られますが、物理モデルや結果の精度には一切影響しません。

3.2 ソルバーの収束状況の確認方法

A plot of solver error on a logarithmic vertical axis against iteration number on the horizontal axis. From 0 to about 1000 iterations the error wanders down from around 1e15 to a noisy floor near 1e-7: this is the initial equilibrium solve. After that the trace shows a series of sharp spikes that rise and then decay, each spike corresponding to a voltage step where the error jumps and is then driven back down. — 反復計算回数に対するソルバーの誤差 \( f(x) \) の変化を示す。左側の緩やかな下降部分は初期平衡状態の求解過程であり、右側の各ピークは電圧ステップに対応する。電圧変更により方程式の解が乱れるため、ソルバーが誤差を再び低減させているのである。

ソルバーの調整を行うには、まずその動作状況を確認する必要があります。OghmaNanoでは、各反復計算におけるニュートン法による誤差（電流誤差）を出力したり、電圧ステップごとに単一行のみを表示するのではなく、収束履歴をファイルに書き出すことも可能です。ニュートン法誤差出力機能および収束ログ機能を有効化した上で、シミュレーションを再実行してください。

最初のシミュレーションステップは特別なものです。モデルが起動した時点ではデバイス状態に関する適切な推定値が存在せず、電位やフェルミ準位もまだ判明していないため、初期状態は解からかなり離れた位置となります。??においては、誤差が非常に大きな値（\( 10^{15} \)程度）から始まり、多数の反復計算を経て徐々に減少し、最終的には\( 10^{-7} \)から\( 10^{-9} \)付近の値に落ち着きます。この時点で方程式は事実上解かれたことになります。この最初の計算が解からかなり遠い位置にあるため、2つの対策が意図的に講じられます。まず、ニュートン法によるステップサイズを厳しく制限することで（後述）、ソルバーが不安定な状態に陥ることなく安全かつ緩やかに収束するようにします。また、十分に安定した解に到達できるよう、通常は約1000回程度の反復計算が許可されます。

平衡状態での計算が完了した後、モデルは指定された電圧値を順次適用していく。各電圧ステップにおいてデバイスは前回の解からわずかにずれるため、誤差が一時的に増大するが、すぐに再び減少していく。これがプロットの右側に見られる誤差のピーク群である。各ステップは直前の解の近傍から開始されるため、計算は迅速に収束する。

3.3 ソルバーの調整

この収束状況のプロットは、第2節で述べた誤差関数 \( f(x) \) そのものであり、チューニングの目的はこの曲線をより短くすること、つまりより少ない反復計算回数で収束させることです。ソルバー設定における主な制御パラメータは以下の通りです：

最初の反復計算における最大反復回数。例えば400回の反復計算を経て最初の計算がすでに収束状態に達している場合、1000回まで計算を続ける意味はない。この上限値を1000から600といった具合に引き下げることで、すでに収束している計算部分を省略できる。こうして削減された反復回数分だけ実際の計算時間が節約されるのである。
最初の反復計算における目標誤差値。通常、最初の計算では\( 10^{-9} \)のような極めて小さな誤差に到達するよう設定される。しかし結果が\( 10^{-7} \)時点で既に十分に収束している場合は、目標値を緩和すべきである。つまり\( 10^{-9} \)を目指す代わりに\( 10^{-7} \)で計算を打ち切ることで、最も計算負荷が高く進捗の遅い段階を省略できるのだ。
クランピング処理 クランプ値とは、完全なニュートンステップのうち実際に適用される割合のことである \( [\Delta x] \) ソルバーが移動を許可される範囲である。完全なステップ幅はヤコビアン行列を用いて算出されるが、クランプ処理により意図的にその一部のみが適用される。例えばニュートン法が変数の移動量を0.5と求めた場合でも、クランプ値により実際の移動量は0.1に制限されることもある。クランプ値が小さいほど収束は遅くなるがソルバーが異常終了するリスクは著しく低減する。逆にクランプ値が大きいほど収束は速まるが失敗の可能性も高まる。最初のステップでは通常0.1程度の保守的なクランプ値が用いられ、構造が安定している場合は0.2や0.3などに引き上げることで初期段階の収束速度が顕著に向上する。その後のステップでもデフォルト値からクランプ値を増加させることは可能だが、過度な値設定は収束不能を招く恐れがある。
ニュートン法による反復計算の最小回数。これにより、ソルバーが終了する前に少なくとも指定された回数だけ計算を実行することが強制され、そうすることで局所的な極小点に陥ってしまうのを回避できる。数回程度の小さな値を設定しておくことで安全策となり得るが、解がすでに求まった後でこの値を過大に設定すると、ノイズレベル付近で計算が進むだけで時間の無駄になる。
最小限のソルバー誤差および最大反復回数（その他のステップは標準設定）。これらは各電圧値における計算の目標精度と反復回数の上限を定めるパラメータである。これらも最初のステップ用パラメータと同様の挙動を示し、デバイス構造の収束速度が把握でき次第、より厳格な値に設定することも可能だ。

処理手順としては、各設定変更の後で収束状況を記録したファイルを再描画し、グラフがより急峻になる様子を確認する。具体的には、最初の計算における平坦な部分を削減し、初期段階での目標誤差を緩和した上で、クリッピング値を増大させることで収束速度を高めるのだ。本例のデバイスにおいては、これらの変更を施した結果、計算時間が約15秒から8秒程度へと短縮されたが、得られる結果自体は変わらなかった。モデルに同梱されている入力ファイルでは控えめな値が用いられているため、ユーザー自身の構造に対してもこれらの値を引き上げる余地が通常は存在する。

節約できた時間は精度向上に活用してください。ソルバーの最適化が完了すれば、得られた時間をより難易度の高い問題の解決に再投資できます。例えばメッシュ点の数を増やすといった具合です。適切にチューニングされたニュートン法の設定こそが、より大規模で細密な、あるいは強く結合した半導体デバイスを現実的な時間内で解析可能にするのです。

よくある落とし穴

クランプ値を急激に引き上げると、ソルバーの処理速度は向上するものの、クラッシュしたり収束しなくなる可能性が高まります。クランプ値は段階的に増加させ、シミュレーションが正常に完了することを確認してください。
最小反復回数の設定値を過度に高くすると、すでに収束しているにもかかわらずソルバーが引き続き計算を続けてしまい、節約しようとした時間が無駄になってしまいます。
目標とする誤差を過度に緩和しすぎると、実際には十分に収束していない状態のまま計算が終了してしまう可能性があります。誤差を緩和する際は、収束プロットが明らかに数値ノイズの領域へと平坦化しているレベルまでに留めるべきです。
各構造ごとに最適な設定値は異なります。あるデバイスで有効な値が別のデバイスでは過度に厳格すぎる可能性があるため、構造を変更する際には必ず収束プロットを再確認してください。

概要と次の手順

OghmaNanoは、メッシュ上でドリフト拡散方程式、ポアソン方程式、および完全なSRHモデルに基づくキャリアの捕捉・放出方程式を解くことにより、エネルギー空間および位置空間におけるキャリア分布を求める。
これらの連立方程式はニュートン法を用いて解かれる。これは反復計算手法であり、誤差関数 \( f(x) \) をゼロに近づけることで超線形収束を実現する。その行列形式は \( [J][\Delta x] = -[f(x)] \) に示されている。
収束状況のグラフでは、まず長時間にわたる初期平衡状態の計算が行われ、その後各電圧ステップにおいて迅速な計算が実施される様子が示されている。
高速化するためには、ディスクへの書き込み量を減らし、初期反復計算の回数制限を引き下げて目標誤差を緩和し、クランピング値を慎重に増加させ、さらに誤った極小値に陥るのを防ぐために最小反復回数を少なめに設定する必要がある。
パラメータ調整によって得られた時間を、より細かいメッシュの生成やより複雑な問題の解析に再利用することが可能だ。

次に、これらの方程式に用いられる材料パラメータの定義方法については材料データベースをご覧ください。