伊辛模型大都会算法

Question

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

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()

你会发现 space 被分成了两个区域。或者，您可以如下初始化网格而不是随机数。

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

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

伊辛模型大都会算法

Ising model metropolis algorithm

python

physics

montecarlo