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強磁性OPアンプPICCPAecaljモンテカルロ解析常微分方程式odeトランジスタ状態密度DOSインターフェース定電流PDS5022スイッチング回路半導体シェルスクリプト乱数レベルシフト分散関係HP6632AI2C可変抵抗トランジスタ技術ブレッドボード温度解析R6452A反強磁性確率論バンドギャップセミナー数値積分熱設計非線形方程式ソルババンド構造絶縁偏微分方程式ISO-I2CLM358マフィンティン半径フォトカプラシュミットトリガカオスLED三端子レギュレータGW近似A/Dコンバータ発振回路PC817C直流動作点解析USBTL431数値微分アナログスイッチカレントミラー74HC4053サーボ量子力学単振り子チョッパアンプ補間2ちゃんねる開発環境bzqltyFFT電子負荷LDAイジング模型BSch基本並進ベクトルブラべ格子パラメトリック解析標準ロジックアセンブラ繰り返し六方最密充填構造SMPコバルトewidthFET仮想結晶近似QSGW不規則合金VCAMaximaGGA熱伝導cygwinスレーターポーリング曲線キュリー温度スイッチト・キャパシタ失敗談ランダムウォークgfortran抵抗相対論位相図スピン軌道相互作用VESTA状態方程式TLP621ラプラス方程式TLP552条件分岐NE555LM555TLP521マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測QNAPUPSWriter509ダイヤモンドデータロガー格子比熱熱力学起電力awkブラウン運動スーパーセルUbuntu差し込みグラフ第一原理計算フェルミ面fsolveCIFxcrysden最大値最小値ubuntu最適化平均場近似OpenMP井戸型ポテンシャル固有値問題シュレディンガー方程式TeX2SC1815結晶磁気異方性OPA2277フラクタルFSM固定スピンモーメントc/a非線型方程式ソルバgnuplot全エネルギーfcc初期値マンデルブロ集合縮退正規分布interp1ウィグナーザイツ胞L10構造multiplotフィルタ面心立方構造PGAハーフメタル二相共存ZnOウルツ鉱構造BaOSIC重積分磁気モーメント電荷密度化学反応クーロン散乱岩塩構造CapSenseノコギリ波デバイ模型キーボード半金属フォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金Realforce軸ラベルグラフの分割凡例線種シンボルMAS830LCK1026LMC662PIC16F785トランス関数フィッティングトラックボールPC等価回路モデルヒストグラムパラメータ・モデル不規則局所モーメント最小二乗法TS-112TS-110直流解析ExcelGimp円周率片対数グラフ両対数グラフspecx.f疎行列三次元ifort文字列不純物問題P-10等高線ジバニャン方程式ヒストグラム確率論マテリアルデザイン入出力境界条件陰解法AACircuit熱拡散方程式HiLAPWMBEEAGLE連立一次方程式ナイキスト線図負帰還安定性Crank-Nicolson法日本語

最新コメント
リンク

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