Scilabで二重積分

Scilabを利用すると1変数の数値積分が簡単に計算できます。

\begin{equation}
\int_{x_0}^{x_1}f(x)\mathrm{d}x
\end{equation}

このブログでも数値積分タグにいくつかの例を見つけることができます。しかしながら、2変数の数値積分はこれまで行ってきませんでした。

\begin{equation}
\int\int f(x,y) \mathrm{d}x\mathrm{d}y
\end{equation}

Scilabには二重積分を計算することが可能な int2d が存在します。今回は高校数学の美しい物語で解析的に解かれている二重積分を数値的に計算してみます。


積分範囲が長方形領域の場合


積分範囲が長方形の領域の場合、すなわち以下のような式で表すことができる場合は、簡単に数値積分できます。

\begin{equation}
\int_{x_0}^{x_1}\int_{y_0}^{y_1}f(x,y)\mathrm{d}x\mathrm{d}y
\end{equation}

Scilabの int2d では長方形領域を2つの三角形のパッチワークとして与えます。
積分範囲を int2d に渡すために行列XとYを用意します。それぞれ2つの三角形の頂点のx座標とy座標を与えます。

\begin{equation}
X =
\begin{pmatrix}
x_{0} & x_{0} \\
x_{1} & x_{1} \\
x_{1} & x_{0} \\
\end{pmatrix},
Y =
\begin{pmatrix}
y_{0} & y_{0} \\
y_{0} & y_{1} \\
y_{1} & y_{1} \\
\end{pmatrix}
\end{equation}

001_20170423145553b2b.png
Fig.1: Scilabのint2dへの積分範囲の与え方


実際に以下の積分を計算して見ます。

\begin{equation}
\int_{0}^{\pi}\int_{0}^{R}x^4 \sin(y)\mathrm{d}x\mathrm{d}y
\end{equation}

clear;

r = 1;

// *** 積分する関数の定義 ***
function z = func(x,y)
z = (x .^ 4) * sin(y)
endfunction
// 積分範囲
x0 = 0; x1 = r;
y0 = 0; y1 = %pi;

// *** 二重積分 ***
X = [x0, x0;
x1, x1;
x1, x0];
Y = [y0, y0;
y0, y1;
y1, y1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = 2*(r^5)/5


数値化解と解析解が同じ値になることが確認できます。

積分範囲が三角形の組み合わせで表せる場合


積分範囲が長方形の場合は2つの三角形の組み合わせで表現されますが、より複雑な形状の場合も任意の個数の三角形の組み合わせで表現できるはずです。今回は逆に簡単になってしまいますが、1個の三角形で表現できる例を計算します。

\begin{equation}
\int \int_D xy^2 \mathrm{d}x\mathrm{y}
\end{equation}

jusekibun.png
Fig.2: 積分領域Dが三角形ひとつ分の例


積分領域が三角形ひとつ分なので、与える行列は3行1列になります。

\begin{equation}
X =
\begin{pmatrix}
x_{0} \\
x_{1} \\
x_{1} \\
\end{pmatrix},
Y =
\begin{pmatrix}
y_{1} \\
y_{0} \\
y_{1} \\
\end{pmatrix}
\end{equation}

この計算を行うScilabスクリプトは以下のようになります。

clear;

// *** 積分する関数の定義 ***
function z = func(x,y)
z = x .* (y .^ 2)
endfunction

// *** 二重積分 ***
X = [0;
1;
1];
Y = [1;
0;
1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = 3/20


このスクリプトも数値解と解析解が同じに値になることが分かります。

同様にしてN個の三角形の組み合わせで表現される積分範囲の場合3行N列の行列で指定することができます。

更に複雑な積分領域の場合


どんなに複雑な積分領域の形状であっても三角形のパッチワークで表現できるはずですが、現実的には大変です。そこでOctaveの精義―フリーの高機能数値計算ツールを使いこなすで紹介されている方法を試してみましたが、現状うまく行っていません。上手く行っていませんがとりあえず方法だけは紹介します。
具体的にはScilabの論理演算で条件分岐の考え方を使って積分領域外では値がゼロになるように被積分関数の定義を行います。

\begin{equation}
\int\int_D -\frac{1}{(2x + y + 1)^2}\mathrm{d}x\mathrm{d}y
\end{equation}

jusekibun2.png
Fig.3: 複雑な積分領域の例


clear;

// *** 積分する関数の定義 ***
function z = func(x,y)
region = y >= x .^ 2
z = - 1 ./ ((2 * x + y + 1) .^ 2) .* region
endfunction

// *** 二重積分 ***
X = [0;
1;
1];
Y = [0;
0;
1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = (1/3) * log(4) - 1/2


原理的にはこのスクリプトでよいはずですが、実際には正しく計算してくれません。Scilab 6.0ではエラーで停止します。Scilab 5.5.2ではそれっぽい値を返しますが、解析解の値とはかなりずれた値となっており、不正確です。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 数値積分 重積分 

デバイ模型

デバイ模型やデバイ温度といった言葉は、固体物理学の分野でたびたび目にします。
近似自体は、非常に単純なものなのですが、あらゆる物性に顔を出す反面、初歩的な教科書には比熱のことしか書いてなかったりするので、ときどき良く分からなくなったりします。
デバイ模型は、そもそも、フォノンの分散関係を近似的に表すための模型であり、そこから(格子比熱を含む)色々な物性を計算することができるようになります。
今回は、デバイ模型の仮定している近似についてまとめると共に、デバイ温度が現れる物性の計算式を列挙します。

001_201511080539137af.png

Fig.1: Quantum ESPRESSOで計算したダイヤモンドのフォノンの分散関係



フォノンの分散関係と物性値


結晶の格子振動(フォノン)にだけに依存する物性値は、フォノンの分散関係が分かればすべて計算できるということになっています。デバイ模型(Debye model)やアインシュタイン模型(Einstein model)は、単純に格子比熱を計算するためだけのモデルだと見做される事もありますが、本来はこの分散関係を近似的にあらわす模型です。したがって、デバイ模型を使えば格子比熱以外の物性も計算することができます。

フォノンの分散関係は、以下の式で表され、一般には複雑な形状をしています。

\begin{equation}
\omega = \omega(\vec{q})
\end{equation}

ここで$\omega$は角振動数で$\vec{q}$は波数ベクトルです。Fig.1に示したのは、Quantum ESPRESSOで計算したダイヤモンドのフォノンの分散関係です。音響フォノンと光学フォノンの両方を持っていること、異方性を持っていることなどが分かります。

デバイ模型の2つの仮定


そこで、デバイ模型では「分散関係をシンプルな形で表す」「異方性をなくす」という二つの方針で簡単な式に置き換えます。

仮定(1) 角振動数$\omega$は波数$q$に比例し、その比例定数を音速$v$とする
\begin{equation}
\omega = v q
\end{equation}

仮定(2) ブリュアンゾーン(BZ)を同じ体積の球で置き換える
\begin{equation}
q_D = \left( \frac{6 \pi^2 N}{\Omega} \right)^{\frac{1}{3}}
\end{equation}
ここで$q_D$はデバイ波数と呼ばれ、体積Ωとその中の原子の数Nのみで表すことができます。

この波数$q$は、$0 \leqq q \leqq q_D$の範囲に限られます。この$q_D$は、体積Ωと原子数Nだけで決まるので、デバイ模型を特徴付けるパラメータは、音速$v$のみということになります。しかしながら、音速$v$はこの用途にはあまり使われず、代わりにデバイ温度$\Theta_D$がよく使われます。音速$v$とデバイ温度$\Theta_D$は、いったんエネルギーの次元で考えた後、以下のように換算できます。
\begin{equation}
k_B \Theta_D = \hbar \omega_D = \hbar v q_D \\
\therefore v = \frac{k_B \Theta_D}{\hbar q_D}
\end{equation}

以降では、デバイ模型のパラメータであるデバイ温度 ΘD が出てくる物性値について列挙します。

熱力学


ヘルムホルツの自由エネルギー(ヘルムホルツエネルギー)は、以下のような式で、エネルギー E の項とエントロピー S の項から表されます。
\begin{equation}
F=E-TS
\end{equation}
エネルギーにせよエントロピーにせよ、色々な原因による項(磁性、配置のエントロピー、電子比熱によるエントロピーなど)がありますが、格子系のエネルギーとエントロピーをデバイ模型を用いて表します。(参考: AkaiKKRで金属の熱物性)

格子系のエネルギー
\begin{equation}
E_D = 3 N k_B T \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D / T} \frac{z^3}{\exp(z)-1} \mathrm{d}z + E_0
\end{equation}

ゼロ点エネルギー
\begin{equation}
E_0 = \frac{9}{8} N k_B \Theta_D
\end{equation}

格子系のエントロピー
\begin{equation}
S_D = 3 N k_B \left[ \frac{4}{3} \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D / T} \frac{z^3}{\exp(z) - 1} \mathrm{d}z \\ - \ln \left\{ 1 - \exp \left(- \frac{\Theta_D}{T} \right) \right\} \right]
\end{equation}

定積比熱は格子系のエネルギーを、温度で微分したものです。
\begin{equation}
C_V = \frac{\partial E_D}{\partial T} = 9 N k_B \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D / T} \frac{z^4 \exp(z)}{(\exp(z)-1)^2} \mathrm{d} z
\end{equation}

格子振動の大きさ


デバイワラー因子
\begin{equation}
M_D = \frac{6 h^2 T}{M k_B \Theta_D^2} \left( \frac{T}{\Theta_D} \int_0^{\Theta_D / T} \frac{z}{\exp(z)-1} \mathrm{d} z + \frac{\Theta_D}{4T} \right) \left( \frac{\sin \theta}{\lambda} \right)^2
\end{equation}

平均二乗変位
\begin{equation}
u_D^2 = \frac{3 \hbar T}{M k_B \Theta_D^2} \left( \frac{T}{\Theta_D} \int_0^{\Theta_D/T} \frac{z}{\exp(z)-1}\mathrm{d}z + \frac{\Theta_D}{4T} \right)
\end{equation}

原子の振幅が一定以上大きくなると、融解が起こるのではないかという考えから、振幅と融点の関係式が提案されました。当然ながら、かなりラフな推定ですが、少なくとも定性的には、デバイ温度の高い物質の融点が高いという関係が見られるようです。(参考: ザイマン 固体物性論の基礎)

リンデマンの融解公式
\begin{equation}
T_m = \frac{x_m^2}{9 \hbar} M k_B \Theta_D r_s^2
\end{equation}

電気抵抗率


金属の電気伝導は、格子系だけでは決まらないので、電子と格子の相互作用に起因するパラメータも必要となります。

電気抵抗率
\begin{equation}
\rho_D = A \left( \frac{T}{\Theta_D} \right)^5 \int_0^{\Theta_D / T} \frac{z^5}{(\exp(z)-1)(1-\exp(-z))} \mathrm{d} z
\end{equation}

BCS理論による超伝導転移温度
\begin{equation}
T_c = \Theta_D \exp \left( - \frac{1}{N(0)V} \right)
\end{equation}

つまりデバイ模型/デバイ温度とは何なのか?


結局のところ、デバイ温度とかデバイ模型とか言われるものは、いったい何なのでしょうか?
結晶の比熱を計算するための近似がデバイ模型であり、デバイ温度はその中のパラメータである、というのはひとつの答えです。
実際、ねがてぃぶろぐでもScilabで数値積分: 固体の比熱のエントリで比熱の計算をしています。

001_20130522024102.png
Fig.2: 銅の比熱の温度依存性(ΘD = 343.5 K)。青破線がアインシュタインモデル、赤実線がデバイモデルによる計算。


しかしながら、上記のように、比熱以外の色々な物理量の計算式に顔を出します。また、直接数式の中に出てこない場合でも、例えば「デバイ温度よりも充分高温のときだけ成り立つ」みたいな文脈で現れることもあります。あるいは「ダイヤモンドは固いのでデバイ温度が高い」のような記述が見つかるかもしれません。

温度がデバイ温度よりも高いかどうかを議論しているときは、格子振動を古典力学的に扱ってよいのか、量子力学的な効果を考慮しなければならないのかを問題にしている場合が多いです。例えばFig.2の銅の比熱の温度依存性では、低い温度では比熱が急激に変化していますが、高温ではほとんど変化しなくなります。高温での比熱の温度依存性がなくなることは、古典力学的から導かれ、デュロン・プティの法則として知られています。その一方で、低温の急激な温度依存性は、量子力学を考えなければ説明できません。
また、物質の硬さを議論しているときには、直感的に言えば、デバイ温度を原子同士をつなぐ「ばね定数」のようなものだと考えています。

そんなわけで、最初に書いたとおり「デバイ模型とは、結晶のフォノンの分散関係に対する最も簡単な近似のひとつ」で「デバイ温度は、デバイ模型の中で、その物質の個性を表す最も重要なパラメータ」が答えとなるわけですが、その近似が単純であるにもかかわらず、色々な応用が考えられるので、出くわすたびに、それまでもっていた理解とは、ほんの少し違う側面が見えることがあります。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: デバイ模型 フォノン Scilab quantumESPRESSO 

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 

Scilabで自発磁化の温度依存性

磁性入門―スピンから磁石まで (材料学シリーズ)では、ワイスの分子場近似を用いた磁化の温度依存性を解析的に表す式が以下のように書かれています。

_eq_htcs.png

今回は、これをScilabの非線型方程式ソルバfsolveを用いて数値的に解きます。(参考:Scilabで非線形方程式ソルバ その1)

001_2015081913143435a.png
Fig.1: 遷移金属の磁気モーメントの温度依存性



自発磁化の温度依存性


磁性入門―スピンから磁石まで (材料学シリーズ)によると、スピン角運動量 S=1/2 の場合、自発磁化 M(T) は

M(T) = N \mu_B \tanh \left( \frac{\alpha \mu_B M(T)}{k_B T} \right)

ここで T=0 (K) のときの磁化 M0=M(0) は

M_0 = N \mu_B

であり、キュリー温度TC

T_C = \frac{\alpha N \mu_B^2}{k_B}

したがって

\frac{M(T)}{M_0} = \tanh \left(\frac{M(T)}{M_0} \frac{T_C}{T} \right)

となるのでScilabなどの非線型方程式ソルバ(fsolve)で解く事が可能です。
したがって基底状態の自発磁化M0とキュリー温度TCの二つのパラメータから磁化温度曲線を内挿することができます。

数値計算


実際にScilabを使って計算した結果が冒頭のFig.1です。
赤、緑、青の点はCrangle and Goodman (1971)により報告されている、鉄・ニッケル・コバルトの実験値です。
大雑把な仮定と簡単な式から計算したことを考えると、そんなに悪くない結果ではないかと思います。

なお、M0とTCはどちらもAkaiKKR(machikaneyam)を用いて計算できます。しかしながら、計算されるキュリー温度はかなりの過大評価です。例えば、体心立方構造(bcc)鉄の場合は、TC=1043 Kの実験値に対して、AkaiKKRの計算では TC=1616 Kと約55%も過大評価となります。(AkaiKKRで鉄のキュリー温度の計算のときはTC=1229 Kでしたがbzqltyを上げてpbeで計算するとTC=1616 Kとなりました。いずれにせよ過大評価です。)

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 非線型方程式ソルバ 強磁性 キュリー温度 平均場近似 

AkaiKKRでBain機構 その1

体心立方構造(bcc)と面心立方構造(fcc)は、一見すると全く別の構造のように見えますが、体心正方構造(bct)のc/aを変化させたものであると考えると、実は非常によく似た構造であることがわかります。

001_201506220522417ce.jpg

Fig.1: bcc/fcc構造のBainの対応関係(Xie et al. (2008)より)


今回は、AkaiKKRを用いてbct鉄の全エネルギーを格子体積とc/aを変化させながら計算することを前提に、入力ファイルのためのマフィンティン半径の設定について考えます。
AkaiKKRでコバルトのc/a その2のときと同様に「格子体積に対するマフィンティン球体積の比が一定になる範囲で、マフィンティン半径を最大にとる」ならばbcc構造においてマフィンティン球同士がタッチングする半径ですべての計算を行えばいいことがわかります。


Bainの機構


体心立方構造(body-centered cubic; bcc)と面心立方構造(face-centered cubic; fcc)は、どちらも体心正方構造(body-centered tetragonal; bct)のc/aが特別な値の場合と考えることができます。すなわちc/a=1のときbcc構造となり、c/a=√2≒1.414のときにfcc構造となります。

見方・考え方 合金状態図によると、これはBainが鋼のマルテンサイト変態を説明しようとしたことに端を発するとのことです。ただし、実際のマルテンサイト変態はBainの考えとは異なるメカニズムで起こっているとのことですが。

今回は、AkaiKKR(machikaneyama)このBainの対応関係をみるために、体心正方構造(bct)を計算セルとして格子体積とc/aを変化させながら全エネルギーの変化を確認します。その1である今回は、マフィンティン半径をどのようにするかについて確認します。

マフィンティン半径


AkaiKKRはマフィンティン近似を用いています。マフィンティン近似では、マフィンティン球の半径を計算パラメータとして与えなければなりません。そして、同じ結晶を計算している場合であっても、得られる全エネルギーはマフィンティン半径によって変化してしまいます。したがって体積やc/aを最適化するために全エネルギーを比較するときは、マフィンティン半径の設定が食い違わないようにする必要があります。

どういう設定にするのがベストなのかは、難しい問題らしいのですが、今回は差し当たりAkaiKKRでコバルトのc/a その1その2のときと同じ条件を考えます。
すなわち「格子体積に対するマフィンティン球体積の比が一定になる範囲で、マフィンティン半径を最大にとる」とします。

AkaiKKRでコバルトのc/a その2のときと同様にScilabを用いて、体積を固定しながらc/aを変化させたときに、マフィンティン球同士が接触する半径を計算しました。

clear;

// 体心正方(bct)の格子体積を固定
v = 100;
// bctのc/aを変化
// c/a = 1: bcc
// c/a = √2: fcc
//eta = linspace(1, sqrt(2));
eta = linspace(1, 1.5);

// 格子定数aの計算
a = (v ./ eta) .^ (1/3);

// ab面内で頂点同士が触れる時
rmt_ab = a ./ 2;
// // ac面内で頂点同士が触れる時
// rmt_ac = eta .* a ./ 2;
// 頂点と体心で触れる時
rmt_diag = sqrt(2 .* (a .^ 2) + (eta .^ 2) .* (a .^ 2)) ./ 4;

// *** グラフの描画 ***
// ab面内
plot(eta, rmt_ab, '-b');
// // ac面内
// plot(eta, rmt_ac, '-g');
// 頂点と体心
plot(eta, rmt_diag, '-r');
// グラフの装飾
xlabel("c/a");
ylabel("muffin-tin radius");


002_20150622052240f6f.png

Fig.: bct構造の格子体積をΩ=100 Bohr3とした場合のタッチング時のMT半径のc/a依存性。赤線は頂点位置にある原子と体心位置にある原子が先に接すると仮定した場合、青線はa-b軸方向の原子が先に触れると仮定した場合。これら二つの値はc/a=√2(fcc構造のとき)で一致する、c/a=1(bcc構造)のときとることのできるMT半径が最小になる。


体心立方構造(bcc)の充填率は0.68で面心立方構造(fcc)の充填率は0.74なので、fcc構造の方が密な構造であると言えます(参考: Wikipedia: 空間充填率)。マフィンティン半径もこのことを反映して、c/a=1(bcc構造に相当)のときに最小値を取ることがわかります。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR マフィンティン半径 シェルスクリプト Scilab 

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

最新コメント
リンク

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