Lecture 17: February 23

The Microcanonical Ensemble

The following applet is generated by the program Ideal.java, which is a translation of PROGRAM ideal on pages 547-548.

Here is the code which generates this applet:

// Ideal.java

import comphys.graphics.*;
import comphys.Easy;

public class Ideal extends Animation {

    int N = 40;                 // number of particles
    double[] v;                 // velocities of particles
    double E = 10;              // total energy
    double E0 = E;              // initial total energy
    double Ed;                  // demon energy

    int mcs;                    // Monte Carlo steps per particle
    double dvmax = 2;           // maximum change in velocity
    double Ecum;                // sum to compute average energy
    double Edcum;               // sum to compute average demon energy
    int accept;                 // number of acceptances

    int nHist = 100;            // number of bins in v histogram
    int[] vHist;                // v histogram
    double vHistMax;            // max value of v in histogram

    void initial () {
        if (v == null || v.length != N)
            v = new double[N];
        Ed = 0;
        // divide energy equally among particles
        double vinitial = Math.sqrt(2 * E / N);
        // all particles have same initial velocities
        for (int i = 0; i < N; i++)
            v[i] = vinitial;
        // initialize sums
        Ecum = 0;
        Edcum = 0;
        accept = 0;
        E0 = E;
        mcs = 0;

        if (vHist == null || vHist.length != nHist)
            vHist = new int[nHist];
        for (int n = 0; n < nHist; n++)
            vHist[n] = 0;
        vHistMax = 3 * Math.sqrt(2 * E / N);
    }

    void change () {
        ++mcs;
        for (int i = 0; i < N; i++) {
            // trial change in velocity
            double dv = (2 * Math.random() - 1) * dvmax;
            // select a random particle
            int ipart = (int) (N * Math.random());
            // trial velocity
            double vtrial = v[ipart] + dv;
            // trial energy change
            double de = 0.5 * (vtrial * vtrial - v[ipart] * v[ipart]);
            if (de < Ed) {
                v[ipart] = vtrial;
                ++accept;
                Ed -= de;
                E += de;
            }
        }
        // accumulate data after each Monte Carlo step per particle
        Ecum += E;
        Edcum += Ed;
        // add velocities to histogram
        for (int n = 0; n < N; n++)
            if (Math.abs(v[n]) < vHistMax)
                ++vHist[(int)(nHist * (v[n] + vHistMax) / (2 * vHistMax))];
    }

    double Esave;               // mean energy per system particle
    double Edave;               // mean demon energy
    double acceptProb;          // acceptance probability

    void averages () {
        if (mcs == 0)
            return;
        double norm = 1.0 / mcs;
        Edave = Edcum * norm;
        norm /= N;
        acceptProb = accept * norm * 100;
        Esave = Ecum * norm;
    }

    int stepsToGraph = 100;

    class Graph extends Plot {
        
        int number = 0;
        int current = 0;
        double[] values = new double[stepsToGraph];
        double vMin = 0;
        double vMax = +1;
        String name;
        
        Graph (String name) {
            setSize(450, 150);
            this.name = name;
        }

        void reset () {
            current = number = 0;
            for (int v = 1; v < values.length; v++)
                values[v] = 0;
        }

        void add (double value) {
            if (values.length < stepsToGraph) {
                double[] temp = values;
                values = new double[stepsToGraph];
                for (int v = 0; v < temp.length; v++)
                    values[v] = temp[v];
            }
            if (number < stepsToGraph - 1)
                ++number;
            if (++current == stepsToGraph)
                current = 0;
            values[current] = value;
        }

        public void paint () {
            clear();
            setColor("black");
            drawAxes(0, stepsToGraph, vMin, vMax);
            setColor("red");
            plotStringCenter(name, 0.5 * stepsToGraph, vMax);
            double xOld = 0;
            double yOld = 0;
            setColor("blue");
            for (int n = 0; n < number; n++) {
                double x = n;
                int iy = current - number + n;
                if (iy < 0)
                    iy += stepsToGraph;
                double y = values[iy];
                if (n == 0)
                    plotLine(x, y, x, y);
                else
                    plotLine(x, y, xOld, yOld);
                xOld = x;
                yOld = y;
            }
        }
    }

    class Histogram extends Plot {
        Histogram () {
            setSize(450, 150);
        }

        public void paint () {
            clear();
            double hMax = 0;
            for (int n = 0; n < nHist; n++)
                if (vHist[n] > hMax)
                    hMax = vHist[n];
            if (hMax == 0)
                hMax = 1;
            setColor("black");
            double xMax = nHist / 2.0;
            drawAxes(-xMax, xMax, 0, hMax);
            setColor("red");
            plotStringCenter("Velocity Histogram", 0, hMax);
            plotOverlay();
            setColor("green");
            for (int n = 0; n < nHist; n++) {
                double x = -xMax - 0.5 + n;
                boxArea(x, x + 1, 0, vHist[n]);
            }
        }
    }

    class Output extends Plot {
        Output () {
            setSize(450, 30);
        }

        public void paint () {
            clear();
            setWindow(0, 20, 0, 3);
            setColor("red");
            plotString("=" + Easy.format(Edave, 4), 0, 1);
            plotStringCenter("=" + Easy.format(Esave, 4), 7, 1);
            plotStringCenter("MC Step: " + mcs, 13, 1);
            plotStringLeft("Accept=" + Easy.format(acceptProb, 3)
                           + " %", 20, 1);
        }
    }

    Graph EGraph, EdGraph;
    Histogram histogram;
    Output output;
    Reader NReader, EReader, dvmaxReader;

    void updateGraphData () {
        EGraph.add(E / N);
        EdGraph.add(Ed);
    }

    void repaintGraphs () {
        EGraph.vMax = 1.5 * E0 / N;
        EdGraph.vMax = 4 * E0 / N;
        EGraph.repaint();
        EdGraph.repaint();
        histogram.repaint();
        averages();
        output.repaint();
    }
            
    public void init () {
        initial();
        add(EGraph = new Graph("Energy Per Particle E/N"));
        add(EdGraph = new Graph("Demon Energy E_d"));
        add(histogram = new Histogram());
        add(output = new Output());
        add(NReader = new Reader("N = ", N));
        add(EReader = new Reader("E = ", E, 4));
        add(dvmaxReader = new Reader("dvmax", dvmax, 4));
        addControlPanel();
        repaintGraphs();
    }

    public void step () {
        change();
        updateGraphData();
        repaintGraphs();
    }

    public void reset () {
        N = NReader.readInt();
        E = EReader.readDouble();
        dvmax = dvmaxReader.readDouble();
        Esave = Edave = acceptProb = 0;
        initial();
        EGraph.reset();
        EdGraph.reset();
        repaintGraphs();
    }

    public static void main (String[] args) {
        Ideal ideal = new Ideal();
        ideal.frame("Demon Algorithm for 1-D Ideal Gas", 500, 650);
    }
}