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()