<p>在这里,我改变了你的代码,以我一贯的方式解决了这个问题</p>
<p>请仔细检查代码和公式</p>
<pre><code>#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
</code></pre>
<p>在模拟结束时,变量磁化和能量是包含每个MC步骤的测量值的阵列。
您可以直接使用这些数组来计算直方图和绘图</p>
<p>注:能量阵列是系统的总能量,而不是能量/L</p>