固体系

ここではTiC固体を例として、固体系における内殻準位の絶対束縛エネルギーの計算方法を説明します。

始状態の計算は以下のように実行します。

    % mpirun -np 112 ./openmx TiC216.dat | tee TiC216.std
 

ここで、特別なキーワードは指定しませんが、 「scf.SpinPolarization=on」としてスピン分極計算が実行されます。ここでのスピン分極計算は 終状態の計算と整合性を取るために実施されています。ただし計算結果は非磁性となります。 ここで用いた入力ファイル「TiC216.dat」はディレクトリ「work」に収容されています。

終状態の計算は以下のように実行します。

    % mpirun -np 112 ./openmx TiC216-CH3.dat | tee TiC216-CH3.std
 

ここで用いた入力ファイル「TiC216-CH3.dat」はディレクトリ「work」に収容されています。 この入力ファイル中で原子種は以下のように定義されています。

   <Definition.of.Atomic.Species
     Ti  Ti7.0-s3p2d2        Ti_PBE13
     C   C6.0_1s-s3p2d1      C_PBE17_1s
     C1  C6.0_1s_CH-s3p2d1   C_PBE17_1s
   Definition.of.Atomic.Species>

そして、以下のように原子種「C1」が原子5に割り当てられます。

   Atoms.Number         216
   Atoms.SpeciesAndCoordinates.Unit   Ang # Ang|AU
   <Atoms.SpeciesAndCoordinates           
      1  Ti   0.000000000000  0.000000000000  0.000000000000   6.0 6.0
      2  Ti   2.163500000000  2.163500000000  0.000000000000   6.0 6.0
      3  Ti   0.000000000000  2.163500000000  2.163500000000   6.0 6.0
      4  Ti   2.163500000000  0.000000000000  2.163500000000   6.0 6.0
      5  C1   2.163500000000  0.000000000000  0.000000000000   3.0 3.0
      6  C    0.000000000000  2.163500000000  0.000000000000   3.0 3.0
      7  C    0.000000000000  0.000000000000  2.163500000000   3.0 3.0
      8  C    2.163500000000  2.163500000000  2.163500000000   3.0 3.0
      ....
      ..
   Atoms.SpeciesAndCoordinates>

以下のキーワードを指定することで、原子5の$1s$状態に内殻ホールが生成されます。

    scf.restart                 on
    scf.restart.filename        TiC216
    scf.coulomb.cutoff          on
    scf.core.hole               on
    scf.system.charge          0.0       # default=0.0

    <core.hole.state  
      5 s 1
    core.hole.state>

終状態の計算においてHartreeポテンシャル$V_{H}$ は二つの寄与から成ります[88]。 次式で定義される周期部分 $V_{H}^{\rm (P)}$ と非周期部分 $V_{H}^{\rm (NP)}$ です。

$$V_{H}({\bf r}) = V_{H}^{\rm (P)}({\bf r}) + V_{H}^{\rm (NP)}({\bf r}),$$ (15)

周期部分 $V_{H}^{\rm (P)}$は 始状態計算で得られた電荷密度を用いて周期境界条件でポアソン方程式を介して計算されます。 キーワード「scf.restart」と「scf.restart.filename」により始状態計算で得られた電荷密度が指定されます。 非周期部分 $V_{H}^{\rm (NP)}$は差電荷密度 $\Delta \rho({\bf r}) = \rho_{\rm f}({\bf r}) - \rho_{\rm i}({\bf r})$に対して 厳密クーロンカットオフ法[91]を用いて計算されます。 ここで $\rho_{\rm f}({\bf r})$ と $\rho_{\rm i}({\bf r})$はそれぞれ終状態と始状態の電荷密度です。 クーロン相互作用のカットオフ半径は最も短い格子ベクトルの長さの半分に設定されます。 厳密クーロンカットオフ法を有効にするにはキーワード「scf.coulomb.cutoff」を「on」に設定します。 内殻ホールはキーワード「scf.core.hole」と「core.hole.state」により生成されます。 今回の場合では、原子5の1s状態に内殻ホールが生成されます。 TiC固体は金属であるためキーワード「scf.system.charge」は「0.0」に設定されていることに注意してください。 エネルギーギャップのある系を扱う際には「scf.system.charge」を「1」に設定しなければなりません。

始状態と終状態の計算を終えると、出力ファイルから全エネルギーが得られます。

     Initial state: -10499.900104007471 (Hartree)
     Final state:   -10489.553360141708 (Hartree)

そして、式 (14)を用いて、束縛エネルギーが以下のとおり計算できます。

$$E^{\rm (metal)}_{\rm b} = E^{(0)}_{\rm f}(N-1) - E^{(0)}_{\rm i}(N),$$

$$= -10489.553360141708 - (-10499.900104007471) = 10.346743865763 ({\rm Hartree}),$$

$$\approx 281.55 ({\rm eV}).$$

得られた281.55eVの値は実験値の281.5eV [92]と比較して良く一致していることが分かります。

Figure 79: (a) ギャップのある系、(b) 半金属(グラフェン)と金属での 内殻ホール間の距離の関数として、 最も収束した結果からの相対値として計算された結合エネルギー。 (a) と(b) の基準の結合エネルギーは各系の最大の単位胞に対して、 それぞれ式(13) と(14) から計算。図中でEq. (B) とEq. (M) は それぞれ式(13) と式(14) に対応。 (c)シリコンでの2p状態での内殻ホールの生成により誘起された差電荷密度。 単位胞は1000原子を含み、内殻ホール間距離は27.15 Å。

ギャップのある系の他の例として、固体シリコンの計算を説明します。 始状態と終状態の計算を以下のように実行します。

    % mpirun -np 256 ./openmx Si-4-SOI.dat | tee Si-4-SOI.std
    % mpirun -np 256 ./openmx Si-4-CH-SOI1.dat | tee Si-4-CH-SOI1.std
    % mpirun -np 256 ./openmx Si-4-CH-SOI6.dat | tee Si-4-CH-SOI6.std
 

「Si-4-SOI.dat」は始状態計算に対する入力ファイル、また 「Si-4-CH-SOI1.dat」 と 「Si-4-CH-SOI6.dat」は終状態計算に対する入力ファイルで、 それぞれ $2p_{3/2} (J=3/2, M=3/2)$ と $2p_{1/2} (J=1/2, M=-1/2)$ で指定される 内殻準位にホールを導入しています。 Si原子の$2p$状態のスピン軌道相互作用を考慮するため、 以下のキーワードを指定し、ノンコリニア計算を実行します。

    scf.SpinPolarization        nc         # On|Off|NC
    scf.SpinOrbit.Coupling      on         # On|Off, default=off

始状態と終状態の計算を終えると、出力ファイルから全エネルギーが得られます。

     Initial state:          -34820.483255130872 (Hartree)
     Final state for SOI1:   -34816.628201335407 (Hartree)
     Final state for SOI6:   -34816.601864921540 (Hartree)

化学ポテンシャルは始状態計算から得られます。 式 (13)を用いて、「Si-4-CH-SOI1.dat」と「Si-4-CH-SOI6.dat」の束縛エネルギーが 以下のとおり計算できます。

$${\rm For}~2p_{3/2}$$

$$E^{\rm (bulk)}_{\rm b} = E^{(0)}_{\rm f}(N-1) - E^{(0)}_{\rm i}(N-1) + \mu_0,$$

$$= -34816.628201335- (-34820.483255131) -0.201641583 = 3.65341221 ({\rm Hartree}),$$

$$\approx 99.41 ({\rm eV}),$$

$${\rm For}~2p_{1/2}$$

$$E^{\rm (bulk)}_{\rm b} = E^{(0)}_{\rm f}(N-1) - E^{(0)}_{\rm i}(N-1) + \mu_0,$$

$$= -34816.601864922 - (-34820.483255131) -0.201641583 = 3.6797486 ({\rm Hartree}),$$

$$\approx 100.13 ({\rm eV}).$$

$2p_{3/2}$ と $2p_{1/2}$ に対する99.41 と 100.13 eVの束縛エネルギー値は 実験値の99.2と99.8eV [92]と比較して良い一致が見られます。 ここで、$2p_{3/2}$ と $2p_{1/2}$ の縮重度は表12に記載されるようにそれぞれ4と2です。 ギャップのある系では単位胞サイズに対する絶対束縛エネルギーの収束性は遅いです。 図78(a)と(b)ではそれぞれギャップのある系と金属での相対束縛エネルギーを 内殻ホール間距離の関数として示しています。 ギャップのある系では収束が遅く、収束値を得るために大きなスーパーセルが必要であることを意味しています。 一方、金属では小さな内殻ホール間距離において束縛エネルギーが速やかに収束していることが分かります。 図 78(c)に$2p$状態への内殻ホールの生成により誘起された固体シリコンでの差電荷密度を示します。 内殻ホールにより生成されたポテンシャルを遮蔽するため、内殻ホールから約$7$Åの周辺まで電荷の再配置が発生します。 ギャップのある系で見られた遅い収束は電荷の遠距離に渡る電荷の再配置が原因となっています。 他方、金属では電荷の再配置は比較的、狭い範囲に限定され、 このため図 78(b)に示したように早い収束を達成されています。

計算の他の例題と計算に用いた入力ファイルはウェブサイト
https://t-ozaki.issp.u-tokyo.ac.jp/vps_pao_core2019/ で公開されていますので、参考にして下さい。 また、本手法の応用例（シリコン、ボロフェン、グラフェン上に分散されたプラチナ単原子）が 文献[93,94,95] で議論されています。