Python 中的 Ising 模型

Question

我目前正在使用 Python3 为 Ising 模型编写代码。我对编码还是很陌生。我有工作代码，但输出结果与预期不符，我似乎找不到错误。这是我的代码：

import numpy as np
import random


def init_spin_array(rows, cols):
    return np.random.choice((-1, 1), size=(rows, cols))


def find_neighbors(spin_array, lattice, x, y):
    left = (x , y - 1)
    right = (x, y + 1 if y + 1 < (lattice - 1) else 0)
    top = (x - 1, y)
    bottom = (x + 1 if x + 1 < (lattice - 1) else 0, y)

    return [spin_array[left[0], left[1]],
            spin_array[right[0], right[1]],
            spin_array[top[0], top[1]],
            spin_array[bottom[0], bottom[1]]]

def energy(spin_array, lattice, x ,y):
    return -1 * spin_array[x, y] * sum(find_neighbors(spin_array, lattice, x, y))


def main():
    lattice = eval(input("Enter lattice size: "))
    temperature = eval(input("Enter the temperature: "))
    sweeps = eval(input("Enter the number of Monte Carlo Sweeps: "))
    spin_array = init_spin_array(lattice, lattice)
    print("Original System: \n", spin_array)
    # the Monte Carlo follows below
    for sweep in range(sweeps):
        for i in range(lattice):
            for j in range(lattice):
                e = energy(spin_array, lattice, i, j)
                if e <= 0:
                spin_array[i, j] *= -1
            elif np.exp(-1 * e/temperature) > random.randint(0, 1):
                spin_array[i, j] *= -1
            else:
                continue
print("Modified System: \n", spin_array)

main()

我认为错误出在蒙特卡洛循环中，但我不确定。该系统应在低温下高度有序，并在超过 2.27 的临界温度时变得无序。换句话说，系统的随机性应该随着 T 接近 2.27 而增加。例如，在 T=.1 时，我们应该看到对齐的大片自旋，即 -1 和 1 的片。 2.27以后系统应该是乱的，我们应该看不到这些补丁。

Answer 1

如果您将系统大小、扫描次数和平均磁化率包括在内，您的问题会更有意义。中间配置有多少是有序的，多少是无序的？ MC 是一种采样技术 - 单独的配置没有任何意义，在低温下可能（并且将会）出现无序状态，而在高 T 下可能会出现有序状态。有意义的是组装特性（平均磁化强度）。

无论如何，您的代码中存在三个错误：小错误、中错误和非常严重的错误。

一个小问题是你在find_neighbors中搜索邻居时忽略了整行整列：

right = (x, y + 1 if y + 1 < (lattice - 1) else 0)

应该是：

right = (x, y + 1 if y + 1 < lattice else 0)

甚至更好：

right = (x, (y + 1) % lattice)

同样适用于bottom。

中间的一个是你对能量差的计算是两倍：

def energy(spin_array, lattice, x ,y):
   return -1 * spin_array[x, y] * sum(find_neighbors(spin_array, lattice, x, y))
          ^^

因子实际上是2*J，其中J是耦合常数，因此有-1意味着：

1. 您的临界温度减半，更重要的是...
2. 你有反铁磁自旋相互作用（J

然而，最严重的错误是使用random.randint() 进行拒绝抽样：

elif np.exp(-1 * e/temperature) > random.randint(0, 1):
    spin_array[i, j] *= -1

您应该改用random.random()，否则转换概率将始终为 50%。

这是对您的程序的修改，它会自动扫描从 0.1 到 5.0 的温度区域：

import numpy as np
import random


def init_spin_array(rows, cols):
    return np.ones((rows, cols))


def find_neighbors(spin_array, lattice, x, y):
    left   = (x, y - 1)
    right  = (x, (y + 1) % lattice)
    top    = (x - 1, y)
    bottom = ((x + 1) % lattice, y)

    return [spin_array[left[0], left[1]],
            spin_array[right[0], right[1]],
            spin_array[top[0], top[1]],
            spin_array[bottom[0], bottom[1]]]


def energy(spin_array, lattice, x ,y):
    return 2 * spin_array[x, y] * sum(find_neighbors(spin_array, lattice, x, y))


def main():
    RELAX_SWEEPS = 50
    lattice = eval(input("Enter lattice size: "))
    sweeps = eval(input("Enter the number of Monte Carlo Sweeps: "))
    for temperature in np.arange(0.1, 5.0, 0.1):
        spin_array = init_spin_array(lattice, lattice)
        # the Monte Carlo follows below
        mag = np.zeros(sweeps + RELAX_SWEEPS)
        for sweep in range(sweeps + RELAX_SWEEPS):
            for i in range(lattice):
                for j in range(lattice):
                    e = energy(spin_array, lattice, i, j)
                    if e <= 0:
                        spin_array[i, j] *= -1
                    elif np.exp((-1.0 * e)/temperature) > random.random():
                        spin_array[i, j] *= -1
            mag[sweep] = abs(sum(sum(spin_array))) / (lattice ** 2)
        print(temperature, sum(mag[RELAX_SWEEPS:]) / sweeps)


main()

20x20 和 100x100 格子和 100 次扫描的结果：

起始配置是完全有序的配置，以防止形成在低温下非常稳定的畴壁。此外，最初执行 30 次额外扫描以使系统热化（在接近临界温度时还不够，但 Metropolis-Hastings 算法无论如何都无法正确处理那里的临界减速）。