VASPのサンプル(cd Si)やってみた

VASPのサンプルやってみたシリーズ。今回はdiamond Si。

http://cms.mpi.univie.ac.at/wiki/index.php/VASP_workshop_2003
http://cms.mpi.univie.ac.at/wiki/index.php/Cd_Si

やってることはfcc Siと同じ。
新しいことは特にないので入力ファイルの説明もなしで。

結果

結果を早速。

横軸格子定数、縦軸Etotなのだが、気になるのはEtotの単位。eV/atomになってるんだろうか?それともeV/cellで2 atomの計算だから2で割らないといけないのだろうか?

それから、プログラムによって異なる単位にも気をつけること。ここで、単位とは2つの意味があり、
 1) いわゆるエネルギーの単位 (→換算)
 2) 原子(イオン)個数の規格化
を考えなければならない。よくマニュアルを確認しておきましょう。
 1)は、ちなみにWien2kの場合 Ry(リュードベリ)系、VASPの場合 eV系である。
 2)は、時々間違う人が居るので注意。ちょっと詳しく述べてみよう。
固体に関するエネルギー表記は、慣習的に 〇〇kJ/molと書かれているけど、これでは1molの定義づけが非常にあいまいである。たとえばNaClのばあい普通Na 1個とCl 1個を合わせて1molと定義するので、正確には 〇〇kJ/mol per NaClなのだろう。ところが、第一原理計算をするときに、Fm3m(岩塩型のF格子)で座標を入力したとなると、現実の計算はNa4Cl4の8原子で計算したことになる。つまり出力は〇〇kJ/mol per Na4Cl4となるので注意。この場合は値を4で割ってあげなければならない。ところが、これは計算するソフトで違ってくる。Wien2kの場合は、入力ファイル作成プログラムがFm3mのシンメトリーを考慮してくれるので、実際の計算はPrimitive Cellつまり単位格子1個NaClで計算してくれる(設定による)。一方、VASPの場合はそういった補助プログラムがないので8原子の結果を与える。ということなので、とっても気をつけてください。

参考:中山将伸のホームページより http://mmnakayama.jimdo.com/study/%E7%AC%AC%E4%B8%80%E5%8E%9F%E7%90%86%E8%A8%88%E7%AE%97%EF%BC%92/

というわけで、多分fcc Siと比較したいとかの場合は2で割る必要がありそうだ(上のグラフは2で割っている)。

構造緩和

ついでに構造緩和のサンプル (cd Si volume relaxation) もしておく。

fcc Siでも既に入力ファイルについては説明済みだが、2つほど新しいタグが出てきたので書きとどめておく。

INCAR

System = diamond Si
ISMEAR = 0
SIGMA  = 0.1
ENMAX  =  240
IBRION = 2 
ISIF   = 3
NSW    = 15
EDIFF  = 0.1E-04
EDIFFG = -0.01

EDIFF: 電子計算のセルフコンシステントループを終了する際のエネルギー差。イタレーション毎のエネルギー差がこの値以下になったら終了。デフォルト値は1.0E-4。

EDIFFG: 緩和計算を終了する際のエネルギー差。イタレーション毎のエネルギー差がこの値以下になったら終了。この値が0以下の場合、全てのforcesが|EDIFFG|以下になったら計算を終了する(理由は分からないが、”This is usually a more convinient setting.”らしい。)。デフォルト値はEDIFF * 10。

結果

結果は、格子定数が5.465Åの時に、E0 = -.10827169E+02であった。グラフにのせると以下。

fcc Siと比較すると、最も低いエネルギーが
fcc Si: -4.878 eV/atom
diamond Si: -5.414 eV/atom

なので、diamond Siの方が安定な結晶構造ということが分かる。ただ、これは0K (T = 0)、0気圧 (P = 0)の時の結果である。(室温、1気圧の時はどうしたら求まるのか?)

DOS

構造緩和後のデータでDOS計算をした結果は以下。

INCARはfcc Si DOSと同じ。(構造緩和で出力されたDOSCARは役に立たないので、正確なDOSを求めたい場合は、構造緩和した後のCONTCARをPOSCARにコピーして、staticな計算(ISTART=1; NSW=0)を再度する必要がある。)

横軸E-EF (eV)、縦軸DOSであるが、E-EF = 0付近にDOS = 0の部分が見られる。これがバンドギャップでfcc Siと違い半導体であることを示している。ただし、DFTはこのバンドギャップの幅を正確に計算できず過小評価する傾向にある。逆にHFは過大評価するらしい。
ちなみにバンド図は以下 (INCARはfcc Si bandstructureと同じ)。

コメントする

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