AkaiKKRでコバルトの格子定数

AkaiKKR(Machikaneyama)を用いて、全エネルギー最小化の条件からhcp構造をもつコバルトの格子定数(a, c/a)を計算しました。
この際、入力ファイルがたくさんになるため、シェルスクリプトを用いてテンプレートとなる入力ファイルの中の格子定数を順次置き換えるようにしました。
計算結果に対して、Scilabでスプライン補間を行い、カラーコンターでプロット、及び全エネルギーが最小になるときの平衡格子定数を計算しました。その結果は、文献値と良い一致を示しました。

000_20140323150909952.png

Fig.1: AkaiKKRによるコバルトの平衡格子定数。第一原理計算でエネルギーが最小となる格子定数(白丸)と文献値(赤クロス: 鮫島研究室電気電子材料講義ノート(PDF), 黄クロス: GitHub Co-Cobalt.cif)



平衡格子定数


AkaiKKRで鉄の安定相と格子定数では、AkaiKKR(Machikaneyama)を用いて、体心立方構造(bcc)や面心立方構造(fcc)の鉄の格子定数aを変化させる計算を行いました。
この計算から求められた全エネルギーが最も小さくなる格子定数が、ゼロ気圧(≒大気圧)におけるその物質の格子定数になります。

鉄はbcc構造、ニッケルはfcc構造を取りますが、これらの他にコバルトなどが取る六方最密充填構造(hcp)も重要な結晶構造です。しかしながら、前者二つが格子定数としてaだけを可変させれば計算できるのに対して、hcp構造ではaとcの二つの格子定数を持つので作らなければならない入力ファイルの数がたくさんになります。

前回の第24回CMDワークショップで、私はアドバンストコースにて、小口先生のHiLAPWの演習と赤井先生のAkaiKKR(Machikaneyama)の演習を受講させていただきました。
その際に、小口先生にシェルスクリプトを利用して入力ファイルの格子定数を置き換える方法を教えていただきました。今回はこのテクニックを応用して、AkaiKKRでhcp構造を持つコバルトの格子定数を計算してみます。
(なお、シェルスクリプトの入門書に関しては坂本文著続・たのしいUNIXを紹介していただきました。)

入力ファイル


入力ファイルのテンプレートはAkaiKKRでニッケル・鉄・コバルトをベースに以下の様にしました。(Co0.in)

 go  data/co_4.75_1.633
hcp 4.75, 1.633, , , , ,
0.001 1.0 sra gga91 mag 2nd
update 4 200 0.024
1
Co 1 0 0 2 27 100

2
0.33333a 0.66667b 0.25c Co
0.66667a 0.33333b 0.75c Co


シェルスクリプト


上記のテンプレートの中の格子定数をシェルスクリプトを使って順次置き換えながらAkaiKKRに第一原理計算を行わせます。
第一原理計算を実行した後のアウトプットファイルに色々な細かい情報が書き出されるわけですが、その中で特に重要な格子定数aと全エネルギーと磁気モーメントの3つは、ポテンシャルファイルと同じdata/フォルダの中にco_4.700_1.600.infoと言ったように拡張子に.infoが付くファイルとして保存されます。
そこで第一原理計算が終わったら、最後の計算結果を表示するようにしました。(なお、状態密度の計算dosでも.infoに追記がされるようですが、その値は読んではいけない値との事。)

#!/bin/csh -f
set path=($path ~/kkr/cpa2002v009c )

set A_LIST=( 4.700 4.705 4.710 4.715 4.720 4.725 4.730 4.735 4.740 4.745 4.750 4.755 4.760 4.765 4.770 4.775 4.780 4.785 4.790 )
set ETA_LIST=( 1.600 1.605 1.610 1.615 1.620 1.625 1.630 1.635 1.640 1.645 1.650 1.655 1.660 1.665 1.670 1.675 1.680 1.685 1.690 )

set A0=4.75
set ETA0=1.633

specx < in/Co.in
specx < in/Co-DOS.in >out/Co-DOS.out

foreach A ( ${A_LIST} )
foreach ETA ( ${ETA_LIST} )
echo " A,ETA= "${A}" "${ETA}
if ( ! -e data/co_${A}_${ETA}.info ) then
cp data/co data/co_${A}_${ETA}
sed '1,6s/'${A0}'/'${A}'/g' in/Co0.in | sed '1,6s/'${ETA0}'/'${ETA}'/g' > in/Co_${A}_${ETA}.in
specx < in/Co_${A}_${ETA}.in > out/Co_${A}_${ETA}.out
endif
tail -n 1 data/co_${A}_${ETA}.info
end
end


これらのファイルを置くためのディレクトリ構成は、以下の様にしました。
Co/─┬─in/─┬─Co.in
    │     ├─Co-DOS.in
    │     └─Co0.in
    ├─out/
    ├─data/
    └─Co-AE.sh


平衡格子定数の計算


第一原理計算が出来たら、その結果をScilabで処理して、全エネルギーが最小になるときの格子定数を決定します。
第一原理計算の出力を多少整形しScilabで読みやすいマトリクスの形にしたのがEnergy.txtです。

このデータを読み込んで、カラーコンターとしてプロットするScilabスクリプトが下記のenergy.sceで、その計算結果が冒頭のFig.1です。

clear;
format('v',16)
bohr = 0.5291772108 // 1 bohr = 0.529 A

// *** データの読み込み ***
e = fscanfMat("Energy.txt")';
// 第一原理計算の格子定数
a = [4.70:0.005:4.79];
c = [1.60:0.005:1.69];

// *** スプライン補間 ***
E = splin2d(a,c,e);
aa = [4.70:0.001:4.79];
cc = [1.60:0.001:1.69];
[AA,CC] = ndgrid(aa,cc);
ee = interp2d(AA, CC, a, c, E);

// *** カラーコンターのプロット ***
// グラフを正方形に
isoview(4.70,4.79,1.60,1.69);
// 色の設定
xset("colormap",jetcolormap(64))
// カラコンターのプロット
Sgrayplot(aa, cc, ee - min(ee));
// カラーバーの表示
colorbar(0, max(ee) - min(ee));
// グラフの軸ラベル
xlabel("a (bohr)");
ylabel("c/a");

// *** 平衡格子定数の計算 ***
// 補間データから最小値を探す
[emin,ind] = min(ee);
// 平衡格子定数のプロット(白丸)
plot(aa(ind(1)),cc(ind(2)),'ow');
// 文献値(実験データ)のプロット
// http://www.tuat.ac.jp/~sameken/lecture/2012raremetal.pdf
plot(2.514 / bohr, 4.105 / 2.514,'xr'); // 赤クロス
// https://github.com/cryos/avogadro/blob/master/crystals/elements/Co-Cobalt.cif
plot(2.5071 / bohr, 4.0686 / 2.5071,'xy'); // 黄クロス


カラーコンターは、全エネルギーが小さいほど青くなり、その格子定数が安定であることを意味しています。
第一原理計算の計算点を補完するために二次元のスプライン補間を行いました(参考: splin2d)。

こうして得られた平衡格子定数が白丸としてプロットしてあります。
これに対して、コバルトの平衡格子定数の文献値は、赤と黄色のクロスシンボルで示してあります。第一原理計算により得られた平衡格子定数は、これらの値と良い一致を示しています。

関連エントリ




参考URL




付録


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


参考文献/使用機器




注意事項


私は、AkaiKKRの開発者である赤井先生、及び、大阪大学の固体電子論グループとは何の関係も無い一個人であり、自分の固体物理の勉強のためにAkaiKKRを利用させていただこうとしている身です。更に言うなら、バンド計算屋さんではなくただの電子工作趣味人です。従って、このブログの内容の正しさに関しては一切保障できませんので、真似をされる方は自己責任でお願いします。おそらくブログ内の記述には、間違っている点もたくさんあると思います。間違いを見つけられた方は、このブログのコメント欄にてお知らせください。

フィードバック



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

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


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


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


tag: Scilab AkaiKKR machikaneyama KKR シェルスクリプト 第一原理計算 

comment

Secret

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

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

最新コメント
リンク

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