cshスクリプトで配列の最大値と最小値

AkaiKKRでコバルトのc/a その2ではcshのシェルスクリプトの中で、配列の中の最大値と最小値を探す方法が分からないので、常に昇順または降順に並んでいると仮定してスクリプトを書きました。

今回は、その点の補足として awk を使った配列の中で最大値と最小値を探すスクリプトを書きました。

#!/bin/csh -f

set LIST=( 1.2 3.00 4.3 2.1 1 3.2 )
set MAX=${LIST[1]}
set MIN=${LIST[1]}

foreach NUM ( ${LIST} )
set MAX=`echo ${NUM} ${MAX} | awk '{if ($1>=$2) print $1; else print $2}'`
set MIN=`echo ${NUM} ${MIN} | awk '{if ($1<=$2) print $1; else print $2}'`
end

echo ${MAX} ${MIN}


...たぶんもっとシンプルな回答があるのだろうとは思いますが。

関連エントリ




参考文献/使用機器




フィードバック



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

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


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


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

tag: シェルスクリプト awk 

AkaiKKRで黄銅の規則・不規則転移

AkaiKKR(machikaneyama)はコヒーレントポテンシャル近似(CPA)を用いて不規則合金の電子状態の計算を行うことができます。今回は、秩序パラメータの異なる立方晶の黄銅(B2型構造とbcc構造の間)の計算を行いました。
更に、その全エネルギーと配置のエントロピーからギブスエネルギーを計算し、相の安定性を議論しました。
その結果、規則・不規則転移温度は430℃となり、この値は実験によって得られている460℃に極めて近い値となっています。

CuZnCIF.png

Fig.1: CsCl型(B2)規則構造のβ-黄銅(CuZn)。Crystallography Open DatabaseのCuZnのcifデータをVESTAで描画した。



β黄銅の規則(B2)相・不規則(bcc)相転移


物質の科学・反応と物性8.2.相転移の種類では二次の相転移の一例として、黄銅(CuZn)の規則不規則転移が紹介されています。

黄銅は、低温ではFig.1に示すようなCsCl型(B2)の規則構造をとります。これを充分高温にすると、体心立方構造(bcc)の不規則構造になります。Fig.1のB2構造は、立方晶の頂点位置に銅原子があり、体心位置に亜鉛原子があります。規則・不規則転移では、温度が上がるに従って、これらの原子がランダムに入れ替わっていき、転移温度(黄銅の場合は460℃)以上では、各サイトのそれぞれの原子の占有率が50%となる、完全に不規則なbcc構造になります。

このように連続的に状態が変わっていくような相転移を二次の相転移と呼びます(ちゃんとした定義はギブスエネルギーのn次微分が不連続になる相転移をn次相転移ということです)。状態の変化を表すパラメータとして秩序パラメータ(規則度)が使われます。いま、体心位置のサイトを亜鉛が占有する割合を r とすると、秩序パラメータは、以下のように表すことができます。
\begin{equation}
\eta = 2r - 1
\end{equation}
秩序パラメータは、完全な規則状態で η=1 完全な不規則状態で η=0 となります。
秩序パラメータは温度に対して連続的に変化しますが(実験的には)比熱の不連続から転移温度を決定することができます。



全エネルギー計算


AkaiKKR(machikaneyama)を用いれば、任意の秩序パラメータを持った黄銅の全エネルギーを計算することができます。

入力ファイルは、単純立方格子の頂点位置に銅を、体心位置に亜鉛を置く構造を基本にして、コヒーレントポテンシャル近似(CPA)を使って占有率を変化させます。格子定数は、すべての組成で a=5.565 Bohr (= 2.945 Å)としました。

c------------------------------------------------------------
go data/CuZn_XIA_ABOHR
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
sc ABOHR , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 1.0 sra pbe nmag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 200 0.01
c------------------------------------------------------------
c ntyp
2
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Cu 2 1 0.0 2
29 XIA
30 XIB
Zn 2 1 0.0 2
30 XIA
29 XIB
c------------------------------------------------------------
c natm
2
c------------------------------------------------------------
c atmicx atmtyp
0 0 0 Cu
1/2 1/2 1/2 Zn
c------------------------------------------------------------


上記が、入力ファイルのテンプレートです。
このテンプレートを元にCシェルスクリプトを利用して、入力ファイルを作成、全エネルギーを計算します。

全エネルギー計算の結果


全エネルギーはB2型規則相で最も低くなりました。この事は、低温では規則構造が現れることと調和的です。

CuZn.png

Fig.3: サイトの占有率と全エネルギー。占有率100%で完全規則状態(B2構造)、50%で完全不規則状態(bcc構造)。


配置のエントロピー


次は、配置のエントロピーを考慮した有限温度での構造について議論します。
相の安定性は、ギブスエネルギーの大小関係から判断することができます(参考:全エネルギーって何だよ?)。
\begin{equation}
G = E + pV - TS
\end{equation}
ここでEは全エネルギー、pは圧力、Vは体積、Tは絶対温度、Sはエントロピーです。大気圧(p≒0)を考えるとpVの項が消えるので、残るのはエントロピーSです。

エントロピーは、振動のエントロピー(参考: AkaiKKRで金属の熱物性)と配置(混合)のエントロピーの2つの項が主要です。今回は、振動のエントロピーは、秩序パラメータにかかわらず一定と仮定して、配置のエントロピーだけを考えることにします。

配置のエントロピーは、スターリングの近似式を用いて、以下のように表すことができます。
\begin{equation}
S=-k_B N \left( \frac{1+\eta}{2}\ln (1+\eta) + \frac{1-\eta}{2}\ln(1-\eta) -\ln 2 \right)
\end{equation}

単位には注意が必要です。
まず、ギブスエネルギーの各項の次元がエネルギーの次元になっていることを確認してください。
次にエネルギーの単位として、何を選ぶかを決めます。ジュール(J)でも、リュードベリ原子単位系(Ry)でも、エレクトロンボルト(eV)でも構いませんが、今回はリュードベリ原子単位系(Ry)とします。するとボルツマン定数も kB = 6.33367 * 10-6 (Ry/K) となります。
さらに原子の数Nもそろえる必要があります。1 molでも 1 原子でも構いませんが、今回は 1 化学式、すなわち N=2 とします。

001_20151021013246770.png

Fig.4: 温度と秩序パラメータの関係


以下のScilabスクリプトを用いて計算した温度と秩序パラメータの関係がFig.4です。
得られた転移温度は430℃となり実験値の460℃と近い値です。また、広い温度範囲にわたって秩序パラメータが変化する二次の相転移の特徴も確認できます。

clear;

// *** 入力パラメータと定数 ***
// ボルツマン定数 6.33367E-6 (Ry/K)
kB = 1.380658E-23/2.17987E-18;
// 原子数
N = 2;

// *** 全エネルギーの読み込み ***
M = fscanfMat("CuZn.txt");
X = M(:,1);
E = M(:,2);
// 全エネルギーの内挿
Xp = linspace(50,99.999,1000);
Ep = interp1(X, E, Xp, 'spline');
// 秩序パラメータ
ETA = 2 .* (Xp ./ 100) - 1;

// *** エントロピー計算 ***
// 配置のエントロピー
Sm = -1 * kB * N * (((1+ETA)./2).*log(1+ETA)+((1-ETA)./2).*log(1-ETA)-log(2));
// 温度ごとにギブスエネルギーを計算
Temp = [0:0.1:600];
Xmin = [];
ETAmin = [];
for i=1:length(Temp) do
T = Temp(i) + 273.15;
Gp = Ep - T .* Sm; // ギブスエネルギーの計算
[gmin, k] = min(Gp); // ギブスエネルギーの最小を探す
// ギブスエネルギーが最小となる占有率と秩序パラメータを保存
Xmin = [Xmin, Xp(k)];
ETAmin = [ETAmin, ETA(k)];
end

// *** グラフのプロット ***
plot(Temp,ETAmin);
xlabel("Temperature (deg C)")
ylabel("Order parameter");


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR Scilab machikaneyama KKR CPA 

AkaiKKRでLDAとGGA その2

AkaiKKRでLDAとGGA その1と同様に、強磁性の体心立方構造(bcc)鉄の全エネルギーとスピン磁気モーメントを格子定数の関数としてAkaiKKR(machikaneyama)で計算しました。

FeMult.png

Fig.1: 鉄の格子定数と全エネルギー、スピン磁気モーメントの関係


前回と同様に、LDA(mjw)は格子定数を過小評価、GGA(pbe)は過大評価しました。


強磁性bcc鉄


AkaiKKRでLDAとGGA その1では非磁性の面心立方構造(fcc)の銅の格子定数をLDA(mjw)とGGA(pbe)で、それぞれ計算しました。その結果、LDAは格子定数を過小評価し、GGAは過大評価することが確認できました。

今回は、スピン分極を考慮した強磁性の体心立方構造(bcc)の鉄に対して、同様の計算を行います。

結果


結果をFig.1に示します。AkaiKKRでLDAとGGA その1と同様に、格子定数に関してはLDA(mjw)が過小評価、GGA(pbe)が過大評価となりました。実験値との比較では、GGA(pbe)の方が近い値となっていることが分かります。

同時にスピン磁気モーメントもプロットしてあります。磁性入門―スピンから磁石まで (材料学シリーズ)に書いてある鉄の磁気モーメントの実験値は M = 2.218 μB です。

LDA(mjw)によって得られた格子定数のときのLDAの磁気モーメントは、かなり実験値に近い値ですで、ほんのわずかに過小評価してるようです。LDAの磁気モーメントは、格子定数が大きくなるほど大きくなっています。GGA(pbe)によって得られた格子定数でのスピン磁気モーメントは、実験値に対して過大評価で、そのずれはLDAによるものよりも大きくなりました。

そんなわけで、強磁性bcc鉄のスピンモーメントの計算に、LDAの代わりにGGAを使うと、一見すると実験との一致が悪くなるという見方もできます。しかし、GGAは一般的にLDAよりも磁性体の計算を改善する傾向にあるということが、知られているようです。
第一原理計算の結果の正しさを検証する上で、実験値のと比較を行うというのは、最も王道的な方法ではありますが、どの程度の確度で議論ができるのかはよく考えないといけないかもしれません。いずれにせよ、今回の計算ではLDA、GGAどちらも格子定数、磁気モーメントをそれなりに良い確度で再現できていると考えられると思っています。

関連エントリ




参考URL




付録


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


参考文献/使用機器







フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR LDA GGA 

AkaiKKRでLDAとGGA その1

密度汎関数理論(DFT)を用いた第一原理計算において交換相関ポテンシャルに何を使うかは、最も経験的に決めなければならないパラメータのひとつです。AkaiKKR(machikaneyama)でもLDA/GGAのなかからいくつかの選択肢があります。今回はmjw型のLDAとpbe型のGGAを使って面心立方構造(fcc)の銅の格子定数を計算してみました。

Cu.png

Fig.1: 銅の格子定数と全エネルギーの関係


その結果、LDAは格子定数を過大評価し、GGAは過小評価する傾向にあることが確認できました。


交換相関ポテンシャル


密度汎関数理論(DFT)を利用した第一原理計算ソフトでは、その交換相関項を表すために何らかの近似を行います。非常に良く使われるのが局所密度近似(LDA)と一般勾配近似(GGA)です。

ありていに言うとGGAの方がLDAよりも一歩進んだ「えらい」近似であるということになっています。しかしながら、実際に物性値を計算してみると必ずしもGGAのほうが良い結果になるとも限らないようです。更に言うならば、GW近似などの、更に進んだ計算手法と比べるとLDAもGGAも五十歩百歩としてどちらも一まとめにLDAと扱われることもあるようです。

AkaiKKR(machikaneyama)も交換相関項にLDAとGGAのどちらかを選ぶことができます。LDA/GGAの中にもいくつかの定式化があるようでAkaiKKRでは以下の選択肢があります。(詳しくは source/potenv.f のコメントを参照)

LDA: vbh, mjw, vwn
GGA: gga91, pbe, ev

いずれにせよ、交換相関ポテンシャルの項に対して、LDA/GGAといった近似を持ち込むことが、DFT計算において最も程度の大きい近似であるので、これらの項にどの近似を使うことによって、結果にどのような影響があるのかは知っておいたほうが良いです。

今回は面心立方構造(fcc)の銅の格子定数を全エネルギーの最小化から求めました。LDAの代表として mjw を選び、GGAの代表としては pbe を選びました。

計算結果


Fig.1に銅の格子定数と全エネルギーの関係をグラフに示します。mjw型のLDAでは格子定数が a=6.75 Bohr 程度になることが分かります。
これに対してpbe型のGGAを用いると a=6.95 Bohr 付近に全エネルギーの最小が存在します。グラフの中に垂直の破線で示したところが、実験から得られている格子定数です。
LDAは格子定数を過小評価、GGAは過大評価していることが分かります。

この傾向は、銅に限った話ではなくLDAは密度の高い相を導きやすくGGAはより隙間の大きい構造を好む様です。(参考: 物質の電子状態〈下〉 (SPRINGER UNIVERSITY TEXTBOOKS))

なおグラフの縦軸は任意目盛りです。そして、LDAとGGAでも別の軸を取っているので、相対的な比較はできません。それぞれの曲線の中で全エネルギーが最小となる格子定数がどの程度になるのかだけ見てください。

関連エントリ




参考URL




付録


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


参考文献/使用機器






フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR LDA GGA 

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

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

最新コメント
リンク

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