Home
|
Course Outline
|
Lectures
|
Homework
|
Files
The following applet, which is generated by the program Autocorrelation.java, measures correlation times in the energy and magnetization:
Here is the code which generates this applet:
// Autocorrelation.java
import comphys.graphics.*;
import comphys.Easy;
public class Autocorrelation extends Animation {
// compute energy and magnetization autocorrelations
// in the 2-d Ising model using the Metropolis method
int L = 16; // linear size of lattice
int N = L * L; // number of spins
int[][] spin; // spin values
double T = 3; // absolute temperature
double H; // magnetic field
int M; // magnetization
int SS; // sum s_i s_j
double e; // energy / spin
double m; // magnetization / spin
double[][] w; // Boltzmann factors
int nequil = 20; // equilibration steps
int pass; // Monte Carlo steps after nequil
int nsave = 10; // values to save for correlations
double[] esave; // saved E values
double[] msave; // saved M values
double[] Ce; // energy correlation sums
double[] Cm; // magnetization correlation sums
double[] ecum; // to accumulate E and E*E
double[] mcum; // to accumulate M and M*M
void initial () {
N = L * L;
if (spin == null || spin.length < L)
spin = new int[L][L];
// initialize all spins up
for (int x = 0; x < L; x++)
for (int y = 0; y < L; y++)
spin[x][y] = +1;
M = N;
SS = 2 * N;
// Boltzmann factors
if (w == null)
w = new double[17][3];
for (int i = -8; i <= 8; i += 4) {
w[i + 8][0] = Math.exp( - (i + 2 * H) / T);
w[i + 8][2] = Math.exp( - (i - 2 * H) / T);
}
mcs = 0;
// initialize correlation data
pass = 0;
if (esave == null || esave.length < nsave) {
esave = new double[nsave];
msave = new double[nsave];
Ce = new double[nsave];
Cm = new double[nsave];
}
for (int n = 0; n < nsave; n++)
esave[n] = msave[n] = Ce[n] = Cm[n] = 0;
if (ecum == null) {
ecum = new double[2];
mcum = new double[2];
}
for (int i = 0; i < 2; i++)
ecum[i] = mcum[i] = 0;
}
void correl () {
++pass;
// accumulate data for time autocorrelation functions
// index0 = array index of earliest saved time
int index0 = (pass - 1) % nsave;
if (pass > nsave) {
// compute Ce and Cm after nsave values are saved
int index = index0;
for (int t = nsave - 1; t >= 0; t--) {
Ce[t] += e * esave[index];
Cm[t] += m * msave[index];
if (++index == nsave)
index = 0;
}
}
// save latest value in array position of earliest value
esave[index0] = e;
msave[index0] = m;
// accumulate for averages
ecum[0] += e;
mcum[0] += m;
ecum[1] += e * e;
mcum[1] += m * m;
}
int sumOfNeighbors (int x, int y) {
// periodic boundary conditions
int left = 0;
if (x == 0)
left = spin[L - 1][y];
else left = spin[x - 1][y];
int right = 0;
if (x == L - 1)
right = spin[0][y];
else right = spin[x + 1][y];
int down = 0;
if (y == 0)
down = spin[x][L - 1];
else down = spin[x][y - 1];
int up = 0;
if (y == L - 1)
up = spin[x][0];
else up = spin[x][y + 1];
return left + right + up + down;
}
int mcs; // Monte Carlo steps per spin
int accept; // accepted Metropolis steps
void Metropolis () {
// one Monte Carlo step per spin
for (int ispin = 0; ispin < N; ispin++) {
int x = (int) (L * Math.random());
if (x == L)
--x;
int y = (int) (L * Math.random());
if (y == L)
--y;
int dSS = 2 * spin[x][y] * sumOfNeighbors(x, y);
if (Math.random() < w[dSS + 8][spin[x][y] + 1]) {
spin[x][y] = -spin[x][y]; // flip spin
++accept;
M += 2 * spin[x][y];
SS -= dSS;
}
}
e = - SS - H * M;
e /= N;
m = M;
if (m < 0)
m = -m;
m /= N;
++mcs;
if (mcs > nequil)
correl();
}
class Lattice extends Plot {
Lattice () {
setSize(200, 200);
}
public void paint () {
clear();
setWindow(0, L, 0, L);
for (int x = 0; x < L; x++) {
for (int y = 0; y < L; y++) {
if (spin[x][y] == 1)
setColor("red");
else
setColor("green");
boxArea(x, x + 1.1, y, y + 1.1);
}
}
}
}
class Graph extends Plot {
String name;
double[] Acum;
double[] A;
Graph (String name, double[] Acum, double[] A) {
this.name = name;
this.Acum = Acum;
this.A = A;
setSize(200, 200);
}
public void paint () {
clear();
setColor("black");
drawAxes(0, nsave, -0.2, 1);
setColor("blue");
plotStringCenter(name, nsave / 2, 1.02);
if (pass <= nsave)
return;
double Aav = Acum[0] / pass;
double A2av = Acum[1] / pass;
double yPrevious = 1;
double xSum = 0, ySum = 0, xxSum = 0, xySum = 0;
boolean negative = false;
setColor("red");
for (int t = 0; t < nsave; t++) {
double y = A[t] / (pass - nsave) - Aav * Aav;
y /= (A2av - Aav * Aav);
plotLine(t, yPrevious, t + 1, y);
yPrevious = y;
// accumulate sums to compute slope
if (y > 0 && !negative) {
double x = t + 1;
xSum += x;
xxSum += x * x;
double logy = Math.log(y);
ySum += logy;
xySum += x * logy;
} else
negative = true;
}
double tau = ySum * xSum - xySum;
if (tau != 0) {
tau = (xxSum - xSum * xSum) / tau;
setColor("magenta");
plotStringLeft("tau = " + Easy.format(tau, 3), nsave, -0.2);
}
}
}
class Output extends Plot {
Output () {
setSize(600, 30);
}
public void paint () {
clear();
setWindow(0, 100, 0, 3);
setColor("blue");
boxArea(0, 100, 0, 3);
setColor("white");
plotString("MC Step = " + mcs, 5, 1);
if (pass == 0) {
if (mcs > 0)
plotString("Thermalizing ...", 30, 1);
return;
} else if (pass <= nsave) {
plotString("Initializing correlation data ...", 30, 1);
return;
}
double eav = ecum[0] / pass;
double mav = mcum[0] / pass;
plotString(" / N = " + Easy.format(eav, 6), 40, 1);
plotString(" / N = " + Easy.format(mav, 6), 70, 1);
}
}
Lattice lattice;
Graph eGraph, mGraph;
Output output;
Reader LReader, TReader, nsaveReader, nequilReader, skipReader;
int skip;
public void init () {
initial();
add(lattice = new Lattice());
add(eGraph = new Graph("Energy Autocorrelation", ecum, Ce));
add(mGraph = new Graph("Magnetization Autocorrelation", mcum, Cm));
add(output = new Output());
add(LReader = new Reader("L = ", L));
add(TReader = new Reader("T = ", T, 4));
add(nsaveReader = new Reader("Corr steps = ", nsave));
add(nequilReader = new Reader("Therm steps = ", nequil));
add(skipReader = new Reader("Skip steps = ", skip));
addControlPanel();
}
public void step () {
for (int step = 0; step < skip; step++)
Metropolis();
Metropolis();
lattice.repaint();
output.repaint();
if (mcs > nequil) {
eGraph.repaint();
mGraph.repaint();
}
}
public void reset () {
L = LReader.readInt();
T = TReader.readDouble();
nsave = nsaveReader.readInt();
nequil = nequilReader.readInt();
skip = skipReader.readInt();
initial();
eGraph.Acum = ecum;
eGraph.A = Ce;
mGraph.Acum = mcum;
mGraph.A = Cm;
lattice.repaint();
eGraph.repaint();
mGraph.repaint();
output.repaint();
}
public static void main (String[] args) {
Autocorrelation autocorrelation = new Autocorrelation();
autocorrelation.frame("Ising Model Autocorrelations", 650, 420);
}
}
Questions or comments:
phygons@acsu.buffalo.edu