ねがてぃぶろぐ

構造緩和計算

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

最初にMg₂Siの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

関連エントリ

• UbuntuにBoltzTraPをインストール
• PWscfの入力作成補助

参考URL

• BoltzTraP
• Levi Lentz's blog
• BoltzTraP
• Materials project
• WinmostarTM チュートリアル BoltzTraP Quantum ESPRESSOとの連携
• 計算科学基礎講座：熱電研究のための第一原理計算入門 第3回 Boltzmann 輸送方程式を用いた熱電特性計算

参考文献/使用機器

フィードバック

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



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

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

以下のようにシリコンの入力ファイルを作ります。
何らかの作成補助を使うとラクかもしれません(参考: 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 を小さくすると計算できるようになります。

関連エントリ

• PWscfのインストール
• PWscfの入力作成補助
• PWscfでZnOの構造緩和

参考URL

• 3. グラフェンのバンド構造と状態密度
• [Pw_forum] Too many bands are not converged from nscf calculation

参考文献/使用機器

フィードバック

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



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

BoltzTraPのインストール

基本的にはLevi Lentz's blogに書かれている通りに従います。

This is not a guide on how to install BoltzTraP, as that documentation is actually pretty straight forward. You can find the software here: BoltzTrap Website. The main thing I will mention is that it will generate two executables: x_trans and BoltzTraP. You can leave x_trans anywhere in your path, but it appears to need BoltzTraP to be placed in /usr/local/bin. I could not tell you why this is, but it was the only way I could get it installed.

まずBolzTraPのウエブページへ行き、最新版のv1.2.5(2019年1月1日現在)をダウンロードします。

次に下記のコマンドで展開します。これには少し時間がかかります。

tar xf BoltzTraP.tar.bz2

boltztrap-1.2.5/src/ に存在する BoltzTraP と x_trans を /usr/local/bin/ にコピーします。

sudo cp boltztrap-1.2.5/src/BoltzTraP /usr/local/bin/
sudo cp boltztrap-1.2.5/src/BoltzTraP /usr/local/bin/

更に qe2bolz.py をダウンロードして実行権限を付けます。

sudo wget -O /usr/local/bin/qe2boltz.py http://levilentz.com/BoltzTraP/qe2boltz.py
sudo chmod +x /usr/local/bin/qe2boltz.py

これでインストールは完了のはずですが Ubuntu 16 以降(?)ではエラーが出るかもしれません。

BoltzTraPのコンパイル

Levi Lentz's blogを参考にチュートリアルを実行すると liblapack.so.3gf が無いというエラーが出て失敗するかもしれません。
これはecaljのインストール(Ubuntu + gfortran)の2018年4月7日追記のところに書いたのと同じ原因なので、同じ対処をすれば解決できます。自分でコンパイルをするので gfortran などがインストールされていることが前提です。

具体的には、まず src/ にある Makefile を下記のように編集します。

#LIBS = -llapack -lblas
LIBS = /usr/lib/liblapack.so.3 /usr/lib/libblas.so.3

これを make した後 /usr/local/bin/ へコピーすれば解決です。

make
sudo cp BoltzTraP /usr/local/bin/

関連エントリ

• ecaljのインストール(Ubuntu + gfortran)

参考URL

• BolzTraP
• Levi Lentz's blog

フィードバック

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



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