アルミニウム結晶 †
このチュートリアルでは金属系に関するk点およびスメアリング関数を用いてフェルミ面を取り扱う際のぼかし幅について収束性のテストに関する説明を行います。
バンド構造および状態密度の計算方法も説明します。
収束性のテスト †
金属系では一般的に半導体・絶縁体に比べて多くのk点が必要になります。
さらにSTATEではフェルミ面を取り扱うために様々な関数系によるスメアリング法、あるいはテトラへドロン法の使用が可能です。
k点に関する収束性の調査 †
ここではMethfessel-Paxton法とテトラへドロン法を用いて全エネルギーのk点に関する収束を調べます。
スメアリング法を使用するためには負のぼかし幅 (WIDTH)を指定します。ぼかし幅の単位はRyです。
正のWIDTHを使用すると多項式のスメアリング関数が使用され、電子のエントロピーの項は考慮されません。
テトラへドロン法を使用する場合にはWIDTHを-10よりも小さい任意の数値にします。
- スメアリング関数(Methfessel-Paxton法)を使用する場合の入力ファイル
0 0 0 0 0 0 : I_CTRL(1:6) (DUMMY)
4.00 8.00 1 1 1 : GMAX GMAXP NTYP NATM NATM2
221 2 : NUM_SPACE_GROUP TYPE
7.5967 7.5967 7.5967 90.0 90.0 90.0 : A B C ALPHA BETA GAMMA
06 06 06 1 1 1 : N1 N2 N3 M1 M2 M3
0 0 : NCORD, NINV
0.00 0.00 0.00 1 0 1 : CPS(1,1:3) IWEI IMDTYP ITYP
13 0.50 26.98 6 1 0.2 : IATOMN ALFA AMION ILOC IVAN
0 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
0 1 : IPRE IPRI
30 30 0 84200.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP
6 1 : WAY_MIX MIX_WHAT
0 20 0.60 : ITER_START KBXMIX MIX_ALPHA
0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM
300.00 4 1 0.50D-09 : DTIO IMDALG IEXPL EDELTA
-0.0020 0.50D+03 0 : WIDTH FORCCR ISTRESS
ggapbe 1 : XCTYPE KSPIN
2.00 : DESTM
101 : NBZTYPE
4 4 4 : NKX NKY NKZ (DUMMY)
4 4 4 : NKX2 NKY2 NKZ2 (DUMMY)
6 : KEG
1 : NEXTST
0 : (DUMMY)
2 : IMSD
0 : EVALUATE_EKO_DIFF
0 : NPDOSAO
0 0.000 : SM_N DOPPING (DUMMY)
- テトラへドロン法を用いる場合の入力ファイル
0 0 0 0 0 0 : I_CTRL(1:6) (DUMMY)
4.00 8.00 1 1 1 : GMAX GMAXP NTYP NATM NATM2
221 2 : NUM_SPACE_GROUP TYPE
7.5967 7.5967 7.5967 90.0 90.0 90.0 : A B C ALPHA BETA GAMMA
06 06 06 1 1 1 : N1 N2 N3 M1 M2 M3
0 0 : NCORD, NINV
0.00 0.00 0.00 1 0 1 : CPS(1,1:3) IWEI IMDTYP ITYP
13 0.50 26.98 6 1 0.2 : IATOMN ALFA AMION ILOC IVAN
0 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
0 1 : IPRE IPRI
30 30 0 84200.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP
6 1 : WAY_MIX MIX_WHAT
0 20 0.60 : ITER_START KBXMIX MIX_ALPHA
0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM
300.00 4 1 0.50D-09 : DTIO IMDALG IEXPL EDELTA
-10.02 0 0.50D+03 0 : WIDTH FORCCR ISTRESS
ggapbe 1 : XCTYPE KSPIN
2.00 : DESTM
101 : NBZTYPE
4 4 4 : NKX NKY NKZ (DUMMY)
4 4 4 : NKX2 NKY2 NKZ2 (DUMMY)
6 : KEG
1 : NEXTST
0 : (DUMMY)
2 : IMSD
0 : EVALUATE_EKO_DIFF
0 : NPDOSAO
0 0.000 : SM_N DOPPING (DUMMY)
上記のファイルの違いはWIDTHの値だけです。
k点メッシュを変えながら計算を実行すると、k点の関数として全エネルギーを以下のように計算できます。
スメアリング幅に関する収束性のテスト †
金属系の全エネルギー(電子のエントロピーを含む自由エネルギー)はスメアリング関数のぼかし幅に敏感です。
特にガウス関数とフェルミ-ディラック関数の場合は顕著です。
以下ではde Gironcoliによって最初に示された全エネルギーのぼかし幅
依存性をアルミニウムについて確認します[Phys. Rev. B51, 6773(R) (1995)]。Marzariによる記事も大変参考になります。
デフォルトのMethfessel-Paxton (一次のHermite-Gaussian関数) 以外のスメアリング関数を使用する場合は&OTHERS ... &ENDというブロックを入力ファイルに追加し、使用したいスメアリング関数名を指定します。
フェルミ-ディラック分布関数の場合
&OTHERS
FERMI_DIRAC
&END
ガウシアン関数の場合
&OTHERS
GAUSSIAN
&END
MarzariとVanderbiltのcold smearingを使用する場合
&OTHERS
COLD
&END
とします。
全エネルギーをぼかし幅の関数として計算すると以下が得られます。
先に述べたようにフェルミ-ディラック関数とガウス関数を使用した場合、全エネルギーのぼかし幅依存性は著しく、Hermite-Gaussianあるいはcold smearingを使用した場合、全エネルギーはぼかし幅にほぼ依存性しないことが分かります。
一般に我々の計算は0Kで行われており、ぼかし幅0の極限での全エネルギーが我々の求めたい値と考えられます。
そのような目的のもとでHermite-Gaussianあるいはcold smearingが適切な選択肢と考えられます。
また一般にぼかし幅が大きければk点サンプリング数は少なくて済みます(k点数に関する収束性のテストは常に必要)。
ぼかし幅とk点数のバランスをうまくとり、効率的に計算を進めることができるパラメーターを選ぶ必要があります。
バンド計算 †
ここではバンド構造の計算方法を説明します。
バンド計算を行うためには、まず最初に電子密度をSCF計算により求める必要があります。
以下はテスト用に比較的少ないk点数を用いたSCF計算用の入力ファイルの例を示します。
- 入力ファイル (nfinp_scf)
0 0 0 0 0 0 : I_CTRL(1:6) (DUMMY)
4.00 8.00 1 1 1 : GMAX GMAXP NTYP NATM NATM2
221 2 : NUM_SPACE_GROUP TYPE
7.5967 7.5967 7.5967 90.0 90.0 90.0 : A B C ALPHA BETA GAMMA
06 06 06 1 1 1 : N1 N2 N3 M1 M2 M3
0 0 : NCORD, NINV
0.00 0.00 0.00 1 0 1 : CPS(1,1:3) IWEI IMDTYP ITYP
13 0.50 26.98 6 1 0.2 : IATOMN ALFA AMION ILOC IVAN
0 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
0 1 : IPRE IPRI
30 30 0 84200.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP
6 1 : WAY_MIX MIX_WHAT
0 20 0.60 : ITER_START KBXMIX MIX_ALPHA
0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM
300.00 4 1 0.50D-09 : DTIO IMDALG IEXPL EDELTA
-0.0020 0.50D+03 0 : WIDTH FORCCR ISTRESS
ggapbe 1 : XCTYPE KSPIN
2.00 : DESTM
101 : NBZTYPE
4 4 4 : NKX NKY NKZ (DUMMY)
4 4 4 : NKX2 NKY2 NKZ2 (DUMMY)
6 : KEG
1 : NEXTST
0 : (DUMMY)
2 : IMSD
0 : EVALUATE_EKO_DIFF
0 : NPDOSAO
0 0.000 : SM_N DOPPING (DUMMY)
次にSCF計算により得られた電荷密度を用いて、(電荷密度を固定した)非SCF (non-SCF) 計算を実行します。
バンド計算においてはブリルアンゾーン内の高対称点をつなぐ線上のk点を選びKohn-Sham方程式を解きます。
以下はnon-SCFバンド計算用の入力ファイルの例です。
- 入力ファイル (nfinp_band)
0 0 0 0 0 0 : I_CTRL(1:6) (DUMMY)
4.00 8.00 1 1 1 : GMAX GMAXP NTYP NATM NATM2
221 2 : NUM_SPACE_GROUP TYPE
7.5967 7.5967 7.5967 90.0 90.0 90.0 : A B C ALPHA BETA GAMMA
06 06 06 1 1 1 : N1 N2 N3 M1 M2 M3
0 0 : NCORD, NINV
0.00 0.00 0.00 1 0 1 : CPS(1,1:3) IWEI IMDTYP ITYP
13 0.50 26.98 6 1 0.2 : IATOMN ALFA AMION ILOC IVAN
22 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
0 1 : IPRE IPRI
30 30 0 84200.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP
6 1 : WAY_MIX MIX_WHAT
0 20 0.60 : ITER_START KBXMIX MIX_ALPHA
0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM
300.00 4 1 0.50D-09 : DTIO IMDALG IEXPL EDELTA
-0.0020 0.50D+03 0 : WIDTH FORCCR ISTRESS
ggapbe 1 : XCTYPE KSPIN
2.00 : DESTM
101 : NBZTYPE
4 4 4 : NKX NKY NKZ (DUMMY)
4 4 4 : NKX2 NKY2 NKZ2 (DUMMY)
6 : KEG
1 : NEXTST
0 : (DUMMY)
2 : IMSD
0 : EVALUATE_EKO_DIFF
0 : NPDOSAO
0 0.000 : SM_N DOPPING (DUMMY)
&KPOINTS_BAND
NKSEG 4
KMESH 40 20 20 20
KPOINTS
0.000 0.000 0.000
0.000 0.500 0.500
0.250 0.500 0.750
0.500 0.500 0.500
0.000 0.000 0.000
&END
バンド計算を実行する時は以下のようにICONDを22に指定します。
22 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
ICONDに加えて以下のように&KPOINTS_BAND ... &ENDというブロックを追加し、バンドを計算するブリルアンゾーン内の線を定義する高対称点を指定します。
&KPOINTS_BAND
NKSEG 4
KMESH 40 20 20 20
KPOINTS
0.000 0.000 0.000
0.000 0.500 0.500
0.250 0.500 0.750
0.500 0.500 0.500
0.000 0.000 0.000
&END
ここで
NKSEG
はk点のセグメント(区分、線)の数を指定します。
各セグメントの刻み幅は
KMESH
でセグメントの数だけ指定します。
高対称点は
KPOINTS
というキーワードに続いて指定します。
指定するk点の数はNKSEG+1となります。
この入力ファイルを用いてnon-SCFバンド計算を実行すると、出力ファイルに加えてk点と固有値が"energy.data"というファイルに出力されます。このデータは直接プロットすることができません。
バンドをプロットするために"energy2band"というプログラムが用意されています。
"energy2band"は以下のように実行します。
$ energy2band
計算に用いられたバンドの数、プロットするバンドの数、バンド計算の際に使用されたk点の数、エネルギーの基準(Hartree単位、一般にはFermi levelを指定)が聞かれるのでそれらを正しく入力すると"band.data"というファイルが生成されます。
これはgnuplotやxmgraceでプロット可能です。
以下はgnuplotスクリプトの例です。
set terminal postscript eps 'Helvetica' 22
set output 'band.eps'
emin=-12.5
emax=12.5
xmin=0.0
xmax=2.54176932
G=0.0
X=0.82709404
W=1.24064106
L=1.82548487
G2=2.54176932
offset=0.75
set xrange [xmin:xmax]
set yrange [emin:emax]
unset key
unset xtics
unset xlabel
set ylabel 'E-E_F (eV)'
set xzeroaxis
set arrow from X,emin to X,emax nohead
set arrow from W,emin to W,emax nohead
set arrow from L,emin to L,emax nohead
set label '{/Symbol G}' at G,emin-offset center
set label 'X' at X,emin-offset center
set label 'W' at W,emin-offset center
set label 'L' at L,emin-offset center
set label '{/Symbol G}' at G2,emin-offset center
plot 'band.data' using 1:2 with lines lw 3
これを"band.gp"というファイルに保存して
$ gnuplot band.gp
を実行すると"band.eps"というファイルが得られます。
練習問題 †
- 格子定数のk点およびスメアリング幅依存性を調べてみましょう。またスメアリング法とテトラへドロン法を用いて得られた格子定数を比較してみましょう。