MCMC方法一维铁磁伊辛模型

Question

我的问题与使用马尔可夫链蒙特卡罗方法（MCMC）对一维伊辛模型进行Python编码有关

我有下面的哈密顿量

$$H = - \sum_{i=1}^{L-1}\sigma_{i}sigma_{i+1} - B\sum_{i=1}^{L}\sigma_{i}$$

我想编写一个python函数，生成一个马尔可夫链，在每个步骤中，它计算并保存磁化（每个站点）和能量

能量是（=哈密顿量），我将磁化定义为：

$$\frac{1}{L}\sum_{i}\sigma_{i}$$

我的概率分布是：

$$p(x) = e^{-H\beta}$$ where, $T^{-1} = \beta$

对于马尔可夫链，我将实现Metropolis-Hastings算法

if $$\frac{P(\sigma')}{P(\sigma)} = e^{(H(\sigma)-H(\sigma'))\beta}$$

我的想法是在必要时接受转换

$$H(\sigma') < H(\sigma)$$

并且只接受转换

$$H(\sigma') > H(\sigma)$$

有可能

$$P = e^{(H(\sigma)-H(\sigma'))\beta}$$

因此，让我设置一些参数，例如：

$L=20$ - Lattice Size

$T=2$ - Temperature

$B=0$ - Magnetic Field

在计算之后，我需要绘制一个磁化强度和能量与步长的直方图。我对这部分没有异议

我的python知识不是很好，但我已经包括了我的粗略（未完成）草稿。我认为我没有取得多大进展。任何帮助都会很好

#Coding attempt MCMC 1-Dimensional Ising Model
import numpy as np
import matplotlib.pyplot as plt

#Shape of Lattice L 
L = 20
Shape = (20,20)

#Spin Configuration
spins = np.random.choice([-1,1],Shape)

#Magnetic moment
moment = 1

#External magnetic field
field = np.full(Shape, 0)

#Temperature 
Temperature = 2
Beta = Temperature**(-1)

#Interaction (ferromagnetic if positive, antiferromagnetic if negative)
interaction = 1

#Using Probability Distribution given
def get_probability(Energy1, Energy2, Temperature):
  return np.exp((Energy1 - Energy2) / Temperature)

def get_energy(spins):
    return -np.sum(
    interaction * spins * np.roll(spins, 1, axis=0) +
    interaction * spins * np.roll(spins, -1, axis=0) +
    interaction * spins * np.roll(spins, 1, axis=1) +
    interaction * spins * np.roll(spins, -1, axis=1)
  )/2 - moment * np.sum(field * spins)

#Introducing Metropolis Hastings Algorithim
x_now = np.random.uniform(-1, 1)    #initial value
d = 10**(-1)                               #delta
y = []
for i in range(L-1):
      #generating next value
      x_proposed = np.random.uniform(x_now - d, x_now + d)
      #accepting or rejecting the value
      if np.random.rand() <  np.exp(-np.abs(x_proposed))/(np.exp(-np.abs(x_now))):
          x_now = x_proposed
      if i % 100 == 0:
          y.append(x_proposed)

Answer 1

在这里，我改变了你的代码，以我一贯的方式解决了这个问题

请仔细检查代码和公式

#Coding attempt MCMC 1-Dimensional Ising Model
import numpy as np
import matplotlib.pyplot as plt

#Shape of Lattice L 
L = 20
#Shape = (20)

#Number of Monte Carlo samples

MC_samples=1000

#Spin Configuration
spins = np.random.choice([-1,1],L)
print(spins)
#Magnetic moment
moment = 1

#External magnetic field
field = 0

#Temperature 
Temperature = 2
Beta = Temperature**(-1)

#Interaction (ferromagnetic if positive, antiferromagnetic if negative)
interaction = 1

#Using Probability Distribution given
def get_probability(delta_energy, Temperature):
  return np.exp(-delta_energy / Temperature)

def get_energy(spins):
    energy=0
    for i in range(L):
        energy=energy+interaction*spins[i-1]*spins[i]
    energy= energy-field*sum(spins)
    return energy
def delta_energy(spins,random_spin):
    #If you do flip one random spin, the change in energy is:
    #(By using a reduced formula that only involves the spin
    # and its neighbours)
    if random_spin==L:
        PBC=0
    else:
        PBC=random_spin+1
    
    return -2*interaction*(spins[random_spin-1]*spins[random_spin]+
              spins[random_spin]*spins[PBC]+field*spins[random_spin])

#Introducing Metropolis Hastings Algorithim
#x_now = np.random.uniform(-1, 1)    #initial value
#d = 10**(-1)                               #delta
#y = []

magnetization=[]
energy=[]
for i in range(MC_samples):
    #Each Monte Carlo step consists in L random spin moves
    for j in range(L):
        #Choosing a random spin
        random_spin=np.random.randint(L-1,size=(1))
        #Compuing the change in energy of this spin flip
        delta=delta_energy(spins,random_spin)
        
        #Metropolis accept-rejection:
        if delta<0:
            #Accept the move if its negative
            spins[random_spin]=-spins[random_spin]
        else:
            #If its positive, we compute the probability
            probability=get_probability(delta,Temperature)
            random=np.random.rand()
            if random<=probability:
                #Accept de move
                spins[random_spin]=-spins[random_spin]
                
                
    #Afer the MC step, we measure the system
    magnetization.append(sum(spins)/L)
    energy.append(get_energy(spins))
print(magnetization,energy)

#Do histograms and plots

在模拟结束时，变量磁化和能量是包含每个MC步骤的测量值的阵列。 您可以直接使用这些数组来计算直方图和绘图

注：能量阵列是系统的总能量，而不是能量/L

Answer 2

我正在寻找一个1D伊辛模型的简单实现，偶然发现了这篇文章。虽然我不是这个领域的专家，但我确实写了一篇相关主题的硕士论文。我在Oriol Cabanas Tirapu的答案中实现了代码，并发现了一些bug（我想）

下面是我的改编版哦他们的代码。希望它对某些人有用

#Coding attempt MCMC 1-Dimensional Ising Model
import numpy as np
import matplotlib.pyplot as plt

#Using Probability Distribution given
def get_probability(delta_energy, Temperature):
    return np.exp(-delta_energy / Temperature)

def get_energy(spins):
    energy=0
    for i in range(len(spins)):
        energy=energy+interaction*spins[i-1]*spins[i]
    energy= energy-field*sum(spins)
    return energy

def delta_energy(spins,random_spin):
    #If you do flip one random spin, the change in energy is:
    #(By using a reduced formula that only involves the spin
    # and its neighbours)
    if random_spin==L-1:
        PBC=0
    else:
        PBC=random_spin+1
        
    old = -interaction*(spins[random_spin-1]*spins[random_spin] + spins[random_spin]*spins[PBC]) - field*spins[random_spin]
    new = interaction*(spins[random_spin-1]*spins[random_spin] + spins[random_spin]*spins[PBC]) + field*spins[random_spin]
    
    return new-old

def metropolis(L = 100, MC_samples=1000, Temperature = 1, interaction = 1, field = 0):
    
    # intializing
    #Spin Configuration
    spins = np.random.choice([-1,1],L)
    
    Beta = Temperature**(-1)
    
    #Introducing Metropolis Hastings Algorithim
    data = []
    magnetization=[]
    energy=[]
    for i in range(MC_samples):
        #Each Monte Carlo step consists in L random spin moves
        for j in range(L):
            #Choosing a random spin
            random_spin=np.random.randint(0,L,size=(1))
            #Compuing the change in energy of this spin flip
            delta=delta_energy(spins,random_spin)

            #Metropolis accept-rejection:
            if delta<0:
                #Accept the move if its negative
                spins[random_spin]=-spins[random_spin]
                #print('change')
            else:
                #If its positive, we compute the probability
                probability=get_probability(delta,Temperature)
                random=np.random.rand()
                if random<=probability:
                    #Accept de move
                    spins[random_spin]=-spins[random_spin]

        data.append(list(spins))

        #Afer the MC step, we measure the system
        magnetization.append(sum(spins)/L)
        energy.append(get_energy(spins))
        
    return data,magnetization,energy
    
    
def record_state_statistics(data,n=4):
    ixs = tuple()
    
    sub_sample = [[d[i] for i in range(n)] for d in data]
    
    # get state number
    state_nums = [int(sum([((j+1)/2)*2**i for j,i in zip(reversed(d),range(len(d)))])) for d in sub_sample]
    
    return state_nums



# setting up problem
L = 200 # size of system
MC_samples = 1000 # number of samples
Temperature = 1 # "temperature" parameter
interaction = 1 # Strength of interaction between nearest neighbours
field = 0 # external field

# running MCMC
data = metropolis(L = L, MC_samples = MC_samples, Temperature = Temperature, interaction = interaction, field = field)
results = record_state_statistics(data[0],n=4) # I was also interested in the probability of each micro-state in a sub-section of the system

# Plotting
plt.figure(figsize=(20,10))

plt.subplot(2,1,1)
plt.imshow(np.transpose(data[0]))
plt.xticks([])
plt.yticks([])
plt.axis('tight')
plt.ylabel('Space',fontdict={'size':20})
plt.title('Critical dynamics in a 1-D Ising model',fontdict={'size':20})


plt.subplot(2,1,2)
plt.plot(data[2],'r')
plt.xlim((0,MC_samples))
plt.xticks([])
plt.yticks([])
plt.ylabel('Energy',fontdict={'size':20})
plt.xlabel('Time',fontdict={'size':20})