The Photon Gas

The Photon Gas#

Note: This exercise was updates compared to previous years.

This exercise gives you an example of a Metropolis Monte Carlo algorithm, where you will apply this type of algorithm to calculate the state occupancy of vibrational states of a \(H_2\) gas. In this exercise, we will approximate the vibrations with the harmanic oscillator. The resulting system is a set of harmonic oscillators in different energy levels.

Ensemble Averages from the Metropolis Monte Carlo Algorithm#

The ensemble average of the state occupancy \(\left<N_i\right>\) of state \(i\) of an ensemble of harmonic oscillators can be calculated analytically. In the canonical ensemble, the ensemble average occupancy of state \(i\) is expressed as:

\[ \begin{align} \langle N_i\rangle = P_i * N = \frac{e^{-E_i\beta}}{Z}N \end{align} \]

where \(Z\) is the partition function of an ensemble of harmonic oscillators.

We can express the energy of a set of \(N\) harmonic oscillators with the same frequency \(\omega\) as:

\[ \begin{align} E_{\{N_i\}} = \sum_i^N N_i \hbar \omega (n_i+\frac12) \end{align} \]

We will use this expression as part of the Metropolis Monte Carlo algorithm to obtain the state occupancy of our ensemble of harmoinic oscillators.

The scheme you will employ is as follows:

Start with an arbitrary set of \(\{N_i\}\), here we will start with all of them in the ground state.
Decide to perform a trial move for one random harminic oscillator to randomly move it to a higher or lower energy level (we limit the move to \(\pm 1\) energy level).
Accept the trial move with probability (i.e. perform a Metropolis Monte Carlo move):

\[ \begin{aligned} P_{acc}(o \rightarrow n)= \min \left(1, e^{-\beta(E(n)-E(o))}\right),\end{aligned} \]

where \(E(n)\) and \(E(o)\) are the energies of the new and old ensemble of harmonic oscillators respectively.
Update averages regardless of acceptance or rejection.
Iterate from 2).

Note that you will need to find the value of \(\hbar \omega\) for a \(H_2\) gas.

Implementing a Metropolis Monte Carlo Algorithm#

import random as r
import math as m
import matplotlib.pyplot as plt
import numpy as np 

#from ipywidgets import interact, interactive, fixed, interact_manual
#import ipywidgets as widgets

import sys
sys.path.append("..")
#import helpers
#helpers.set_style()

# Random Seed
r.seed(42)

Exercise 3

To express the occupation of our ensemble of harmonic ascillators, we will use a list where each item is the number of particle in the corresponding level. For example occupancy[0]=10 means that there are 10 particles in the ground state (level 0) Create a function that calculates the energy of an occupancy list using the hints provided.

def energy(occupancy,hbaromega):
    # Loop over the energy levels
    # For each energy level, calculate the energy using the expression provided in the introduction (you will use hbaromega and n, the number of particles in that energy level)
    # Add it to the total energy
    return total_energy

Exercise 4

Make modifications in the code, right after the commented section MODIFICATION ... END MODIFICATION. Include the entire code within your report and comment upon the part that you wrote.

Exercise 5

How can this scheme retain detailed balance when \(N_i = 0\)? Note that \(N_i\) cannot be negative.

# This function is given to you and we will use to randomly pick a level from a list of all levels available.
def select_level(occupancy):
    n = np.sum(occupancy)
    k = r.randint(1,n)
    cumulative_occupancy = np.cumsum(occupancy)
    level = np.argmax(cumulative_occupancy >= k)
    return level

# Input the correct values here, for a di-hydrogen gas at room temperature.
numberOfIterations = 1000
T = # FILL IN K
omega = # FILL IN cm-1
levels=100
n=40

# Constants
kb=1.4*10**-23
hbar = 1.05*10**-34

beta = 1/(kb*T)
hbaromega = hbar*omega*3*10**10 # We convert cm-1 to SI units by multiplying by the speed of ligth in cm/s

def calculateOccupancy(beta,hbaromega,numberOfIterations,levels,n):

    # We create an occupancy vector to keep track of how many particles are in each level.
    occupancy = []
    for level in range(levels):
        occupancy.append(0)
    # We put all the particles in the ground state
    occupancy[0] = n

    # This occupancy_sum array is there to collect the statistics of the occupation
    occupancy_sum = np.zeros(len(occupancy))

    for i in range(numberOfIterations):
        # Create new_occupancy, which is a copy of the current occupancy, to store the new occupancy after performing a trial move.
        new_occupancy = occupancy.copy()

        # We select from which level to move a particle
        initial_level = select_level(occupancy)

        """  MODIFICATION
    Metropolis algorithm implementation to calculate {N_i}
    Tasks:
    1) Loop from int i = 0 to numberOfiterations (this was done above)
    2) Randomy select from which level (initial_level) to move a particle by using the select_level function provided to you. (this was done above)
    (modifications start here)
    3) Call r.randint(0,1) to perform a trial move to move the particle from a higher or lower level (final_level = initial_level +/- 1). 
        Note: randint(0,1) returns random integers from 0 to 1 (i.e. 0 OR 1)
               but you need to extract -1 OR 1
        Hint:   what happens to the extraction if you multiply r.randint(0,1)*A, where A is an integer?
                what happens to the extraction if you sum r.randint(0,1)+B, where B is an integer (positive or negative)?
                The above hints shall suggest how to obtain the desired numbers from the extraction of 0 or 1
    4) Test if final_level < 0, if it is, force it to be 0
    5) Test if final_level > levels-1, if it is, force it to be levels-1: by adding these lines:
            if final_level > levels-1:
                final_level = levels-1
    5) Update the initial_level and final_level of new_occupancy so that the number of particles in initial_level is decreased by 1, 
        and the number of particles in final_level is increased by 1.
    5) Accept the trial move with probability defined in the Theory section
        Note: Accepting the trial move means updating the occupancy (old_occupancy).
        with the new occupancy (new_occupancy);
    6) sum occupancy to occupancy_sum to kkeep track of the statistics
    7) Average the occupancy_sum by the number of iterations to compute the average occupancy after numberOfiterations iterations

    END MODIFICATION
    """

        # Call r.randint(0,1) to perform a trial move to move the particle from a higher or lower level

        # check that the level to move the particle to is allowed by our problem

        if final_level > levels-1:
            final_level = levels-1

        # update the occupancy the inital and final level

        # check the Metropolis acceptance criterion
        
        # sum occupancy

    # average sum_occupancy by number of iterations

    return avg_occupancy

# perform a single calculation
estimatedOccupancy = calculateOccupancy(beta,hbaromega,numberOfIterations,levels,n)

Exercise 6 - Bonus

Bonus: Why do we add these lines in our code: if final_level > levels-1: final_level = levels-1 How is this realted to the number of levels?

Exercise 7

Using your code, plot the average occupancy. Do it for two different temperatures and comment on your observations. Calculate the anyltical solution with the expression in the theory above and the partition function for a harmonic oscillator.

\[ \begin{align} Z = \frac{e^{\frac12 \beta \hbar \omega}}{e^{\beta \hbar \omega}-1} \end{align} \]

Plot your calculated values versus those from the analytical solution (two curves in the same plot) and include your curve in your report (you will include 2 plts, one for each temperature chosen).

Now, do the same process with 3 or more values for numberOfIterations. What is the influence of the number of MC iterations on the estimated result vs the analytical one? Why? Justify in words and include supporting plots.

Note: this question requires more time, you need to clearly justify your statements.

# define the 2 temperatures
T1 = # FILL IN K
T2 = #FILL IN K
beta_1 = 1/(kb*T1)
beta_2 = 1/(kb*T2)
numerOfIterations= # fill here with different values to study the influence of this variable on the results.

# run the Metropolis algorithm for each T
estimatedOccupancy_1 = calculateOccupancy(beta_1,hbaromega,numberOfIterations,levels,n)
estimatedOccupancy_2 = calculateOccupancy(beta_2,hbaromega,numberOfIterations,levels,n)

# analytical function

def calculate_analytic_occupancy(beta,hbaromega,levels,n):
    # calculate the partition function Z
    Z = (np.exp(0.5*beta*hbaromega))/(np.exp(beta*hbaromega)-1)
    # calculate the occupancy of each level using the expression provided in the introduction and append it to a list
    analytic_occupancy = []
    for i in range(levels):
        analytic_occupancy.append(n*np.exp(-hbaromega*(i + 0.5)*beta)/Z)
    return analytic_occupancy

# Calculate analytical occupancy for each T
analytical_occupancy_1=calculate_analytic_occupancy(beta_1,hbaromega,levels,n)
analytical_occupancy_2=calculate_analytic_occupancy(beta_2,hbaromega,levels,n)

plt.plot(range(levels), analytical_occupancy_1, label='analytical T1')
plt.plot(range(levels), estimatedOccupancy_1, label='estimated T1')
plt.plot(range(levels), analytical_occupancy_2, label='analytical T2')
plt.plot(range(levels), estimatedOccupancy_2, label='estimated T2')
plt.xlabel('Level')
plt.ylabel('Occupancy')
plt.xlim(0,10)
plt.legend()
plt.show()

Exercise 8

Why does ignoring rejected moves lead to erroneous results? Hint: define \(P'(o \rightarrow o)\) (i.e the probability that you stay in the old configuration) and recall that the transition probability \(P'\) is normalised.

Exercise 9 - Bonus

Based on your comprehension of the exercise, what do you think is the influence of the number of partciples?

The Photon Gas

Contents

The Photon Gas#

Ensemble Averages from the Metropolis Monte Carlo Algorithm#

Implementing a Metropolis Monte Carlo Algorithm#