VASPのサンプルやってみたシリーズ。今回はbeta-tin Si。
Si has the diamond structure under air pressure. However, it has beta-tin structure under more than about 100 kbar. Can DFT support this phenomenon?
Siは大気圧ではダイアモンド構造をとるが、100kbar付近の圧力でβスズ構造に変化する。これを第一原理で表せるのか？という話。
特に入力ファイルについて追加で説明することはない。緩和計算の後、fcc Siのようにバッチファイルを作って格子パラメータを変えた時のエネルギーも計算した。
Contents
結果
diamond Siも一緒に載せた結果が以下。（EもVも2 (atoms)で割ってある。）
横軸が格子パラメータで縦軸がEtot (eV/cell)。
横軸を体積Vに取るE-V曲線というのがよく使われるようだ。なので、いまさらながら横軸をVにしてグラフにしてみる。エネルギーと体積は、”vasprun.xml”もしくは”OUTCAR”から抽出した。
熱力学第一法則より内部エネルギー変化
は、系に取り込まれた熱
から系のした仕事
をひいて、
(1) 
これは可逆な断熱過程では
および
であるから、
(2) 
よって圧力
は、
(3) 
となる。これはE-V曲線の傾きが圧力であることを示している。つまり、2本のE-V曲線があるとき、その共通接線の傾きは相転移時の圧力と等しい。共通接線の単位はeV/Åである。これをGPaに変換するには、
(4) ![\begin{equation*} 1 [GPa] = 0.00624150948 [eV/\AA] \end{equation*}](https://blog.kiyohiroyabuuchi.com/wp-content/ql-cache/quicklatex.com-e4f8430913cb46920b4b653d1b43351c_l3.png "Rendered by QuickLaTeX.com")
を用いる。
ということなので、共通接線をひいて傾きを求めてみた。とりあえず、各EV曲線を2次方程式で近似して（Excel等で近似曲線が求められる。）共通接線を求め（Maximaで連立方程式を解いた。他にもっと楽なやり方があったら教えて。）、傾きから上記の換算式で圧力に変換すると、およそ10 GPa。ものの本によると100 kbar（10 GPa）あたりで相変態が起こるらしいので、ほぼ一致している。（すごい。）
Birch-Murnaghanの状態方程式
上で2次式で近似して求めたが、これは2次のテイラー展開が有効な範囲、つまり平衡位置の周りのごく狭い範囲でのみ正しいといえる。実際に平衡位置を中心に2次方程式で近似したグラフを書くとほとんど合わない（下図）。
エネルギーと体積（格子定数）の関係を表す式として、等方性固体についてのBirch–Murnaghanの状態方程式というのがあり、
(5) ![\begin{equation*} E (V) = E_0 + \frac{9 V_0 B_0}{16} \left\{ \left[ \left( \frac{V_0}{V} \right)^\frac{2}{3} - 1 \right]^3 B'_0 + \left[ \left( \frac{V_0}{V} \right)^\frac{2}{3} - 1 \right]^2 \left[ 6 - 4 \left( \frac{V_0}{V} \right)^\frac{2}{3} \right] \right\} \end{equation*}](https://blog.kiyohiroyabuuchi.com/wp-content/ql-cache/quicklatex.com-eeb62616bc0f638ebbe6d3f2c22e71a1_l3.png "Rendered by QuickLaTeX.com")
で表される。この式を用いてフィッティングを行うとびっくりするくらい合う（下図）。
