python - 一维 Ising 模型图和配置

Question

所以我对我编写的这个 Python 代码有疑问。它看起来不错，我没有错误并且它运行。但是，它并没有停止运行，这让我认为它必须卡在某种无限循环上。当我停止运行时，我每次都会收到最新的调用 elif rand() < np.exp(-cost*beta):。所以我认为这就是问题所在。但是，我找不到任何问题。任何帮助将不胜感激，我的代码正在模拟和绘制量子和统计力学中的一维伊辛链模型。

import matplotlib.pyplot as plt
import numpy as np
from numpy.random import rand

#The 1D Ising Model
#MC (or mc) is just short for Monte Carlo throughout the code
#This code will look the same as the code for the 2D Ising model code, just changing variables to get a result that is 2-dimensional
def initialstate(N): #This will generate a random spin configuration for initial condition as mentioned in section 4.11
    state = 2*np.random.randint(2, size=(N))-1
    return state

def mcmove(config, beta, N): #This is the Monte Carlo move using Metropolis algorithm talked about in the textbook section 4.11
    for i in range(N):
        for j in range(N):
            a = np.random.randint(0, N)
            s = config[a]
            nb = config[(a+1)%N] + config[(a-1)%N] #%N allows for toroidal boundary conditions
            cost = 2*s*nb
            if cost < 0:
                s *= -1
            elif rand() < np.exp(-cost*beta):
                config[a] = s
    return config
    print(config)


def calculateEnergy(config, N): #This is the energy of a given configuration
    energy = 0
    for i in range(len(config)):
        for j in range(len(config)):
            S = config[i]
            nb = config[(i+1)%N] + config[(i-1)%N]
            energy += -nb * S
    return energy/2 #divides the resultinng energy by 4

def calculateMagnetization(config): #This is the magnetization of the given configuration
    magnetization = np.sum(config) #The magnetization is just a sum
    return magnetization

def main():

#The following will be the parameters
    N_T = 100 #This is the number of temperature points
    N = 200 #This will be the size of the lattice
    equilsteps = 1000 #Number of Monte Carlo sweeps for equilibration
    mcsteps = 1000 #Number of Monte Carlo sweeps for calculation

    T = np.linspace(1.01, 3.28, N_T);
    E,M,C,X = np.zeros(N_T), np.zeros(N_T), np.zeros(N_T), np.zeros(N_T) #Energy (E), Magnetization (M), Specific Heat (C), and Susceptibility (X)
    n_1, n_2 = (1.0)/(mcsteps*N), (1.0)/(mcsteps*mcsteps*N)

    for tt in range(N_T):
        E_1 = M_1 = E_2 = M_2 = 0
        config = initialstate(N)
        i_T = 1.0/T[tt]; i_T_2 = i_T*i_T;

        for i in range(equilsteps): #This is the equilibrate
            mcmove(config, i_T, N)
            Energy = calculateEnergy(config, N) #defined above, which will calculate the energy
            Magnet = calculateMagnetization(config) #also defined above, will calculate the magnetization

            M_1 = M_1 + Magnet
            E_1 = E_1 + Energy
            M_2 = M_2 + Magnet*Magnet
            E_2 = E_2 + Energy*Energy

        E[tt] = n_1*E_1 #Energy
        M[tt] = n_1*M_1 #Magnetization
        C[tt] = (n_1*E_2 - n_2*E_1*E_1)*i_T_2 #Specific Heat
        X[tt] = (n_1*M_2 - n_2*M_1*M_1)*i_T #Susceptibility

    f = plt.figure(figsize = (18, 10)); #This will plot the calculated values

#This will be for the Energy
    f.add_subplot(2, 2, 1);
    plt.scatter(T, E, s=50, marker = 'o', color = 'g') #This will plot the energy values vs temperature on a scatter plot
    plt.xlabel("Temperature (kT/J)", fontsize = 15); #This labels the x-axis
    plt.ylabel("Energy (E/NJ)", fontsize = 15); #This labels the y-axis
    plt.axis('tight');

#This will be for the Magnetization
    f.add_subplot(2, 2, 2);
    plt.scatter(T, M, s=50, marker = 'o', color = 'r') #This will plot the magnetization values vs temperature on a scatter plot
    plt.xlabel("Temperature (kT/J)", fontsize = 15); #This labels the x-axis
    plt.ylabel("Magnetization (M(T))", fontsize = 15); #This labels the y-axis
    plt.axis('tight');

#This will be for the Specific Heat
    f.add_subplot(2, 2, 3);
    plt.scatter(T, C, s=50, marker = 'o', color = 'g') #This will plot the specific heat values vs temperature on a scatter plot
    plt.xlabel("Temperature (kT/J)", fontsize = 15); #This labels the x-axis
    plt.ylabel("Specific Heat (C/Nk)", fontsize = 15); #This labels the y-axis
    plt.axis('tight');

#This will be for the Susceptibility
    f.add_subplot(2, 2, 4);
    plt.scatter(T, X, s=50, marker = 'o', color = 'r') #This will plot the susceptibility values vs temperature on a scatter plot
    plt.xlabel("Temperature (kT/J)", fontsize = 15); #This labels the x-axis
    plt.ylabel("Susceptibility (X)", fontsize = 15); #This labels the y-axis
    plt.axis('tight');

    plt.show()

main()

标签： pythonnumpy