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 固体電子論 

Scilabでマクスウェルの速度分布則

金持徹 編著 固体電子論の§2.2ガス模型で紹介されているマクスウェルの速度分布をScilabで計算してみました。

N2.png

Fig.1: N2分子のマクスウェルの速度分布。赤が-187℃、青が20℃。



マクスウェルの速度分布


マクスウェルの速度分布則によると気体分子全体の速度分布は、以下の式で表すことができます。
\[
f(x) = \frac{4}{\sqrt{\pi}}x^2 \exp (- x^2)
\]
ここで $x = V / V_m$ は、規格化された粒子の速度です。$V_m$ は「最大確率速度」と呼ばれる代表的な速度で、以下のように温度 $T$ と 粒子の質量 $m$ から計算することができます。
\[
V_m = \left(\frac{2 k_B T}{m} \right)^{1/2}
\]
ここで $k_B$ はボルツマン定数です。

Scilabスクリプト


マクスウェルの速度分布則の式の横軸は $V_m$ で規格化されています。グラフで外形を確認するだけなら $0 \leqq x \leqq 3$ 位を計算すればよいでしょう。縦軸は $x$ を無限大まで積分すると面積が1になるように規格化されています。したがって、グラフを描画するときには、横軸を $V_m$ 倍して、縦軸を $V_m$ で割れば良いことになります。

clear;

// *** 物理定数 ***
// ボルツマン定数 (m^2 kg s^-2 K^-1)
kB = 1.38064852E-23;
// アボガドロ定数 (/mol)
na = 6.02214076E23;
// 窒素の原子量 (g/mol)
mass_N = 14.0067;
// 粒子の質量 (kg)
mass = 2 * mass_N / na * 1E-3; // N2 の質量 (kg)
//mass = 9.10938356E-31; // 電子の質量 (kg)

// *** Maxwellの速度分布則 ***
// 無次元化速度
X = linspace(0, 3);
// Maxwellの速度分布則
N = 4 ./ sqrt(%pi) .* (X .* X) .* exp(- X .* X);

// *** 最大確率速度の計算 ***
// 温度 (K)
t = -187 + 273.15; // 低温 -187 ℃
// 最大確率速度
vm = sqrt(2 * kB * t / mass)
// グラフの描画
plot(vm * X, (1 / vm) * N, "-b");

// 温度 (K)
t = 20 + 273.15; // 常温 20 ℃
// 最大確率速度
vm = sqrt(2 * kB * t / mass)
// グラフの描画
plot(vm * X, (1 / vm) * N, "-r");

// *** グラフの装飾 ***
xlabel("Speed (m/s)");
ylabel("Number of particle");
xgrid(color("gray"));


ガス模型による電子の速度


粒子の質量 $m$ を、窒素分子の質量から、電子の質量へ変更することで、ガス模型から計算した電子の速度分布を計算できます。

electron.png
Fig.2: 電子のマクスウェルの速度分布。赤が-187℃、青が20℃。


ところが実験によれば, 実際の金属結晶内では, 0 Kに近づいても自由電子は活発に飛び回っているものとみられ, このガスモデルの示す状態とはひどく違っているのである (§3.6参照)


実際、金属は温度を下げていくほど、電気抵抗が減っていきます。この温度依存性は、温度を下げるほど電子の速度が下がっていくというガス模型の結果とは一致しません。

誤植?


ところでP17の最後の段落に以下の様に書いてあります。
すなわち, 固体結晶内で同一方向に同一速度で飛ぶ電子は, 電子の自転を考慮しても2個以上は存在しないと考える.


これは
すなわち, 固体結晶内で同一方向に同一速度で飛ぶ電子は, 電子の自転を考慮しても2個以上はまでしか存在しないと考える.

の誤植じゃないでしょうか…?

付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 固体電子論 

AkaiKKRのpmix

AkaiKKRの入力ファイルにおけるpmixは、セルフコンシステント計算における収束の速さを決めるパラメータです。大きい値を使うほど速く解に近づきますが、振動して収束しない可能性が出てきます。逆に小さい値を使うほど収束しやすくなりますが、収束に時間がかかるようになります。(たぶん、小さい値を使う方が準安定な解に収束しやすくもなる?)

最近(といっても1年ぐらい前ですが)、収束させる方法としてBroyden法が実装されました。これに伴って、今まで通りpmixに数値だけを与えた場合のデフォルトの動作は、最初に、以前から使われていたTchebyshev法を使い、適当なところでBroyden法に切り替わる挙動になったとの事です。

以前のようにTchebyshev法だけを使いたい場合は pmix のところに数値だけではなく 0.035tch の様に書きます。逆にBroyden法だけをつかいたい場合は 0.035bry と書けば良い様です。

普通は特に何も考えずに数値だけを書いておいて、両方の方法で収束させるのが良いと思います。


tch → bry で収束させる


AkaiKKRに付属しているbcc Feの入力ファイルをそのまま実行するとTchebyshev法とBroyden法の両方が利用されます。

c----------------------Fe------------------------------------
go data/fe
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
bcc 5.27 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 1.0 nrl mjwasa mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 8 100 0.035
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Fe 1 0 0.0 2 26 100
c------------------------------------------------------------
c natm
1
c------------------------------------------------------------
c atmicx(in the unit of a) atmtyp
0 0 0 Fe
c------------------------------------------------------------


pmixに単純に0.035と入力した場合、以下のような出力が得られます。
   maxitr=100  pmix= 0.03500  mixtyp=tchb-brydn

mixtyptchb-brydn となっていることが分かります。

tch のみで収束させる


他の部分は同じなので、変更した pmix の部分だけを示します。今回はTchebyshev法だけを使う例です。
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 8 100 0.035tch
c------------------------------------------------------------

出力は以下のようになります。
   maxitr=100  pmix= 0.03500  mixtyp=tchebyshev

確かに mixtyptchebyshev となっています。

bry のみで収束させる


Broyden法だけを使う場合は、以下のようになります。
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 8 100 0.035bry
c------------------------------------------------------------

出力は以下の通りです。
   maxitr=100  pmix= 0.35000  mixtyp=broyden

mixtypbroyden となっている事の他に pmix の値が10倍になっています。source/spmain.f を見てみると mixtypbroyden の場合は、値を10倍にするような処理が入っているようです。なので、たぶんこれは正しい挙動なのだと思います…。

比較と使い分け


今回の bcc Fe の計算では、どの方法を使っても計算は収束しました。収束に必要だったiterationの回数は tch+bry で35回、tch のみなら57回、 bry のみなら43回でした。全エネルギーも以下の通り、ほとんど同じ値になりました。
total energy= -2522.835012453 (tch+bry)
total energy= -2522.835012462 (tch)
total energy= -2522.835012456 (bry)

第37回CMDワークショップにて、これらの手法の使い分けについて赤井先生に尋ねました。AkaiKKRでは結晶の計算の前に、原子の計算を行うのですが、 bry だけでは原子の計算に失敗することがあり、この場合は tch を使う必要があり、このため tch+bry がデフォルトの動作になっているとの事。また、デフォルトの tch+bry で上手く収束させられない場合は tch だけにする(つまり以前の挙動に戻す)と収束しやすくなることがあるかもしれないとも。ただし、上手く収束させられないときは ewidth とか pmix の値そのものとかを変えてみるのが先決という感じでした。

関連エントリ




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama pmix Tchebyshev Broyden 

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

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

最新コメント
リンク

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