Алгоритм Метрополиса модели Изинга: решетка не уравновешивается

У меня есть некоторый код для модели Изинга в python (2d), и решетка не достигнет равновесия. Здесь код выводит число вращений, которые переворачиваются для каждой развертки Монте-Карло, и то же самое число переворачивается для каждой развертки. Если я прав, то число, которое переворачивается, должно уменьшаться с каждой разверткой, поскольку решетка достигает равновесия. Может ли кто-нибудь увидеть какие-либо ошибки в коде?



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]
516   1  

1 ответ:

Я думаю, что петли 

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

Не годятся: вы должны выбрать спин, который хотите протестировать случайным образом.

T и другие переменные целочисленны, и у вас будет целочисленное деление в экспоненциальном np.exp(-1 * delta_E[i, j] / T), что неверно.

С этим кодом у меня есть насыщение (заметьте, я беру T=1. вместо 5. и метет = 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()
1  
2017-04-04 10:30:15

Comments

    Ничего не найдено.
Click here to cancel reply.