import math as math
import random as random
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
m = 3.2 # mass is in unit of Kb
Kb = 1.0
N = 100
u0 = 150.00
dv = 0.01
def starter(key):
u = [0.0 for k in range(N)]
if key == 0:return u
else:
for k in range(N):
u[k] = u0
return u
def hamiltonian(u):
H = 0.0
for k in range(N):
H = H + 0.5*u[k]**2
return H
def pick_random_particle():
n = random.randint(0, N-1)
return n
def thermalize(u,T,nruns):
irun = 0
h_stor = [0.0 for k in range(nrun)]
while irun < nrun:
L = len(u)
h_old = hamiltonian(u)
n = pick_random_particle()
ov = u[n]
du = random.uniform(-dv,dv)
u[n] = u[n] + du
h_new = hamiltonian(u)
dh = 0.5*m*(u[n]**2-ov**2)
if dh < 0:
h_stor[irun] = h_new
else:
frac = math.exp(-dh/(Kb*T))
b = random.uniform(0.0,1.0)
if b < frac:
h_stor[irun] = h_new
else:
u[n] = u[n] - du
h_stor[irun] = h_old
if u[n] != ov :
irun = irun +1
return h_stor,u
N = 10
u = starter(1)
T = 300.0
nrun = 100
H,u = thermalize(u,T,nrun)
N = 1000
u = starter(1)
T = 300.0
nrun = 10000
H,u = thermalize(u,T,nrun)
X = np.arange(0,len(H),1)
plt.figure(figsize = [15,4])
plt.grid(True)
plt.plot(X,H,"-")
plt.show()
num_bins = 20
plt.figure(figsize = [10,8])
plt.hist(u,num_bins, density= 1.0, facecolor='green', alpha = 0.5)
plt.grid(True)
plt.show()
num_bins = 100
plt.figure(figsize = [10,8])
plt.hist(u,num_bins, density = 1.0, facecolor='green', alpha = 0.5)
plt.grid(True)
plt.show()