Ising模型都会算法

Question

我正在尝试使用都会算法来模拟伊辛模型。我遇到的问题是代码不会一直稳定。对于高 beta 值，能量偏好应该是自旋都在同一方向，但是在我的代码中，即使它在大部分时间都有效，偶尔会出现水平或垂直的向上旋转带，然后向下旋转网格（网格2）。然后，这使得平均旋转计数稳定在一个不是 1 或 -1 的值，而高 Beta 值应为 1 或 -1。我知道问题出在哪里 - 当这种情况发生时，我在边界处的 delta 能量是 4，因此在高 beta 值下，这使得自旋不太可能翻转，因此最终会被卡住。我在想，也许我运行的时间不够长，但查看变量资源管理器上的average_spin_count，它似乎在最多 200 圈后稳定下来。我也在想也许使用对角线，但其他代码和模拟似乎没有使用这个想法，所以我希望有人能指出我哪里出错了。 （您必须多次运行代码，直到遇到并非所有旋转都向上/向下的情况）。谢谢

import numpy as np
import numpy.random as rnd

N = 20
#grid size
random_number = np.random.randint(0,1000)
prob = 0.5 
#probability of starting with spin down.

np.random.seed(random_number)
grid =np.random.choice(list(range(-1, 0)) + list(range(1, 2)), size=(N, N), p=[prob,1-prob])



def find_energy(grid, r, c):

    if r < N-1 and c < N-1:
        energy = -((grid[r,c]*grid[r+1,c])+(grid[r,c]*grid[r-1,c])+
                    (grid[r,c]*grid[r,c+1])+(grid[r,c]*grid[r,c-1]))

    if r == (N-1) and c < N-1:
        energy = -((grid[r,c]*grid[0,c])+(grid[r,c]*grid[r-1,c])+
                    (grid[r,c]*grid[r,c+1])+(grid[r,c]*grid[r,c-1]))
    

    if c == (N-1) and r < N-1:
            energy = -((grid[r,c]*grid[r+1,c])+(grid[r,c]*grid[r-1,c])+
                      (grid[r,c]*grid[r,0])+(grid[r,c]*grid[r,c-1]))
    if c== (N-1) and r==(N-1):
            energy = -((grid[r,c]*grid[0,c])+(grid[r,c]*grid[r-1,c])+
                    (grid[r,c]*grid[r,0])+(grid[r,c]*grid[r,c-1]))
    
    return energy


def get_coordinate():
    """Returns a random coordinate as a tuple"""
    return tuple(np.random.randint(0, high=N, size=(1, 2))[0])


def probability( delta_energy, beta):

    probability = np.exp(-1*beta*delta_energy)

    return probability

def compare_energies(n, beta):

    total_number_spin_up=np.array([])
    total_number_spin_down=np.array([])
    for i in range(n):
        original_point = get_coordinate()
        original_energy = find_energy(grid, original_point[0], original_point[1])
        new_energy = original_energy * -1
        delta_energy = new_energy-original_energy
        x = (original_point[0])
        y = (original_point[1])
        if delta_energy <= 0:
              grid[x,y] = grid[x, y]*-1
              #flips the spin
        elif delta_energy > 0:
            if probability(delta_energy, beta) > rnd.uniform(0,1):
                grid[x,y] = (grid[x, y])*-1  
        number_of_spin_up = np.count_nonzero(grid == 1)
        total_number_spin_up = np.append(total_number_spin_up,number_of_spin_up)
        total_number_spin_down = np.append(total_number_spin_down, pow(N,2)-number_of_spin_up)
    spin_up = total_number_spin_up[-1]
    spin_down = total_number_spin_down[-1]
    average_spin = (spin_up-spin_down)/N**2

    return grid, average_spin



beta = 10

#runs the code to settle it.
average_spin_count = np.array([])
for i in range(300):
    grid_2, average_spin = compare_energies(N**2, beta)
    average_spin_count = np.append(average_spin_count, average_spin)
    y = average_spin_count[-1]

Answer 1

如果您绘制这种情况的最终网格，您将看到会发生什么，例如使用matplotlib，如下所示：

import matplotlib.pyplot as plt
plt.imshow(grid)
plt.show()

您会发现空间被分成两个区域。或者，您可以按如下方式初始化网格，而不是随机数。

grid = np.zeros((N, N))
grid[:N//2] = 1

您会看到类似的行为。解决方案非常简单。如果计算边界处自旋翻转的能量和相应概率，您会发现能量值为 -3，导致概率为 9.76e-27，基本上是 0< /code> 对于您给定的测试版。因此，这种状态非常稳定，在 300 次迭代内您可能不会看到任何变化。