🗂
2次元イジングモデルのモンテカルロシミュレーション

2025/10/18に公開
 はじめにイジングモデルは、統計力学における基本的なモデルの一つであり、磁性体の性質を理解するために広く用いられています。本記事では、2 次元イジングモデルのモンテカルロシミュレーションを Python で実装し、その結果を解析します。

 モデル
 2 次元イジングモデルのハミルトニアンイジングモデルのハミルトニアンは以下のように定義されます。

H = -J \sum_{\langle i,j \rangle} S_i S_j - h \sum_i S_i
ここで、S_iはスピン変数であり、S_i = \pm 1を取ります。Jはスピン間の相互作用定数、hは外部磁場を表します。\langle i,j \rangleは隣接するスピンのペアを示します。

 平衡状態平衡系の統計力学の結果から、イジングモデルが有限温度Tの状態において、その状態の確率は以下の分布に従うことが知られています。
確率分布:

  P(\{S\}_{i}) = \frac{1}{Z} e^{-\beta H(\{S\}_{i})}
  
分配関数:

  Z = \sum_{\{S\}_{i}} e^{-\beta H(\{S\}_{i})}
  ここで\beta = \frac{1}{k_B T}は逆温度、k_Bはボルツマン定数です。
この確率分布はカノニカル分布と呼ばれます。

 物理量の定義イジングモデルの物理量は以下で定義されます。
磁化 M: 全スピンの和で定義されます。

  M(\{S\}_{i}) = \sum_i S_i
  
エネルギー E: ハミルトニアンの値で定義されます。

  E(\{S\}_{i}) = -J \sum_{\langle i,j \rangle} S_i S_j - h \sum_i S_i
  
上記で定義される物理量は、スピンのゆらぎにより変化します。そのため、これらの物理量の期待値を計算することが重要です。

 物理量の期待値物理量の期待値は、物理量の確率分布に対する平均として定義されます。平衡状態における物理量の期待値は以下のように計算されます。
磁化の期待値 \langle M \rangle:

  \langle M \rangle = \frac{1}{Z} \sum_{\{S\}_{i}} M(\{S\}_{i}) e^{-\beta E(\{S\}_{i})}
  
エネルギーの期待値 \langle E \rangle:

  \langle E \rangle = \frac{1}{Z} \sum_{\{S\}_{i}} E(\{S\}_{i}) e^{-\beta E(\{S\}_{i})}
  
また比熱や磁化率の物理量は以下のように計算されます。
比熱 C: エネルギーの温度に対する変化率で定義されます。

  C = \frac{\partial E}{ \partial T}
  平衡状態ではエネルギーの分散を用いて計算出来ます。

  C = \frac{1}{k_B T^2} \left( \langle E^2 \rangle - \langle E \rangle^2 \right)
  
磁化率 \chi: 磁化の外部磁場に対する変化率で定義されます。

  \chi = \frac{\partial \langle M \rangle}{\partial h}
  平衡状態では磁化の分散を用いて計算出来ます。

  \chi = \frac{1}{k_B T} \left( \langle M^2 \rangle - \langle M \rangle^2 \right)
  

 厳密計算上記の物理量を計算してみましょう。
コードfrom numba import njit
import numpy as np
from typing import Dict, List
from main import plot_results


@njit()
def energy(spins, L, J, h):
    E = 0.0
    for idx in range(L * L):
        i = idx // L
        j = idx % L

        s = spins[idx]

        right_idx = i * L + ((j + 1) % L)
        E += s * spins[right_idx]

        down_idx = ((i + 1) % L) * L + j
        E += s * spins[down_idx]
    E = -J * E - h * np.sum(spins)
    return E


@njit()
def magnetization(spins):
    return np.sum(spins)


@njit()
def magnetization_abs(spins):
    return np.abs(np.sum(spins))


@njit()
def exact_calculation_ising(L, J, h, T):
    N = L * L
    total_states = 2 ** N

    beta = 1.0 / T

    Z = 0.0
    M_abs_sum = 0.0
    M2_sum = 0.0
    E_sum = 0.0
    E2_sum = 0.0

    for state_i in range(total_states):
        spins = np.empty(N, dtype=np.float64)
        for j in range(N):
            spins[j] = 1.0 if (state_i >> j) & 1 else -1.0

        E = energy(spins, L, J, h)
        M = magnetization(spins)
        M_abs = magnetization_abs(spins)

        boltzmann = np.exp(-beta * E)

        Z += boltzmann
        E_sum += E * boltzmann
        E2_sum += E * E * boltzmann
        M_abs_sum += M_abs * boltzmann
        M2_sum += M * M * boltzmann

    avg_E = E_sum / Z
    avg_E2 = E2_sum / Z
    avg_M_abs = M_abs_sum / Z
    avg_M2 = M2_sum / Z

    energy_mean = avg_E / N
    magnetization_mean = avg_M_abs / N
    specific_heat = (avg_E2 - avg_E * avg_E) / (T * T * N)
    susceptibility = (avg_M2 - avg_M_abs * avg_M_abs) / (T * N)

    return energy_mean, magnetization_mean, specific_heat, susceptibility


def temperature_scan(L: int, J: float, h: float, T_values: List[float]):
    results = {}
    for T in T_values:
        energy, magnetization, specific_heat, susceptibility = exact_calculation_ising(
            L, J, h, T)
        results[T] = {
            'energy': energy,
            'magnetization': magnetization,
            'specific_heat': specific_heat,
            'susceptibility': susceptibility
        }
    return results


def main():
    L = 4
    J = 1.0
    h = 0.0
    T_values = np.linspace(1.0, 4.0, 10)

    results = temperature_scan(L, J, h, T_values)
    plot_results(results)


if __name__ == '__main__':
    main()
直接計算は、L \times L格子の場合2^{L \times L}通りの状態を列挙し期待値を計算する必要がある。そのため、格子サイズが大きくなると計算量が指数関数的に増加し、実用的ではなくなる。

 モンテカルロシミュレーション
 考え方モンテカルロシミュレーションは、厳密計算の代替手法として一部の状態を選択的にサンプリングし、物理量の期待値を近似的に求める方法です。これにより、計算量を大幅に削減できます。具体的には、カノニカル分布P(\{S\}_{i})に従う状態をサンプリングし、物理量の期待値を以下のように近似します。

\langle M \rangle \approx \frac{1}{N_s} \sum_{k=1}^{N_s} M(\{S\}_{i}^{(k)})
ここで、N_sはサンプリングした状態の数、\{S\}_{i}^{(k)}はサンプリングされた状態を表します。
しかし、カノニカル分布に従って状態を直接サンプリングすることは困難です。なぜなら、確率分布を計算するには分配関数Zの計算が必要であり、これは全ての状態を列挙する必要があるためです。
そこで、どのようにしてカノニカル分布に従う状態をサンプリングするかが問題となります。

 マルコフ連鎖モンテカルロ法そこで、カノニカル分布に従う状態を生成するために、以下の戦略を取ります。

マルコフ連鎖を構築: 状態から状態への遷移過程を定義

遷移確率を設計: 連鎖の遷移確率をカノニカル分布を定常分布とするように設計

サンプリング: 十分に長い時間連鎖を進めた後の状態をサンプリング

 マルコフ連鎖とはマルコフ連鎖とは、次の状態が現在の状態のみに依存する確率過程です。

 定義状態\{S\}_iから状態\{S\}_jへの遷移確率W(\{S\}_i \to \{S\}_j)を定義します。

マルコフ性: 次の状態は現在の状態のみに依存します。

P(\{S\}_{n+1} | \{S\}_n, \{S\}_{n-1}, \ldots) = P(\{S\}_{n+1} | \{S\}_n)

つまり、過去の履歴\{S\}_{n-1}, \{S\}_{n-2}, \ldotsは次の状態に影響しません。

正規化条件: 全ての遷移確率の和は 1 になります。

W(\{S\}\_i \to \{S\}\_j) \geq 0, \quad \sum_j W(\{S\}\_i \to \{S\}\_j) = 1

 定常分布マルコフ連鎖を十分に長く進めると、状態の確率分布がある特定の分布に収束する場合があります。この分布を定常分布と呼びます。定常分布\pi(\{S\}_i)は以下の条件を満たします。

P(\{S\}_j) = \sum_i P(\{S\}_i) W(\{S\}_i \to \{S\}_j)
この条件は、現在の状態の分布がPのとき、マルコフ連鎖を 1 ステップ進めた後の状態の分布もPであることを意味します。
この式を満たす遷移確率W(\{S\}_i \to \{S\}_j)を設計することで、特定の定常分布を持つマルコフ連鎖を構築できます。

 詳細釣り合い条件定常分布の条件を満たすための十分条件として、以下の条件を考えます。

P(\{S\}_i) W(\{S\}_i \to \{S\}_j) = P(\{S\}_j) W(\{S\}_j \to \{S\}_i)
この条件を詳細釣り合い条件と呼びます。直感的には、状態\{S\}_iから状態\{S\}_jへの遷移の確率と、状態\{S\}_jから状態\{S\}_iへの遷移の確率が等しいことを意味します。この条件を満たすとき、定常分布の条件も自動的に満たされます。
具体的に、カノニカル分布を定常分布とするマルコフ連鎖を構築するために、以下のように遷移確率を設計します。

\frac{e^{-\beta E(\{S\}_i)}}{Z} W(\{S\}_i \to \{S\}_j) = \frac{e^{-\beta E(\{S\}_j)}}{Z} W(\{S\}_j \to \{S\}_i)
整理すると、

\frac{W(\{S\}_i \to \{S\}_j)}{W(\{S\}_j \to \{S\}_i)} = e^{-\beta (E(\{S\}_j) - E(\{S\}_i))}
となります。

 メトロポリス法詳細釣り合い条件から、以下の関係を満たす遷移確率Wを設計する必要があります。

\frac{W(\{S\}_i \to \{S\}_j)}{W(\{S\}_j \to \{S\}_i)} = e^{-\beta \Delta E}
ここで、\Delta E_{ij} = E(\{S\}_j) - E(\{S\}_i)です。
遷移確率Wを以下のように分割して考えます。

W(\{S\}_i \to \{S\}_j) = g(\{S\}_i \to \{S\}_j) A(\{S\}_i \to \{S\}_j)
ここで提案確率g(\{S\}_i \to \{S\}_j)は状態\{S\}_iにおいて状態\{S\}_jへの遷移を提案する確率を表し、 受理確率A(\{S\}_i \to \{S\}_j)は提案された候補\{S\}_jを受け入れる確率を表します。
今回はランダムに 1 つのスピンを選び反転させる提案を行います。

このときの提案確率は、

g(\{S\}_i \to \{S\}_j) = \left\{\begin{array}{ll}
\frac{1}{N} & \text{if }\{S\}_j \text{ is obtained by flipping one spin of } \{S\}_i \\
0 & \text{otherwise}
\end{array}\right.
となります。ここで、Nはスピンの総数です。この提案確率は対称であるため、
g(\{S\}_i \to \{S\}_j) = g(\{S\}_j \to \{S\}_i)が成り立ちます。
また、受理確率Aは以下のように定義します。

A(\{S\}_i \to \{S\}_j) = \left\{\begin{array}{ll}
1 & \text{if } \Delta E \leq 0 \\
e^{-\beta \Delta E} & \text{if } \Delta E > 0
\end{array}\right.
このように選択することで、詳細釣り合い条件を満たすことができます。

 アルゴリズムメトロポリス法を用いたモンテカルロシミュレーションのアルゴリズムは以下の通りです。
スピンの選択:

系からランダムに 1 つのスピンを選択します。
エネルギー変化の計算:

選択したスピンを反転させたときのエネルギー変化\Delta Eを計算します。

\Delta E = 2 J S_i \sum_{j \in \text{nn}(i)} S_j + 2 h S_i

ここで、\text{nn}(i)はスピンiの隣接スピンの集合を表します。系全体のエネルギーを再計算する必要はなく、選択したスピンとその隣接スピンのみを考慮すればよいです。
受理判定:

エネルギー変化\Delta Eに基づいて、スピンの反転を受理するかどうかを決定します。

\Delta E \leq 0の場合、スピンの反転を無条件に受理します。

\Delta E > 0の場合、確率e^{-\beta \Delta E}でスピンの反転を受理します。具体的には、[0,1)の一様乱数rを生成し、r < e^{-\beta \Delta E}ならば受理します。

スピンの更新:

受理された場合、選択したスピンを反転させます。受理されなかった場合は、スピンをそのままにします。

 初期化と平衡化初期状態はカノニカル分布に従わないため、シミュレーションの最初の数ステップは平衡化のために使用します。平衡化が完了した後、物理量のサンプリングを開始します。具体的には、以下の手順で進めます。
初期状態の設定: ランダムにスピンを配置します。
平衡化ステップ: 一定数のモンテカルロステップを実行し、系を平衡状態に近づけます。
サンプリングステップ: 平衡化後、物理量のサンプリングを行います。
物理量の計算: 各サンプリングステップでエネルギーと磁化を計算し、記録します。

 実装例上記のアルゴリズムを実装したコード例を以下に示します。
コードfrom typing import Dict, List

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from numba import int64, float32
from numba.experimental import jitclass
import time

sns.set_style("darkgrid")

"""
Simple implementation of a 2D Ising model.

Hamiltonian:
H = -J * sum(s_i * s_j) - h * sum(s_i)
where:
- J is the interaction strength between neighboring spins
- h is the external magnetic field
- s_i are the spin variables (+1 or -1)
- The first sum is over all neighboring pairs of spins
"""

spec = [
    ('L', int64),
    ('J', float32),
    ('h', float32),
    ('T', float32),
    ('state', float32[:, :])
]


@jitclass(spec)
class Model:
    def __init__(self, L, J, h, T):
        self.L = L
        self.J = J
        self.h = h
        self.T = T
        self.state = (2 * np.random.randint(0, 2, (L, L)) - 1).astype(
            np.float32)

    def energy(self):
        energy = 0.0
        for i in range(self.L):
            for j in range(self.L):
                s = self.state[i, j]
                energy += s * self.state[i, (j + 1) % self.L]
                energy += s * self.state[(i + 1) % self.L, j]
        return -self.J * energy - self.h * np.sum(self.state)

    def magnetization(self):
        return np.sum(self.state)

    def delta_energy(self, i, j):
        delta_energy_first_term = 2 * self.J * self.state[i, j] * (
                self.state[(i + 1) % self.L, j] +
                self.state[(i - 1) % self.L, j] +
                self.state[i, (j + 1) % self.L] +
                self.state[i, (j - 1) % self.L]
        )
        delta_energy_second_term = 2 * self.h * self.state[i, j]
        return delta_energy_first_term + delta_energy_second_term

    def step(self):
        for _ in range(self.L * self.L):
            i = np.random.randint(0, self.L)
            j = np.random.randint(0, self.L)
            delta_E = self.delta_energy(i, j)
            if delta_E < 0 or np.random.rand() < np.exp(-delta_E / self.T):
                self.state[i, j] *= -1

    def run(self, n_steps):
        for _ in range(n_steps):
            self.step()

    def thermalize(self):
        L = self.L
        self.run(100 * L * L)

    def sampling_observables(self, n_samples, n_steps_between):
        energy = np.empty(n_samples, dtype=np.float32)
        magnetization = np.empty(n_samples, dtype=np.float32)
        for n in range(n_samples):
            self.run(n_steps_between)
            energy[n] = self.energy()
            magnetization[n] = self.magnetization()

        N = self.L * self.L

        energy_mean = np.mean(energy) / N
        magnetization_mean = np.mean(np.abs(magnetization)) / N

        specific_heat = np.var(energy) / (self.T * self.T * N)
        susceptibility = np.var(magnetization) / (self.T * N)

        return energy_mean, magnetization_mean, specific_heat, susceptibility


def plot_results(results: Dict[float, Dict[str, np.ndarray]]):
    T_values = sorted(results.keys())
    avg_energy = [np.mean(results[T]['energy']) for T in T_values]
    avg_magnetization = [np.mean(np.abs(results[T]['magnetization'])) for T in
                         T_values]
    avg_specific_heat = [results[T]['specific_heat'] for T in T_values]
    avg_susceptibility = [results[T]['susceptibility'] for T in T_values]

    plt.figure(figsize=(12, 10))

    plt.subplot(2, 2, 1)
    plt.plot(T_values, avg_energy, marker='o')
    plt.title('Average Energy vs Temperature')
    plt.xlabel('Temperature (T)')
    plt.ylabel('Average Energy')

    plt.subplot(2, 2, 2)
    plt.plot(T_values, avg_magnetization, marker='o', color='orange')
    plt.title('Average Magnetization vs Temperature')
    plt.xlabel('Temperature (T)')
    plt.ylabel('Average Magnetization')

    plt.subplot(2, 2, 3)
    plt.plot(T_values, avg_specific_heat, marker='o', color='green')
    plt.title('Specific Heat vs Temperature')
    plt.xlabel('Temperature (T)')
    plt.ylabel('Specific Heat')

    plt.subplot(2, 2, 4)
    plt.plot(T_values, avg_susceptibility, marker='o', color='red')
    plt.title('Susceptibility vs Temperature')
    plt.xlabel('Temperature (T)')
    plt.ylabel('Susceptibility')

    plt.tight_layout()
    plt.show()


def temperature_scan(L: int, J: float, h: float, T_values: List[float],
                     n_samples: int, n_steps_between: int) -> Dict[
    float, Dict[str, np.ndarray]]:
    results = {}
    for T in T_values:
        model = Model(L=L, J=J, h=h, T=T)
        model.thermalize()

        energy, magnetization, specific_heat, susceptibility = model.sampling_observables(
            n_samples, n_steps_between)
        results[T] = {
            'energy': energy,
            'magnetization': magnetization,
            'specific_heat': specific_heat,
            'susceptibility': susceptibility
        }
    return results


def main(name):
    print(f'Hello, {name}!')
    L = 4
    J = 1.0
    h = 0.0
    T_values = np.linspace(1.0, 4.0, 10)
    n_samples = 100
    n_steps_between = 10

    start = time.perf_counter()
    results = temperature_scan(L=L, J=J, h=h, T_values=T_values,
                               n_samples=n_samples,
                               n_steps_between=n_steps_between)
    end = time.perf_counter()
    print(f'Simulation completed in {end - start:.2f} seconds.')
    plot_results(results)


if __name__ == '__main__':
    main('PyCharm')

 まとめ本記事では、2 次元イジングモデルのモンテカルロシミュレーションを Python で実装しました。メトロポリス法を用いてカノニカル分布に従う状態をサンプリングし、物理量の期待値を近似的に計算しました。シミュレーション結果は厳密計算と比較して良好な一致を示しました。
はじめに

モデル

2 次元イジングモデルのハミルトニアン

平衡状態

物理量の定義

物理量の期待値

厳密計算

モンテカルロシミュレーション

考え方

マルコフ連鎖モンテカルロ法

マルコフ連鎖とは

定義

定常分布

詳細釣り合い条件

メトロポリス法

アルゴリズム

初期化と平衡化

実装例

まとめ

Discussion
