====== 概要 ======
一般にはDFT+Uを使ってバンドギャップの精度を増すことができるが、場合によってはそのまま計算しただけでは基底状態にならないこともある。
ここでは例として反強磁性体FeOが単純にはうまく計算できないことと、その解決方法を示す。
このページは以下のスライドを参考に作成したので、詳しく知りたい場合はそちらを参照してください。
* http://www.iiserpune.ac.in/~smr2626/hands_on/week1/july4/hands-on_dft_pu_last.pdf
なお、FeOのバンドギャップは光学吸収測定により2.4eVと見積もられている。
* H. K. Bowen, D. Adler, and B. H. Auker, J. Solid State Chem. 12, 355 (1975).
====== DFT+U ======
まずは普通に計算してみよう。
&control
calculation='scf',
prefix='FeO',
pseudo_dir = './',
outdir='./tmp/'
/
&system
ibrav=0,
celldm(1)=8.19,
nat=4,
ntyp= 3,
ecutwfc = 30.0, ecutrho = 240.0,
occupations='smearing',
smearing='mp',
degauss=0.02,
nspin=2,
starting_magnetization(2)= 0.5,
starting_magnetization(3)=-0.5,
lda_plus_u = .true.,
Hubbard_U(2) = 4.0,
Hubbard_U(3) = 4.0,
/
&electrons
/
CELL_PARAMETERS
0.50 0.50 1.00
0.50 1.00 0.50
1.00 0.50 0.50
ATOMIC_SPECIES
O 16.00 O.pbe-rrkjus.UPF
Fe1 55.85 Fe.pbe-nd-rrkjus.UPF
Fe2 55.85 Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
O 0.25 0.25 0.25
O 0.75 0.75 0.75
Fe1 0.0 0.0 0.0
Fe2 0.5 0.5 0.5
K_POINTS automatic
4 4 4 1 1 1
&control
calculation='nscf',
prefix='FeO',
pseudo_dir = './',
outdir='./tmp/'
/
&system
ibrav=0,
celldm(1)=8.19,
nat=4,
ntyp= 3,
ecutwfc = 30.0, ecutrho = 240.0,
occupations = 'tetrahedra',
nspin=2,
starting_magnetization(2)= 0.5,
starting_magnetization(3)=-0.5,
lda_plus_u = .true.,
Hubbard_U(2) = 4.0,
Hubbard_U(3) = 4.0,
/
&electrons
/
CELL_PARAMETERS
0.50 0.50 1.00
0.50 1.00 0.50
1.00 0.50 0.50
ATOMIC_SPECIES
O 16.00 O.pbe-rrkjus.UPF
Fe1 55.85 Fe.pbe-nd-rrkjus.UPF
Fe2 55.85 Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
O 0.25 0.25 0.25
O 0.75 0.75 0.75
Fe1 0.0 0.0 0.0
Fe2 0.5 0.5 0.5
K_POINTS automatic
8 8 8 0 0 0
&dos
prefix='FeO',
outdir='./tmp/',
DeltaE=0.01
/
$ pw.x < FeO.scf.in > FeO.scf.out
$ pw.x < FeO.nscf.in > FeO.nscf.out
$ dos.x < FeO.dos.in > FeO.dos.out
状態密度を描いてみればわかるように金属となってしまう。
実は、基底状態ではない準安定なところに収束してしまったため、このような結果になってしまう。
====== 問題の検討 ======
FeO.scf.outで一回目のSCF計算が終わった部分を見てみると、次のようになっている。
(略)
iteration # 1 ecut= 30.00 Ry beta=0.70
Davidson diagonalization with overlap
ethr = 1.00E-02, avg # of iterations = 2.2
--- enter write_ns ---
LDA+U parameters:
U( 2) = 4.00000000
alpha( 2) = 0.00000000
U( 3) = 4.00000000
alpha( 3) = 0.00000000
atom 3 Tr[ns(na)] (up, down, total) = 4.99738 1.04703 6.04441
spin 1
eigenvalues:
0.999 0.999 0.999 1.000 1.000
eigenvectors:
0.000 0.018 0.016 0.527 0.438
0.333 0.589 0.054 0.021 0.002
0.333 0.033 0.610 0.001 0.022
0.000 0.016 0.018 0.438 0.527
0.333 0.343 0.301 0.013 0.010
occupations:
1.000 -0.000 -0.000 -0.000 -0.000
-0.000 0.999 -0.000 -0.000 0.000
-0.000 -0.000 0.999 0.000 0.000
-0.000 -0.000 0.000 1.000 -0.000
-0.000 0.000 0.000 -0.000 0.999
spin 2
eigenvalues:
0.166 0.166 0.204 0.256 0.256
eigenvectors:
0.513 0.459 0.000 0.001 0.026
0.017 0.001 0.333 0.576 0.073
0.001 0.017 0.333 0.376 0.273
0.459 0.513 0.000 0.026 0.001
0.010 0.009 0.333 0.021 0.627
occupations:
0.169 -0.006 -0.006 -0.000 -0.012
-0.006 0.237 -0.016 -0.010 0.016
-0.006 -0.016 0.237 0.010 0.016
-0.000 -0.010 0.010 0.169 -0.000
-0.012 0.016 0.016 -0.000 0.237
atomic mag. moment = 3.950343
(略)
ここで注目したいのが、atom 3(Feのアップスピンサイト)のspin 2(ダウンスピン)のeigenvaluesである。
spin 2
eigenvalues:
0.166 0.166 0.204 0.256 0.256
固有値(eigenvalue)が大きいほどエネルギー準位が低い。
もともとの格子は立方対称なので準位の低い三重縮退のt2g軌道と、高い二重縮退のeg軌道に分裂している。
反強磁性秩序によってt2g軌道の対称性が壊れて二重縮退の準位(0.256)と縮退のない準位(0.204)に分裂している。
Feは価電子を6つ持っているので、フントの規則によりegとt2gのすべてにアップスピンの電子がすべて詰まったのち、残りの1つがダウンスピンとしてどこかの準位に入る。
すなわち、もしも系が絶縁体なのであれば、このダウンスピンの準位は縮退のない準位に入らなければならない(縮退のあるバンドに電子が1つだけ入れば、バンドが中途半端に満たされるので金属となる)ので、縮退のない準位はもっとも大きい固有値にならなければならない。
しかし、FeO.scf.outの結果を見てわかるように、縮退のない準位は固有値がもっとも小さくなってしまっている。
(略)
End of self-consistent calculation
--- enter write_ns ---
LDA+U parameters:
U( 2) = 4.00000000
alpha( 2) = 0.00000000
U( 3) = 4.00000000
alpha( 3) = 0.00000000
atom 3 Tr[ns(na)] (up, down, total) = 4.98974 1.73816 6.72790
spin 1
eigenvalues:
0.996 0.996 0.999 1.000 1.000
eigenvectors:
0.680 0.304 0.000 0.004 0.012
0.009 0.002 0.333 0.656 0.000
0.000 0.011 0.333 0.164 0.492
0.304 0.680 0.000 0.012 0.004
0.007 0.003 0.333 0.164 0.492
occupations:
0.996 -0.000 -0.000 -0.000 -0.000
-0.000 0.999 -0.000 -0.000 0.000
-0.000 -0.000 0.999 0.000 0.000
-0.000 -0.000 0.000 0.996 -0.000
-0.000 0.000 0.000 -0.000 0.999
spin 2
eigenvalues:
0.099 0.221 0.221 0.599 0.599
eigenvectors:
0.000 0.544 0.342 0.008 0.105
0.333 0.066 0.010 0.555 0.036
0.333 0.002 0.074 0.287 0.303
0.000 0.342 0.544 0.105 0.008
0.333 0.047 0.029 0.044 0.547
occupations:
0.264 -0.049 -0.049 -0.000 -0.098
-0.049 0.403 -0.152 -0.085 0.152
-0.049 -0.152 0.403 0.085 0.152
-0.000 -0.085 0.085 0.264 -0.000
-0.098 0.152 0.152 -0.000 0.403
atomic mag. moment = 3.251583
(略)
このように、本来ほしい状態と違う状態に収束してしまうのを避けるため、starting_ns_eigenvalueというものを使う。
====== starting_ns_eigenvalue ======
FeO.scf.inをコピーし、$systemにstarting_ns_eigenvalueを追加した以下の入力ファイルを作成する。
&control
calculation='scf',
prefix='FeO',
pseudo_dir = './',
outdir='./tmp/'
/
&system
ibrav=0,
celldm(1)=8.19,
nat=4,
ntyp= 3,
ecutwfc = 30.0, ecutrho = 240.0,
occupations='smearing',
smearing='mp',
degauss=0.02,
nspin=2,
starting_magnetization(2)= 0.5,
starting_magnetization(3)=-0.5,
lda_plus_u = .true.,
Hubbard_U(2) = 4.0,
Hubbard_U(3) = 4.0,
starting_ns_eigenvalue(3,2,2)=1.0,
starting_ns_eigenvalue(3,1,3)=1.0,
/
&electrons
/
CELL_PARAMETERS
0.50 0.50 1.00
0.50 1.00 0.50
1.00 0.50 0.50
ATOMIC_SPECIES
O 16.00 O.pbe-rrkjus.UPF
Fe1 55.85 Fe.pbe-nd-rrkjus.UPF
Fe2 55.85 Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
O 0.25 0.25 0.25
O 0.75 0.75 0.75
Fe1 0.0 0.0 0.0
Fe2 0.5 0.5 0.5
K_POINTS automatic
4 4 4 1 1 1
ここで、starting_ns_eigenvalueの意味は次のとおりである。
^変数^説明^
|starting_ns_eigenvalue(m,i,n)|n番目の原子のi番目のスピンの、m番目の固有値にある電子数の初期値|
ここでm=3としているのは、FeO.scf.outで1回目のSCF計算で縮退していない固有値が左から3番目だったからである。
縮退のない準位に電子が入ってほしかったことを思い出そう。
$ pw.x < FeO_2.scf.in > FeO_2.scf.out
さきほどの全エネルギーと比較すると、今回のほうが安定になっていることがわかる。
! total energy = -175.09491797 Ry
! total energy = -175.17890387 Ry
状態密度の計算はさきほどと同様に行う。
ギャップが開いていることがわかる。
{{:quantumespresso:feo_u_dos.png?400|}}