Python 中的伊辛模型

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=0.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 < 0），所以即使在非常低的温度下也没有有序状态。

然而，最严重的错误是使用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 算法无论如何都无法正确处理那里的关键减速）。