mhjensen
46 lines
1.5 KiB
Plaintext
46 lines
1.5 KiB
Plaintext
# Program to test the Metropolis algorithm with one particle at given temp in
|
|
# one dimension
|
|
#!/usr/bin/env python
|
|
import numpy as np
|
|
import matplotlib.mlab as mlab
|
|
import matplotlib.pyplot as plt
|
|
import random
|
|
from math import sqrt, exp, log
|
|
# initialize the rng with a seed
|
|
random.seed()
|
|
# Hard coding of input parameters
|
|
MCcycles = 100000
|
|
Temperature = 2.0
|
|
beta = 1./Temperature
|
|
InitialVelocity = -2.0
|
|
CurrentVelocity = InitialVelocity
|
|
Energy = 0.5*InitialVelocity*InitialVelocity
|
|
VelocityRange = 10*sqrt(Temperature)
|
|
VelocityStep = 2*VelocityRange/10.
|
|
AverageEnergy = Energy
|
|
AverageEnergy2 = Energy*Energy
|
|
VelocityValues = np.zeros(MCcycles)
|
|
# The Monte Carlo sampling with Metropolis starts here
|
|
for i in range (1, MCcycles, 1):
|
|
TrialVelocity = CurrentVelocity + (2.0*random.random() - 1.0)*VelocityStep
|
|
EnergyChange = 0.5*(TrialVelocity*TrialVelocity -CurrentVelocity*CurrentVelocity);
|
|
if random.random() <= exp(-beta*EnergyChange):
|
|
CurrentVelocity = TrialVelocity
|
|
Energy += EnergyChange
|
|
VelocityValues[i] = CurrentVelocity
|
|
AverageEnergy += Energy
|
|
AverageEnergy2 += Energy*Energy
|
|
#Final averages
|
|
AverageEnergy = AverageEnergy/MCcycles
|
|
AverageEnergy2 = AverageEnergy2/MCcycles
|
|
Variance = AverageEnergy2 - AverageEnergy*AverageEnergy
|
|
print(AverageEnergy, Variance)
|
|
n, bins, patches = plt.hist(VelocityValues, 400, facecolor='green')
|
|
|
|
plt.xlabel('$v$')
|
|
plt.ylabel('Velocity distribution P(v)')
|
|
plt.title(r'Velocity histogram at $k_BT=2$')
|
|
plt.axis([-5, 5, 0, 600])
|
|
plt.grid(True)
|
|
plt.show()
|