(cache) ねがてぃぶろぐ　状態方程式

AkaiKKRで鉄の安定相と格子定数

AkaiKKRでダイヤモンド型構造半導体やAkaiKKRでニッケル・鉄・コバルトでは、色々な結晶構造を持つ固体の電子構造の計算を行いましたが、その格子定数はあらかじめ分かっているものとして入力ファイルを作成しました。

今回は、全エネルギー最小化という条件から、結晶構造や格子定数を第一原理的に決定することを目的にα, γ, δ-Feの計算をAkaiKKRで行い、計算結果をMurnaghanの状態方程式へフィッティングしました。

その結果α-Feが最も安定で、その格子定数は a=2.81(Å) となる事がわかりました。これは、実測値である a=2.87 (Å) と極めて近い値です。

第一原理計算と格子定数

『第一原理計算とは』といった説明は、調べてみると色々な方が説明されていますが、私にとって直感的に分かりやすかったのは第一原理計算（基礎編）の下記のものです。

第一原理(英語はfirst-principlesまたはab initio(英語ではない))とは、なんら実験データや経験パラメーターを使わないで理論計算をする方法の総称。でも、この業界では電子状態計算、平たく言うとシュレディンガー方程式を解く計算のことを暗に指すことが多い。

第一原理計算というために大切なことは、実験データを見てパラメータをいじるということをしていないということです。

とは言うものの実際に第一原理計算パッケージを使うときには、これまでAkaiKKRを使ったエントリにも当てはまることですが、結晶構造やその格子定数は計算の前に入力する必要があります。

例えば鉄の場合について考えます。

鉄は常温では体心立方(body-centered cubic; bcc)構造となり(フェライト; α-Fe)強磁性を持ちます。温度を上げていくと911℃で面心立方(face-centered cubic; fcc)構造に相転移します(オーステナイト; γ-Fe)。さらに温度を上げていくと1392℃で再びbcc構造へ戻りますが(デルタフェライト; δ-Fe)、このときは既にキュリー温度を超えているため常磁性です。

AkaiKKRの入力ファイルを作成するとき、常温(というよりは絶対零度)でα, γ, δのどの鉄が安定なのか、また、その格子定数がいくつなのかは実験データを見ないと分からないことになります。

第一原理計算が実験データを使わないとは言っても、入力する結晶構造と格子定数だけは多くの場合大目に見てくれるというのが暗黙の了解の様ですが、厳しい見方をするなら第一原理でないということも出来ます。

この問題は、いろいろな結晶構造でいろいろな格子定数を使って計算をした後、全エネルギーを比較することによって解決することが出来ます。
実例として、アドバンスソフト社が発売しているAdvance/PHASEで行われた計算例が公開されています。

Fig.1: Advance/PHASEによるFeの全エネルギーと体積の関係

今回のエントリではこの計算をAkaiKKRで行います。

specx.fの設定

今回の計算では、specx.fを下記の設定にしてmakeしました。

     & (natmmx=4, ncmpmx=4, msizmx=198, mxlmx=3, nk1x=2200, nk3x=2688,
     & msex=201, ngmx=15, nrpmx=650, ngpmx=650, npmx=350, msr=400)

入力ファイル

磁性を考慮した(magtyp=mag)bcc鉄(alpha_in.txt)、及び考慮しない(magtyp=nmag)bcc鉄(delta_in.txt)とfcc鉄(gamma_in.txt)の3つの入力ファイルを用意し、それぞれ1原子あたりの体積が8-15(Å3)となるようにしました。

Murnaghanの状態方程式へのフィッティング

計算結果の体積とエネルギーの対応関係(Energy_dat.txt)から、最もエネルギーが低いときの体積を求めるため、gnuplotを用いてMurnaghanの状態方程式へ最小二乗法フィッティングを行います。

Murnaghanの状態方程式は下記の式で表されます。

$E_{tot}(V) = \frac{BV}{B^{'}(B^{'}-1)}\left[B^{'}\left(1-\frac{V_0}{V}\right)+\left(\frac{V_0}{V}\right)^{B^{'}}-1\right]+E_0$

ここでEtotが全エネルギー、Bが体積弾性率、B'が体積弾性率の圧力微分、Vが体積です。
フィッティングとグラフの描画を行うgnuplotのスクリプトはEnergy_plt.txtです。

結果

結果をFig.2に示します。Advance/PHASEで行われた計算結果と比較すると(定量的なものはともかく)強磁性bcc構造のα-Feが1原子辺りの体積が11(Å3)のあたりで最もエネルギーが低く安定であることが再現できています。

Fig.2: AkaiKKRによるFeの全エネルギーと体積の関係

実際に測定された鉄の格子定数、Advance/PHASEによる計算結果と今回AkaiKKRで計算したものをTable.1にまとめました。

table.1: 強磁性bcc鉄の格子定数と体積弾性率
	実測	PHASE	AkaiKKR
a (Å)	2.87	2.845	2.81
B (GPa)	168	177.237	190

実測とAdvance/PHASEの値は第一原理シミュレータ入門 PHASE&CIAOから引用しました。

ひょっとしたらたまたまなのかも知れませんが、格子定数に関しては非常によく一致しています。

Appendix: β-Fe？

今回のエントリでは「強磁性bcc鉄=α-Fe」「常磁性bcc鉄=δ-Fe」という様な書き方をしてしまいましたが、これは厳密には間違いです。鉄のキュリー温度は770℃なので、常温付近では常磁性だったbcc鉄はfcc(γ-Fe)へ相転移する911℃よりも低温で常磁性になります。
この770-911℃の間で安定な常磁性bcc鉄のことを以前はβ-Feと呼んでいたようですが、現在ではα-Feの一部として扱われています。

強磁性bcc鉄(α-Fe)
　｜
　｜ 770℃
　↓
常磁性bcc鉄(これもα-Fe,以前はβ-Feと呼ばれていた。)
　｜
　｜ 911℃
　↓
fcc鉄(γ-Fe)
　｜
　｜ 1392℃
　↓
常磁性bcc鉄(δ-Fe)

AkaiKKRでは絶対零度の計算しか出来ないため『もしも絶対零度のとき常磁性bcc構造が安定だったとしたらどうなるか』という計算を(本来高温の相である)δ-Feと呼ぶのはおかしいのですが、表記をややこしくしないためこう書きました。

なお、第一原理シミュレータ入門 PHASE&CIAOの表3.4にはα,β,γと書いてありますが、これは単なるα,γ,δの誤植だと思われます。