AkaiKKRで結晶の回転

AkaiKKRでは、直交座標系(x, y, z)のなかに結晶の軸(a, b, c)を置いて第一原理計算を行っています。この際、計算の種類によっては直交座標系のz方向が特別扱いされることがあります。例えばスピン軌道相互作用を入れた計算では、結晶のスピンがz軸方向に揃うとして計算します。

AkaiKKRでコバルトの結晶磁気異方性では、六方最密充填のc軸方向をz軸方向に置いた場合(デフォルト設定)とx軸周りに徐々に回転させて行ったときで全エネルギーがどのように変化するかを計算しました。

この様な回転をさせる際に、結晶構造をブラべ格子ではなく、基本ベクトルでしていすると便利だと書きました。実は、その後の2018年2月15日のアップデートで、より簡単に結晶を回転させる方法がAkaiKKRに追加されています。今回は、その方法を使ってコバルトを回転させてみます。


最後に、結晶を傾けることができるようになった点。

Also, it is now possible to tilet the crystal axis. It is
done by using "tlt" keyword attached after bravai type
or "aux" ("prv"), like

fcctlt 0 90 0 6.65

or

auxtlt -0.5 0.5 0.5 0.5 -0.5 0.5 0.5 0.5 -0.5 0 90 0 5.27

In both cases crystal is tilted by Euler angle (0, 90, 0).

AkaiKKRでは磁化を(数値計算上の)z軸方向にとるので、磁気異方性を計算するときなどは結晶を回転させなければいけなかったのですが、これが簡単に出来るようになったということですね。


入力ファイル


通常のhcp Coの入力ファイルは以下のようになります。

c----------------------Co------------------------------------
go data/Co
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
hcp 4.74 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 0.8 srals pbe mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 12 200 0.023
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Co 1 1 0.0 2
27 100
c------------------------------------------------------------
c natm
2
c------------------------------------------------------------
c atmicx atmtyp
0a 0b 0c Co
1/3a 2/3b 1/2c Co
c------------------------------------------------------------


これに対して、結晶を回転させるためにtltとオイラー角を追加したものが以下になります。
とりあえず、全く回転させない場合を見てみます。

c----------------------Co------------------------------------
go data/Co_tlt_0
c------------------------------------------------------------
c brvtyp angle a c/a b/a alpha beta gamma
hcptlt 0 0 0 4.74 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 0.8 srals pbe mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 12 200 0.023
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Co 1 1 0.0 2
27 100
c------------------------------------------------------------
c natm
2
c------------------------------------------------------------
c atmicx atmtyp
0a 0b 0c Co
1/3a 2/3b 1/2c Co
c------------------------------------------------------------


brvtyphcpからhcptltに変わり、その後にオイラー角の数値3つが追加されています。上の例では全く回転させないので0 0 0です。

次にx軸周りに90度回転させて、結晶のc軸をxy平面の中に持ってきます。x軸周りに回転させる場合、オイラー角の2番目の数値を変更します。

c----------------------Co------------------------------------
go data/Co_tlt_90
c------------------------------------------------------------
c brvtyp angle a c/a b/a alpha beta gamma
hcptlt 0 90 0 4.74 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 0.8 srals pbe mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 12 200 0.023
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Co 1 1 0.0 2
27 100
c------------------------------------------------------------
c natm
2
c------------------------------------------------------------
c atmicx atmtyp
0a 0b 0c Co
1/3a 2/3b 1/2c Co
c------------------------------------------------------------


関連エントリ




参考URL




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR 強磁性 コバルト 結晶磁気異方性 基本並進ベクトル 相対論 スピン軌道相互作用 六方最密充填構造 

質量パーセントと原子パーセント その2

最近はハイエントロピー合金と呼ばれる5種類以上の元素が同程度ずつ含まれた合金の研究が流行っています。実験の分野では、合金の濃度は質量パーセント(wt.%)で表される事が多いですが、AkaiKKRなどの計算の分野では、原子パーセント(at.%)で表されている方が便利です。
質量パーセントと原子パーセントでは、二元合金の濃度の質量パーセント(wt.%)と原子パーセント(at.%)の換算を行いました。今回は、多元合金の間で同様の換算を行います。


計算方法


$n$元合金の$i$番目の元素の原子量を$m_i$、その原子パーセントでの濃度を$100 \times x_i$ (at.%)、質量パーセントでの濃度を$100 \times y_i$ (wt.%)と表すことにします。

原子パーセントから質量パーセントへの変換


全体の質量$M$は
\[
M = \sum_{i=1}^{n} m_i x_i
\]

元素$i$の質量比$y_i$は
\[
y_i = \frac{m_i x_i}{M}
\]

例としてCo(30at.%), Cr(25at.%), Fe(20at.%), Mn(15at.%), Ni(10at.%)合金の濃度を質量パーセントに換算するScilabスクリプトを以下に示します。
結果はCo(31.595034wt.%), Cr(23.229955wt.%), Fe(19.959602wt.%), Mn(14.726585wt.%), Ni(10.488825wt.%)となりました。

clear;

// *** 原子量 ***
mCo = 58.933195;
mCr = 51.9961;
mFe = 55.845;
mMn = 54.938044;
mNi = 58.6934;
M_i = [mCo; mCr; mFe; mMn; mNi];

// *** 原子比 x ***
xCo = 0.3;
xCr = 0.25;
xFe = 0.2;
xMn = 0.15;
xNi = 0.1;
X_i = [xCo; xCr; xFe; xMn; xNi];

// *** 質量の計算 ***
m_tot = sum(M_i .* X_i);

// *** 質量比の計算 ***
Y_i = (M_i .* X_i) ./ m_tot;
// 質量パーセント(wt.%)
100 * Y_i


質量パーセントから原子パーセントへの変換


全体の原子数$N$は
\[
N = \sum_{i=1}^{n} \frac{y_i}{m_i}
\]

元素$i$の原子比$x_i$は
\[
x_i = \frac{\frac{y_i}{m_i}}{N}
\]

例としてCo(30wt.%), Cr(25wt.%), Fe(20wt.%), Mn(15wt.%), Ni(10wt.%)合金の濃度を原子パーセントに換算するScilabスクリプトを以下に示します。
結果はCo(28.416343at.%), Cr(26.839607at.%), Fe(19.991833at.%), Mn(15.241404at.%), Ni(9.5108131at.%)となりました。

clear;

// *** 原子量 ***
mCo = 58.933195;
mCr = 51.9961;
mFe = 55.845;
mMn = 54.938044;
mNi = 58.6934;
M_i = [mCo; mCr; mFe; mMn; mNi];

// *** 質量比 y ***
yCo = 0.3;
yCr = 0.25;
yFe = 0.2;
yMn = 0.15;
yNi = 0.1;
Y_i = [yCo; yCr; yFe; yMn; yNi];

// *** 原子数の計算 ***
n_tot = sum(Y_i ./ M_i);

// *** 原子比の計算 ***
X_i = (Y_i ./ M_i) ./ n_tot;
// 原子パーセント(at.%)
100 * X_i


関連エントリ




フィードバック



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

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


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


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

tag: Scilab 

Scilabで一次元周期的境界条件の分散関係

Scilabで一次元井戸型ポテンシャルの分散関係では金持徹 編著 固体電子論の§2.5で紹介されている一次元井戸型ポテンシャルの$E-k$関係をプロットしました。今回は一次元の周期境界条件の$E-k$関係をプロットします。


Periodic.png
Fig.1 周期的境界条件の$E-k$関係



Scilabスクリプト


Scilabスクリプトは以下のようになりました。
clear;

// *** リュードベリ(Rydberg)原子単位系 ***
hbar = 1; // 換算プランク定数
h = 2 * %pi * hbar; // プランク定数
m = 0.5; // 電子の質量
a0 = 1; // ボーア半径

// *** 系の設定 ***
// 原子の数
//n_atom = 20;
n_atom = 100;
// 原子間距離
r0 = 1.0 * a0;
// 結晶の長さ
l = n_atom * r0;
// 価数と電子の総数
n_valence = 1;
n_electron = n_valence * n_atom;

// *** 分散関係の計算 ***
g = 0;
// 波数ベクトル (式 3.11)
kg = 2 * %pi / l * g;
Kg = [kg];
// エネルギー (式 3.13)
eg = hbar * hbar / (2 * m) * (2 * %pi / l * g) * (2 * %pi / l * g);
Eg = [eg];
for g=1:int((n_electron - 2) / 4)
// 波数ベクトル (式 3.11)
kg = 2 * %pi / l * g;
Kg = [-kg, Kg, kg];
// エネルギー (式 3.13)
eg = hbar * hbar / (2 * m) * (2 * %pi / l * g) * (2 * %pi / l * g);
Eg = [eg, Eg, eg];
end
// フェルミ波数ベクトル (1/Bohr)
kf = kg
// フェルミエネルギー (Ry)
ef = eg
// gの最大値
gmax = g;

// *** グラフの描画 ***
plot(0, 0); // グラフ描画範囲のためのダミー
plot(Kg, Eg, ".b");

// *** グラフの装飾 ***
xgrid(color("gray"));
xlabel("Wave vector (1/Bohr)");
ylabel("Energy (Ry)");
//zoom_rect([-kf, 0, kf, ef]);


原子数Nの影響


原子数 N = 100 個で計算してみると、フェルミエネルギー$E_f$やフェルミ波数ベクトル$k_f$がScilabで一次元井戸型ポテンシャルの分散関係のときとわずかに異なりました。そこで原子数を変えながらフェルミエネルギーがどの様に変わるかを確認しました。


Periodic_g.png
Fig.2 原子の数とフェルミエネルギーの関係。破線で示したのが井戸型ポテンシャルのフェルミエネルギー。


周期的境界条件の計算では、十分に多数の原子が存在すると仮定していますが、原子の数が千個から一万個程度は無いと正しいフェルミエネルギーが得られないようです。

関連エントリ




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 固体電子論  

Scilabで複素数の進行波

Scilabで定在波と進行波では金持徹 編著 固体電子論の§3.1 数学的準備(周期関数)で紹介されている定在波と進行波を位置と時間の実数関数としてプロットしてみました。
これに対して電子の波動関数は、複素数として表してある方が便利です。今回は式(3.8)で示されている複素数の進行波(の実部)をScilabでプロットしました。

TravelingWave.png

Fig.1 進行波のプロット。実関数としてプロットした場合と同じ結果が得られる(ので画像は使いまわし)。



複素数の進行波


位置のみに依存する波は以下の式で表されます。
\[
\psi(x) = A \exp ikx
\]
時間のみに依存する波は、以下の式で表されます。
\[
\phi(t) = B \exp(- i 2 \pi f t)
\]
進行波はこれらの積で表されます。
\[
\Psi(x, t) = \psi(x)\phi(t) = C \exp i (kx - 2 \pi f t)
\]

Scilabスクリプト


Scilabスクリプトは以下のようになりました。
clear;

// *** パラメータ ***
a = 1;
b = 1;
c = a * b;
lambda = 1;
k = 2 * %pi / lambda;
f = 1;

// *** 位置に関する周期関数 ***
function phi_x = phi_x(x)
// phi_x = a .* cos(2 * %pi / lambda .* x)
phi_x = a .* exp(%i * k .* x)
endfunction

// *** 時間に関する周期関数 ***
function phi_t = phi_t(t)
// phi_t = b .* cos(2 * %pi * f .* t)
phi_t = b .* exp(- %i * 2 * %pi * f .* t)
endfunction

// *** 進行波 ***
function phi = traveling(x, t)
phi = phi_x(x) .* phi_t(t)
endfunction

// *** 進行波の計算 ***
X = linspace(0,1);
T = linspace(0,1);
[X2D, T2D] = ndgrid(X, T);
PHI_travering = traveling(X2D, T2D);

// *** 進行波のプロット ***
xlabel("position, x");
ylabel("time, t");
zlabel("PHI(x,t)");
mesh(X2D, T2D, real(PHI_travering));
//surf(X2D, T2D, real(PHI_travering));


関連エントリ




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 固体電子論  

Scilabで定在波と進行波

金持徹 編著 固体電子論の§3.1 数学的準備(周期関数)で紹介されている定在波と進行波を位置と時間の関数としてプロットしてみました。

余談ですが、§3.1は数学的準備という名前がついていながら、§3の中に波動関数を扱うパートはありません。実をいうと§2で扱った井戸型ポテンシャルの中の電子は定在波であるのに対して、§3.2で扱う周期的境界条件の中での電子は進行波です。著者は、数学的な準備のためというよりは、扱う波の種類の違いをイメージしやすいように§3.1でこれらの波の違いに言及したのかもしれません。


定在波


定在波の式は、位置$x$と時間$t$の関数として、以下のようにあらわされます(式3.3)。
\[
\Psi(x, t) = C \left(\cos \frac{2 \pi x}{\lambda} \right) \cdot (\cos 2\pi f t)
\]
これをScilabで描画すると以下のようになりました。

StandingWave.png
Fig.1 定在波のプロット。波が時間経過に伴って振動しているが、波の腹になっている位置、節になっている位置は変化していない。


進行波


進行波の式は、以下のようにあらわされます(式3.4)。
\[
\Psi(x, t) = A \cos 2 \pi \left( \frac{x}{\lambda} - f t \right)
\]
Scilabのプロットは以下の通りです。

TravelingWave.png

Fig.2 進行波のプロット。時間経過ともに波の頂点が移動していっていることが分かる。


Scilabスクリプト


Scilabスクリプトは以下のようになりました。
clear;

// *** パラメータ ***
a = 1;
b = 1;
c = a * b;
lambda = 1;
f = 1;

// *** 位置に関する周期関数 ***
function phi_x = phi_x(x)
phi_x = a .* cos(2 * %pi / lambda .* x)
endfunction

// *** 時間に関する周期関数 ***
function phi_t = phi_t(t)
phi_t = b .* cos(2 * %pi * f .* t)
endfunction

// *** 定在波 ***
function phi = standing(x, t)
phi = phi_x(x) .* phi_t(t)
endfunction

// *** 進行波 ***
function phi = traveling(x, t)
phi = c * cos(2 * %pi * (1 / lambda .* x - f .* t))
endfunction

// *** 位置と時間に依存する波の計算 ***
X = linspace(0,1);
T = linspace(0,1);
[X2D, T2D] = ndgrid(X, T);
// 定在波の計算
PHI_standing = standing(X2D, T2D);
// 進行波の計算
PHI_travering = traveling(X2D, T2D);

// *** 定在波のプロット ***
scf(0);
xlabel("position, x");
ylabel("time, t");
zlabel("PHI(x,t)");
mesh(X2D, T2D, PHI_standing);
//surf(X2D, T2D, PHI_standing);

// *** 進行波のプロット ***
scf(1);
xlabel("position, x");
ylabel("time, t");
zlabel("PHI(x,t)");
mesh(X2D, T2D, PHI_travering);
//surf(X2D, T2D, PHI_travering);


関連エントリ




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 固体電子論  

Scilabで一次元井戸型ポテンシャルの分散関係

金持徹 編著 固体電子論の§2.5で紹介されている一次元井戸型ポテンシャルの$E-k$関係の図(図2.8)をScilabで計算してみました。

Wall100.png
Fig.1: $E-k$関係



一次元井戸型ポテンシャル


今回の一次元井戸型ポテンシャルのイメージとしては、太さが1原子分の針金です。長さ方向には、たくさん原子がありますが、太さ方向には隣り合う原子が存在しない状況です。針金の両端はどこにも接続されていないので、電子が外へ出ていくことができず、ポテンシャルが無限大になっています。

Scilabスクリプト


電子の計算を行うときには、数値が非常に小さいものになりがちです。それを避けるために、原子単位系を用います。AkaiKKRやecaljで慣れているので、今回はRydberg(リュードベリ)原子単位系を使いました。

clear;

// *** リュードベリ(Rydberg)原子単位系 ***
hbar = 1; // 換算プランク定数
h = 2 * %pi * hbar; // プランク定数
m = 0.5; // 電子の質量
a0 = 1; // ボーア半径

// *** 系の設定 ***
// 原子の数
//n_atom = 20;
n_atom = 100;
// 原子間距離
r0 = 1.0 * a0;
// 結晶の長さ
l = n_atom * r0;
// 価数と電子の総数
n_valence = 1;
n_electron = n_valence * n_atom;

// *** 分散関係の計算 ***
// 占有されている準位の計算
Kg = [];
Eg = [];
for g=1:int(n_electron / 2)
// 波数ベクトル (式 2.8)
kg = %pi / l * g;
Kg = [Kg, kg];
// エネルギー (式 2.18)
eg = h * h / (8 * m * l * l) * g * g;
Eg = [Eg, eg];
end
// フェルミ波数ベクトル (1/Bohr)
kf = kg
// フェルミエネルギー (Ry)
ef = eg

// *** グラフの描画 ***
plot(0, 0); // グラフ描画範囲のためのダミー
plot(Kg, Eg, ".b");

// *** グラフの装飾 ***
xgrid(color("gray"));
xlabel("Wave vector (1/Bohr)");
ylabel("Energy (Ry)");
//zoom_rect([0, 0, kf, ef]);


高圧や不純物の効果


上記のスクリプトでは、原子間隔をボーア半径の1倍、電子の個数を1原子当たり1個としています。これらの値を変更することで、高圧をかけた時にどの様に分散関係が変化するか、原子の一部を価数の異なるものに変えた時にどうなるかを確認できます。

関連エントリ




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 固体電子論 

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

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAOPアンプPICecalj状態密度常微分方程式モンテカルロ解析トランジスタodeDOSインターフェースPDS5022定電流スイッチング回路分散関係半導体確率論シェルスクリプトレベルシフト乱数HP6632Aバンド構造温度解析可変抵抗I2CR6452AブレッドボードPWscfバンドギャップトランジスタ技術セミナー数値積分反強磁性絶縁熱設計偏微分方程式非線形方程式ソルバフォトカプラシュミットトリガA/DコンバータLM358LED三端子レギュレータカオスISO-I2CGW近似順列・組み合わせマフィンティン半径発振回路74HC4053カレントミラーアナログスイッチPC817CTL431直流動作点解析USBサーボ数値微分固体電子論基本並進ベクトルQuantumESPRESSOスピン軌道相互作用相対論補間FFT開発環境チョッパアンプ単振り子電子負荷BSch量子力学アセンブラパラメトリック解析標準ロジック2ちゃんねるトレーナーバトルポケモンGOスーパーリーグbzqlty状態方程式コバルトLDAイジング模型VESTAブラべ格子六方最密充填構造位相図繰り返し熱伝導Maximaキュリー温度TLP621FETラプラス方程式抵抗失敗談SMPスイッチト・キャパシタgfortranスレーターポーリング曲線ewidth最適化Quantum_ESPRESSOVCAGGA仮想結晶近似ランダムウォークcygwin不規則合金QSGWダイヤモンドQNAPデータロガーマントルシュレディンガー方程式固有値問題自動計測条件分岐格子比熱井戸型ポテンシャルUPSFXA-7020ZRTLP521LM555NE555TLP552Writer509テスタMCU詰め回路三角波過渡解析ガイガー管OpenMPZnOフォノンハーフメタルubuntu最小値最大値熱力学xcrysdenCIF結晶磁気異方性ゼーベック係数スーパーセルUbuntu第一原理計算平均場近似起電力差し込みグラフfsolveブラウン運動フェルミ面awk疎行列縮退文字列不規則局所モーメント入出力化学反応クーロン散乱ヒストグラムRealforceマンデルブロ集合フラクタルキーボード三次元凡例HiLAPW両対数グラフ関数フィッティング熱拡散方程式陰解法片対数グラフグラフの分割シンボルGimp線種Crank-Nicolson法軸ラベル円周率MAS830L負帰還安定性ナイキスト線図EAGLEMBEAACircuitP-10フィルタノコギリ波CapSense2SC1815OPA2277PGALMC662PIC16F785TS-112TS-110等価回路モデルパラメータ・モデル日本語Excel直流解析CK1026トランス連立一次方程式トラックボールPC最小二乗法fcccif2cellPWgui状態図ハイパーリーグPvP擬ポテンシャル不純物問題SIC二相共存重積分電荷密度磁気モーメントEMTO-CPASPR-KKRBroydenTchebyshevルチル構造gnuplot面心立方構造pmixPbTeLmtARTFPLOAMULETsedOpenPBSウルツ鉱構造BaOinterp1初期値ウィグナーザイツ胞L10構造非線型方程式ソルバ正規分布等高線specx.fifortマテリアルデザインヒストグラム確率論ジバニャン方程式TeXFSMedeltquantumESPRESSOリジッドバンド模型スワップ領域岩塩構造デバイ模型半金属全エネルギー固定スピンモーメントc/amultiplot合金境界条件

最新コメント
リンク

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