スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。


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

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。