PWscfでSiのフォノンバンド

5月14日(木)にCCMS Webハンズオン: MateriApps LIVE!講習会を受講させていただきました。この講習会の自由演習で、シリコンのフォノン分散の計算方法を習ったので、今回はその復習をします。

001_20200622164420840.png

Fig.1: Siのフォノン分散



Hands on phonoのエクセサイズ


実習の内容はQuantum EspressoでSiのフォノンバンド図を作成するに書かれているとおりです。この内容自体も2014年のSummer Schoolの資料を基にしています。そのHands on Phononには4つのエクセサイズが収められています。
  1. SiのΓ点フォノン
  2. Siのフォノン分散
  3. AlAsのΓ点フォノン
  4. Alのフォノン分散

なおHands_on_phonons.tar.gzの中のPDFファイルに解説があるので自習が可能です。CCMS Webハンズオンで行ったのはエクセサイズ2のSiのフォノン分散計算です。

Hands-on-Phononの準備


Hands-on-Phonon全体の準備として、ファイルのダウンロードと展開を行います。計算に必要な擬ポテンシャルファイルもついでにダウンロードしておきます。
これらを行うコマンドは以下の通りです。

mkdir qe-ph-test
cd qe-ph-test
wget http://www.iiserpune.ac.in/~smr2626/hands_on/week2/july9/Hands_on_phonons.tar.gz
tar xzvf Hands_on_phonons.tar.gz
mkdir pseudo
wget -O ./pseudo/Si.pbe-rrkj.UPF http://www.quantum-espresso.org/wp-content/uploads/upf_files/Si.pbe-rrkj.UPF
wget -O ./pseudo/Al.pz-vbc.UPF https://www.quantum-espresso.org/upf_files/Al.pz-vbc.UPF
wget -O ./pseudo/As.pz-bhs.UPF https://www.quantum-espresso.org/upf_files/As.pz-bhs.UPF
wget -O ./pseudo/Al.pbe-rrkj.UPF https://www.quantum-espresso.org/upf_files/Al.pbe-rrkj.UPF
cd Hands_on_phonons


Siのフォノンバンド計算の概要


今回のエントリの本編であるSiのフォノンバンド計算のコマンドをいかに示します。

cd exercise-2/
mpirun -np 4 -x OMP_NUM_THREADS=1 pw.x < si_scf.in | tee si_scf.out
sed -i "1itestsi" si_ph.in
mpirun -np 4 -x OMP_NUM_THREADS=1 ph.x < si_ph.in | tee si_ph.out
mpirun -np 4 -x OMP_NUM_THREADS=1 q2r.x < q2r.in | tee q2r.out
mpirun -np 4 -x OMP_NUM_THREADS=1 matdyn.x < matdyn_Si_disp.in |tee matdyn_Si_disp.out
plotband.x
si.freq
0 500
siph.plot
siph.ps
0
100 0

plotband.xは計算したフォノンバンドの数値データから対話的にグラフを作成するプログラムです。plotband.x以降の行は、このプログラムに対話的に与えるパラメータを示しています。

上記の手順で実際に行っているのは
  1. セルフコンシステント計算(電子状態計算)
  2. フォノン計算
  3. q→r計算(動的マトリクスのフーリエ変換)
  4. フォノン分散計算
  5. グラフの作成

です。通常の電子のバンド構造計算や電子の状態密度計算では、セルフコンシステント計算を行ったのちポストスクリプト処理としてバンド構造や状態密度の計算を行います。フォノンの計算では1-3の計算を行ったのちに、フォノンのバンド構造計算やフォノンの状態密度計算などのポストスクリプト処理を行うイメージです。

1. セルフコンシステント計算(電子状態計算)


セルフコンシスト計算は、通常のものです。フォノンの計算であっても電子状態計算は必要です。
以下がセルフコンシステント計算用の入力ファイル si_scf.in です。
&control
calculation = 'scf'
restart_mode='from_scratch',
prefix='si',
tstress = .true.
tprnfor = .true.
pseudo_dir = '../../pseudo',
outdir='tmp'
/
&system
ibrav= 2, nat= 2, ntyp= 1,celldm(1)=10.3463,
ecutwfc =20.0, ecutrho=80.0,nbnd=8
occupations='fixed',
/
&electrons
mixing_mode = 'plain'
mixing_beta = 0.5
conv_thr = 1.0d-10
/
&ions
/
&cell
/

ATOMIC_SPECIES
Si 28.086 Si.pbe-rrkj.UPF

ATOMIC_POSITIONS {crystal}
Si 0.00000 0.00000 0.00000
Si 0.25000 0.25000 0.25000


K_POINTS {automatic}
6 6 6 1 1 1


実行ファイル pw.x で計算を行います。
mpirun -np 4 -x OMP_NUM_THREADS=1 pw.x < si_scf.in | tee si_scf.out


2. フォノン計算


以下がフォノン計算用の入力ファイル si_ph.in です。
nq1=4,nq2=4,nq3=4 がq点の分割数です。fildyn='si.dynmat',は計算結果の保存先ファイル名で、実際にはこの文字列の後に連番のついたファイルが生成されるようです(si.dynmat0, sidynmat1...など)。これらのファイルを次のq→r計算で利用するようです。tr2_ph=1.0d-14,が収束条件で 1.0d-14 より小さくする方がよいとの事。

testsi
&inputph
outdir='./tmp/',
prefix='si',
ldisp=.true.
nq1=4,nq2=4,nq3=4
fildyn='si.dynmat',
tr2_ph=1.0d-14,
amass(1)=28.086
/


上記の入力ファイルには既に書き込まれていますが Hands_on_phonons.tar.gz の入力ファイルでは、最初の行にタイトルが書かれていません。テキストエディタで行頭に適当な文字列(今回の場合は"testsi")を書き込む必要があります。今回はテキストエディタではなく sed コマンドを用いてファイルの先頭行に文字列を挿入してみました(参考:sedコマンドでファイルの先頭行に文字列を挿入する)。フォノン計算の実行ファイルは、電子のセルフコンシステント計算のときとは異なり ph.x を用います。
sed -i "1itestsi" si_ph.in
mpirun -np 4 -x OMP_NUM_THREADS=1 ph.x < si_ph.in | tee si_ph.out


3. q→r計算(動的マトリクスのフーリエ変換)


q空間(波数空間)からr空間(実空間)への変換...との事です。

以下が入力ファイル。短いですが理解しておくべきことは多いです。
まず fildyn='si.dynmat', が入力ファイルです。フォノン計算で指定した出力ファイルと同じ文字列を書きます。次に flfrc='si444.fc' が出力ファイルで、次のフォノン分散計算やフォノン状態密度計算の入力ファイルとして使われます。zasr='simple' が Acoustic sum rule (ASR) の設定です。今回は simple を採用していますが、この設定だとΓ点の音響フォノンが(厳密には)ゼロになりません。普通の三次元結晶の場合は zasr='crystal' とする方がよいとされています(参考: Quantum ESPRESSO入力ファイル作成手順14.Phonon計算(バンド構造))。
&input
fildyn='si.dynmat',
zasr='simple'
flfrc='si444.fc'
/


実行ファイルは q2r.x です。
mpirun -np 4 -x OMP_NUM_THREADS=1 q2r.x < q2r.in | tee q2r.out


4. フォノン分散計算


ポストスクリプト処理として、ついにフォノンの分散関係(バンド構造)の計算を行います。

まず前述の ASR の設定が asr='simple', です。q→rの計算で crystal を指定した場合は、こちらもそろえた方がよいでしょう。flfrc='si444.fc' が、q→rの計算結果の入力ファイル。flfrq='si.freq'が出力ファイル名の指定です。ファイル後半の数字の羅列は、計算する q点の経路です。5 と書いてあるのが、通る特徴的な q点の数(つまりこの後に何行あるか)で、それに続く各行が q点の座標とその分割数です。下記の例の場合 Γ→X→W→L→Γ という経路をそれぞれ10分割するという意味です。
 &input
asr='simple',
amass(1)=28.0855,
flfrc='si444.fc', flfrq='si.freq'
q_in_band_form=.true.
/
5
0.0000000 0.0000000 0.0000000 10
0.0000000 0.0000000 1.0000000 10
0.5000000 0.0000000 1.0000000 10
0.5000000 0.5000000 0.5000000 10
0.0000000 0.0000000 0.0000000 1


実行ファイルは matdyn.x です。
mpirun -np 4 -x OMP_NUM_THREADS=1 matdyn.x <  matdyn_Si_disp.in |tee matdyn_Si_disp.out

実行すると si.freqsi.freq.gp という二つのファイルが出力されます。後者の gp ファイルの方が gnuplot などでプロットするのに向いている数値ファイルです。

5. グラフの描画


ここまででフォノンの分散関係の計算結果が数値データとして得られています。この si.freq.gp の冒頭部分を以下に示します。(話がそれますがΓ点の音響フォノンが厳密にゼロになっていますね...たぶん私がASRにcrystalを指定して計算したのだと思います...このブログ通りに計算するとちょっと違った値になるかもしれません。)
  0.000000      -0.0000   -0.0000    0.0000  499.9411  499.9411  499.9411
0.100000 30.1467 30.1467 50.2804 497.9306 497.9306 499.4990
0.200000 59.2966 59.2966 99.5601 491.9361 491.9361 497.8325
0.300000 86.2601 86.2601 147.0152 482.3522 482.3522 494.2054
0.400000 109.5528 109.5528 192.0843 470.4851 470.4851 488.0585
0.500000 127.5939 127.5939 234.4387 458.6896 458.6896 479.3392
0.600000 139.4024 139.4024 273.8876 449.5627 449.5627 468.3264
0.700000 145.4470 145.4470 310.2766 444.5044 444.5044 455.1737
0.800000 147.6191 147.6191 343.4150 439.6134 442.9610 442.9610
0.900000 148.0519 148.0519 373.0548 421.0602 443.1545 443.1545
1.000000 148.0663 148.0663 398.9590 398.9590 443.4316 443.4316
1.050000 149.2370 149.2370 397.5904 397.5904 443.8141 443.8141
1.100000 152.6482 152.6482 393.6993 393.6993 444.7982 444.7982
1.150000 158.0147 158.0147 387.7863 387.7863 446.0406 446.0406
1.200000 164.9082 164.9082 380.4216 380.4216 447.2434 447.2434
1.250000 172.7949 172.7949 372.1743 372.1743 448.2391 448.2391

1列目はq点の経路の道のりで 2-6列目がフォノンの振動数(cm-1)です。基本的にはこれをそのままプロットすればよいのですが、きれいな図にするためには q点の道のりと特徴的なq点の記号の対応関係を知る必要があります。別に何も難しくないのですが、わざわざそれ用の plotband.x という対話式の実行ファイルが用意されているのでそれを利用します。以下のように順次入力していきます。
plotband.x
si.freq
0 500
siph.plot
siph.ps
0
100 0

ここで重要なのが0 500を入力した後に出力される以下の部分です。
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 0.0000
high-symmetry point: 0.0000 0.0000 1.0000 x coordinate 1.0000
high-symmetry point: 0.5000 0.0000 1.0000 x coordinate 1.5000
high-symmetry point: 0.5000 0.5000 0.5000 x coordinate 2.2071
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 3.0731

この x coordinate の後に出力されるのが、その前に出力される特徴的な q点に対応する経路の道のりです。

これを踏まえて、フォノンバンド分散をプロットする gnuplot スクリプトを以下のように作成しました。
set terminal pngcairo size 520,390
set output "siph.png"

set xrange [0:3.073132]
set yrange [0:500]

set grid x

set ylabel "Frequency (cm^{-1})"
set xtics ("{/Symbol G}" 0.0000, "{X}" 1.0000, "{W}" 1.5000, "{L}" 2.2071, "{/Symbol G}" 3.0731)

# plot "si.freq.gp" u 1:2 w l lc 1 not,\
# "si.freq.gp" u 1:3 w l lc 1 not,\
# "si.freq.gp" u 1:4 w l lc 1 not,\
# "si.freq.gp" u 1:5 w l lc 1 not,\
# "si.freq.gp" u 1:6 w l lc 1 not,\
# "si.freq.gp" u 1:6 w l lc 1 not,\
# "si.freq.gp" u 1:7 w l lc 1 not
plot "siph.plot" u 1:2 w l lc 1 not


6. フォノン状態密度計算


ついでにフォノンの状態密度(DOS)の計算方法について書きます。Hands_on_phonons.tar.gzの中に入力ファイルが収められているにもかかわらず付属のPDFなどでも触れられていませんが、簡単です。以下が入力ファイルです。nk1=4, nk2=4, nk3=4,がq点の分割数なのですが 4*4*4 だと心もとない雰囲気を感じたので、私は 16*16*16 で計算しました。
 &input
asr='simple',
amass(1)=28.0855,
flfrc='si444.fc', flfrq='si.freq.dos'
dos=.true.,
fldos='si.dos'
deltaE=1.d0,
nk1=4, nk2=4, nk3=4,
/

実行ファイルはフォノンの分散関係で用いたのと同じ matdyn.x です。
mpirun -np 4 -x OMP_NUM_THREADS=1 matdyn.x <  matdyn_Si_dos.in | tee matdyn_Si_dos.out

プログラムを実行すると si.dos というファイルが出力されます。これが状態密度の数値データを格納したファイルです。
これをプロットするために以下のような gnuplot のスクリプトを書きました。
set terminal pngcairo size 520,390
set output "si.phdos.png"

set xrange [0:500]
set xlabel "Frequency (cm^{-1})"

set yrange [0:0.1]
set ylabel "DOS"

plot "si.dos" u 1:2 w l lc 1 not

するとFig.2のような画像が得られます。

002_20200622164421f48.png
Fig.2: フォノンの状態密度(DOS)


参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: PWscf フォノン 

PWscfでMg2Siのゼーベック係数

UbuntuにBoltzTraPをインストールではUbuntuでPWscf(Quantum ESPRESSO)と連携することを前提にBoltzTraPをインストールしました。今回は Mg2Si を例にゼーベック係数を計算してみます。

300K.png
Fig.1: Mg2Si の 300 K でのゼーベック係数のエネルギー依存性



構造緩和計算


PWscfの入力作成補助にあるように xtl2pw.py がインストールされているという前提でMg2Si の CIF ファイルから PWscf (Quantum ESPRESSO) の入力ファイルを作成します。

最初にMg2SiのCIFファイルを入手します。今回はMaterials projectのページから Primitive Cell のモノをダウンロードしました。このCIFファイルを VESTA で開いて xtl 形式でエクスポートします(Mg2Si.xtl)。

DFT計算用のマシンに Mg2Si/ ディレクトリを作成し、そこに Mg2Si.xtl を置きます。これを xtl2pw.py に渡して入力ファイルのひな型を作ります。

xtl2pw.py Mg2Si.xtl vc-relax 1


作成された入力ファイルのひな型を基に、構造緩和計算用の入力ファイルを作成します。

mv Mg2Si.in Mg2Si.relax.in
emacs -nw Mg2Si.relax.in


擬ポテンシャルファイルやカットオフを正しく設定することに加えて system に以下を追記します。

    occupations = 'smearing' ,
degauss = 0.005 ,


私の環境では、以下のような入力ファイルになりました。

&control
calculation='vc-relax' ,
restart_mode='from_scratch' ,
prefix='Mg2Si' ,
outdir = './Mg2Si/' ,
wfcdir = './Mg2Si/' ,
pseudo_dir = './' ,
disk_io='default' ,
forc_conv_thr= 0.001 ,
verbosity = 'default' ,
nstep = 100 ,
/
&system
ibrav= 5 ,
celldm(1) = 8.50304569 ,
celldm(4) = 0.50000000 ,
nat = 3 ,
ntyp = 2 ,
ecutwfc = 50.0 ,
ecutrho = 250.0 ,
occupations = 'smearing' ,
degauss = 0.005 ,
/
&electrons
electron_maxstep = 100 ,
mixing_beta = 0.7 ,
! use smaller conv_thr for better results ,
conv_thr = 1.0d-12 ,
/
&ions
ion_dynamics='bfgs' ,
/
&CELL
cell_dynamics = 'bfgs' ,
press = 0.001,
press_conv_thr = 0.05 ,
! cell_dofree = 'xyz' ,
/
ATOMIC_SPECIES
Mg 24.305 Mg.pbe-n-kjpaw_psl.0.3.0.upf
Si 28.086 Si.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
Mg 0.750000 0.750000 0.750000 0 0 0
Mg 0.250000 0.250000 0.250000 0 0 0
Si 0.000000 0.000000 0.000000 0 0 0
K_POINTS automatic
4 4 4 0 0 0


これを pw.x に渡して構造緩和計算を行います。

mpirun -np 4 pw.x < Mg2Si.relax.in > Mg2Si.relax.out


nscf計算


次にnscf計算を行います。まずは構造緩和計算の入力ファイルを基にnscf計算の入力ファイルを作成します。

cp Mg2Si.relax.in Mg2Si.nscf.in
emacs -nw Mg2Si.nscf.in


変更点は以下です。
    calculation='nscf' ,
verbosity = 'high' ,
K_POINTS automatic
10 10 10 0 0 0


k点はもっと細かく採った方がいいかもしれませんが、今回は練習なのでこの程度にしておきます。

&control
calculation='nscf' ,
restart_mode='from_scratch' ,
prefix='Mg2Si' ,
outdir = './Mg2Si/' ,
wfcdir = './Mg2Si/' ,
pseudo_dir = './' ,
disk_io='default' ,
forc_conv_thr= 0.001 ,
verbosity = 'high' ,
nstep = 100 ,
/
&system
ibrav= 5 ,
celldm(1) = 8.50304569 ,
celldm(4) = 0.50000000 ,
nat = 3 ,
ntyp = 2 ,
ecutwfc = 50.0 ,
ecutrho = 250.0 ,
occupations = 'smearing' ,
degauss = 0.005 ,
/
&electrons
electron_maxstep = 100 ,
mixing_beta = 0.7 ,
! use smaller conv_thr for better results ,
conv_thr = 1.0d-12 ,
/
&ions
ion_dynamics='bfgs' ,
/
&CELL
cell_dynamics = 'bfgs' ,
press = 0.001,
press_conv_thr = 0.05 ,
! cell_dofree = 'xyz' ,
/
ATOMIC_SPECIES
Mg 24.305 Mg.pbe-n-kjpaw_psl.0.3.0.upf
Si 28.086 Si.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
Mg 0.750000 0.750000 0.750000 0 0 0
Mg 0.250000 0.250000 0.250000 0 0 0
Si 0.000000 0.000000 0.000000 0 0 0
K_POINTS automatic
10 10 10 0 0 0


これを実行します。

mpirun -np 4 pw.x < Mg2Si.nscf.in > Mg2Si.nscf.out


BoltzTraPの実行


BoltzTraP に渡すためのファイルを作成します。最初にnscf計算の結果からフェルミエネルギーを読み出します。

grep "Fermi" Mg2Si.nscf.out


つぎに qe2boltz.py を用いて BoltzTraP に与えるバンド構造の数値データや計算の設定ファイルを作成します。

qe2boltz.py Mg2Si pw 4.5310 0


4.5310の部分は、先ほど読みだしたフェルミエネルギーの値です。
これを走らせると Mg2Si.intrans や Mg2Si.energy が出力されます。
Mg2Si.intrans が BoltzTraP に与える計算の設定が書かれたファイルなので、適宜編集します。

emacs -nw Mg2Si.intrans


編集したら以下のように BoltzTraP を実行します。

x_trans BoltzTraP


300Kのデータのプロット


色々なファイルが出力されますが、もっとも重要なのは Mg2Si.trace です。
開いてみると1列目にエネルギー(フェルミエネルギーと書いてありますが、ただのエネルギーだと思います、状態密度の横軸とおなじ意味での)、2列目に温度、3列目以降に色々な物性が出力されており、今回の目的のゼーベック係数は5列目にあります。
今回は300Kにおけるゼーベック係数をエネルギーの関数としてプロットしてみます。そのままだとプロットしにくいので300Kのデータだけ抜き出します。

awk '{if($2 == 300){print $0}}' Mg2Si.trace > 300K.txt


こうして得られたデータをプロットするためのgnuplot用のpltファイルを下記のように準備しました。

set terminal pngcairo size 520,390
set output "300K.png"

set ylabel "Seebeck coefficient (uV/K)"
set title "Seebeck coefficient of Mg2Si at 300 K"

set xlabel "Energy (eV)"
set xrange [-2:2]
plot "300K.txt" u (($1-0.334102802189)*13.6058):($5*1E6) w l not


関連エントリ




参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: PWscf QuantumESPRESSO ゼーベック係数 

PWscfでシリコンのバンド構造

3. グラフェンのバンド構造と状態密度を参考に Quantum ESPRESSO(PWscf)を用いてシリコンのバンド構造を計算してみました。

Si-band.png
Fig.1: シリコンのバンド構造



セルフコンシステント計算


以下のようにシリコンの入力ファイルを作ります。
何らかの作成補助を使うとラクかもしれません(参考: PWscfの入力作成補助)。
今回は構造緩和を行うため vc-relax を利用しました。

xtl2pw.pw Si.xtl vc-relax 1


このままだと単位胞に原子を8個含む単純立方格子なので、ブラベ格子を fcc (ibrav=2)にして ATOMIC_POSITIONS もそれに倣うように編集し、これを走らせます。

mpirun -np 4 pw.x < Si.relax.in | tee Si.relax.out


&control
calculation='vc-relax' ,
restart_mode='from_scratch' ,
prefix='Si' ,
outdir = './Si/' ,
wfcdir = './Si/' ,
pseudo_dir = './' ,
disk_io='default' ,
forc_conv_thr= 0.001 ,
verbosity = 'default' ,
nstep = 100 ,
/
&system
ibrav= 2 ,
celldm(1) = 10.26221441
nat = 2 ,
ntyp = 1 ,
ecutwfc = 50.0 ,
ecutrho = 250.0 ,
/
&electrons
electron_maxstep = 100 ,
mixing_beta = 0.7 ,
! use smaller conv_thr for better results ,
conv_thr = 1.0d-12 ,
/
&ions
ion_dynamics='bfgs' ,
/
&CELL
cell_dynamics = 'bfgs' ,
press = 0.001,
press_conv_thr = 0.05 ,
! cell_dofree = 'xyz' ,
/
ATOMIC_SPECIES
Si 28.086 Si.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
Si 0.000000 0.000000 0.000000 0 0 0
Si 0.250000 0.250000 0.250000 0 0 0
K_POINTS automatic
4 4 4 0 0 0


nscf計算


セルフコンシステント計算の入力ファイルをコピーして nscf 計算の入力ファイルを作成します。

cp Si.relax.in Si.nscf.in


その後、以下の点を編集します。

  • calculation='bands'
  • nbnd=12
  • K_POINTS {crystal_b}


nbnd はバンドの数です。省略すると価電子数の半分×原子の個数(Siの場合はnbnd=4)となります。そうすると、半導体の場合は、価電子帯のみが計算されることになるので、少し大きめに取っておくほうがいいと思います。今回は nbnd=12 としました。

K_POINTS {crystal_b} の後の L 20 とかは、バンド構造をプロットする特徴的なk点のパスです。
最初の数字が特徴的なk点の数(だと思う)で、文字の後の数字が分割数(だと思う)です。LはL点、gGはΓ点で、この標記で使える点の名前は Doc/brillouin_zones.pdf を参照との事です。

mpirun -np 4 pw.x < Si.nscf.in | tee Si.nscf.out


なお nscf 計算用の入力ファイルの格子定数などは、構造緩和をする前の値のままですが、実際の nscf 計算では、ちゃんと構造緩和後の値を使ってくれます。

&control
calculation='bands' ,
restart_mode='from_scratch' ,
prefix='Si' ,
outdir = './Si/' ,
wfcdir = './Si/' ,
pseudo_dir = './' ,
disk_io='default' ,
forc_conv_thr= 0.001 ,
verbosity = 'default' ,
nstep = 100 ,
/
&system
ibrav= 2 ,
celldm(1) = 10.26221441
nat = 2 ,
ntyp = 1 ,
ecutwfc = 50.0 ,
ecutrho = 250.0 ,
nbnd = 12 ,
/
&electrons
electron_maxstep = 100 ,
mixing_beta = 0.7 ,
! use smaller conv_thr for better results ,
conv_thr = 1.0d-12 ,
/
&ions
ion_dynamics='bfgs' ,
/
&CELL
cell_dynamics = 'bfgs' ,
press = 0.001,
press_conv_thr = 0.05 ,
! cell_dofree = 'xyz' ,
/
ATOMIC_SPECIES
Si 28.086 Si.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
Si 0.000000 0.000000 0.000000 0 0 0
Si 0.250000 0.250000 0.250000 0 0 0
K_POINTS {crystal_b}
5
L 20
gG 20
X 20
W 20
gG 20


プロット用データの出力


このまででバンド計算は完了しています。次に gnuplot でプロットしやすい形式に直します。
以下のように Si.band.in を用意します。

&bands
outdir = './Si/' ,
prefix='Si' ,
filband='Si.band' ,
lsym=.true.
/


実行ファイルは bands.x です。

mpirun -np 4 bands.x < Si.band.in | tee Si.band.out


実行後に作成された Si.band.gnu が gnuplot 用の数値データです。

プロット


gnuplot でプロットするための plt ファイルとして以下のようなものを準備しました。

set terminal pngcairo size 520,390
set output "Si-band.png"

## *** Plot range ***
x1=0.8585
x2=1.8499
x3=2.3455
xmax=3.4539
set xrange [0:xmax]
set yrange [-15:10]

ef=5.9827

set xzeroaxis
set grid x

set ylabel "Energy (eV)"
set xtics ("{L}" 0, "{/Symbol G}" x1, "{X}" x2, "{W}" x3, "{/Symbol G}" xmax)
set x2tics ("{L}" 0, "{/Symbol G}" x1, "{X}" x2, "{W}" x3, "{/Symbol G}" xmax)

plot "Si.band.gnu" u 1:($2-ef) w lp not


特徴的なk点の座標は Si.band.out の中を見ると、以下のように書かれています。

     high-symmetry point:  0.4957 0.4957 0.4957   x coordinate   0.0000
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 0.8585
high-symmetry point: 0.0000 0.9913 0.0000 x coordinate 1.8499
high-symmetry point: 0.4957 0.9913 0.0000 x coordinate 2.3455
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 3.4539


エネルギーの単位は eV のようです。

Too many bands... と言われる場合


バンド計算をすると Too many bands are not converged from nscf calculation というメッセージが出て計算が止まることがあります。Google 検索をすると対策として下記の[Pw_forum] Too many bands are not converged from nscf calculationがヒットします。

Hi,Zhiting Tian
"Too many bands are not converged" can be solved for two way:
1.increase ecutwfc
2.decrease conv_thr
or both do them.

Best
Regards
--
Yun Song,Kang
Physical Science and Technology of Inner Mongolia University.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.democritos.it/pipermail/pw_forum/attachments/20110926/8f5f3309/attachment.htm


ここに書かれている通り ecutwf を大きくするか、conv_thr を小さくすると計算できるようになります。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: バンド構造 QuantumESPRESSO PWscf 分散関係 

PWscfでZnOの構造緩和

PWscf(Quantum ESPRESSO)を用いてウルツ鉱構造のZnOの高圧力下での構造最適化を行いました。

003_201807261443018d6.png

Fig.1: ウルツ鉱構造ZnOの格子定数の圧力変化



入力ファイルのテンプレート


まずPWscfの入力作成補助で紹介したxlt2pw.pyを利用して、ZnOの入力ファイルの雛形を作ります。入力の元になるxtlファイルは、CIFからVESTAを通じて出力します。CIFはDFT計算用データベース MatNaviで紹介した無機材料データベース(AtomWork)などからダウンロードします。

xtl2pw.py ZnO.xtl vc-relax 1
cp ZnO.in ZnO.relax.in


擬ポテンシャルとカットオフエネルギー


PWscfの擬ポテンシャルを参考に入力ファイルを編集します。
擬ポテンシャルの種類には PAW ポテンシャルを選びました。

ATOMIC_SPECIES
Zn 65.38 Zn.pbe-dnl-kjpaw_psl.1.0.0.UPF
O 15.999 O.pbe-n-kjpaw_psl.0.1.UPF


カットオフエネルギーは、擬ポテンシャルファイルに書かれている推奨値を参考にします。波動関数のカットオフの推奨値は酸素の方が大きく、電荷密度の推奨値は亜鉛の方が大きいことに注意が必要です。それぞれ大きいほうの1.5倍ぐらいを選びました。

    ecutwfc = 75.0 ,
ecutrho = 420.0 ,


圧力


圧力下の構造最適化を行うこともできます。その場合は press に圧力を設定します。単位は kbar です(1 GPa = 10 kbar)。

原子位置の制約


xtl2pw.py の最後の入力パラメータは、原子位置の制約の有無に関連しています。
0 を入力すると全ての原子について位置の最適化を行います。
1 を入力すると結晶学的に動かなさそうな原子を固定して計算します。

下記は xtl2pw.py の最後のパラメータに 1 をしていた場合の入力ファイルです。原子の座標の後ろに更に3つの数字が続いていますが、これが原子を動かすか否かのフラグになっていて 0 なら固定、1 なら動かせることを意味しています。今回の例では、酸素原子のz方向だけ原子を最適化するようになっていることが分かります。

ATOMIC_POSITIONS crystal
Zn 0.333333 0.666667 0.000000 0 0 0
Zn 0.666667 0.333333 0.500000 0 0 0
O 0.333333 0.666667 0.389010 0 0 1
O 0.666667 0.333333 0.889010 0 0 1


入力ファイル


結局、入力ファイルは以下のようにしました。

&control
calculation='vc-relax' ,
restart_mode='from_scratch' ,
prefix='ZnO' ,
outdir = './ZnO/' ,
wfcdir = './ZnO/' ,
pseudo_dir = './' ,
disk_io='default' ,
forc_conv_thr= 0.001 ,
verbosity = 'default' ,
nstep = 100 ,
/
&system
ibrav= 4 ,
celldm(1) = 6.14709011 ,
celldm(3) = 1.60158627686 ,
nat = 4 ,
ntyp = 2 ,
ecutwfc = 75.0 ,
ecutrho = 420.0 ,
/
&electrons
electron_maxstep = 100 ,
mixing_beta = 0.7 ,
! use smaller conv_thr for better results ,
conv_thr = 1.0d-12 ,
/
&ions
ion_dynamics='bfgs' ,
/
&CELL
cell_dynamics = 'bfgs' ,
press = 0.001,
press_conv_thr = 0.05 ,
! cell_dofree = 'xyz' ,
/
ATOMIC_SPECIES
Zn 65.38 Zn.pbe-dnl-kjpaw_psl.1.0.0.UPF
O 15.999 O.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
Zn 0.333333 0.666667 0.000000 0 0 0
Zn 0.666667 0.333333 0.500000 0 0 0
O 0.333333 0.666667 0.389010 0 0 1
O 0.666667 0.333333 0.889010 0 0 1
K_POINTS automatic
6 6 4 0 0 0


結果


圧力を変化させながら格子定数がどのように変化するかを調べた結果が Fig.1 です。40 GPaまでは正常に圧縮されていっていますが 50 GPa で異常が見られます。
pwout2xtl.pyをもちいて出力したxtlファイルをVESTAで描画してみると、酸素原子のz位置が大きく動いて、結晶全体がc軸方向につぶれたことが良く分かります。

001_20180726144258855.jpg
002_20180726144259aba.jpg

Fig.2-3: ZnOの結晶構造の変化


関連エントリ




参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: PWscf QuantumESPRESSO 最適化 ZnO 

PWscfの入力作成補助

Quantum ESPRESSO(PWscf)の入力ファイル作成支援に使えるソフトウエアとして、以下の3つがあげられます。


cif2cell


cif2cellCIFからecalj入力の作成
でCIFファイルからecaljの入力ファイルを作成するために利用しました。インストール方法はそちらのエントリに書いてあります。

cif2cellがインストールしてあれば、格子定数や原子位置が記入された入力ファイルの雛形を作成することが出来ます。
フォーマットは以下の通りです。
cif2cell -p プログラム -f CIFファイル -o PWscf入力ファイルの雛形


例えばシリコンのCIFファイル Si.cif の場合は、以下のようにします。

cif2cell -p pwscf -f Si.cif -o Si.in


他にも cif2cell -h とタイプすることによってヘルプを表示することが出来ます。

とはいえ、擬ポテンシャルの設定などの自明でないパラメータは一切出力されないため、追加の編集が必要になります。私は、追加の編集をPWguiで行いたいと考えたのですが、非自明なパラメータが一切出力されないため、そのまま PWgui で開くことが出来ないようです。

xtl2pw.py


他の入力支援用のPythonスクリプトがQuantum-ESPRESSOとVestaの連携QuantumESPRESSO_空間群入力で公開されています。

ダウンロードした後、適当にPATHの通ったディレクトリへ置いて実行権限を付けます。

wget http://www.misasa.okayama-u.ac.jp/~masami/xtl2pw.py
mv xtl2pw.py ~/bin/
chmod +x ~/bin/xtl2pw.py


wget http://www.misasa.okayama-u.ac.jp/~masami/pwout2xtl.py
mv pwout2xtl.py ~/bin/
chmod +x ~/bin/pwout2xtl.py


wget http://www.misasa.okayama-u.ac.jp/~masami/vesta2pw.py
mv vesta2pw.py ~/bin/
chmod +x ~/bin/vesta2pw.py


xtl2pw.py は cif2cell と異なり、色々なパラメータをとりあえず埋めてくれるため、直接 PWgui で開くことが出来ます。ただし変換するのは cif ファイルからではなく xtl ファイルからなので、 VESTA で xtl ファイルをあらかじめ出力しておく必要があります。
具体的には、以下の手順になると思います(シリコンの例)。

  1. VESTA で cif ファイルを xtl ファイルへ変換する。
    File → Export data → Fractional Coordinate (*.xtl) のファイル形式で保存
  2. xtl2pw.py Si.xtl scf 1
    2つ目の引数は scf, relax, vc-relax が選べる
    3つ目の引数は原子位置の制約
     0:制約を点けない
     1:制約をつける
  3. 出力された入力ファイル Si.in をpwgui で編集


Ubuntu への PWgui のインストール


PWguiのページから Self-contained standalone executables つまりコンパイル済みのバイナリをダウンロードします。
そのまえに iwidgets4 をインストールしておく必要があるかもしれません。

sudo apt-get update
sudo apt-get install iwidgets4
wget http://www-k3.ijs.si/kokalj/pwgui/download/pwgui-6.1-linux-x86_64.tgz
tar xzvf pwgui-6.1-linux-x86_64.tgz


必要に応じて pwgui を PATH に追加します。

Windows への PWgui のインストール


Linuxサーバー上で動作する PWgui の GUI を転送して Windows クライアントで動作させるのは重いので、Windows 上で PWgui が動作するとありがたいのですが、ネットで探してみても最新版の PWgui に関しては失敗報告しか見かけません。古いバージョンの2.1.3ならWindows版がPWgui の旧バージョンのダウンロードページから入手できますが、最新版は無いようです。

恐ろしく回りくどい方法ですが、現状では VirtualBox 上で Ubuntu を動かして、そこで PWgui を動かすのが一番マシかもしれません。

関連エントリ




参考URL




参考文献/使用機器





フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Quantum_ESPRESSO PWscf VESTA PWgui cif2cell 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj状態密度モンテカルロ解析常微分方程式odeトランジスタインターフェースDOSPDS5022定電流スイッチング回路確率論半導体分散関係シェルスクリプト乱数レベルシフトHP6632Aトランジスタ技術温度解析可変抵抗I2CブレッドボードR6452A反強磁性数値積分バンド構造バンドギャップセミナー絶縁偏微分方程式非線形方程式ソルバPWscf熱設計シュミットトリガLED三端子レギュレータ順列・組み合わせLM358GW近似カオスマフィンティン半径ISO-I2CフォトカプラA/Dコンバータ発振回路74HC4053数値微分直流動作点解析サーボPC817CアナログスイッチUSB補間TL431カレントミラーbzqltyVESTA電子負荷イジング模型LDA開発環境ブラべ格子FFT量子力学2ちゃんねるチョッパアンプ単振り子ポケモンGOスーパーリーグ標準ロジックQuantumESPRESSO基本並進ベクトルパラメトリック解析アセンブラBSchトレーナーバトル抵抗Maximaラプラス方程式失敗談状態方程式SMPキュリー温度スイッチト・キャパシタ位相図繰り返し熱伝導gfortranコバルトewidthTLP621不規則合金ランダムウォーク六方最密充填構造FET最適化相対論スピン軌道相互作用QSGWQuantum_ESPRESSOGGAVCA仮想結晶近似スレーターポーリング曲線cygwinZnOシュレディンガー方程式フォノンNE555詰め回路条件分岐固有値問題最大値ダイヤモンドガイガー管TLP552マントル自動計測データロガーQNAPUPSCIF井戸型ポテンシャルMCUxcrysdenゼーベック係数格子比熱最小値LM555フェルミ面fsolve過渡解析差し込みグラフ三角波起電力スーパーセル第一原理計算ブラウン運動FXA-7020ZROpenMPTLP521Ubuntuハーフメタル熱力学Writer509ubuntu平均場近似テスタawkLMC662フィルタMAS830LCK1026トランスPIC16F785AACircuit負帰還安定性ハイパーリーグCapSenseナイキスト線図ノコギリ波2SC1815EAGLEPvPP-10OPA2277MBEPGA入出力固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型磁気モーメント電荷密度重積分不純物問題擬ポテンシャル状態図cif2cellPWguiSIC二相共存リジッドバンド模型edeltquantumESPRESSOスワップ領域ルチル構造ウルツ鉱構造BaO岩塩構造ヒストグラム確率論マテリアルデザインフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデル文字列不規則局所モーメント陰解法熱拡散方程式HiLAPWCrank-Nicolson法連立一次方程式specx.fifort境界条件両対数グラフ片対数グラフGimp円周率ヒストグラムシンボル線種グラフの分割軸ラベル凡例トラックボール

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ