Нахождение физических величин с помощью Монте-Карло. Модель Изинга

Модель Изинга - это решётка взаимодействующих спинов NxN c периодическими граничными условиями.

\[ E(S)=-J \sum_{i \sim j} S_i S_j \]

Алгоритм Метрополиса в данном случае позволяет генерировать элементы из NVT ансамбля согласно распределению Гиббса.

С помощью этого можем получить средние величины, просто усредняя их по получившейся выборке.

import numpy as np
import matplotlib.pyplot as plt
from numba import njit, prange

@njit()
def E_calc(system: np.array, N: int) -> float: # Подсчёт энергии с учётом ПГУ
    
    E = 0
    
    for i in range(N):
        for j in range(N):
            
            if (i != 0) and (j != 0):
                E -= (system[i-1, j]*system[i, j] + system[i, j]*system[i, j-1])
                continue
            if (i == 0) and (j != 0):
                E -= (system[N-1, j]*system[i, j] + system[i, j]*system[i, j-1])
                continue
            if (i != 0) and (j == 0):
                E -= (system[i-1, j]*system[i, j] + system[i, j]*system[i, N-1])
                continue
            
    return E

@njit()
def izing(N: int,T: float, steps=5000000):
    
    system = np.ones((N, N)) # Инициализируем систему с сонаправленными спинами. Наш x_0
     
    E = E_calc(system, N) # Считаем энергию
    
    for i in range(steps): # Алгоритм Метрополиса (версия Гиббса)
        
        # выбираем случайный элемент
        i1 = np.random.randint(0, N) 
        i2 = np.random.randint(0, N)
        
        # определяем индексы соседей с учётом ПГУ
        i_up = i1 - 1 
        i_down = i1 + 1
        i_left = i2 - 1
        i_right = i2 + 1
        
        if i_up == -1:
            i_up = N-1
        
        if i_down == N:
            i_down = 0
        
        if i_left == -1:
            i_left = N-1
            
        if i_right == N:
            i_right = 0
        
        # смотрим на изменение энергии при изменении значения спина i1, i2
        deltaE = 2*system[i1, i2]*(system[i_up, i2] + system[i_down, i2] + system[i1, i_left] + system[i1, i_right])
            
        if deltaE < 0: # Если энергия стала ниже, то всегда принимаем новую конфигурацию
            system[i1, i2] *= -1 # измененный system наше x_{t+1} состояние
            E += deltaE # отслеживаем значение энергии
        else:
            if np.random.uniform(0,1) < np.exp(-deltaE/T): # Если нет, то используем рпспределение Гиббса как P(x)
                system[i1, i2] *= -1 # измененный system наше x_{t+1} состояние
                E += deltaE
                
    return E, system.mean() # выводим энергию на последнем получившемся x_t и средний спин на состоянии x_t

from tqdm import tqdm

N = 100
T = 10

Ms = []
Es = []
Ts = np.linspace(1, 3, 50)

#avg E from T
for T in tqdm(Ts):
    
    E, M = 0, 0
    
    # Усредняем 100 раз по разным экспериментам, чтобы избавиться от корреляции
    for i in range(100): 
        deltaE, deltaM = izing(N,T, steps=1000000)
        E += deltaE
        M += deltaM
    E /= 100
    M /= 100
    
    Es.append(E)
    Ms.append(M)

100%|███████████████████████████████████████████████████████████████████████████████████| 50/50 [04:29<00:00,  5.39s/it]

# Теоретические выкладки (ищите сами, откуда :))

from scipy.optimize import fsolve
from scipy.special import ellipk
from scipy.special import ellipe

def T_crit(T_cand):
    return np.sinh(2/T_cand) - 1
T_c_theory = fsolve(T_crit, 2)[0]

def E_ud_theory(T):
    
    kappa = 2*np.sinh(2./T)/(np.cosh(2./T)**2)
    E = 0
    E += -2*np.tanh(2./T)
    E += -(np.sinh(2./T)**2 - 1)/np.sinh(2./T)/np.cosh(2./T)*(2./np.pi*ellipk(kappa**2) - 1)

    return E

def C_theory(T):
    
    kappa = 2*np.sinh(2./T)/(np.cosh(2./T)**2)
    
    C = 4./np.pi
    
    C *= (1./T/np.tanh(2./T))**2
    
    C *= (ellipk(kappa**2) - ellipe(kappa**2) - (1 - np.tanh(2./T)**2)*(np.pi/2 + (2*np.tanh(2./T)**2 - 1)*ellipk(kappa**2)))
    return C

Es = np.array(Es)

Es_theory = E_ud_theory(Ts)

plt.plot(Ts, Es/N**2, label = "exp")
plt.plot(Ts, Es_theory, label = "theory")
plt.legend()
plt.title("E_уд(T)")

Text(0.5, 1.0, 'E_уд(T)')

plt.plot(Ts, Ms)
plt.title("M(T)")

Text(0.5, 1.0, 'M(T)')

dT = Ts[1] - Ts[0]
Cvs = (Es[2:]-Es[0:-2])/2/dT/N**2

plt.plot(Ts[1:-1], Cvs)
plt.plot(Ts[1:-1], C_theory(Ts[1:-1]))
plt.title("C(T)")
print(f"Tc_exp = {Ts[Cvs.argmax()]}")
print(f"Tc_theory = {T_c_theory}")

Tc_exp = 2.3061224489795915
Tc_theory = 2.269185314213022