Lecture 18: February 26

The One-Dimensional Ising Model

The following applet is generated by the program Demon.java, which is a translation of PROGRAM Demon on pages 553-555.

Here is the code which generates this applet:

// Demon.java

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

public class Demon extends Animation {

    int N = 200;                // number of spins
    int[] spin;                 // Ising spins in 1-d
    int Etot = -120;            // desired total energy
    int E;                      // energy of spins
    int M;                      // net magnetization
    int Ed;                     // energy of demon

    int mcs;                    // Monte Carlo steps per spin
    double Ecum;                // accumulate energy and other
    double Edcum;               // observables to compute
    double Mcum;                // ensemble averages
    double M2cum;
    int accept;

    double Eave;                // average energy
    double Mave;
    double M2ave;
    double Edave;               // average demon energy
    double T;                   // corresponding absolute temperature
    double acceptProb;
    
    void initial () {
        if (spin == null || spin.length != N)
            spin = new int[N];
        Etot = 4 * (Etot / 4);
        for (int n = 0; n < N; n++)
            spin[n] = 1;
        M = N;
        //compute initial system energy
        E = -N;
        Ed = Etot - E;
        mcs = 0;
        // initialize sums
        Ecum = Edcum = Mcum = M2cum = 0;
        accept = 0;
    }

    void change () {
        for (int n = 0; n < N; n++) {
            // choose a random spin
            int isite = (int) (N * Math.random());
            // neighbors with periodic boundary conditions
            int left = spin[N - 1];
            if (isite > 0)
                left = spin[isite - 1];
            int right = spin[0];
            if (isite < N - 1)
                right = spin[isite + 1];
            // trial energy change
            int de = 2 * spin[isite] * (left + right);
            if (de <= Ed) {
                // spin flip dynamics
                spin[isite] = -spin[isite];
                M += 2 * spin[isite];
                ++accept;
                Ed -= de;
                E += de;
            }
        }
        ++mcs;
    }

    void data () {
        // accumulate data
        Ecum += E;
        Edcum += Ed;
        Mcum += M;
        M2cum += M * M;
    }

    void averages () {
        double norm = 1;
        if (mcs > 0)
            norm /= mcs;
        Edave = Edcum * norm;
        T = 4 / Math.log(1 + 4 / Edave);
        Eave = Ecum * norm;
        Mave = Mcum * norm;
        M2ave = M2cum * norm;
        acceptProb = accept * norm / (double) N * 100.0;
    }

    class Lattice extends Plot {

        public void paint () {
            clear();
            setWindow(-2, 2, -2, 2);
            for (int n = 0; n < N; n++) {
                double theta = n * 2 * Math.PI / (double) N;
                double r = 1.2;
                double x = r * Math.cos(theta);
                double y = r * Math.sin(theta);
                double dr = 1;
                if (spin[n] == 1) {
                    setColor("red");
                    dr = 1.2;
                } else {
                    setColor("cyan");
                    dr = 0.8;
                }
                double x1 = dr * x;
                double y1 = dr * y;
                plotLine(x, y, x1, y1);
            }
            setColor("black");
            double x = 0;
            double y = 0;
            plotStringCenter("MC step = " + mcs, x, y);
            setColor("magenta");
            y = 1.8;
            x = -1.8;
            plotString(" = " + Easy.format(Edave,5), x, y);
            x = 0;
            plotStringCenter(" = " + Easy.format(Eave,6), x, y);
            x = 1.8;
            plotStringLeft(" = " + Easy.format(T, 4), x, y);
            setColor("blue");
            y = -1.8;
            x = -1.8;
            plotString(" = " + Easy.format(Mave, 4), x, y);
            x = 0;
            plotStringCenter(" = " + Easy.format(M2ave, 4), x, y);
            x = 1.8;
            plotStringLeft("Accept = " + Easy.format(acceptProb, 3)
                           + " %", x, y);
        }
    }

    Lattice lattice;
    Reader NReader, EtotReader;

    public void init () {
        initial();
        data();
        averages();
        add(lattice = new Lattice());
        add(NReader = new Reader("N = ", N));
        add(EtotReader = new Reader("E_tot = ", Etot));
        addControlPanel();
    }

    public void step () {
        change();
        data();
        averages();
        lattice.repaint();
    }

    public void reset () {
        N = NReader.readInt();
        Etot = EtotReader.readInt();
        initial();
        lattice.repaint();
    }

    public static void main (String[] args) {
        Demon demon = new Demon();
        demon.frame("Demon algorithm: 1-d Ising Model", 450, 550);
    }
}