はじめに

こんにちは．東京大学理学部物理学科B2の灰陶です．Physics Lab. では統計物理班と数理物理班に参加しています． 本記事は，東京大学理学部物理学科の有志による 理物 Advent Calender 2022 の17日目の記事です．

https://qiita.com/advent-calendar/2022/ributu

本記事では，プログラミング言語 Julia について紹介し，Julia を用いて 2次元 Ising モデルのシミュレーションを行ってみます．

What is Julia?

Julia は科学技術計算向けの言語として開発されました．2018年に安定版のバージョン1がリリースされたばかりの，まだまだ成長途中の言語です．

Julia には，以下のような特徴があります．

1. Python のように平易な文法
2. 動的型付け言語でありながら，C や Fortran などの静的型付け言語に迫る計算速度
3. MATLAB と同等の線型代数計算の記述
4. Unicode 文字対応に起因する高い可読性

これらの特徴についての詳しい説明はここでは避けます．詳しく知りたいという方は，文献1を読んでみてください．

シミュレーション概要

シミュレーションに必要な Ising モデルの知識とシミュレーションの方法について説明します．

Ising モデルミニマム

格子状の2次元 Ising モデルの Hamiltonian は次式で与えられます：

ここで，第1項の和は次のように書き換えられます：

外部磁場のない（\(h=0\)）2次元 Ising モデルについては厳密解が見つかっています．それによると，2次元 Ising モデルでは，

今回のシミュレーションでは，\(L_x\times L_y\hspace{-0.2em}\hspace{0.2em}\) 格子状2次元 Ising モデルのスピン1つ当たりのエネルギー，磁化，比熱の温度依存性をメトロポリスモンテカルロ法により再現し，相転移が起こっていることを確認します．

メトロポリスモンテカルロ法ミニマム

モンテカルロ法とは，乱数を用いた人工的なダイナミクスを導入することで物理系のシミュレーションを行う方法です．その中でも，系の分布の列を Markov 連鎖として生成する方法をメトロポリスモンテカルロ法と呼び，多粒子系の静的な性質を研究するために用いられます．

今回のシミュレーションでは，メトロポリスモンテカルロ法に従って，あるスピン配位 \(\hspace{-0.2em}\hspace{0.2em}C_k\hspace{-0.2em}\hspace{0.2em}\) が与えられたときに，新たなスピン配位 \(\hspace{-0.2em}\hspace{0.2em}C_{k+1}\hspace{-0.2em}\hspace{0.2em}\) を定めます．その概略は以下のようになります：

1. スピン配位 \(\hspace{-0.2em}\hspace{0.2em}C_k\hspace{-0.2em}\hspace{0.2em}\) に対し，無作為にスピンを選び，そのスピンを反転させた試行配位 \(\hspace{-0.2em}\hspace{0.2em}C_k'\hspace{-0.2em}\hspace{0.2em}\) を作成する．
2. 古いスピン配位と試行配位のエネルギー差 \(\hspace{-0.2em}\hspace{0.2em}\Delta E(C_k\to C_k')\hspace{-0.2em}\hspace{0.2em}\) を計算する．
3. エネルギーが増加する，即ち \(\hspace{-0.2em}\hspace{0.2em}\Delta E(C_k\to C_k') > 0\hspace{-0.2em}\hspace{0.2em}\) の場合，確率 \(\hspace{-0.2em}\hspace{0.2em}\exp[-\beta\Delta E(C_k\to C_k')]\hspace{-0.2em}\hspace{0.2em}\) で \(\hspace{-0.2em}\hspace{0.2em}C_{k+1}=C_k'\hspace{-0.2em}\hspace{0.2em}\) とし（試行配位の採択），確率 \(\hspace{-0.2em}\hspace{0.2em}1-\exp[-\beta\Delta E(C_k\to C_k')]\hspace{-0.2em}\hspace{0.2em}\) で \(\hspace{-0.2em}\hspace{0.2em}C_{k+1}=C_k\hspace{-0.2em}\hspace{0.2em}\) とする（試行配位の棄却）．
4. エネルギーが減少するなら，常に \(\hspace{-0.2em}\hspace{0.2em}C_{k+1}=C_k'\hspace{-0.2em}\hspace{0.2em}\) とする．

この過程を繰り返すことで系のダイナミクスを再現し，系の平衡状態の性質を調べていきます．なお，メトロポリスモンテカルロ法を用いるためには，生成したい状態の集合が次の条件を満たしていなければなりません：

• 任意の状態が，他の任意の状態から高々有限回の過程の繰り返しによって到達できる．
• 周期性がない．

以上の条件を満たす Marcov 連鎖はエルゴード的であると言います．

シミュレーションの設定

今回は，以下の設定の下で Ising モデルのシミュレーションを行います：

• 格子のサイズは \(\hspace{-0.2em}\hspace{0.2em}50\times50\)．
• \(\hspace{-0.2em}J=1,~h=0\)．
• 格子の境界においては，周期的境界条件を課して計算を行う．
• メトロポリスモンテカルロ法の総ステップ数を \(\hspace{-0.2em}\hspace{0.2em}60000\)，熱平衡に達するまでのステップ数を \(\hspace{-0.2em}\hspace{0.2em}10000\hspace{-0.2em}\hspace{0.2em}\) とし，物理量は \(\hspace{-0.2em}\hspace{0.2em}10001\hspace{-0.2em}\hspace{0.2em}\) ステップ目から \(\hspace{-0.2em}\hspace{0.2em}10\hspace{-0.2em}\hspace{0.2em}\) ステップ間隔で値を記録し，平均をとる．
• 各種物理量を温度 \(\hspace{-0.2em}\hspace{0.2em}k_{\text{B}}T/J=0\hspace{-0.2em}\hspace{0.2em}\) から \(\hspace{-0.2em}\hspace{0.2em}k_{\text{B}}T/J=5\hspace{-0.2em}\hspace{0.2em}\) で計算．

コード

今回のシミュレーションのコードは以下のようになります．なお，このコードは文献2を参考に，一部変更を加えて作成しました．Julia1.8.2 を使用し，Jupyter Notebook 上で動作することを確認しています．

using Random
using Plots
using LaTeXStrings

# スピン配位の初期化
function spins_initialize(Lx, Ly, rng)
    return rand(rng, [-1, 1], Lx, Ly)
end

# Siの計算
function Si_calculation(Ck, ix, iy, Lx, Ly)
    jx = ix + 1
    if jx > Lx
        jx -= Lx
    end
    jy = iy
    Si = Ck[jx, jy]

    jx = ix - 1
    if jx < 1
        jx += Lx
    end
    jy = iy
    Si += Ck[jx, jy]

    jx = ix
    jy = iy + 1
    if jy > Ly
        jy -= Ly
    end
    Si += Ck[jx, jy]

    jx = ix
    jy = iy - 1
    if jy < 1
        jy += Ly
    end
    Si += Ck[jx, jy]

    return Si
end

# エネルギーの計算
function E_measure(Ck, J, h, Lx, Ly)
    E = 0
    for iy = 1:Ly
        for ix = 1:Lx
            Si = Si_calculation(Ck, ix, iy, Lx, Ly)
            σi = Ck[ix, iy]
            E += -(J/2)*σi*Si - h*σi
        end
    end
    return E
end

# エネルギー差の計算
function ΔE_calculation(Ck, J, h, ix, iy, Lx, Ly)
    Si = Si_calculation(Ck, ix, iy, Lx, Ly)
    return 2J*Ck[ix, iy]*Si + 2h*Ck[ix, iy]
end

# メトロポリスモンテカルロ法によるスピンの更新
function spins_update(Ck, ix, iy, J, h, T, Lx, Ly, rng)
    ΔE = ΔE_calculation(Ck, J, h, ix, iy, Lx, Ly)
    σi = Ck[ix, iy]
    σ_new = ifelse(rand(rng)≤exp(-ΔE/T), -σi, σi)
    return σ_new
end


# メトロポリスモンテカルロ法の実行，各種物理量の計算
function Ising_Montecarlo(step_equilibration, step_MC, measure_interval, J, h, T, Lx, Ly)
    rng = MersenneTwister(123)  # 乱数の種
    count_measure = 0  # 測定回数
    E_mean = 0  # エネルギーの時間平均
    E2_mean = 0  # エネルギーの2乗の時間平均
    m_mean = 0  # 磁化の時間平均
    absm_mean = 0  # 磁化の絶対値の時間平均

    Ck = spins_initialize(Lx, Ly, rng)

    for trial = 1:step_equilibration+step_MC
        for ix = 1:Lx
            for iy = 1:Ly
                Ck[ix, iy] = spins_update(Ck, ix, iy, J, h, T, Lx, Ly, rng)
            end
        end

        if trial > step_equilibration && trial % measure_interval == 0
            count_measure += 1
            E = E_measure(Ck, J, h, Lx, Ly) / (Lx*Ly)
            E_mean += E
            E2_mean += E^2
            m = sum(Ck) / (Lx*Ly)  # 磁化の計算
            m_mean += m
            absm_mean += abs(m)
        end
    end

    E_mean /= count_measure
    E2_mean /= count_measure
    m_mean /= count_measure
    absm_mean /= count_measure
    C_mean = (E2_mean-E_mean^2) / T^2

    return E_mean, m_mean, absm_mean, C_mean
end

# 上記の設定でシミュレーションを実行
function test()
    step_equilibration = 10000  # 熱平衡に達するまでのステップ数
    step_MC = 50000  # 時間平均計算に用いるステップ数
    measure_interval = 10  # 数値の測定間隔
    J = 1  # 交換相互作用定数
    h = 0  # 外部磁場
    Ts = range(0.01, 5, length=200)  # 温度
    Lx = 50  # 格子のサイズ
    Ly = 50
    E_T = []  # 各温度におけるエネルギー
    m_T = []  # 各温度における磁化
    absm_T = []  # 各温度における磁化の絶対値
    C_T = []  # 各温度における比熱

    # 各温度で各種物理量を計算，記録
    for T in Ts
        E_mean, m_mean, absm_mean, C_mean = Ising_Montecarlo(step_equilibration, step_MC, measure_interval, J, h, T, Lx, Ly)
        push!(E_T, E_mean)
        push!(m_T, m_mean)
        push!(absm_T, absm_mean)
        push!(C_T, C_mean)
    end

    # グラフ描画
    plot(Ts, E_T, label="", xlabel=L"k_{\textrm{B}}T", ylabel="Energy", title="Temperature Dependence of Energy")
    savefig("E_Tdep.png")
    plot(Ts, m_T, label="", xlabel=L"k_{\textrm{B}}T", ylabel="Magnetization", title="Temperature Dependence of Magnetization")
    savefig("m_Tdep.png")
    plot(Ts, absm_T, label="", xlabel=L"k_{\textrm{B}}T", ylabel="Absolute Value of Magnetization", title="Temperature Dependence of Magnetization")
    savefig("absm_Tdep.png")
    plot(Ts, C_T, label="", xlabel=L"k_{\textrm{B}}T", ylabel="Specific Heat", title="Temperature Dependence of Specific Heat")
    savefig("C_Tdep.png")
end

コードを読んでいただければわかると思いますが，文法がかなり Python に似ており，直感的に理解しやすいですよね．

結果

シミュレーションによって得られた結果を示します．

エネルギーの温度依存性

上図はエネルギーの温度依存性を示したグラフです．高温になるにつれて隣同士のスピンが反対向きになり，エネルギーが大きくなっていく様子がみられます．

磁化の温度依存性

磁化の温度依存性は上図のようになります．低温側ではスピンが揃っているために磁化の絶対値は \(\hspace{-0.2em}\hspace{0.2em}\mu_0\hspace{-0.2em}\hspace{0.2em}\) または \(\hspace{-0.2em}\hspace{0.2em}-\mu_0\hspace{-0.2em}\hspace{0.2em}\) となっています．今回は外部磁場をゼロとしているために，スピンの揃う向きが「上向き」でも「下向き」でもエネルギーの値は変わりません．そのため，初期配位によってはどちらの向きにも揃い得ます． 磁化の絶対値の温度依存性を示したこちらの図の方が分かりやすいですね．臨界温度付近で磁化が急激にゼロに近づいていく様子が見て取れます．

比熱の温度依存性

比熱の温度依存性は上図のグラフの通りです．臨界温度付近で値が急激に増加していますね．相転移の兆候が現れていると言って良さそうです．

まとめ

いかがだったでしょうか．Julia という言語を紹介しようと思い書き始めたものの，どうせならちょっと背伸びした題材にしようとした結果，当初の目的からは少しズレた中途半端な記事になってしまいました．

余談ながら，Julia のアドベントカレンダーも存在するようなのでこちらで紹介しておきます．

https://qiita.com/advent-calendar/2022/julia

皆さんも Julia と共に良きシミュレーションライフを！ ありがとうございました．

参考文献

Julia については文献1, 2, 3を，Ising モデルについては文献4, 5, 6を，メトロポリスモンテカルロ法については文献7を参考にしました．

1. Julia 1.8 Documentation https://docs.julialang.org/en/v1/
2. Bogumił Kamiński, Przemysław Szufe（2019）「Julia プログラミングクックブック」オーライリー・ジャパン
3. 永井佑紀（2022）「1週間で学べる！　Julia 数値計算プログラミング」講談社
4. 田崎晴明（2008）「統計力学1」（新物理学シリーズ37）培風館
5. 田崎晴明（2008）「統計力学2」（新物理学シリーズ38）培風館
6. 藤田重次（1990）「量子統計物理学」裳華房
7. Jos. M. Thijssen（2003）「計算物理学」シュプリンガー・フェアラーク東京

特に，文献2は物理系の学生・研究者が Julia に入門する際にうってつけの1冊です．この記事で Julia に興味を持ってくださった方，ぜひ読んでみてください．