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アンプPICCPA強磁性常微分方程式モンテカルロ解析odeトランジスタ状態密度インターフェーススイッチング回路ecaljPDS5022DOS定電流半導体シェルスクリプト乱数レベルシフトHP6632Aブレッドボード分散関係温度解析R6452Aトランジスタ技術I2C可変抵抗反強磁性セミナー数値積分確率論偏微分方程式バンド構造非線形方程式ソルババンドギャップ絶縁熱設計シュミットトリガLEDA/Dコンバータ三端子レギュレータLM358ISO-I2CGW近似カオスフォトカプラマフィンティン半径TL431数値微分PC817Cアナログスイッチ直流動作点解析発振回路USBサーボカレントミラー74HC4053パラメトリック解析LDAbzqltyチョッパアンプ量子力学FFT2ちゃんねるアセンブラBSch開発環境電子負荷ブラべ格子イジング模型補間基本並進ベクトル標準ロジック単振り子キュリー温度繰り返しMaxima状態方程式失敗談相対論スピン軌道相互作用FETランダムウォーク熱伝導六方最密充填構造コバルトewidthTLP621GGAQSGW不規則合金位相図抵抗SMPcygwinラプラス方程式スレーターポーリング曲線gfortranスイッチト・キャパシタ詰め回路TLP552三角波格子比熱TLP521条件分岐LM555MCUNE555QNAPマントルテスタ過渡解析FXA-7020ZRダイヤモンドデータロガーガイガー管自動計測Writer509UPSシュレディンガー方程式ブラウン運動awk差し込みグラフ熱力学平均場近似仮想結晶近似VCAfsolve井戸型ポテンシャルVESTA起電力スーパーセルOpenMP第一原理計算ubuntu固有値問題L10構造OPA2277interp12SC1815fccウィグナーザイツ胞面心立方構造フィルタジバニャン方程式ヒストグラム確率論マテリアルデザインspecx.f等高線正規分布PGAフェルミ面非線型方程式ソルバ初期値固定スピンモーメントスワップ領域ルチル構造リジッドバンド模型edeltquantumESPRESSO岩塩構造BaOSIC二相共存ZnOウルツ鉱構造フォノンデバイ模型c/aノコギリ波全エネルギーFSMTeXgnuplotmultiplotハーフメタルCapSense半金属合金結晶磁気異方性Ubuntu文字列入出力TS-110TS-112疎行列Excel直流解析ヒストグラム円周率不規則局所モーメントトラックボールPC等価回路モデルパラメータ・モデルキーボードRealforce三次元マンデルブロ集合フラクタル化学反応重積分縮退日本語最小二乗法関数フィッティングGimpMAS830LHiLAPW熱拡散方程式両対数グラフナイキスト線図負帰還安定性陰解法Crank-Nicolson法P-10クーロン散乱境界条件連立一次方程式片対数グラフEAGLEPIC16F785LMC662トランスシンボルCK1026線種凡例MBEAACircuitグラフの分割軸ラベルifort

最新コメント
リンク

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