import numpy as np
import matplotlib.pyplot as plt
beta = 1 # defines the temperature
N = 1000 # number of particles
e = np.zeros(N) # size N array, Nth element is the energy of Nth particle
avg_e = [] # array to track the average energy throughout the simulation
for i in range(100000):
p = np.random.randint(0,N) # select random particle of interest
move = np.random.randint(0,2) # trial move, coin toss: if 1, consider adding 1 energy. if 0, consider adding -1 energy
prob = np.exp(-beta) # boltzmann factor, probability of successfully raising the energy
if move == 1:
r = np.random.random() # random number, r, to determine whether the move to add energy is successful
if r > prob:
e[p] = e[p] + 1 # accept the move if the probability is above e^(-b*delta*E)
else:
if e[p] > 0:
e[p] -= 1 # accept the move if it lowers the particle's energy, but not below zero
avg_e.append(np.mean(e)) # update the average energy array
# plot histogram of energies
plt.subplot(2,2,1)
plt.hist(e,bins=50)
plt.title('histogram of energies')
plt.xlabel('energy')
# plot average energy throughout the simulation
plt.subplot(2,2,2)
plt.plot(range(len(avg_e)), avg_e)
plt.title('average energy over iterations')
plt.xlabel('iterations')
plt.ylabel('average energy')
#plt.figure(figsize=(10,10))
#plt.figure(figsize=(20,15))
plt.tight_layout()
plt.show()