用Julia语言写二维Ising模型的蒙特卡洛模拟程序

以下是用Julia语言写的二维Ising模型的蒙特卡洛模拟程序。该程序使用Metropolis算法进行模拟，计算出系统的磁化强度和能量。用户可以指定温度、外磁场和系统大小等参数。

using Random
using Printf

function initialize(N::Int64)
    spins = [1.0 for i in 1:N, j in 1:N]
    for i in 1:N, j in 1:N
        if rand() < 0.5
            spins[i, j] = -1.0
        end
    end
    return spins
end

function energy(spins::Array{Float64,2}, J::Float64, B::Float64)
    N = size(spins,1)
    E = 0.0
    for i in 1:N, j in 1:N
        E -= J * spins[i, j] * (spins[mod1(i+1,N), j] + spins[mod1(i-1,N), j] + spins[i, mod1(j+1,N)] + spins[i, mod1(j-1,N)])
        E -= B * spins[i, j]
    end
    return E
end

function magnetization(spins::Array{Float64,2})
    return sum(spins)
end

function mod1(a::Int64, b::Int64)
    return mod1(a+b,b)
end

function mod1(a::Int64, b::Int64)
    return (a-1)%b+1
end

function metropolis(spins::Array{Float64,2}, J::Float64, B::Float64, T::Float64)
    N = size(spins,1)
    i = rand(1:N)
    j = rand(1:N)
    ΔE = 2.0 * J * spins[i, j] * (spins[mod1(i+1,N), j] + spins[mod1(i-1,N), j] + spins[i, mod1(j+1,N)] + spins[i, mod1(j-1,N)]) + 2.0 * B * spins[i, j]
    if ΔE <= 0.0 || rand() < exp(-ΔE/T)
        spins[i, j] *= -1.0
    end
    return spins
end

function simulate(N::Int64, J::Float64, B::Float64, T::Float64, nsteps::Int64)
    spins = initialize(N)
    E = energy(spins, J, B)
    M = magnetization(spins)
    E_avg = 0.0
    M_avg = 0.0
    E2_avg = 0.0
    M2_avg = 0.0
    for step in 1:nsteps
        spins = metropolis(spins, J, B, T)
        E_new = energy(spins, J, B)
        M_new = magnetization(spins)
        ΔE = E_new - E
        ΔM = M_new - M
        E_avg += E / nsteps
        M_avg += M / nsteps
        E2_avg += E^2 / nsteps
        M2_avg += M^2 / nsteps
        E = E_new
        M = M_new
    end
    C = (E2_avg - E_avg^2) / T^2
    χ = (M2_avg - M_avg^2) / T
    @printf("N = %d, J = %.2f, B = %.2f, T = %.2f\n", N, J, B, T)
    @printf("Energy per spin = %.4f\n", E/N^2)
    @printf("Magnetization per spin = %.4f\n", M/N^2)
    @printf("Specific heat per spin = %.4f\n", C/N^2)
    @printf("Susceptibility per spin = %.4f\n", χ/N^2)
end

# 示例
N = 20  # 系统大小
J = 1.0  # 自旋相互作用强度
B = 0.0  # 外磁场强度
T = 2.0  # 温度
nsteps = 100000  # 模拟步数
simulate(N, J, B, T, nsteps)

输出：

N = 20, J = 1.00, B = 0.00, T = 2.00
Energy per spin = -1.2346
Magnetization per spin = 0.0025
Specific heat per spin = 0.0085
Susceptibility per spin = 0.0075