====== 概要 ======

  * 酸化ニッケル NiO を例に反強磁性相の計算方法を解説します。
  * スピン配置を指定する際にVESTAを使用します。


====== 手順 ======

===== 準備 =====

事前に NiO.struct を作成します。

  * 結晶構造は以下のサイトのCIFファイルを使用しました。
  * http://www.crystallography.net/cod/1010093.html

===== 超格子の作成 =====

（反強磁性相転移をしても単位胞が拡大しない場合はこの操作は不要です）

==== WIEN2kを使う場合 ====

  - supercellというプログラムを使って超格子を作成します。ここでは単位胞を各方向に2倍して、面心立方格子から単純立方格子に変更しています<code>
$ x_lapw supercell 
 Program generates supercell from a WIEN struct file.
 
 Filename of struct file: 
NiO.struct
 
 Number of cells in x direction: 
2
 Number of cells in y direction: 
2
 Number of cells in z direction: 
2
 Optional shift all atoms by the same amount (fractional coordinates). 
 Please enter x shift: 
0
 Please enter y shift: 
0
 Please enter z shift: 
0
 
 Current structure has lattice type F   
 Enter your target lattice type: (P,B,F)
P
 Target lattice type will be P   
 
 Add vacuum in x-direction for surface-slab [bohr]:
0
 Add vacuum in y-direction for surface-slab [bohr]:
0
 Add vacuum in z-direction for surface slab [bohr]:
0
 
 Supercell generated sucessfully.
 Stored in struct file: NiO_AFM_man_super.struct                                
               
 You may need to replace an atom by an impurity or distort the positions, ....
</code>
    * 超格子のstructファイルである NiO_super.struct ができます。 


==== VESTAを使う場合 ====

  - structファイルをVESTAで開きます
  - メニューから [Edit] -> [Edit Data] -> [Unit cell] -> [Transform]
  - 開いたウィンドウの右に書いてある説明を参考に、変換行列を指定します
    * 今回の場合は対角成分をすべて2にしてください
  - [Search Atoms ...] を選択します 

===== スピン配置の指定 =====

  - VESTAを起動して NiO_super.struct を開きます。
  - 反強磁性のスピン配置を指定します。ダウンスピンサイトとなるNi原子を別の元素（例えばFeなど）に変更します
    * [Edit]->[Edit Data]->[Structure Parameters]で変更できます。Symbolだけ変更してください。
    * 原子をダブルクリックすると、それが何番目の原子か表示されます。 
    * NiOのスピン配置は、例えば https://arxiv.org/abs/1806.07968 のFig.1を見てください
  - [File]->[Export]で NiO_super.cif という名前で保存してください。

下の図ができあがった超格子です。灰色がNi（アップスピンサイト）、茶色がFe（ダウンスピンサイト）、赤がOです。

{{ :wien2k:nio_super.png?400 |}}
===== 反強磁性のSCF計算 =====

  - 新たにプロジェクト NiO_AFM 作成します。この名前のフォルダを作成し、その中にさきほどの NiO_super.cif をコピーしてください。
  - このCIFファイルからstructファイルを作成します<code>
$ cif2struct NiO_super.cif
$ mv NiO_super.struct NiO_AFM.struct
</code>
  - 空間群の最適化を行います<code>
$ x_lapw sgroup -f NiO_AFM
$ mv -f NiO_AFM.struct_sgroup NiO.struct
</code>
    * structファイルを開いて空間群が166 R-3mになっていることを確認してください。なっていない場合は設定したスピン配置が間違っています。
  - ダウンスピンサイトの原子番号を正しく設定します
    * 例えばダウンスピンサイトにFeを使った場合はATOMを Fe2 -> Ni2 に、Zを 26 -> 28 にします
  - マフィンティン半径を設定します<code>
$ setrmt_lapw NiO
$ mv -f NiO.struct_setrmt NiO.struct
</code>
  - SCF計算のスピンの初期配置を指定します。Ni1をアップスピンサイト、Ni2をダウンスピンサイト、Oを非磁性サイトとします<code>
$ instgen -ask
 3 Atoms found:  with labels Ni1  Ni2  O 1
y
generate atomic configuration for atom 1 : Ni1
select spinpolarization up, dn or non-magnetic ( u, d, n )
u
generate atomic configuration for atom 2 : Ni2
select spinpolarization up, dn or non-magnetic ( u, d, n )
d
generate atomic configuration for atom 3 : O 1
select spinpolarization up, dn or non-magnetic ( u, d, n )
n
</code>
  - SCF計算の初期設定を行います。オプション-spでスピン分極を許可します<code>
$ init_lapw -b -numk 1000 -rkmax 8.0 -sp
</code>
    * 通常は反強磁性体の計算ではユニットセルが二倍になってスピンごとに計算する必要も出てくるので計算が非常に重くなります。最初は numk を10くらいから、結果が収束するまで少しずつ増やしていってください。
  - スピン分極ありのSCF計算を実行します。<code>
$ runsp_lapw -cc 0.0001 -ec 0.00001 -i 100 -p
</code>
    * 磁気モーメントやエネルギーギャップは次のようにしてチェックできます。磁気モーメントの単位はμBです<code>
$ grep ":MMTOT:" NiO.scf      # 全磁気モーメント。今は反強磁性なのでゼロ。
$ grep ":MMI001:" NiO.scf     # 1番目の原子（Niのアップスピンサイト）のモーメント
$ grep ":MMI002:" NiO.scf     # 2番目の原子（Niのダウンスピンサイト）のモーメント
$ grep ":MMI003:" NiO.scf     # 3番目の原子（O）のモーメント。単位はμB。NiOの場合はゼロ。
$ grep ":GAP:" NiO.scf        # エネルギーギャップ。NiOは反強磁性秩序によってギャップが開きます。
</code>

===== 物理量の計算 =====

スピン分極を許した計算なので、例えば次に状態密度を計算する場合などで x_lapw を実行する際には、アップ（ダウン）スピンのそれぞれで -up (-dn) を必ずつけてください。
<code>
$ x_lapw lapw1 -up
</code>

====== LDA+U ======

結果はアップスピンのモーメントが1.38μB、ギャップが0.95eVくらいになると思います。
これは実験で得られているもの（手法にもよりますが4eV前後）よりもだいぶ小さいです。
  * 参考: https://doi.org/10.1103/PhysRevB.38.3449
  * このページの例だとギャップの大きさはk点数に対して収束していないです。どのみち正確ではないので増やしてもあまり意味がありませんが...

実は第一原理計算は電子相関を無視した計算になっているため、電子相関の強いNiOではこのように実験に合わない結果となります。
これを補正するため、原子にオンサイトのクーロン相互作用Uの補正を入れるLDA+Uという方法を使います。

以下ではNiのd軌道に+Uの補正を入れてみましょう。


===== 手順 =====

事前にスピン分極のSCF計算を済ませてください。

  - LDA+U計算用の入力ファイルを作成します。今回はNiのd軌道(l=2)に U = 6 eV (6/13.606 = 0.4409...) のポテンシャルを入れます。Jはゼロにしました。<code>
$ init_orb_lapw -orb
 You have the following atoms:
1 : Ni1
2 : Ni2
3 : O
3 : O
Enter the name, l, U(Ry) and J of the Atoms (eg. Fe 2 0.3 0.0; exit with RETURN):Ni 2 0.441 0
Enter the name, l, U(Ry) and J of the Atoms (eg. Fe 2 0.3 0.0; exit with RETURN):
</code>
  - 入力ファイルが表示されるので、正しく設定されたか確認します
    * NiO_AFM.indm
      * 2行目: +Uを入れる原子の数
      * 3,4行: +Uを入れる原子のインデックス（NiO_AFM.structで確認。Ni1の場合は1、Ni2の場合は2になっているはず）、+Uをいれる軌道の数（今回はd軌道の1つ）、+Uを入れる軌道の種類（d軌道ならl=2）
      * 5行目: LDA+Uは 0 0 で固定。詳しくはUser's Guideの7.10節<file - NiO_AFM.indm>
-12.0         Emin cutoff
2               number of atoms
1  1 2      index of atom, number of l, l
2  1 2      index of atom, number of l, l
 0 0          r-index,(l,s) index
</file>
    * NiO_AFM.inorb
      * 1行目: nmodはLDA+U計算なら1、natorbは+Uする原子の数、 iprは1にすれば出力がより詳細になる
      * 2行目: 変更不可
      * 3, 4行目: iatomは+Uを入れる原子のインデックス、nlorbはその原子で+Uを入れる軌道の数、lorbはその軌道の軌道角運動量（今の場合はd電子なのでl=2）
      * 5行目: ポテンシャルの入れ方。詳しくはUser's guideの7.4節参照
        * AMF(0) ... 金属などの弱相関系に使用。
        * SIC(1) ... モット絶縁体などの強相関電子系に使用。
        * HMF(2) ... 平均場的にハバードUを入れる。現在は非推奨。
      * 6, 7行目: クーロン相互作用Uの値とフント結合Jの値。単位はRyなので注意<file - NiO_AFM.inorb>
 1 2 0        nmod, natorb, ipt
PRATT   1.0  
1  1 2      index of atom, number of l, l
2  1 2      index of atom, number of l, l
  1               nsic 0..AMF, 1..SIC, 2..HFM
 0.441   0      U J (Ry)
 0.441   0      U J (Ry)
</file>
    * https://wien.zeus.theochem.tuwien.ac.narkive.com/dV0E1fCH/how-to-include-d-and-f-states-in-calculation
  * ブロイデン法の履歴を削除して、オプション-orbをつけて実行します<code>
$ rm *.broyd*
$ runsp_lapw -cc 0.0001 -ec 0.00001 -i 100 -orb -p -NI
</code>
    * k点数やUをかえて再度実行するときは、ブロイデン法の履歴を削除しなくてもオプション -NI をつけるだけで十分です。
    * Uを変更するときは NiO_AFM.inorb を直接編集してください。

今回の計算結果はバンドギャップが約3.29eV、モーメントが約1.72μBになりました。

===== 物理量の計算 =====

+UでSCF計算が終わったあと、状態密度の計算などで x_lapw lapw1 を実行するときは、オプション-orbをつけてください。
<code>
$ x_lapw lapw1 -up -orb
</code>