MCMC方法1D铁磁Ising模型

问题描述

我的问题与使用马尔可夫链蒙特卡洛方法（MCMC）的一维Ising模型的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=1) +
    interaction * spins * np.roll(spins,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)

解决方法

在这里，我更改了您的代码以一如既往地解决问题。

请仔细检查代码和公式

#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。

,

我一直在寻找1D Ising模型的简单实现，并且遇到了这篇文章。虽然我不是该领域的专家，但是我确实写了一个相关主题的硕士。我在Oriol Cabanas Tirapu的答案中实现了代码，发现了一些错误（我认为）。

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

const testString = 'This #here is a #test.';
console.log(testString.match(/[#_]?\w+|[^\w#_]+/g));