ステップ 2：NEGF計算

A. 左電極$\vert$デバイス$\vert$右電極の設定

図 39 に示した領域$L_0$、$C_0$、$R_0$は次のようにして設定できます。

中心領域$C_0$の幾何学的構造は、キーワード「Atoms.Number」および 「Atoms.SpeciesAndCoordinates」により、以下のように指定します。

     Atoms.Number        18
     <Atoms.SpeciesAndCoordinates
       1  C  3.000  0.000  0.000  2.0 2.0
        .....
      18  C 28.500  0.000  0.000  2.0 2.0
     Atoms.SpeciesAndCoordinates>

左側の電極部分$L_0$の幾何学的構造は、キーワード「LeftLeadAtoms.Number」および
「LeftLeadAtoms.SpeciesAndCoordinates」により、以下のように指定します。

     LeftLeadAtoms.Number         3
     <LeftLeadAtoms.SpeciesAndCoordinates         
       1  C -1.500  0.000  0.000  2.0 2.0
       2  C  0.000  0.000  0.000  2.0 2.0
       3  C  1.500  0.000  0.000  2.0 2.0
     LeftLeadAtoms.SpeciesAndCoordinates>

右側の電極部分$R_0$の幾何学的構造は、キーワード「RightLeadAtoms.Number」および
「RightLeadAtoms.SpeciesAndCoordinates」により、以下のように指定します。

     RightLeadAtoms.Number        3
     <RightLeadAtoms.SpeciesAndCoordinates         
       1  C 30.000  0.000  0.000  2.0 2.0
       2  C 31.500  0.000  0.000  2.0 2.0
       3  C 33.000  0.000  0.000  2.0 2.0
     RightLeadAtoms.SpeciesAndCoordinates>

ここで示した例は、上述の節で取り上げた炭素鎖の場合のものです。 中心領域$C_0$は18個の炭素原子から構成され、左側の領域$L_0$と右側の領域$R_0$は、それぞれ3個の炭素原子から成り、すべての隣接原子間距離は1.5 Å です。 デバイス領域$C_0$と電極$L_0$および$R_0$の幾何学的構造を設定することで、OpenMXは図 39 に示した拡張中心領域 $C (\equiv L_0\vert C_0\vert R_0)$を構成します。 中心領域$C_0$とその外側の$L_0(R_0)$の境界付近での電子構造の緩和を考慮するために、拡張中心領域$C$のグリーン関数が自己無撞着に決定されます。 NEGF法を用いて、拡張中心領域$C$のグリーン関数を計算するために、次の二つの条件を満たすように計算モデルを構築する必要があります [73]。

1. 領域$C_0$における局在基底軌道$\phi$は、領域$L_0$および$R_0$の局在基底軌道と重なるが、領域$L_1$および$R_1$と重なることはない。

2. $L_i$ ($R_i$)領域における局在基底軌道$\phi$は、その最近接セル$L_{i-1}$ ($R_{i-1}$)および$L_{i+1}$ ($R_{i+1}$)より遠方のセル内の 基底軌道と重なることはない。

OpenMXの基底関数は実空間において厳密に局在しているため [41,42]、 特定のカットオフ半径を持つ局在軌道を各領域に対して割り当てると、$L_i$および$R_i$のユニットセルの大きさを調整することにより、 上記の二つの条件を常に満たすことができます。 領域$L_0$、$C_0$、$R_0$の単位胞を指定する必要はありませんが、周期性を暗黙のうちに仮定していることに留意して下さい。 半無限電極の構築は、ステップ 1のバンド構造計算で用いた単位胞を利用して自動的に行われます。 そのための情報はファイル「NEGF.filename.hks」に保存されています。 また、図 39 に示した配置構造のため、左側と右側の電極のbc面内の単位胞ベクトルは一致していなければなりません。 また拡張中心領域$C$に対するbc面内の単位胞ベクトルは、電極部分の単位胞ベクトルと同じであることが 暗黙のうちに仮定されています。ユーザーは、こうした制約の枠内で、幾何構造を設定できます。

原子位置の指定で用いられる単位は、以下の様に指定します。

    Atoms.SpeciesAndCoordinates.Unit   Ang # Ang|AU

ただし、NEGF計算では「Ang」または「AU」のみに対応しています。 「FRAC」での指定はできませんので、注意して下さい。

どの様にOpenMXが指定した幾何学的な配置構造を解析したのか確認するために、標準出力にその情報が 記載されています。上述の炭素鎖の例では以下の情報が出力されています。

   <TRAN_Calc_GridBound>

   *******************************************************
   The extended cell consists of Left0-Center-Right0.
   The cells of left and right reads are connected as.
   ...|Left2|Left1|Left0-Center-Right0|Right1|Right2...

   Each atom in the extended cell is assigned as follows:
   where '12' and '2' mean that they are in 'Left0', and
   '12' has overlap with atoms in the Left1,
   and '13' and '3' mean that they are in 'Right0', and
   '13' has overlap with atoms in the 'Right1', and also
   '1' means atom in the 'Center'.
   ********************************************************

   Atom1  = 12 Atom2  = 12 Atom3  = 12 Atom4  =  1 Atom5  =  1 Atom6  =  1 Atom7  =  1 
   Atom8  =  1 Atom9  =  1 Atom10 =  1 Atom11 =  1 Atom12 =  1 Atom13 =  1 Atom14 =  1 
   Atom15 =  1 Atom16 =  1 Atom17 =  1 Atom18 =  1 Atom19 =  1 Atom20 =  1 Atom21 =  1 
   Atom22 = 13 Atom23 = 13 Atom24 = 13

$L_0\vert C_0\vert R_0$からなる拡張された中心領域の原子には数値が割り当てられています。「12」および「2」は原子が$L_0$に属し、 「12」は$L_1$内の原子と基底関数間の重なりを持つことを意味します。 「13」および「3」は原子が$R_0$に属し、「13」は$R_1$にある原子と基底関数間の重なりを持つことを意味します。 また「1」は原子が$C_0$に属していることを意味します。 この出力を調べることで、計算する系の幾何学的な構造が適切に構築されているか否かを確認することができます。

B. キーワード

ステップ 2のNEGF計算を実行するために、キーワード「scf.EigenvalueSolver」を 次のように指定します。

    scf.EigenvalueSolver       NEGF

NEGF計算に関連するキーワードを以下に列挙します。

    NEGF.filename.hks.l     lead-chain.hks
    NEGF.filename.hks.r     lead-chain.hks

    NEGF.Num.Poles             100       # defalut=150
    NEGF.scf.Kgrid             1 1       # defalut=1 1

    NEGF.bias.voltage          0.0       # default=0.0 (eV)
    NEGF.bias.neq.im.energy    0.01      # default=0.01 (eV)
    NEGF.bias.neq.energy.step  0.02      # default=0.02 (eV)

各キーワードの内容を以下に説明します。

    NEGF.filename.hks.l     lead-chain.hks
    NEGF.filename.hks.r     lead-chain.hks

電極部分の情報を含むファイルは、上記の二つのキーワードにより指定できます。 「NEGF.filename.hks.l」は左側、「NEGF.filename.hks.r」は 右側の電極部分に対応します。

    NEGF.Num.Poles             100       # defalut=150

平衡密度行列は、周回積分法により評価します[73,74]。 この方法で用いる極の数は、キーワード
「NEGF.Num.Poles」により指定します。

    NEGF.scf.Kgrid             1 1       # defalut=1 1

逆格子ベクトル ${\bf\tilde{b}}$および ${\bf\tilde{c}}$を離散化するためのk点数を、 キーワード「NEGF.scf.Kgrid」により指定します。 a軸方向には周期性を持っていないため、a軸に対しては指定する必要がありません。

    NEGF.scf.Iter.Band          6        # defalut=6

SCF計算の最初の数ステップでは、b軸およびc軸と同様にa軸方向にも周期性を仮定することにより、 通常の対角化法を用いた方が最終的な収束が加速されることが分かっています。 自己無撞着なNEGF法においてしばしば問題となるSCF計算での収束困難性は、この方法によってかなり軽減します。 通常の対角化法を使用するSCFステップ数を、キーワード「NEGF.scf.Iter.Band」により指定します。 「NEGF.scf.Iter.Band」で指定したSCFステップ数までは通常の対角化法を用い、それ以降のステップにおいては、 NEGF法が適用されることになります。デフォルト値は6です。

    NEGF.bias.voltage          0.0       # default=0.0 (eV)

キーワード「NEGF.bias.voltage」により、 左右の電極間に印加するソース・ドレインバイアス電圧をeVの単位で指定します。 このeVの単位は電圧に換算するとボルト(V)に対応しています。 ソースとドレイン間の電位差のみが物理的に意味を持つため、電位の差である ソース・ドレインバイアス電圧を本キーワードによって与えることになります。

    NEGF.bias.neq.im.energy    0.01      # default=0.01 (eV)
    NEGF.bias.neq.energy.step  0.02      # default=0.02 (eV)

有限のソース・ドレインバイアス電圧を印加すると、密度行列の一部は非平衡グリーン関数から計算されます。 非平衡グリーン関数は複素平面上で一般に解析的ではありませんので、平衡グリーン関数で用いられた周回積分法は使用できません。 現在の実装では、非平衡グリーン関数は微小虚部を付与した実軸上で単純な矩形求積法を用いて評価しています。 この際、虚数部はキーワード「NEGF.bias.neq.im.energy」により指定します。 またステップ幅はキーワード「NEGF.bias.neq.energy.step」により、eVの単位で与えます。 通常はデフォルト値で十分な精度が確保できますが、収束性の詳細な議論は文献[73]を参照して下さい。 非平衡グリーン関数を評価する実軸上のエネルギー点数は標準出力およびファイル「System.Name.out」で確かめることができます。 「NEGF-Chain.dat」の例では、バイアス電圧が0.5Vの際には120個のエネルギー点が割り当てられており、以下のように 標準出力から確認できます。

   Intrinsic chemical potential (eV) of the leads
     Left lead:  -7.752843837400
     Right lead: -7.752843837400
     add voltage =  0.0000 (eV) to the  left lead: new ChemP (eV):  -7.7528
     add voltage =  0.5000 (eV) to the right lead: new ChemP (eV):  -7.2528

   Parameters for the integration of the non-equilibrium part
     lower bound:           -8.706843837400 (eV)
     upper bound:           -6.298843837400 (eV)
     energy step:            0.020000000000 (eV)
     number of steps:           120

グリーン関数を評価するエネルギー点の総数は、非平衡グリーン関数を評価する実軸上のエネルギー点と、平衡グリーン関数を 評価する極の数の総和で与えられます。計算時間は、エネルギー点の総数に比例することに留意して下さい。

    NEGF.Poisson.Solver       FD     # FD|FFT, default=FD

NEGF法では差電子密度に対する静電ポテンシャルを2次元FFT+1次元有限差分法（FD）[73]もしくは 3次元FFT法（FFT）[75]のどちらかの方法で評価します。 このポアソン・ソルバーの選択はキーワード「NEGF.Poisson.Solver」で行います。 非極性系についてはどちらの方法も同様の静電ポテンシャルを与えますが、極性系については両者の差異は大きくなります。 前者のFDでは、電極部分と中心領域の境界条件が厳密に満されるため、理論的にはより正しい方法です。 ここでの境界条件とは、電極部分と中心領域の界面での静電ポテンシャルがステップ 1の計算におけるポテンシャル値と 同一であるというものです。 SCF収束性の観点からは、FDの収束性が悪く、一方、後者のFFTの方がよりSCF収束が容易であるという傾向が分かっています。 デフォルトはFDです。

C. SCFの収束条件

NEGF法では、キーワード「scf.criterion」により与えられるSCFの収束条件は 残差ノルム「NormRD」に適用されます。NEGF法以外の場合には、SCFの収束条件はdUeleに適用されます。

D. ゲート・バイアス電圧

OpenMX Ver. 3.9の実装では、次式で定義されるポテンシャルを加えることで、ゲート電圧$V_{\rm g}(x)$が取り扱われます。

$$V_{\rm g}(x) = V_{\rm g}^{(0)} \exp\left[-\left(\frac{x-x_{\rm c}}{d}\right)^{8}\right],$$

ここで、 $V_{\rm g}^{(0)}$はゲート電圧に対応する定数であり、 キーワード「NEGF.gate.voltage」により指定できます。

      NEGF.gate.voltage   1.0    # default=0.0 (in eV)

また$x_{c}$は領域$C_0$の中心位置、$d$は領域$C_0$のa軸方向のセルベクトルの長さです。 式の形が示すように、印加されたゲート電圧は、主として拡張中心領域$C$内の領域$C_0$に作用します。 このポテンシャルは、鏡像電荷によって生成されるポテンシャル形状に類似しています [76]。

E. 状態密度（DOS）

NEGF計算において、状態密度 (DOS)は次のキーワードを設定することにより計算できます。

    Dos.fileout                 on              # on|off, default=off
    NEGF.Dos.energyrange     -15.0 25.0 5.0e-3  #default=-10.0 10.0 5.0e-3 (eV)
    NEGF.Dos.energy.div        200              # default=200
    NEGF.Dos.Kgrid             1 1              # default=1 1

NEGF法においてDOSを計算する際には、「Dos.fileout」を「on」にして下さい。 また、DOSを計算するエネルギー範囲は、キーワード「NEGF.Dos.energyrange」により与えます。 ここで、最初と二番目の数値はエネルギーの下限と上限で、三番目の数値はDOSを滑らかにするために用いる虚数値です。 「NEGF.Dos.energyrange」により指定されたエネルギー範囲は、 キーワード「NEGF.Dos.energy.div」で指定した数値で分割され、そのエネルギー点上で DOSが計算されます。 逆格子ベクトル ${\bf\tilde{b}}$および ${\bf\tilde{c}}$を離散化するためのk点のそれぞれの数を、 キーワード「NEGF.Dos.Kgrid」により指定します。 「NEGF.Dos.Kgrid」で与える数を 「NEGF.scf.Kgrid」より大きな値にすることで計算精度が向上します。 NEGF計算が完了すると、「System.Name.Dos.val」と「System.Name.Dos.vec」の二つのファイルを生成します。 その後の解析手順は、一般の場合と同様です。 また、エネルギーの原点は左側の電極の化学ポテンシャルに設定されていることに留意して下さい。