kawasaki_test.py

kawasaki_test.py

import random
from math import exp

def energy(spins):
    if spins == [0,1,1,0] or spins == [1,0,0,1]:
        return 16.0
    else:
        return 0.0

def random_pair(_):
    """
    隣接するとは限らないスピンペアをランダムに選ぶ
    """
    return random.sample([0, 1, 2, 3], 2)

def random_neighbor_pair(_):
    """
    隣接するスピンペアをランダムに選ぶ
    """
    i = random.randint(0,3)
    x = random.choice([1,-1,2,-2])
    j = (i + x + 4) %4
    return i, j

def random_opposite_neighbor_pair(spins):
    """
    隣接するスピンペアのうち、逆向きのものをランダムに選ぶ
    """
    while True:
        i = random.randint(0,3)
        x = random.choice([1,-1,2,-2])
        j = (i + x + 4) %4
        if spins[i] != spins[j]:
            return i, j

def swap(spins, beta, swap_function):
    i, j = swap_function(spins)
    ns = spins.copy()
    ns[i], ns[j] = ns[j], ns[i]
    de = energy(ns) - energy(spins)
    if (de < 0.0 or exp(-de*beta) >= random.random()):
        spins[i], spins[j] = spins[j], spins[i]

def calc_energy(swap_function, beta):
    spins = [0,0,1,1]
    total_step = 1000000
    energy_sum = 0.0
    beta = 0.01
    for _ in range(total_step):
        swap(spins, beta, swap_function)
        energy_sum += energy(spins)

    # エネルギーの推定値
    energy_sum /= total_step
    # 厳密解
    U = 16.0/(1.0 + 2.0*exp(16.0*beta))
    print(f"推定値：{energy_sum} 厳密解{U}")

if __name__ == '__main__':
    beta = 0.1
    print(f"逆温度: {beta}")
    print("隣接するとは限らないスピンペアをランダムに選ぶ")
    calc_energy(random_pair, beta)
    print()
    print("隣接するスピンペアをランダムに選ぶ")
    calc_energy(random_neighbor_pair, beta)
    print()
    print("隣接するスピンペアのうち、逆向きのものをランダムに選ぶ")
    calc_energy(random_opposite_neighbor_pair, beta)
    

実行結果

逆温度: 0.1
隣接するとは限らないスピンペアをランダムに選ぶ
推定値：4.798016 厳密解4.780369305434378

隣接するスピンペアをランダムに選ぶ
推定値：4.77664 厳密解4.780369305434378

隣接するスピンペアのうち、逆向きのものをランダムに選ぶ
推定値：5.41848 厳密解4.780369305434378