python - 伊辛模型大都会算法 : lattice won't equilibrate

Question

我在 python (2d) 中有一些 Ising 模型的代码，晶格不会达到平衡。在这里，代码打印出每次 Monte Carlo 扫描翻转的自旋数，并且每次扫描都翻转相同的数字。如果我是正确的，那么随着晶格达到平衡，翻转的次数应该随着每次扫描而减少。谁能看到代码有任何错误？

import numpy as np
from numpy import random as rn
N=20
a=[1,-1]
b=[N,N]

#first make an array
init_lattice=rn.choice(a,size=(N,N))

#specify how many Monte Carlo sweeps
number_of_sweeps=1000

#main code to flip spins and record how many are flipped
def new_lattice(lattice,T):
    delta_E=np.zeros(b)
    switch1=np.zeros(number_of_sweeps)
    switch=0
    for sweep in range(number_of_sweeps):
        for i in range(N):
            for j in range(N):
                Si=lattice[i,j]
                sum_Sj=lattice[i,(j+1)%N]+lattice[(i+1)%N,j]+lattice[i,(j-1)%N]+lattice[(i-1)%N,j]
                delta_E[i,j]=2*Si*sum_Sj
                if delta_E[i,j]<0:
                    lattice[i,j]*=-1
                    switch+=1
                elif np.exp(-1*delta_E[i,j]/(T))>rn.random():
                    lattice[i,j]*=-1
                    switch+=1
        switch1[sweep]=switch
    return lattice,switch1

#print how many spins have flipped
switch= new_lattice(init_lattice,5)[1]
print switch
for i in range(len(switch)):
    print switch[i+1]-switch[i-1]

Answer 1

我认为循环

    for i in range(N):
        for j in range(N):

不好:您必须随机选择要测试的旋转。

T 和其他变量是整数，您将在指数 np.exp(-1*delta_E[i,j]/T) 中进行整数除法，这是错误的。

使用这段代码我有一个饱和度(注意我取 T=1. 而不是 5. 并且 sweeps = 10000 而不是 1000):

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

latticeDimensions = [20, 20]
nx = latticeDimensions[0]
ny = latticeDimensions[1]
nodesNumber = nx*ny
spinValues = [-1,1]
sweeps = 10000
T = 1.

lattice = init_lattice = rn.choice(spinValues,size=(nx,ny))
#lattice = [ random.choice([-1,1]) for i in range(nodesNumber) ]
t = 0
switch = 0
res = []
while t < sweeps:
    selectedNode = random.choice(nodeIdx)
    i = random.choice(range(nx))
    j = random.choice(range(ny))
    Si = lattice[i,j]
    sum_Sj=lattice[i,(j+1)%ny] + lattice[(i+1)%nx,j] + lattice[i,(j-1)%ny] + lattice[(i-1)%nx,j]
    delta = 2 * Si * sum_Sj
    probaTemp = np.exp(- delta / T)
    if delta > 0. and rn.random() < probaTemp:
        lattice[i,j] *= -1
        switch += 1
    elif delta <= 0.:
        lattice[i,j] *= -1
        switch += 1
    t += 1
    print t
    res.append(switch)

plt.plot(res)
plt.show()