Scilabでブラウン運動 その3

Scilabでブラウン運動 その2では、単位時間あたりに任意の方向に距離1だけ移動するランダムウォークのスクリプトを書きました。これはScilabでブラウン運動 その1では移動できる角度が制限されていたことに対する拡張です。

今回は、さらに発展させて、単位時間あたりに移動する距離もランダムに決定するように拡張を行います。

001_20141109111602617.png
Fig.1: 一次元のランダムウォーク。各ステップでの移動距離を1に固定したもの(青)と移動距離の平均二乗変位平方根が1となるような正規分布から計算したもの(赤)。赤線は、移動距離が1よりも短いこともあれば長いこともある。



一次元のランダムウォーク


二次元から始めてもよいのですが、移動距離を固定したものとの比較が簡単なため、一次元のスクリプトも示します。
各ステップにおける移動量を乱数で与える場合、どのような乱数を使うかの選択肢がいくつかあると思いますが、今回は正規分布(ガウス関数)としました。Scilabにおいてはgrandを利用することによって正規分布に従う乱数を生成する事ができます。

clear;

// *** 計算の設定 ***
u = 1; // 1ステップの間の平均二乗変位平方根
av = 0;
tnum = 20; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
// 移動距離が正規分布に従う場合
xigauss = grand(tnum,1,'nor',av,u);
Sgauss = [0; cumsum(xigauss)];
// 移動距離が一定の場合
xifix = 2 * u * (rand(tnum,1) >= 0.5) - u;
Sfix = [0; cumsum(xifix)];

// *** グラフのプロット ***
plot(t,Sgauss,'-or');
plot(t,Sfix,'-ob');
xlabel("x position");
ylabel("y position");


二次元のランダムウォーク


二次元への拡張を行います。ソースコードを比較するなら今回の一次元のものと比較するよりもScilabでブラウン運動 その2の二次元のものと比較する方がわかりやすいです。ほとんど r=1 を r=abs(grand(1,tnum,'nor',av,u)) にしただけなので。

002_20141109111601aa4.png
Fig.2: 二次元のランダムウォーク


clear;

// *** 計算の設定 ***
u = 1; // 1ステップの間の平均二乗変位平方根
av = 0;
tnum = 10000; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
// 移動する距離
r = abs(grand(1,tnum,'nor',av,u));
// 移動する方向
theta = 2 * %pi * rand(1,tnum);
// 移動量
xi = r .* cos(theta);
yi = r .* sin(theta);
// 最終的な座標
S = [zeros(2,1), [cumsum(xi,'c'); cumsum(yi,'c')]];

// *** グラフのプロット ***
plot(S(1,:),S(2,:),'-b');
xlabel("x position");
ylabel("y position");


三次元のランダムウォーク


三次元への拡張もこれまでと全く同様です。

003_20141109111601d60.png
Fig.3: 三次元のランダムウォーク


clear;

// *** 計算の設定 ***
u = 1; // 1ステップの間の平均二乗変位平方根
av = 0;
tnum = 10000; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
// 移動する距離
r = abs(grand(1,tnum,'nor',av,u));
// 移動する方向
theta = 2 * %pi * rand(1,tnum);
phi = 2 * %pi * rand(1,tnum);
// 移動量
xi = r .* sin(theta) .* cos(phi);
yi = r .* sin(theta) .* sin(phi);
zi = r .* cos(phi);
// 最終的な座標
S = [zeros(3,1), [cumsum(xi,'c'); cumsum(yi,'c'); cumsum(zi,'c')]];

// *** グラフのプロット ***
param3d(S(1,:),S(2,:),S(3,:));
xlabel("x position");
ylabel("y position");
zlabel("z position");


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 乱数 モンテカルロ解析 確率論 ランダムウォーク 正規分布 

Scilabでブラウン運動 その2

Scilabでブラウン運動 その1では二次元のランダムウォークの際に、斜めの方向にしか動けないようなモデル化を行っています。

002_20141109070543ce1.png

Fig.1: Scilabでブラウン運動 その1での二次元のシミュレーションの結果。各ステップにおいては、粒子は斜めの方向にしか移動できない。



二次元への拡張への別アプローチ


Scilabでブラウン運動 その1では、Scilabで楽しむ確率論(PDF)の対称ランダムウォークを三次元まで拡張しました。その際の方針は、Scilabで楽しむ確率論(PDF)で一次元から二次元に拡張した際のものを踏襲しました。

しかしながら一次元のランダムウォークを『単位時間(1ステップ)あたりに絶対値で1だけ移動するが、その方向はランダムである』と解釈するならば、二次元に拡張する際に別の考え方をしなければならなさそうです。

そこで今回はx軸からの角度θが一様分布に従った乱数で与えられるようにして、単位時間後に半径1の円周上のどこかに移動している場合のシミュレーションを行います。移動量はそれぞれ以下のように、極座標で計算することができます。

xi = r cosθ
yi = r sinθ

001_20141109095932f86.png
Fig.2: 二次元のランダムウォーク。任意の角度θへ移動できるようにしたバージョン。


Scilabスクリプトは以下のようになります。

clear;

// *** 計算の設定 ***
r = 1; // 1ステップの間に移動する距離
tnum = 10000; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
theta = 2 * %pi * rand(1,tnum);
xi = r .* cos(theta);
yi = r .* sin(theta);

S = [zeros(2,1), [cumsum(xi,'c'); cumsum(yi,'c')]];

// *** グラフのプロット ***
plot(S(1,:),S(2,:),'-b');
xlabel("x position");
ylabel("y position");


三次元への拡張


三次元への拡張も問題ないと思います。

xi = r sinθcosφ
yi = r sinθsinφ
zi = r cosφ

002_2014110909593264a.png
Fig.3: 三次元のランダムウォーク。三次元の極座標では角度θとφの二つのパラメータが必要。


Scilabスクリプトは以下のようになります。

clear;

// *** 計算の設定 ***
r = 1; // 1ステップの間に移動する距離
tnum = 10000; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
theta = 2 * %pi * rand(1,tnum);
phi = 2 * %pi * rand(1,tnum);
xi = r .* sin(theta) .* cos(phi);
yi = r .* sin(theta) .* sin(phi);
zi = r .* cos(phi);

S = [zeros(3,1), [cumsum(xi,'c'); cumsum(yi,'c'); cumsum(zi,'c')]];

// *** グラフのプロット ***
param3d(S(1,:),S(2,:),S(3,:));
xlabel("x position");
ylabel("y position");
zlabel("z position");


ランダウの計算物理学による分類


ランダウの計算物理学 基礎編には以下のようにあります。

ランダムステップをどのように発生させるかで,異なる結果に至ることもありうる.以下に, 2次元のランダムウォークを発生するためのいくつかの方法をあげた.

  1. 方位角θを[0,2π]の間の乱数として選ぶ.次にΔx=cosθおよびΔy=sinθとする.(こうすると,一様な乱数を三角関数で写像することになり,当然のことだが水平・垂直の格子点上をランダムウォークするのとは異なる振る舞いとなる.)
  2. Δxを[-√2,√2]の範囲の乱数として選ぶ.これと独立にΔyも[-√2,√2]の範囲の乱数として選ぶ.こうすれば,xとyそれぞれの方向について,正負のステップが同じ確率で発生することになる.
  3. Δxを[-1,1]の範囲の乱数として選ぶ.次にΔy=±√(1-Δx2)とする.(符号もランダムに与える.)
  4. ステップの方向として(北, 東, 南, 西)をランダムに選択する(こうすると三角関数が不要になる).4つの方位から一つを選ぶのは[1,4]の整数を選択するのと等価であることに注意せよ.
  5. ステップの方向として(北, 北東, 東, 南東, 南, 南西, 西, 北東)をランダムに選択する(こうすると三角関数が不要になる).8つの方位から一つを選ぶのは[1,8]の整数を選択するのと等価であることに注意せよ.

今回の例は1.の方法を採用したものです。
そして、Scilabでブラウン運動 その1の方法は4.の方法と同じです。(移動する距離は二次元では√2となっていますが。)

関連エントリ




参考URL




付録


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


参考文献/使用機器





フィードバック



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

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


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


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

tag: Scilab 乱数 モンテカルロ解析 確率論 ランダムウォーク ブラウン運動 

Scilabでブラウン運動 その1

Scilabで楽しむ確率論(PDF)では一次元と二次元のランダムウォークをシミュレーションするScilabスクリプトが紹介されています。これらのスクリプトを参考にして三次元の対称ランダムウォークのシミュレーションを行いました。

003_20141109070543e4d.png
Fig.1: 三次元の対称ランダムウォーク



一次元対称ランダムウォーク


一次元の場合、t=0でx=0にあった粒子が1ステップの時間の間に、x軸の方向に半々の確率で+1または-1だけ移動するとします。この場合、Nステップ後に粒子がどこにいるのかをシミュレーションしました。

001_20141109070543015.png
Fig.2: 一次元の対称ランダムウォーク


以下はScilabで楽しむ確率論(PDF)による一次元対称ランダムウォークのスクリプト(に多少の編集とコメントを加えたもの)です。

clear;

// *** 時間ステップ数 ***
tnum = 10000; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
xi = 2 * (rand(tnum,1) >= 0.5) - 1;
S = [0;cumsum(xi)];

// *** グラフのプロット ***
plot(t,S,'-b');
xlabel("Time");
ylabel("Position");


二次元対称ランダムウォーク


二次元の場合t=0で(x,y)=(0,0)にあった粒子が1ステップの間に、x軸の方向には1/2の確率で+1, もう1/2の確率で-1移動し、同様にy軸の方向に対しても半々の確率で正の方向と負の方向にそれぞれ1ずつ移動するとします。この場合、Nステップ後に粒子がどこにいるのかをシミュレーションしました。

002_20141109070543ce1.png
Fig.3: 二次元の対称ランダムウォーク


以下はScilabで楽しむ確率論(PDF)による一次元対称ランダムウォークのスクリプト(に多少の編集とコメントを加えたもの)です。

clear;

// *** 時間 ***
tnum = 10000; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
xi = 2 * (rand(2,tnum) >= 0.5) - 1;
S = [zeros(2,1), cumsum(xi,'c')];

// *** グラフのプロット ***
plot(S(1,:),S(2,:),'-b');
xlabel("x position");
ylabel("y position");


三次元対称ランダムウォーク


一次元から二次元への拡張を見れば、三次元への拡張もほぼ自明です。
なお、三次元空間における軌跡の表示にはparam3dを利用します。

clear;

// *** 時間 ***
tnum = 10000; // 時間ステップ数
t = (0:1:tnum)'; // 時間のベクトル

// *** 位置の計算 ***
xi = 2 * (rand(3,tnum) >= 0.5) - 1;
S = [zeros(3,1), cumsum(xi,'c')];

// *** グラフのプロット ***
param3d(S(1,:),S(2,:),S(3,:));
xlabel("x position");
ylabel("y position");
zlabel("z position");


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 乱数 モンテカルロ解析 確率論 ランダムウォーク ブラウン運動 

n行ごとに改行を挿入するシェルスクリプト

Gnuplotで3次元プロットを行うためには、入力するデータのファイルに改行が入っていないとエラーが出ます。
そこでn行ごとに改行を挿入するシェルスクリプトinsbreak.shを書きました。

#!/bin/bash

tmpfile=$(mktemp)
cat ${1} | awk '{print $0} NR%'${2}'==0 {printf "\n"}' > ${tmpfile}
mv ${tmpfile} ${1}


使い方は以下のようにコマンドラインに打ち込みます。

./insbreak.sh hoge.dat 5



上記の例ではhoge.datにたいして5行ごとに改行が入ります。
仮にhoge.datが以下の内容だったとします。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20


するとhoge.datは以下のように上書きされます。

1
2
3
4
5

6
7
8
9
10

11
12
13
14
15

16
17
18
19
20



本当は安全な一時ファイルの削除の仕方とかもあるようですが、今回はここまで。(どうせこのスクリプトは自分では使わないし)
結果を上書きしてしまうので、もし使うならファイルのバックアップは取っておいた方がいいです。

関連エントリ




参考URL




付録


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


フィードバック



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

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


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


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

tag: cygwin シェルスクリプト awk 

AkaiKKRで反強磁性クロム

AkaiKKR(machikaneyama)を用いて整合反強磁性クロムの計算を行いました。反強磁性の計算の正しいやり方は分かりませんが、格子定数を大きくして無理矢理に強磁性状態にしたクロムのポテンシャルファイルを初期ポテンシャルとして計算を行うと、強磁性状態とも常磁性状態とも異なる反強磁性状態の解が得られませした。
しかしながら、全エネルギーが最小となる格子定数においては、反強磁性状態が安定とはなりませんでした。

CrDOS.gif
Fig.1: 反強磁性クロム(赤)と常磁性クロム(緑)の部分状態密度の格子体積依存性。格子体積が小さくなるとともに反強磁性的な磁気モーメントが消える。



反強磁性クロム


体心立方構造(bcc)のクロムは、反強磁性体の金属として知られています。Wikipediaの反強磁性の項目には以下のように、クロムを広義の反強磁性体として分類しています。

体心立方構造の頂点位置と体心位置のスピンが反対方向を向いているだけでなく、スピンの大きさが単位胞ごとに正弦関数的に変化している。

しかしながら、第一原理計算では、しばしば単純に体心と頂点のスピンが反対を向いているだけの反強磁性クロムの計算がなされます。(例えば第一原理計算ソフトウェア Advance/PHASEの応用機能と解析事例計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)のABCAPによる計算例など。)このような反強磁性を整合(commensurate)反強磁性と呼びます。



今回はAkaiKKR(machikaneyama)を用いて上記の単純な整合反強磁性クロムの計算を行いうことを目標にします。

反強磁性体のポテンシャルファイル


クロムの入力ファイルを作成し、そのままmagで計算を行っても反強磁性状態の解は得られません。

計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)のAkaiKKRの項目には、強磁性状態の鉄のポテンシャルファイルからfmgを用いてスピングラス状態のポテンシャルファイルを作成する方法が書かれています。さらにそのページには反強磁性状態の計算をする場合にも同じようにしてデータファイルを作ることができるとあります。

スピングラス状態の計算では、強磁性状態の計算で得られたポテンシャルファイルのスピンの向きを反転させたポテンシャルファイルを作成し、それらをCPAする事でスピングラス状態を実現していました。いうなればスピンの異なる原子の不規則合金です。ねがてぃぶろぐでもAkaiKKRで鉄のキュリー温度にて計算を行いました。
このことから類推するに、反強磁性状態も、スピンの異なる原子を交互に並べることで実現できそうです。つまりスピンの異なる原子の規則合金という事です。

そのためにはまず、強磁性状態のポテンシャルファイルを作成します。bcc構造のクロムの格子定数はa=2.88Å(a=5.45bohr)程度ですが、この格子定数で磁性状態を含む計算(mag)を行っても、常磁性状態の解が得られるのみです。
そこで意図的に大きな格子定数(例えばa=8.0bohrなど)を指定して強磁性解を得ておきます。

その後、utilにあるfmgを用いて、反転させたポテンシャルファイルを初期ポテンシャルとしてgo計算を行います。
シェルスクリプトやfmgのための入力ファイル、計算のテンプレートなどは以下のようになりました。実行ファイルfmgはパスの通ったところにおいてあるとし、入力ファイルのテンプレートはtemplateという名前のフォルダを作成して入れておきます。

#!/bin/csh -f

set A_LIST=( 5.60 5.55 5.50 5.45 5.40 5.35 5.30 5.25 5.20 )

rm analysis/CrAFM.dat
rm analysis/CrNMG.dat

specx < in/CrFMG.in > out/CrFMG.out
fmg < crafm

foreach A ( ${A_LIST} )
sed 's/'ABOHR'/'${A}'/g' template/CrAFM0.in > in/CrAFM_${A}.in
specx < in/CrAFM_${A}.in > out/CrAFM_${A}.out
sed 's/'ABOHR'/'${A}'/g' template/CrNMG0.in > in/CrNMG_${A}.in
specx < in/CrNMG_${A}.in > out/CrNMG_${A}.out

sed -e '$d' data/crafm.info | sed -n '$p' >> analysis/CrAFM.dat
sed -e '$d' data/crnmg.info | sed -n '$p' >> analysis/CrNMG.dat
end

grep "spin moment" out/CrAFM_*.out | awk 'NR % 4 == 1;NR % 4 == 2'




結果


ある一つのクロム原子に着目した部分状態密度(PDOS)の格子体積依存性がFig.1のgif動画です。格子体積が大きく、反強磁性的な磁気モーメントが存在する間はDOSの形状が非磁性状態のものと異なっていますが、体積が小さくなり磁気モーメントが消えるとともにDOSの形状も非磁性状態のものと同じになります。

DOS.png
Fig.3: 格子定数a=5.45bohrにおける反強磁性クロム(赤)と常磁性クロム(緑)の部分状態密度


Fig.3に示したのが実験により得られた格子定数付近(a=5.45bohr)における反強磁性クロムと常磁性クロムの部分状態密度です。反強磁性状態ではフェルミ準位近傍に特徴的なふくらみを見ることができます。

CrAFM-Energy.png
Fig.4: 全エネルギーの格子体積依存性


Fig.4に示したのが全エネルギーの格子体積依存性です。磁気モーメントが大きい間は、反強磁性相の方がエネルギー的に安定であることがわかります。磁気モーメントが小さくなるにしたがって全エネルギーの差が小さくなることも読み取れます。
しかし、全エネルギーが最小化するのはほとんど磁気モーメントが消失しているa=5.40bohr付近でエネルギー的にはほとんど差が見られません。それどころか、残念なことにわずかながら常磁性状態の方がエネルギー的に安定であるという結果になっています。

これが整合反強磁性を仮定したためなのか、単純にbzqlty(k点の数)が低いためなのか、数値計算の限界なのかは私には判断できません。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR 反強磁性 DOS 状態密度 

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関数フィッティング

最新コメント
リンク

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