3. 磁性を考慮した計算

前章では非磁性の計算を行いましたが，磁性を考慮した計算を行うことも可能です。ここでは体心立方格子の鉄結晶を例に，磁性を考慮した計算を行う方法を紹介します。例題の入力ファイルは wsl/Fe 以下に配置されています。計算の流れはシリコン結晶の場合と同様で，SCF計算を行い，結果得られた電荷密度を入力としてバンド構造や状態密度の計算を行います。同様の手続きでEV曲線の計算なども行うことができます。

3.1. SCF計算

磁性を考慮する場合も基本となるのはSCF計算です。例題の入力ファイルは wsl/Fe/scf の下にあります。

3.1.1. 入力データの確認

SCF計算の入力の内容は下記の通り(前章と同じ部分は一部省略)

accuracy{
  ksampling{
    mesh{ nx = 10, ny = 10, nz = 10 }
  }
}
structure{
  unit_cell{
    #units angstrom
    a = 2.845, b = 2.845, c = 2.845
    alpha = 90, beta = 90, gamma = 90
  }
  symmetry{
    method = automatic
    tspace{
      lattice_system = bodycentered
    }
  }
  magnetic_state = ferro
  atom_list{
    atoms{
      #tag rx ry rz element
      0.000 0.000 0.000 Fe
    }
  }
  element_list{
    #tag element atomicnumber zeta
    Fe 26 0.275
  }
}

原子配置の与え方としてBravais格子を与える方法を採用しています。単位胞としては体心立方格子の基本格子ではなく立方格子のそれをあたえ，体心立方格子であることはsymmetryのtspaceブロックにおけるlattice_system = bodycenteredで指定しています。このように指定すると，プログラムは格子ベクトルを自動的に基本格子のそれに変換します。原子座標の指定もそれに合わせ，原点のみの指定となっています。
structureブロックのmagnetic_stateにferroを指定しています。このキーワードはややミスリーディングですが，強磁性体に限られるわけではなく，スピンを考慮した計算を行う，という設定です。
element_listにおいてzetaの値を指定しています。zetaは“初期スピン分極”の設定です。スピン分極\(\zeta\)とは，アップスピンの電子数を\(n_{\uparrow}\)ダウンスピンの電子数を\(n_{\downarrow}\)とした場合\(\zeta = \frac{n_{\uparrow} - n_{\downarrow}}{n_{\uparrow} + n_{\downarrow}}\)と定義されます。
金属なので，Si結晶の場合と比較するとサンプリングk点を多くしています。

3.1.2. 計算の実行

これまでと同じ手続きで実行することができます。

mpirun -n 2 ~/phase0_2023.01/bin/phase ne=1 nk=2

以下の要因から，原子数自体は鉄結晶よりも多いシリコン結晶の計算より有意に時間のかかる計算となります。

スピンを考慮している
サンプリングk点の数が多い
バンド数（電子数）が多い
d電子を含んでいる

3.1.3. 結果の解析

スピンを考慮したので，スピン分極の値をもとめることができます。以下のようなコマンドによってアップスピンおよびダウンスピンの電子数の履歴を調べることができます。

$ grep NEW output001
!*--- input-file style = NEW
!NEW total charge (UP, DOWN, SUM) = 5.09024391 (+) 2.90975609 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.08351594 (+) 2.91648406 (=) 8.00000000
...
!NEW total charge (UP, DOWN, SUM) = 5.17484294 (+) 2.82515706 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.17483127 (+) 2.82516873 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.17483897 (+) 2.82516103 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.17483935 (+) 2.82516065 (=) 8.00000000

最終的にアップスピンの電子数が5.1748935ダウンスピンの電子数が2.82516065となりました。すなわち，スピン分極は0.294程度ということになります。

3.2. バンド計算

バンド計算は，Siの場合と同様固定電荷計算を用いて行います。例題の入力ファイルは wsl/Fe/band 以下に配置されています。

3.2.1. 入力データの確認

3.2.1.1. nfinp.dataファイル

Siの場合と同様，controlブロックにおいてcondition = fixed_chargeという設定がほどこされています。また，accuracyのksamplingブロックにおいてmethod = fileとすることによってk点生成方法をファイルからの読み込みとしています。

3.2.1.2. kpoint.dataファイル

Siの場合と同様kpoint.dataファイルを生成する必要があります。以下のコマンドを実行してください。

$ ~/phase0_2023.01/bin/band_kpoint.pl bandkpt.in
Distance of 1 = 1
Distance of 2 = 0.707106781186548
Distance of 3 = 0.5
Distance of 4 = 0.866025403784439
Distance of 5 = 0.707106781186548
division numbers = 100 70 50 86 70
i=0 -0.005 0.005 0.005
i=1 0.00714285714285714 -0.00714285714285714 0
i=2 0.005 0.005 -0.005
i=3 -0.00290697674418605 -0.00290697674418605 -0.00290697674418605
i=4 0 0 0.00714285714285714
1 : 0 0 0
    => 0/200 0/200 0/200
...

3.2.2. 計算の実行

固定電荷計算はekcalプログラムで行います。以下の要領で実行してみてください。

mpirun -n 2 ~/phase0_2023.01/bin/ekcal ne=1 nk=2

3.2.3. 計算結果の解析

バンド構造の計算結果はnfenergy.dataファイルに記録されます。Siの場合と違い，アップスピンとダウンスピン状態の固有値が両方とも記録されます。

band.plをたとえば以下のように実行します。

~/phase0_2023.01/bin/band.pl nfenergy.data bandkpt.in -erange=-10,10 -color -with_fermi

鉄の場合はエネルギーの範囲が広いので，-erangeオプションを用いてフェルミエネルギーを基準に-10 eVから10 eVの状態を対象としました。-color, -with_fermiオプションも有効にしました。

3.3. 状態密度計算

Siの場合と同様，固定電荷計算を用いて状態密度の計算を行います。対応する入力ファイルは wsl/Fe/dos 以下に配置されています。

3.3.1. 入力データの確認

3.3.1.1. nfinp.dataファイル

Siの場合と同様，accuracyブロックのksamplingブロックにおいてmethodをmeshとし，さらにsmearingブロックにおいてmethodをtetrahedralに設定しています。また，postprocessingブロックのdosブロックにおいてsw_dosをonとし，methodをtetrahedralとしています。

3.3.2. 計算の実行

固定電荷計算はekcalプログラムで行います。以下の要領で実行してみてください。

mpirun -n 2 ~/phase0_2023.01/bin/ekcal ne=1 nk=2

3.3.3. 計算結果の解析

スピン分極を考慮している場合，dos.dataファイルにはアップスピンの状態密度とダウンスピンの状態密度が記録されます。そのため，考慮していない場合と比較して1行のカラム数が増えます。具体的には，各行次のような並びでデータが記録されます。

エネルギーの識別子/ハートリー単位のエネルギー/ハートリー単位の場合のアップスピン電子の状態密度/ハートリー単位の場合のダウンスピン電子の状態密度/eV単位のエネルギー/eV単位の場合のアップスピン電子の状態密度/eV単位の場合のダウンスピン電子の状態密度/アップスピン電子の積算状態密度/ダウンスピン電子の積算状態密度/全電子の積算状態密度

eV単位の場合0がフェルミエネルギーになるようにシフトされている点はスピンを考慮していない場合と同じです。

dos.plをたとえば以下のように実行します。

~/phase0_2023.01/bin/dos.pl dos.data -with_fermi -color -erange=-10,10

バンド計算の場合と同様 -erangeオプションを用いてフェルミエネルギーを基準に-10 eVから10 eVの状態を対象としました。-color, -with_fermiオプションもシリコン結晶の場合と同様有効にしました。図 3.2 のような状態密度図が得られるはずです。