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で行われた計算例が公開されています。

basis_spin_01.png
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)のあたりで最もエネルギーが低く安定であることが再現できています。


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


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

実測PHASEAkaiKKR
a (Å)2.872.8452.81
B (GPa)168177.237190
table.1: 強磁性bcc鉄の格子定数と体積弾性率


実測と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にはα,β,γと書いてありますが、これは単なるα,γ,δの誤植だと思われます。

関連エントリ




参考URL




付録


このエントリで使用したAkaiKKR用入力ファイルなどを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器





フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。


tag: AkaiKKR machikaneyama KKR 安定相 状態方程式 

comment

Secret

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj常微分方程式モンテカルロ解析状態密度トランジスタodeDOSインターフェース定電流スイッチング回路PDS5022半導体シェルスクリプト分散関係レベルシフト乱数HP6632AR6452A可変抵抗トランジスタ技術温度解析ブレッドボードI2C反強磁性確率論数値積分セミナーバンドギャップバンド構造偏微分方程式非線形方程式ソルバ熱設計絶縁ISO-I2Cカオス三端子レギュレータLM358GW近似マフィンティン半径A/DコンバータフォトカプラシュミットトリガLEDPC817C発振回路数値微分直流動作点解析サーボカレントミラーTL431アナログスイッチUSB74HC4053bzqltyVESTA補間電子負荷アセンブライジング模型BSch量子力学単振り子2ちゃんねるチョッパアンプLDA開発環境基本並進ベクトルFFT標準ロジックブラべ格子パラメトリック解析抵抗SMPMaxima失敗談ラプラス方程式繰り返し位相図スイッチト・キャパシタ熱伝導状態方程式キュリー温度gfortranコバルトTLP621不規則合金Quantum_ESPRESSO六方最密充填構造ランダムウォーク相対論ewidthスピン軌道相互作用FETQSGWVCAcygwinスレーターポーリング曲線GGA仮想結晶近似PWscfシュレディンガー方程式LM555ハーフメタル固有値問題NE555最小値ガイガー管QNAPUPS自動計測ダイヤモンドマントルTLP552格子比熱最適化MCU井戸型ポテンシャル最大値xcrysdenCIF条件分岐詰め回路フェルミ面差し込みグラフスーパーセルfsolveブラウン運動awk過渡解析起電力三角波第一原理計算FXA-7020ZRWriter509Ubuntuテスタ熱力学データロガーTLP521OpenMPubuntu平均場近似MAS830LトランスCK1026PIC16F785PGA2SC1815EAGLEノコギリ波負帰還安定性ナイキスト線図MBEOPA2277P-10フィルタCapSenseAACircuitLMC662文字列固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型電荷密度重積分SIC二相共存磁気モーメント不純物問題PWgui擬ポテンシャルゼーベック係数ZnOウルツ鉱構造edeltquantumESPRESSOフォノンリジッドバンド模型スワップ領域BaO岩塩構造ルチル構造ヒストグラム確率論マテリアルデザインフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデルcif2cell入出力陰解法熱拡散方程式HiLAPW両対数グラフCrank-Nicolson法連立一次方程式specx.fifort境界条件片対数グラフグラフの分割円周率ヒストグラム不規則局所モーメントGimpシンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ