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カウンター
カテゴリ
ユーザータグ

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

最新コメント
リンク

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