Lecture 15: February 19

Kink-jump polymer dynamics

1. Relaxation time from the initial configuration, and

2. The autocorrelation time measured from the time-autocorrelation function.

The following applet, which is generated by the program PolymerDynamics.java, demonstrates the "kink-jump" enrichment algorithm.

Here is the code which generates this applet:

// PolymerDynamics.java

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

public class PolymerDynamics extends Animation {

    int L = 12;                 // size of lattice in each dimension
    int N = 8;                  // number of bonds
    int[] xMonomer, yMonomer;   // positions of monomers

    int polymerNumber;
    double RSquaredSum;
    double DeltaRCMSquaredSum;
    double xCM0, yCM0;          // initial position of center of mass

    double cm (int[] z) {
        double sum = z[0];
        for (int n = 0; n < N; n++)
            sum += z[n + 1];
        return sum /= (N + 1);
    }

    void data () {
        double dx = xMonomer[0] - xMonomer[N];
        double dy = yMonomer[0] - yMonomer[N];
        RSquaredSum += dx * dx + dy * dy;
        dx = cm(xMonomer) - xCM0;
        dy = cm(yMonomer) - yCM0;
        DeltaRCMSquaredSum += dx * dx + dy * dy;
    }

    void initial () {
        if (xMonomer == null || xMonomer.length != N + 1) {
            xMonomer = new int[N + 1];
            yMonomer = new int[N + 1];
        }
        xMonomer[0] = L / 2 - (int) (N / Math.sqrt(2) / 2);
        yMonomer[0] = xMonomer[0];
        for (int i = 0; i < N; i++) {
            xMonomer[i + 1] = xMonomer[i];
            yMonomer[i + 1] = yMonomer[i];
            if (i % 2 == 0)
                ++xMonomer[i + 1];
            else
                ++yMonomer[i + 1];
        }
        polymerNumber = 1;
        xCM0 = cm(xMonomer);
        yCM0 = cm(yMonomer);
        RSquaredSum = DeltaRCMSquaredSum = 0;
        data();
    }

    boolean occupied (int x, int y) {
        for (int i = 0; i <= N; i++)
            if (xMonomer[i] == x && yMonomer[i] == y)
                return true;
        return false;
    }

    void move () {
        // select a monomer at random
        int n = (int) ((N + 1) * Math.random());
        int x = xMonomer[n];
        int y = yMonomer[n];
        if (n != 0 && n != N) {    // not an end site
            if (xMonomer[n - 1] != xMonomer[n])
                x = xMonomer[n - 1];
            if (xMonomer[n + 1] != xMonomer[n])
                x = xMonomer[n + 1];
            if (yMonomer[n - 1] != yMonomer[n])
                y = yMonomer[n - 1];
            if (yMonomer[n + 1] != yMonomer[n])
                y = yMonomer[n + 1];
        } else {                   // an end site
            int i = n == 0 ? 1 : N - 1;
            int plusOrMinus1 = 1;
            if (Math.random() > 0.5)
                plusOrMinus1 = -1;
            if (x == xMonomer[i]) {
                x += plusOrMinus1;
                y = yMonomer[i];
            } else {
                y += plusOrMinus1;
                x = xMonomer[i];
            }
        }
        if (!occupied(x, y)) {
            xMonomer[n] = x;
            yMonomer[n] = y;
        }
    }

    void oneMonteCarloStep () {
        for (int n = 0; n <= N; n++)
            move();
        ++polymerNumber;
        data();
    }

    class Movie extends Plot {
        Movie () {
            setSize(400, 400);
            setBackground("white");
        }

        int pbc (int i) {
            while (i < 0)
                i += L;
            while (i >= L)
                i -= L;
            return i;
        }

        public void paint () {
            clear();
            // the lattice
            setWindow(-0.5, L - 0.5, -0.5, L - 0.5);
            setColor("black");
            for (int i = 0; i < L; i++) {
                for (int j = 0; j < L; j++) {
                    plotPoint(i, j);
                }
            }
            // the polymer
            double r = 0.15;
            int x = pbc(xMonomer[0]);
            int y = pbc(yMonomer[0]);
            for (int i = 0; i < N; i++) {
                int x1 = pbc(xMonomer[i + 1]);
                int y1 = pbc(yMonomer[i + 1]);
                setColor("blue");
                if (Math.abs(x - x1) < 2 && Math.abs(y - y1) < 2)
                    plotLine(x, y, x1, y1);
                setColor("green");
                floodCircle(x - r, x + r, y - r, y + r);
                x = x1;
                y = y1;
            }
            floodCircle(x - r, x + r, y - r, y + r);
        }
    }

    class Output extends Plot {
        Output () {
            setSize(400, 30);
            setBackground("blue");
        }

        public void paint () {
            clear();
            setWindow(0, 100, 0, 3);
            setColor("white");
            plotString("Polymer # " + polymerNumber, 5, 1);
            plotStringCenter(" = " +
                             Easy.format(RSquaredSum / polymerNumber, 4),
                             50, 1);
            plotStringLeft(" = " +
                           Easy.format(DeltaRCMSquaredSum / polymerNumber, 4),
                           95, 1);
        }
    }

    Movie movie;
    Output output;
    Reader NReader, LReader;
        
    public void init () {
        initial();
        add(movie = new Movie());
        add(output = new Output());
        add(NReader = new Reader("N =", N));
        add(LReader = new Reader("L =", L));
        addControlPanel();
    }

    public void step () {
        oneMonteCarloStep();
        movie.repaint();
        output.repaint();
    }

    public void reset () {
        N = NReader.readInt();
        L = LReader.readInt();
        initial();
        movie.repaint();
        output.repaint();
    }
    
    public static void main (String[] args) {
        PolymerDynamics polymerDynamics = new PolymerDynamics();
        polymerDynamics.frame("Polymer Dynamics Simulation", 430, 600);
    }
}