## Harmonic Oscillator using HMC¶

### Import the Library¶

In [7]:
'''import section'''
import numpy as np
import matplotlib.pyplot as plt
import math as  math
import random as random
import seaborn as sns
sns.set()


We want to apply HMC to collection of 100 independent Harmonic Oscillator to get equilibrium configuration.

### Hamiltonian¶

Hamiltonian of Harmonic Oscillator in 1D is:

$H = \frac{1}{2} p^{2} + \frac{1}{2}q^{2}$ with $m = 1,k = 1$

This function calculates the total Hamiltonian of the configuration

In [9]:
def hamiltonian(x,p,np):
'''x,p: x and p are list of position and momentum'''
'''np : number of particles in the system '''
H = 0.0
for k in range(np):
H = H + ((x[k]*x[k])/2.0 + (p[k]*p[k])/2.0 )
return H


### Generating Random Momentum¶

In order to generate random momentum we use "random.gauss"

In [10]:
def drawp(np):
'''this function returns a list of random numbers'''
t = [0.0 for k in range(np)]
for k in range(np):
r = random.gauss(0.0,1.0)
t[k] = r
return(t)


One can check whether the generated numbers are normally distributed or not by doing:

In [13]:
N = 100000
p = [0.0 for k in range(N)]
p = drawp(N)
num_bins = 100
plt.figure(figsize = [10,6])
plt.hist(p,num_bins, density= 1.0, facecolor='green', alpha = 0.5)
plt.show()


### Leap Frog¶

We will use leap frog approximation to evolve the system according to time.

In [14]:
def leap_frog(N,dt,ix,ip,np):

''' N : number of steps to evolve
dt:  fraction of time ie T  = dt*N
ix,ip : initial position and momentum
np : number of the particles in the system
'''
''' Returns
x,p : final position and momentum'''

x = ix
p = ip
k = 0
while k < N:
if k == 0:
for i in range(np):
p[i] = p[i] - ((dt/2.0)*x[i])
elif k > 0 :
if k < N - 1:
for i in range(np):
x[i] = x[i] + (dt*p[i])
p[i] =   p[i] - (dt*x[i])
#S1 = hamiltonian(x,p,np)
#print "k =",k,"S1=",S1

elif k == N - 1:
for i in range(np):
p[i] = (p[i] - (dt/2.0)*x[i])

k = k+1
return x,p


### HMC¶

Here we run the HMC - simulation

In [21]:
def HMC(np,N,dt,steps,x0):

''' np : number of particles in the system
N = number of steps in Leap - Frog
dt = fraction of time in Leap - Frog
steps: total steps in HMC '''

xt = [0.0 for k in range(np)]
pt = [0.0 for k in range(np)]

p0 = drawp(np)
H = [0.0 for k in range(steps)]

S0 = hamiltonian(x0,p0,np)
#print ("=======>", 0,"S0=", S0)

chain = 1
total_frac =  0.0
while chain < steps:
s_stor = [0.0]
xt,pt = leap_frog(N,dt,x0,p0,np)
S1 = hamiltonian(xt,pt,np)
frac = math.exp(-(S1-S0))
#print frac
a = min(1,frac)
b = random.uniform(0.0,1.0)

if b < a:
#print("=======>", chain, "S1=",S1,frac,a,b)
H[chain] = S1
x0 = xt
p0 = drawp(np)
S0 = hamiltonian(x0,p0,np)
else:
H[chain] = S0
p0 = drawp(np)

chain = chain+1

return H


### Seting Constants¶

In [22]:
np = 1000
N = 1000
dt = 0.001
steps = 100


Call HMC

In [23]:
x0 = [1.0 for k in range(np)]
x0 = [random.uniform(0.0,1.0) for k in range(np)]
H = HMC(np,N,dt,steps,x0)


### Plot¶

In [20]:
t = [1.0*k for k in range (steps)]
plt.figure(figsize = [15,4])
plt.scatter(t,H)
plt.show()

In [ ]: