Canonic bosons.py

From Werner KRAUTH

This is the python2 program Canonic_bosons.py useful for the homework session of week 13 of my ICFP Lectures on statistical mechanics. Here, the density of state is 1,3,6,10,15, which corresponds to the three-dimensional isotropic harmonic trap. In the below program, we integrate in the complex plane from (-pi, lambda) to (pi, lambda) (neglecting the integral from (-pi, 0) to (-pi, lambda)), and always get the same result for the partition function. At T=1, for example, we get Z=17.373..., a result we can also obtain from the naive sum over the states, as implemented in Naive_bosons.py.

import math, cmath, numpy, pylab
 
dos = [1, 3, 6, 10, 15]
T = 1.0
beta = 1.0 / T
Zint = complex(0.0, 0.0)
dellambda = 0.01
oldl = complex(-math.pi, dellambda)
complexi = complex(0.0, 1.0)
for RL in numpy.arange(-math.pi, math.pi, 0.00001):
    newl = RL + dellambda * complexi
    integrand = cmath.exp(-complexi * 5 * newl)
    for E in range(5):
        integrand /= (1.0 - cmath.exp( - beta * E + complexi * newl)) ** dos[E]
    Zint += integrand * (newl - oldl) / (2.0 * math.pi)
    oldl = newl
print Zint