Naive bosons.py

From Werner KRAUTH

This is the python program Naive_bosons.py useful for the homework session of week 13 of my ICFP Lectures on statistical mechanics.

import math, numpy, pylab

E = [0] + 3 * [1] + 6 * [2] + 10 * [3] + 15 * [4]
N0vec = []
Tvec = []
for T in numpy.arange(0.1, 1.5, 0.1):
    Tvec.append(T)
    Z = 0.0
    N0mean = 0.0
    beta = 1.0 / T
    for s1 in range(35):
        for s2 in range(s1, 35):
            for s3 in range(s2, 35):
                for s4 in range(s3, 35):
                    for s5 in range(s4, 35):
                        Etot=E[s1] + E[s2] + E[s3] + E[s4] + E[s5]
                        svec = [s1, s2, s3, s4, s5]
                        N0mean += svec.count(0) * math.exp(-beta * Etot)
                        Z += math.exp( -beta * Etot)
    print T, Z, N0mean / Z / 5.0, 'T, Z, N0'
    N0vec.append(N0mean / Z / 5.0)
pylab.plot(Tvec, N0vec, 'r-')
pylab.title('Naive $N=5$ 5-boson trap model')
pylab.xlabel('$T$')
pylab.xlabel('$N_0 / N $')
pylab.show()