VASPで相変態の計算(beta-tin Si)

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のようにバッチファイルを作って格子パラメータを変えた時のエネルギーも計算した。

結果

diamond Siも一緒に載せた結果が以下。(EもVも2 (atoms)で割ってある。)

Fig. 1 E-V curve of beta-tin Si and diamond Si

横軸が格子パラメータで縦軸がEtot (eV/cell)。

横軸を体積Vに取るE-V曲線というのがよく使われるようだ。なので、いまさらながら横軸をVにしてグラフにしてみる。エネルギーと体積は、”vasprun.xml”もしくは”OUTCAR”から抽出した。

熱力学第一法則より内部エネルギー変化dEは、系に取り込まれた熱\delta Qから系のした仕事\delta Wをひいて、

(1)   \begin{equation*} dE = \delta Q - \delta W \end{equation*}

これは可逆な断熱過程では\delta Q = 0および\delta W = PdVであるから、

(2)   \begin{equation*} dE = -PdV \end{equation*}

よって圧力Pは、

(3)   \begin{equation*} P = -\frac{dE}{dV} \end{equation*}

となる。これはE-V曲線の傾きが圧力であることを示している。つまり、2本のE-V曲線があるとき、その共通接線の傾きは相転移時の圧力と等しい。共通接線の単位はeV/Åである。これをGPaに変換するには、

(4)   \begin{equation*} 1 [GPa] = 0.00624150948 [eV/\AA] \end{equation*}

を用いる。

ということなので、共通接線をひいて傾きを求めてみた。とりあえず、各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*}

で表される。この式を用いてフィッティングを行うとびっくりするくらい合う(下図)。

コメントする

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