Scilabでイジング模型 その4

Scilabでイジング模型 その2では一次元のイジング模型の磁区が形成されていく様子をシミュレーションしました。
今回はイジング模型のモンテカルロシミュレーションから巨視的な物理量(エネルギーや磁化など)の計算を行います。
巨視的な系のエネルギー<E>は、状態|αj>におけるエネルギーE(αj)の統計平均として以下のように表されます。

\langle E \rangle = \sum_{\alpha_j}E(\alpha_j)P(\alpha_j)
P(\alpha_{j})=\frac{1}{Z(T)}\exp\left(- \frac{E(\alpha_{j})}{kT} \right)
Z(T)=\sum_{\alpha_j}\exp\left(-\frac{E(\alpha_j)}{kT}\right)

しかしながら、実際には分配関数Z(T)を計算するのが困難であるので、P(αj)に比例した確率で|αj>が出現するようなモンテカルロシミュレーションをm回行い、その平均から

\langle E \rangle = \frac{1}{m} \sum_{t=1}^{m}E_{t}(\alpha)

のように計算を行います。


一次元イジング模型における巨視的な物理量


Scilabでイジング模型 その2では一次元のイジング模型の磁区が形成されていく様子をシミュレーションしました。
今回はイジング模型のモンテカルロシミュレーションから巨視的な物理量(エネルギーや磁化など)の計算を行います。

n粒子スピン系におけるある状態ベクトル|αj>における系のエネルギーは以下のようになることをScilabでイジング模型 その2にも書きました。

E(\alpha_{j}) = - J \sum_{i=1}^{n-1}s_{i}s_{i+1}

よってある温度Tのときの状態|αj>が決まれば、系のエネルギーが決まることになりますが、実際には一意に決まるわけではなく、ボルツマン因子exp(-E(αj)/kT))に比例した確率P(αj)で色々な状態を取り得ます。

P(\alpha_{j})=\frac{1}{Z(T)}\exp\left(- \frac{E(\alpha_{j})}{kT} \right)
Z(T)=\sum_{\alpha_j}\exp\left(-\frac{E(\alpha_j)}{kT}\right)

ここでZ(T)は分配関数と呼ばれすべての状態の和です。

ひとたび確率P(αj)が求まれば、巨視的な物理量<A>は状態|αj>における物理量A(αj)を用いて以下のようにあらわすことができます。

\langle A \rangle = \sum_{\alpha_j}A(\alpha_j)P(\alpha_j)

例えばエネルギーEの場合は以下のようになります。

\langle E \rangle = \sum_{\alpha_j}E(\alpha_j)P(\alpha_j)

しかし実際にZ(T)をすべて計算するのは不可能です。代わりにP(αj)に比例した確率で|αj>を出現させるアルゴリズムがメトロポリスのアルゴリズムでした。

そこで系の巨視的なエネルギー(等の物理量)を計算する際には、メトロポリスのアルゴリズムを複数回繰り返して得られた状態に対するエネルギーの平均を代わりに用います。m回繰り返す場合は

\langle E \rangle = \frac{1}{m} \sum_{t=1}^{m}E_{t}(\alpha)

となります。

他にも微分方程式による物理現象のモデル化(PDF)には磁化M、比熱C、磁化率χの表式が以下のように与えられています。

M(\alpha_j)=\sum_{i=1}^{n}s_i
\langle M \rangle = \frac{1}{m} \sum_{t=1}^{m}M_{t}(\alpha)

C = \frac{\langle E^2 \rangle - \langle E \rangle^2}{kT}

\chi = \frac{\langle M^2 \rangle - \langle M \rangle^2}{kT}

Scilabスクリプト


これらを踏まえたScilabスクリプトを以下に示します。

clear;

// *** 定数の設定 ***
n = 100; // 粒子の数
m = 200; // 熱力学的な平均を取る回数
J = 1; // 交換エネルギー
rand("uniform"); // 乱数は一様乱数とする
tmax = 5 * n; // 時間の最大ステップ
h = 0.0; // 外部磁場
// 温度
ktmin = 0.5; // 最低温度
ktmax = 5.0; // 最高温度
nkt = 19; // 温度の分割数
//T = linspace(ktmin, ktmax, nkt); // 低温から開始
//spin = ones(1,n); // 各粒子におけるスピン(コールドスタート)
T = linspace(ktmax, ktmin, nkt); // 高温から開始
spin = 1 - 2 * round(rand(1,n)); // 各粒子におけるスピン(ランダム)

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

// *** 行列の初期化 ***
E = []; // エネルギーの和
E2 = []; // エネルギーの二乗の和
M = []; // 磁化の和
M2 = []; // 磁化の二乗の和

// *** 温度のループ ***
for kt = 1:nkt do
// エネルギーの初期化
ene1 = 0; // エネルギーの和
ene2 = 0; // エネルギーの二乗和
// 磁化の初期化
mag1 = 0; // 磁化の和
mag2 = 0; // 磁化の二乗和f
// *** 熱力学平均のループ ***
for samp = 1:m do
// *** 時間発展のループ ***
for t = 1:tmax do
oldenergy = energy(spin);
element = ceil(n * rand()); // 粒子を一つ選ぶ
spin(element) = -1 * spin(element); // スピンを反転
newenergy = energy(spin);
if (newenergy > oldenergy) & (exp((- newenergy + oldenergy) / T(kt)) < rand()) then
spin(element) = -1 * spin(element); // 棄却
end
end
ene1 = ene1 + energy(spin); // エネルギーの和
ene2 = ene2 + energy(spin)^2; // エネルギーの二乗の和
mag1 = mag1 + sum(sum(spin)); // 磁化の和
mag2 = mag2 + sum(sum(spin))^2; // 磁化の二乗和
end
E = [E, ene1 / m]; // エネルギーの和
E2 = [E2, ene2 / m]; // エネルギーの二乗の和
M = [M, mag1 / m]; // 磁化の和
M2 = [M2, mag2 / m]; // 磁化の二乗和
end

// *** エネルギーと磁化の揺らぎ ***
C = (E2 - E .^ 2) ./ (n * T .^ 2); // 比熱
X = (M2 - M .^ 2) ./ (n * T); // 磁化率

// *** 厳密解の計算 ***
// 温度ベクトル
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, 1 ./ X, 'or');
plot(Ta,RXa,'--g');
xlabel("kT/J");
ylabel("N/JX");


結果はScilabでイジング模型 その1と同様になります(なので下記の画像は使いまわしです)。コードもほとんど同じです。

001_20141130214554d90.png

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


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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


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

comment

Secret

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

最新コメント
リンク

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