はじめてのVASP (fcc Si)

とりあえずVASPのサンプルやってみる。
http://cms.mpi.univie.ac.at/wiki/index.php/VASP_workshop_2003
http://cms.mpi.univie.ac.at/wiki/index.php/VASP_example_calculations
より、fcc Si。

やってることは、格子定数を変えてエネルギーを計算する(エネルギーの極小値が安定な格子定数を与える)。

INCAR

System = fcc Si # names of System
ISTART = 0 # startjob: 0-new 1-cont 2-samecut
ICHARG = 2 # charge: 1-file 2-atom 10-const
ENCUT = 240 # energy cutoff in eV
ISMEAR = 0 # part. occupancies: -5-Blochl -4-tet -1-fermi 0-gaus 0 MP
SIGMA = 0.1 # broadening in eV -4-tet -1-fermi 0-gaus
ISTARTWAVECAR(波動関数のパラメータ)を読み込んで使用するかどうか。初期設定で計算する場合は0、ファイルを読み込む場合は1。設定しない場合、WAVECARがあれば1、なければ0になる。WAVECARが何かはまだ詳しく見ていないが、波動関数に関する情報が記述されているようだ。
ICHARG最初の電荷密度をどうするか。初期軌道(? initial orbitalsが何を意味するのかわからない。孤立した状態での電荷密度か?)を使う場合は0。1: CHGCARを読み込む。2: スーパーポジション(?)の電荷密度を使う。4: VASP5.1以上で使用。POTファイルからポテンシャルを読み込む(よくわからない。)。10: セルフコンシステントループを使わない(何のための設定いかよく理解出来ていない)。デフォルト値はISTART=0の時に2、それ以外は0。
ENCUTカットオフエネルギーを設定。デフォルト値はPOTCARの最大値を使用。
ISMEAR(Ion Smearing)よく理解していないが、電子のつまり方をどの関数を使って引き延ばす(Smearing)かということらしい(収束性向上のため)。-1: Fermiモデル。0: ガウシアンモデル。1…N: Methfessel-Paxton (MP)モデル。-2: WAVECARを読み込む。-3: INCARで設定された値(SMEARINGS)を読み込む。-4: tetrahedron method without Blochl corrections。-5: tetrahedron method with Blochl corrections。デフォルトは1。total energyとDOS の計算には-5がいい。金属の力や応力テンソルの計算には5~10%の誤差が生じる。フォノン計算には0がおすすめ。MPモデルはSIGMA (the width of smearing)の値を適切に選べばtotal energyも精度よく求められる(マニュアル7.4参照)。smearing-parameter (SIGMAのことか?)が大きすぎるとtotal energyには良くない。小さすぎるとk点メッシュを増やさなければならない。free energyとtotal energyの差(entropy T*S)で出来る限り大きな値を設定すべき。ガウシアンモデルは大抵の場合良い結果をもたらす。
SIGMA(the width of smearing)デフォルト値0.2。ISMEARに-4, -1, 0を設定した時に使用される?

KPOINTS

k-points
0
Monkhorst Pack
11 11 11
0 0 0

1行目: コメント文
2行目: k点の数。0の場合自動生成。
3行目: “G” or “g”で始まっていたらGamma centered grid、”M” or “m”で始まっていたらMonkHorst-Pack scheme。それ以外は素人にはお勧めできないようだ。
4行目: 分割数。
5行目: optional shift of the mesh。普通は0。

POSCAR

fcc Si:
3.9
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5
1
cartesian
0 0 0

1行目: コメント文
2行目: 格子定数
3-5行目: 格子ベクトル
6行目: 各原子数。順番はPOTCAR、INCAR内の記載順と合わせる必要あり。
7行目: 原子座標の指定方法
8行目- : 原子座標

POTCAR

使いたいポテンシャルをコピペ。PAW_PBEを使用。

loop.sh

格子定数を3.5~4.3まで変えて計算してくれるバッチファイル。それぞれの計算の出力ファイルを残すためには書き換える必要あり。もちろんVASPのパスも。

計算結果

計算結果のうちOSZICARから抜き出したtotal energy他がSUMMARY.fccにまとめられている。gnuplotやexcelなどでプロットすれば下記のグラフが得られる。

この結果から格子定数が3.9Aあたりが一番安定なことが分かる。

構造緩和

ちなみに、INCARを

System = fcc Si # names of System
ISTART = 0 # startjob: 0-new 1-cont 2-samecut
ICHARG = 2 # charge: 1-file 2-atom 10-const
ENCUT = 240 # energy cutoff in eV
ISMEAR = 0 # part. occupancies: -5-Blochl -4-tet -1-fermi 0-gaus 0 MP
SIGMA = 0.1 # broadening in eV -4-tet -1-fermi 0-gaus
IBRION = 1 # relaxation method
NSW = 100 # maximum steps of ion
NELMIN = 4 # minimum of self consistent loop
ISIF = 3 # calculation mode for force etc.

のように記載(IBRION以下)すると、構造緩和を行って最も安定な配置とエネルギーを求めてくれる。

IBRION: 構造緩和する方法。デフォルト-1。
NSW: 構造緩和する最大ステップ数。デフォルト0。EDIFF(デフォルト10-4)以下になると終了。
NELMIN: セルフコンシステントループの最小ループ数。IBRION=1(準ニュートン法)の場合、4~8の間の値が良いらしい。
ISIF: 応力テンソル計算や自由度の設定。(詳細はマニュアル参照)

計算結果(OUTCAR)を見ると、格子定数が3.86 Aの時にEtot=-4.878 eVとなる。

k点メッシュ

k点メッシュの数を変えて、計算の収束を見てみる。
loop.shを下記のように書き換えて実行。

! /bin/bash
BIN=vasp_std
rm WAVECAR SUMMARY.fcc
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ; do
cat KPOINTS !
K-Points
0
Monkhorst Pack
$i $i $i
0 0 0
!
echo "a= $i" ; $BIN
E=awk '/F=/ {print $0}' OSZICAR ; echo $i $E SUMMARY.kpoints
done
cat SUMMARY.kpoints

横軸k点メッシュ数、縦軸エネルギーでプロットしたのが下図。

メッシュ数4以上で収束してる感じか。

コメントする

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です