二维Ising模型的蒙特卡洛模拟程序

以下是一个简单的二维Ising模型的蒙特卡洛模拟程序，使用Metropolis算法进行更新。该程序使用周期性边界条件，并计算了磁化、内能和比热的平均值。

import numpy as np
import matplotlib.pyplot as plt

# 定义常量
J = 1  # 交换耦合常数
kB = 1  # 玻尔兹曼常数
T = 1  # 温度
N = 20  # 系统大小
nsteps = 10000  # 迭代次数

# 初始化随机磁化和格点
magnetization = np.random.choice([-1, 1], size=(N, N))
lattice = np.random.choice([-1, 1], size=(N, N))

# 计算能量
def energy(lattice):
    E = 0
    for i in range(N):
        for j in range(N):
            neighbor_sum = lattice[(i+1)%N,j] + lattice[(i-1)%N,j] + \
                           lattice[i,(j+1)%N] + lattice[i,(j-1)%N]
            E += -J * lattice[i,j] * neighbor_sum
    return E

# 定义Metropolis算法
def metropolis(lattice, magnetization):
    i = np.random.randint(N)
    j = np.random.randint(N)
    s = lattice[i, j]
    neighbor_sum = lattice[(i+1)%N,j] + lattice[(i-1)%N,j] + \
                   lattice[i,(j+1)%N] + lattice[i,(j-1)%N]
    delta_E = 2 * J * s * neighbor_sum
    if delta_E <= 0:
        s *= -1
        magnetization += 2 * s
        lattice[i, j] = s
    elif np.exp(-delta_E / (kB * T)) > np.random.rand():
        s *= -1
        magnetization += 2 * s
        lattice[i, j] = s

# 进行模拟
E_list = []
M_list = []
Cv_list = []
for step in range(nsteps):
    metropolis(lattice, magnetization)
    E = energy(lattice)
    M = np.sum(magnetization)
    E_list.append(E)
    M_list.append(M)
    if step % 100 == 0:
        Cv_list.append((np.mean(np.array(E_list[-100:])**2) - np.mean(np.array(E_list[-100:]))**2) / (kB * T**2))

# 计算平均磁化、内能和比热
M_avg = np.mean(np.abs(np.array(M_list)))
E_avg = np.mean(np.array(E_list))
Cv_avg = np.mean(np.array(Cv_list))

# 绘制磁化、内能和比热随迭代次数的变化
plt.figure(figsize=(8, 6))
plt.subplot(3, 1, 1)
plt.plot(M_list)
plt.ylabel('magnetization')
plt.subplot(3, 1, 2)
plt.plot(E_list)
plt.ylabel('energy')
plt.subplot(3, 1, 3)
plt.plot(Cv_list)
plt.ylabel('heat capacity')
plt.xlabel('iteration')
plt.show()

print(f'Magnetization: {M_avg}')
print(f'Energy: {E_avg}')
print(f'Heat capacity: {Cv_avg}')

运行程序后，可以得到磁化、内能和比热随迭代次数的变化图像，以及平均磁化、内能和比热的值。