Scilabでイジング模型 その2

計算物理学 応用編ising.cをScilabへ移植しました。

001_20141201124503785.png

Fig.1: 一次元のイジングモデル。コールドスタートから磁区が発達していく様子が再現されている。(n=100, kt=1.5, J=1, tmax=1000)


その結果、交換エネルギーJの符号によって強磁性状態や反強磁性状態が安定になることが確認できました。


一次元のイジング模型


一次元のリング状に原子(磁気双極子)が等間隔で並んでいる一次元結晶を考えます。それぞれの原子がもつ磁気双極子は、三次元空間では、さまざまな方向を向く可能性がありますが、一次元空間の場合は、方向が正と負しかないので、上向きと下向きの2種類のスピンのみを考えます。するとi番目の粒子がもつスピンは+1と-1のどちらかを取るとして si = ±1 と書くことができます。

リンク状の一次元結晶にN個の原子があるとすると、一次元結晶全体のスピンの配置は

j> = |s1, s2, ..., sN>
= {±1, ±1, ... ±1} (j = 1, 2, 3, ... 2N)

という量子状態ベクトルで表すことができます。

具体的に N=3 の場合に書き下してみます。

1> = {-1, -1, -1}
2> = {-1, -1, +1}
3> = {-1, +1, -1}
4> = {-1, +1, +1}
5> = {+1, -1, -1}
6> = {+1, -1, +1}
7> = {+1, +1, -1}
8> = {+1, +1, +1}

以上のようにN=3のときには23=8種類の状態を取りうることがわかります。
イジング模型では、隣り合う原子だけが相互作用すると考えます。相互作用のパラメータをJとすると、外部磁場がない場合、ある状態αjにおけるエネルギーは

E( \alpha_j ) = - J \sum_{i=1}^{N-1}s_i s_j


この式をすべてのαjに関して計算すればよいことになります。とは言うものの、原子の数がN=3ならば状態の種類が8種類しか存在しないので、すべてを計算することも可能でしょうが、N=100などになってしまうと2100≒1.26*1030となり、とてもすべてを計算するのは現実的ではありません。このため、数値シミュレーションでは乱数を使って計算を行います。

メトロポリスのアルゴリズム


詳しいアルゴリズムの説明は計算物理学 応用編の通りですが、その中のising.cまたはising.fをScilabに移植します。

clear;

// *** 定数の設定 ***
n = 100; // 粒子の数
kt = 1.5; // 温度
J = 1; // 交換エネルギー (1: 強磁性, -1:反強磁性)
rand("uniform"); // 乱数は一様乱数とする
tmax = 1000; // 時間の最大ステップ

// *** 初期化 ***
// 各粒子におけるスピン
spin = ones(1,n); // コールドスタート
//spin = 1 - 2 * round(rand(1,n)); // 乱数スタート

// *** エネルギーの計算関数 ***
function e = energy(spin)
e = - J * sum(spin .* [spin(2:n), spin(1)]);
endfunction

SPIN = [];

// *** 時間発展 ***
for t = 1:tmax do
oldenergy = energy(spin);
element = ceil(n * rand()); // 粒子を一つ選ぶ
spin(element) = -1 * spin(element); // スピンを反転
newenergy = energy(spin);
spin(element) = (- 2 * ((newenergy > oldenergy) & (exp((- newenergy + oldenergy) / kt) <= rand())) + 1) * spin(element);
// 時刻tにおけるスピンを保存
SPIN = [SPIN;spin];
end

// *** スピンの時間変化をプロット ***
Matplot((SPIN' + 1) .* 10);
zoom_rect([0,0,tmax,n]);
xlabel("Time");
ylabel("Position");


反強磁性状態の計算


交換エネルギーを負にとると、反強磁性状態になります。

002_2014120112573532a.png
Fig.2: 反強磁性状態


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 確率論 乱数 イジング模型 モンテカルロ解析 強磁性 反強磁性 

Scilabでイジング模型 その1

微分方程式による物理現象のモデル化(PDF)の熱統計力学の章にある一次元イジングモデルのOctaveのプログラムをScilabへ移植しました。

001_20141130214554d90.png

Fig.1: 一次元のイジングモデル



Scilabでモンテカルロシミュレーション


常微分方程式タグを付けたエントリでは、微分方程式による物理現象のモデル化(PDF)で紹介されているOctaveによる常微分方程式のプログラムをScilabへ移植するという事を行ってきました。Scilabで乱数の生成からのいくつかのエントリでは、微分方程式による物理現象のモデル化(PDF)の熱統計力学の章に入るに際して、Scilabで楽しむ確率論(PDF)Scilabで乱数を扱う方法に触れました。

今回は、また微分方程式による物理現象のモデル化(PDF)に戻り、イジング模型のScilabシミュレーションを行います。
ただし、今回は微分方程式による物理現象のモデル化(PDF)のスクリプトを単純にScilabで実行できるように書き直しただけにしておきます。正直なところ微分方程式による物理現象のモデル化(PDF)の内容だけでは、何の計算をしているのかを理解することは困難です。おそらく元ネタと思われる計算物理学 応用編の併読が必須と思われます。

clear;

N = 200; // 粒子数
J = 1; // 交換相互作用
H = 0.0; // 外部磁界
im = [0:N - 1];
im(1) = N;
ip = [2:N + 1];
ip(N) = 1;
Neq = 5 * N;
Samp = 200;
DIV = 19;

rand("uniform"); // 乱数は一様乱数とする
spin = ones(1,N); // スピン
T = linspace(0.5,5,DIV); // 温度

// *** モンテカルロ計算 ***
for kt = 1:DIV do
// エネルギー
e1 = 0;
e2 = 0;
mag1 = 0;
mag2 = 0;
for ks = 1:Samp do
for kq = 1:Neq do
is = ceil(N * rand());
de = J * 2 * spin(is) * (spin(im(is)) + spin(ip(is)));
if de > 0 & exp(- de / T(kt)) < rand() then
spin(is) = 1 * spin(is);
else
spin(is) = -1 * spin(is);
end
end
mag1 = mag1 + sum(spin);
mag2 = mag2 + sum(spin) ^ 2;
e1 = e1 - J * (spin * [spin(2:$),spin(1)]');
e2 = e2 + (J * (spin * [spin(2:$),spin(1)]')) ^ 2;
end
M(kt) = mag1 / Samp;
M2(kt) = mag2 / Samp;
E(kt) = e1 / Samp;
E2(kt) = e2 / Samp;
X(kt) = (M2(kt) - M(kt) .^ 2) ./ T(kt) / N;
C(kt) = (E2(kt) - E(kt) .^ 2) ./ T(kt) .^ 2 / N;
end
// 磁気感受率の逆数
RX = 1 ./ X;

// *** 厳密解の計算 ***
// 温度ベクトル
Ta = linspace(0.1,5,50);
// 粒子1個あたりの平均エネルギー
Ea = - tanh(J ./ Ta);
// 比熱
Ca = (J ./ Ta) .^ 2 ./ cosh(J ./ Ta) .^ 2;
// 磁化
Ma = sinh(H ./ Ta) ./ sqrt(sinh(H ./ Ta) .^ 2 + exp(-4 * J ./ Ta));
// 磁気感受率
Xa = exp(2 * J ./ Ta) ./ Ta;
RXa = 1 ./ Xa;

// *** グラフのプロット ***
// エネルギー
subplot(2,2,1);
plot(T,E / N,'or');
plot(Ta, Ea, '--g');
xlabel("kT/J");
ylabel("E/NJ");
// 比熱
subplot(2,2,2);
plot(T,C,'or');
plot(Ta,Ca,'--g');
xlabel("kT/J");
ylabel("C/Nk");
// 磁化
subplot(2,2,3);
plot(T,M/N,'or');
plot(Ta,Ma/N,'--g');
xlabel("kT/J");
ylabel("M/N");
// 磁気感受率の逆数
subplot(2,2,4);
plot(T,RX,'or');
plot(Ta,RXa,'--g');
xlabel("kT/J");
ylabel("N/JX");


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 確率論 乱数 イジング模型 モンテカルロ解析 強磁性 

質量パーセントと原子パーセント

A100-xBx二元系合金の質量パーセント(wt.%)と原子パーセント(at.%)を相互に換算すると以下のようになります。
ただしAとBの原子量はそれぞれMAとMBとします。
金属Aに金属Bが原子パーセントで x (at.%)固溶しているときの金属Bの質量パーセント y (wt.%)は以下の式で表すことができます。

y = \frac{x M_{B}}{(100-x)M_{A}+xM_{B}} \times 100

逆に金属Aに金属Bが質量パーセントで y (wt.%)固溶しているときの金属Bの原子パーセント x (at.%)は以下の式で表すことができます。

x = \frac{y M_{A}}{(100-y)M_{B}+yM_{A}} \times 100

簡単な話なのですが、しょっちゅう忘れるのでエントリにまとめておきます。


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


AB二元系合金で、Bの濃度が原子パーセントで x (at.%)のとき、これを質量パーセント y (wt.%)へ変換することを考えます。質量パーセントの定義は (Bの質量)÷(全質量) なので、それぞれの項を計算するために、原子パーセント x (at.%)に原子量を掛け算して質量にすればよい。

A100-xBxの1化学式あたりのBの質量は x * MBなので全質量でわると

y = \frac{x M_{B}}{(100-x)M_{A}+xM_{B}} \times 100

となります。

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


同様に質量パーセント y (wt.%)から原子パーセント x (at.%)への変換も考えます。
もちろん前述の式をxについて解きなおせばよいだけですが、逆の変換のときと同様に定義からいきます。
原子パーセントの定義は (Bの原子数)÷(全原子数) なので、それぞれの項を計算します。(Bの原子数)=(Bの質量)÷(Bの原子量)なので、単位質量あたりのBの原子数は y / MB なので全原子数で割ると

x = \frac{\frac{y}{M_{B}}}{\frac{(100-y)}{M_{A}}+\frac{y}{M_{B}}} \times 100 = \frac{y M_{A}}{(100-y)M_{B}+yM_{A}} \times 100

金銅合金の場合


試しに金と銅の合金 Au100-xCuxについて計算してみます。
ここでそれぞれの原子量は MAu = 196.96657, MCu = 63.546 です。

001_2014113014415168a.png

002_20141130144151953.png

Fig.1-2: 金銅合金の質量パーセントと原子パーセントの相互変換


Scilabのコードは、それぞれ以下のようになりました。

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

clear;

// *** 原子量 ***
Mcu = 63.546;
Mau = 196.96657;

// *** 濃度 ***
// 原子パーセント(at.%)
x = linspace(0,100);

// *** 関数の定義 ***
function y = at2wt(x,Ma,Mb)
y = 100 * x .* Mb ./ ((100 - x) .* Ma + x .* Mb)
endfunction

// *** グラフのプロット ***
plot(x,at2wt(x,Mau,Mcu));
xgrid(color("gray"));
xlabel("Atomic % Cu");
ylabel("Weight % Cu");


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

clear;

// *** 原子量 ***
Mcu = 63.546;
Mau = 196.96657;

// *** 濃度 ***
// 質量パーセント(wt.%)
y = linspace(0,100);

// *** 関数の定義 ***
function x = wt2at(x,Ma,Mb)
x = 100 * y .* Ma ./ ((100-y) .* Mb + y .* Ma)
endfunction

// *** グラフのプロット ***
plot(y,wt2at(y,Mau,Mcu));
xgrid(color("gray"));
xlabel("Weight % Cu");
ylabel("Atomic % Cu");


付録


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


フィードバック



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

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


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


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

tag: Scilab 

AkaiKKRで銅と銅亜鉛合金のフェルミ面

AkaiKKR(Machikaneyama)を用いると状態密度やバンド構造が簡単に描画できます。いくつかの第一原理計算パッケージには、これらに加えてフェルミ面の描画機能があります。しかしながら、AkaiKKRにはおそらくその機能はありません。そこで力任せにBlochスペクトル関数を計算してCuとCu70Zn30のフェルミ面の断面図を作成しました。

その結果は金属電子論〈上〉で紹介されている実験結果をよく再現しました。


金属のフェルミ面


金属の電子構造の個性を表現するために「状態密度(DOS)」や「エネルギー分散(バンド構造)」などがよく図としてプロットされます。「エネルギー分散」は電子の持っている波数とエネルギーの関係をプロットしたもので「状態密度」はエネルギーとそのエネルギーを持つ電子の個数をプロットしたものです。これらは共に、電子の状態をエネルギーの関数として表現しています。

しかしながら、実際の金属の物性は、その多くがフェルミエネルギーの電子の性質だけで決まります。そこで、フェルミ準位だけに限って波数ベクトルを表示した「フェルミ面」も金属の電子状態を表現するために利用されます。

純金属のフェルミ面は、非常によく研究されており、その一覧はウエブ上のデータベースでも見ることができます。(参考: Fermi Surface ExplorerThe Fermi Surface Database)

Cu.jpg
Fig.1: 銅のフェルミ面


Fig.1に示したのはThe Fermi Surface Databaseから引用した銅のフェルミ面です。銅のフェルミ面は、ほとんど自由電子的な球に近い形状をしています。しかしながらブルリアンゾーンのL点の周囲でフェルミ面がブリルアンゾーンに接触しています。

AkaiKKRでフェルミ面の描画


AkaiKKR(Machikaneyama)には標準ではフェルミ面を描画する機能は、ありません(多分)。そもそも不規則合金では、純金属と異なり、フェルミ面という概念自体が必ずしも妥当なものではなくなります。これは電子のエネルギー分散(バンド構造)を考えた際に不規則合金では、各バンドがにじんでしまったことと同じです。(参考: 密度汎関数法の発展 -マテリアルデザインへの応用)

バンド分散に関しては、Blochスペクトル関数を各k点に対して計算したものをプロットすることで表現できました。
同様にして、力まかせに片っ端からBlochスペクトル関数を計算し、フェルミエネルギー(近く)のものだけプロットするという方法でフェルミ面の断面の描画を行う事を考えます。

今回は純金属であるCuと、不規則合金でありながら比較的きれいにフェルミ面の残るCu70Zn30合金のフェルミ面の断面を描いてみることにします。これらの金属のフェルミ面の断面に関する実験的結果は金属電子論〈上〉に紹介されています。

シェルスクリプト


詳しい説明は、別のエントリにで行いたいと思いますが、おおよそ次のような事を行うシェルスクリプトを作成します。

  1. 波数空間上のベクトルPA, PBを2辺とする平行四辺形を考える
  2. 2つのベクトルをそれぞれn分割、m分割した空間メッシュを作成する
  3. spc計算用の入力ファイルのテンプレートから上記のk点を加えた入力ファイルを作成しspecxを起動
  4. 計算結果のspcファイル群からGNUPLOTに適したdatファイルの作成


結果


CuのXΓX断面とKΓX断面を計算したものがFig.2-3です。

XGXCu.png
Fig.2: Cuのフェルミ面のXΓX断面

XGKCu.png
Fig.3: Cuのフェルミ面のKΓX断面


XΓX断面におけるフェルミ面は、自由電子的な真円に近い形状をしています。黒のラインで示したのが第一ブリルアンゾーンの境界で、図中に書き込んではありませんが、斜めになっている角がW点、斜面の真ん中がL点です。

KΓX断面におけるフェルミ面は、前述したとおり円形からひずみ、第一ブルリアンゾーンの境界に触れています。触れている中心がL点で斜めの部分の上の角がU点です。

同様にCu70Zn30の断面図をプロットしたのがFig.4-5です。

XGXCuZn.png
Fig.4: Cu70Zn30のフェルミ面のXΓX断面

XGKCuZn.png
Fig.5: Cu70Zn30のフェルミ面のKΓX断面


金属電子論〈上〉で解説されている通り、Znの合金化によって1原子あたりの平均価電子数e/aが増加します。この結果として、フェルミ面が大きくなっていることがわかります。このことはXΓX断面では円の半径が大きくなっていること、KΓX断面では第一ブリルアンゾーンに触れている領域が広がっていることとして表れています。

関連エントリ






参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR CPA フェルミ面 

AkaiKKRでテスト計算

AkaiKKRのインストールにてAkaiKKR(machikaneyama)のインストールが完了したので、今回はテスト計算として強磁性体心立方構造の鉄のセルフコンシステント計算とそのポテンシャルファイルを利用した状態密度の計算方法について書きました。


AkaiKKRの実行


AkaiKKRのインストールによってcygwinまたはubuntuにAkaiKKRのインストールができたはずなので、次はテスト計算を行います。

インストールが完了してさえいれば、テスト計算の方法はKKR-Green関数法によるバンド計算計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)など色々な文献に書いてある通りなので、簡単なはずです。

AkaiKKRの実行ファイルはspecxです。
specx < in/infile.in > out/outfile.out &
AkaiKKRのインストールの通り実行ファイルの存在するディレクトリにパスが通っていれば、上記のようにして実行できます。

今回は、下記に示すbcc構造強磁性鉄の入力ファイルを例にします。
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 mjw mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 3 50 0.023
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Fe 1 1 0.0 2
26 100
c------------------------------------------------------------
c natm
1
c------------------------------------------------------------
c atmicx(in the unit of a) atmtyp
0 0 0 Fe
c
c------ The following types of inputs and their combination
c are also allowed.
c------ In those cases a, b, b mean primitive unit vectors
c and x, y, z mean shortest unit cell vectors along
c x, y, and z axses.
c 1/2 1/2 1/2 Fe
c 0a 0b 0c Fe
c 0x 0y 0z Fe
c 1/2a 1/2b 1/2c Fe
c------------------------------------------------------------

上記のファイルは、このエントリを書いている時点での最新版(cpa2002v009c September 30, 2014)の場合inフォルダにfeという名前で保存されています。このファイルを(念のため)コピーしてfe.inを作成し、specxに渡します。

cp in/fe in/fe.in
specx < in/fe.in > out/fe.out &

出力ファイル


上記の計算を行うと合計で3つのファイルが作成されます。outフォルダのfe.out、dataフォルダのfeとfe.infoです。dataフォルダのfeがメインの計算結果であるポテンシャルファイルです。しかしポテンシャルファイルはバイナリファイルなので人間には読めません。人間にわかりやすい形式の出力ファイルはoutフォルダのfe.outとdataフォルダのfe.infoです。

out/fe.outの最初の方に計算の設定に関する内容が書き込まれます。設定に問題がなければ、以下に示すようなイテレーションのループに入ります。

   ***** self-consistent iteration starts *****
Fe
itr= 1 neu= -2.3732 moment= 0.1175 te= -2522.45174794 err= 0.633
itr= 2 neu= -1.8783 moment= 0.1674 te= -2523.05786722 err= 0.455
itr= 3 neu= 2.0397 moment= -0.0543 te= -2520.36515324 err= 0.510
itr= 4 neu= 2.0695 moment= 0.5212 te= -2520.72070737 err= 0.468
itr= 5 neu= 0.1956 moment= 0.2876 te= -2522.77976600 err= -0.346
itr= 6 neu= -1.4162 moment= 0.3321 te= -2523.09445307 err= 0.159
itr= 7 neu= -0.8089 moment= 0.9042 te= -2522.80324234 err= -0.427
itr= 8 neu= -0.2096 moment= 0.9756 te= -2522.79459829 err= -0.437
itr= 9 neu= 0.2350 moment= 1.4317 te= -2522.80324021 err= -0.390
(中略)
itr= 45 neu= -0.0001 moment= 2.1755 te= -2522.82229189 err= -2.591
itr= 46 neu= -0.0001 moment= 2.1755 te= -2522.82229254 err= -2.658
itr= 47 neu= -0.0002 moment= 2.1755 te= -2522.82229304 err= -2.732
itr= 48 neu= -0.0002 moment= 2.1755 te= -2522.82229340 err= -2.811
itr= 49 neu= -0.0002 moment= 2.1755 te= -2522.82229364 err= -2.895
itr= 50 neu= -0.0002 moment= 2.1755 te= -2522.82229381 err= -2.982
*** no convergence


イテレーションはerrが-6より小さくなるか、イテレーションの回数が入力ファイルのmaxitrに達するかするまで続けられます。上記の例では、イテレーションの回数が入力ファイルで指定されている50回に達したので収束なし(no convergence)として計算を終了しています。

入力ファイルのrecordを2ndにしておけば、もういちどspecxを実行することにより、続きから計算できます。
specx < in/fe.in > out/fe.out &

逆に、前回までのポテンシャルファイルを無かったことにして最初から計算したい場合は、recordをinitにする、或いは、dataフォルダにあるポテンシャルファイルfeを削除するなどしてから再実行します。

   ***** self-consistent iteration starts *****
Fe
itr= 1 neu= -0.0002 moment= 2.1755 te= -2522.82229381 err= -2.982
itr= 2 neu= -0.0002 moment= 2.1755 te= -2522.82229384 err= -2.998
itr= 3 neu= -0.0002 moment= 2.1755 te= -2522.82229388 err= -3.025
itr= 4 neu= -0.0001 moment= 2.1755 te= -2522.82229394 err= -3.054
itr= 5 neu= -0.0001 moment= 2.1755 te= -2522.82229402 err= -3.089
(中略)
itr= 45 neu= -0.0000 moment= 2.1755 te= -2522.82229503 err= -5.611
itr= 46 neu= -0.0000 moment= 2.1755 te= -2522.82229503 err= -5.696
itr= 47 neu= -0.0000 moment= 2.1755 te= -2522.82229503 err= -5.782
itr= 48 neu= -0.0000 moment= 2.1755 te= -2522.82229503 err= -5.868
itr= 49 neu= -0.0000 moment= 2.1755 te= -2522.82229503 err= -5.956
itr= 50 neu= -0.0000 moment= 2.1755 te= -2522.82229503 err= -6.044


収束したのを確認したら、次は状態密度を表示してみます。

状態密度の計算


作成したポテンシャルファイルから状態密度を計算することができます。
このためには、入力ファイルのgoの部分をdosに変更します。
ついでにbzqltyを大きくするときれいな状態密度を描くことができます。

ただしbzqltyを大きくし過ぎると ***err in bzmesh...too many k-pointのエラーが出てしまいます。このエラーを回避するためにはsource/specx.fのnk1xの値をk点の数(nk)よりも大きくしてmakeすれば良いのですが、bzqltyとnkの対応関係は実際に計算してみないとわかりません。

とりあえず今回はbzqlty=12としました。これはnk=413に対応します。cpa2002v009c (September 30, 2014)ではデフォルトはnk1x=500となっています。

結局、以下の内容のファイルをfedos.inとしてinフォルダに保存します。

c----------------------Fe------------------------------------
dos 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 mjw mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 12 50 0.023
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Fe 1 1 0.0 2
26 100
c------------------------------------------------------------
c natm
1
c------------------------------------------------------------
c atmicx(in the unit of a) atmtyp
0 0 0 Fe
c
c------ The following types of inputs and their combination
c are also allowed.
c------ In those cases a, b, b mean primitive unit vectors
c and x, y, z mean shortest unit cell vectors along
c x, y, and z axses.
c 1/2 1/2 1/2 Fe
c 0a 0b 0c Fe
c 0x 0y 0z Fe
c 1/2a 1/2b 1/2c Fe
c------------------------------------------------------------

specx < in/fedos.in > out/fedos.out &
状態密度の値は、出力ファイルのtotal DOSからintegrated DOSの間に出力されます。これをグラフソフトでプロットします。

001_20141226050748757.png
Fig.1: Excelでプロットした状態密度


スピン分極を考慮した計算では、全状態密度が二組あるので正と負にプロットしました。

最後にdataフォルダのfe.infoについて。

       5.2700       -2522.8222938   2.17548
5.2700 -2522.8222950 2.17547
5.2700 -2517.8475965 0.30790


左の列から、格子定数、全エネルギー、磁気モーメントです。
specxを実行するごとに1行ずつ末尾に値が追加されていきます。上記の例でいえば、1行目が1回目のgo計算でこの段階ではまだ収束していません。2行目は2回目のgo計算で完全に収束した結果です。3行目は状態密度(dos)の計算のときの出力で、この3行目の値は意味のない値です。
したがって上記の例でいえば2行目だけが信用できる値です。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR 強磁性 

AkaiKKRの使い方

AkaiKKR(machikayema)は、東京大学の赤井久純先生が開発された第一原理計算パッケージです。結晶の電子状態を計算することのできる第一原理計算パッケージは、多くの種類が公開されていますが、AkaiKKRはKKR-CPA法を採用することで化学的な不規則(不純物、不規則合金、固溶体)のバンド計算が可能であるというのが最大の特徴です。また、入力ファイルも単純であり、スペックの低い旧式のwindowsノートパソコンでも動作するといったところも優れています。

今回のエントリでは、これまでに書いたAkaiKKR関連のエントリを整理してリンクします。


1. はじめに

この章では、AkaiKKRのインストール方法とテスト計算、及び、使い方の最初の一歩についてまとめます。AkaiKKRは、公式なマニュアルであるKKR-Green関数法によるバンド計算Machikaneyama2000の使用に関するメモ及び、計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)密度汎関数法の発展 -マテリアルデザインへの応用を用いて独習することも可能ですが、やはり年に二回開催されているCMD®ワークショップに参加するのが最大の近道です。

2. 基本的な使い方

第一原理計算の出力結果として基本的なものは、全エネルギー、状態密度、バンド分散です。第2章では、面心立方構造(fcc)、体心立方構造(bcc)、六方最密充填構造(hcp)、ダイヤモンド構造、そして立方晶ペロフスカイト構造に関して状態密度とブロッホスペクトル関数(バンド分散)を出力する入力ファイルの例を示します。

3. 応用的な計算

第2章が基礎編とするならば、第3章は応用編となります。例えば平衡格子定数を計算したい場合は、格子定数を変化させながら、全エネルギーが最小となる点を探します。

4. 設定パラメータ

第一原理計算は、格子定数と原子番号以外には入力パラメータが必要ないということになっていますが、実際には計算の精度に関するパラメータなどを指定しなければなりません。この章ではパラメータの意味と決定方法について書きます。

5. その他

6. 参考URL

6.1 AkaiKKR関連

6.2 第一原理計算全般

7. 参考書籍


8. データベース

9. あとがきと更新履歴

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


フィードバック

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

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

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

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

tag: AkaiKKR machikaneyama KKR CPA  

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

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

最新コメント
リンク

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