Home
|
Course Outline
|
Lectures
|
Homework
|
Files
Here is the code which generates this applet:
// Ising.java
import comphys.graphics.*;
import comphys.Easy;
public class Ising extends Animation {
// 2-D Ising model simulation
int L = 16; // linear dimension of lattice
int N = L * L; // number of spins
int[][] spin; // 2-d array of spin variables
boolean antiFerromagnetic; // J = -1 if true otherwise J = 1
double T = 2; // temperature
double H = 0; // magnetic field
double[][] w; // Boltzmann factors
double E; // total energy
int M; // total magnetic moment
int SM; // staggered magnetization
int SS; // sum s_i * s_j
boolean hotStart = true; // random initial configuration
int thermSteps = 200; // number of thermalization steps
boolean thermalizing; // thermalization steps not done
void initial () {
// allocate spin array if necessary
if (spin == null || spin.length < L)
spin = new int[L][L];
N = L * L;
// initialize spins
M = 0;
SM = 0;
for (int x = 0; x < L; x++) {
for (int y = 0; y < L; y++) {
if (antiFerromagnetic && (x + y) % 2 == 1)
spin[x][y] = -1;
else
spin[x][y] = 1;
if (hotStart)
if (Math.random() > 0.5)
spin[x][y] = -1;
M += spin[x][y];
if ( (x + y) % 2 == 0 )
SM += spin[x][y];
else
SM -= spin[x][y];
}
}
// compute total energy
SS = 0;
for (int x = 0; x < L; x++) {
int right = x + 1;
if (right == L)
right = 0;
for (int y = 0; y < L; y++) {
int up = y + 1;
if (up == L)
up = 0;
int sum = spin[x][up] + spin[right][y];
SS += spin[x][y] * sum;
}
}
int J = 1;
if (antiFerromagnetic)
J = -1;
E = - J * SS - H * M;
computeBoltzmannFactors();
initializeSums();
thermalizing = true;
}
void computeBoltzmannFactors () {
// compute Boltzmann probability ratios
if (w == null)
w = new double[17][3];
int J = 1;
if (antiFerromagnetic)
J = -1;
for (int i = -8; i <= 8; i += 4) {
w[i + 8][0] = Math.exp( - (i * J + 2 * H) / T);
w[i + 8][2] = Math.exp( - (i * J - 2 * H) / T);
}
}
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];
if ( (x + y) % 2 == 0 )
SM += 2 * spin[x][y];
else
SM -= 2 * spin[x][y];
SS -= dSS;
}
}
++mcs;
data();
}
double Esum; // accumulator to compute
double EEsum; // accumulator to compute
double Msum; // accumulator to compute
double MMsum; // accumulator to compute
double absMsum; // accumulator to compute <|M|>
void initializeSums () {
Esum = 0;
EEsum = 0;
Msum = 0;
MMsum = 0;
absMsum = 0;
mcs = 0;
accept = 0;
}
double ePerSpin;
double mPerSpin;
void data () {
int J = 1;
if (antiFerromagnetic)
J = -1;
E = - J * SS - H * M;
ePerSpin = E / N;
Esum += E;
EEsum += E * E;
if (antiFerromagnetic) { // anti-ferromagnetic case
mPerSpin = SM / (double) N;
Msum += SM;
MMsum += SM * (double)SM;
absMsum += Math.abs(SM);
} else {
mPerSpin = M / (double) N;
Msum += M;
MMsum += M * (double)M;
absMsum += Math.abs(M);
}
}
double acceptAverage;
double eAverage;
double e2Average;
double mAverage;
double m2Average;
double abs_mAverage;
double cPerSpin;
double chiPerSpin;
void computeAverages () {
double norm = 1.0 / N;
if (mcs > 0)
norm /= mcs;
acceptAverage = accept * norm;
eAverage = Esum * norm;
e2Average = EEsum * norm;
mAverage = Msum * norm;
m2Average = MMsum * norm;
abs_mAverage = absMsum * norm;
cPerSpin = (e2Average - N * eAverage * eAverage) / (T * T);
chiPerSpin = (m2Average - N * mAverage * mAverage) / T;
}
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);
}
}
if (thermalizing) {
setColor("black");
plotStringCenter("T H E R M A L I Z I N G", L / 2, L / 2);
}
}
}
int stepsToGraph = 200;
class ScrollingGraph extends Plot {
int number = 0;
int current = 0;
double[] values = new double[stepsToGraph];
double vMin = -1;
double vMax = 1;
String legend = new String("");
ScrollingGraph () {
setSize(200, 200);
}
void clearValues () {
current = number = 0;
for (int v = 1; v < values.length; v++)
values[v] = 0;
}
void addValue (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;
}
double vAvg, stdDev;
void computeStatistics () {
vAvg = 0;
double v2Avg = 0;
for (int n = 0; n < number; n++) {
vAvg += values[n];
v2Avg += values[n] * values[n];
}
stdDev = 0;
if (number > 0) {
vAvg /= number;
v2Avg /= number;
stdDev = Math.sqrt(v2Avg - vAvg * vAvg);
}
}
public void paint () {
clear();
setColor("black");
drawAxes(0, stepsToGraph, vMin, vMax);
if (thermalizing)
setColor("red");
else
setColor("blue");
plotStringCenter(legend, stepsToGraph / 2,
vMax + 0.02 * (vMax - vMin));
int j = 0;
for (int n = 1; n < number; n++) {
j = current - number + n;
int j0 = j - 1;
if (j0 < 0)
j0 += stepsToGraph;
if (j < 0)
j += stepsToGraph;
plotLine(n - 1, values[j0], n, values[j]);
}
plotOverlay();
computeStatistics();
setColor("lightGray");
boxArea(0, number, vAvg - stdDev, vAvg + stdDev);
setColor("magenta");
plotLine(0, vAvg, number, vAvg);
}
}
class Output extends Plot {
Output () {
setSize(200, 200);
}
public void paint () {
clear();
setWindow(0, 10, 0, 8.5);
double x1 = 1;
double x2 = 6;
double y = 7.5;
plotString("MC steps", x1, y);
plotString("" + mcs, x2, y);
plotString("", x1, y -= 1);
plotString(Easy.format(acceptAverage, 4), x2, y);
plotString("", x1, y -= 1);
plotString(Easy.format(eAverage, 4), x2, y);
plotString("", x1, y -= 1);
plotString(Easy.format(mAverage, 4), x2, y);
plotString("<|M| / spin>", x1, y -= 1);
plotString(Easy.format(abs_mAverage, 4), x2, y);
plotString("C / spin", x1, y -= 1);
plotString(Easy.format(cPerSpin, 4), x2, y);
plotString("chi / spin", x1, y -= 1);
plotString(Easy.format(chiPerSpin, 4), x2, y);
if (antiFerromagnetic)
plotString("J = -1 staggered M and chi", x1, y -= 1);
}
}
Lattice lattice;
ScrollingGraph mGraph, eGraph;
Output output;
Reader LReader, ferroReader, HReader, TReader, thermStepsReader,
eMinReader, eMaxReader, mMinReader, mMaxReader,
startReader, skipReader, graphReader;
int stepsToSkip;
public void init () {
initial();
add(lattice = new Lattice());
add(mGraph = new ScrollingGraph());
mGraph.legend = "magnetization per spin";
add(eGraph = new ScrollingGraph());
eGraph.vMin = -2;
eGraph.vMax = 0;
eGraph.legend = "energy per spin";
add(output = new Output());
add(LReader = new Reader("L = ", L));
add(HReader = new Reader("H = ", H, 3));
add(TReader = new Reader("T = ", T, 4));
add(thermStepsReader = new Reader("Therm Steps = ", thermSteps));
add(startReader = new Reader("Random Initial Config", hotStart));
add(ferroReader = new Reader("Anti-Ferromagnetic", antiFerromagnetic));
add(eMaxReader = new Reader("e_max = ", eGraph.vMax, 4));
add(mMaxReader = new Reader("m_max = ", mGraph.vMax, 4));
add(eMinReader = new Reader("e_min = ", eGraph.vMin, 4));
add(mMinReader = new Reader("m_min = ", mGraph.vMin, 4));
add(skipReader = new Reader("Skip steps = ", stepsToSkip));
add(graphReader = new Reader("Graph steps = ", stepsToGraph));
addControlPanel();
}
void updateGraphData () {
eGraph.addValue(ePerSpin);
mGraph.addValue(mPerSpin);
}
void updatePlots () {
updateGraphData();
computeAverages();
lattice.repaint();
mGraph.repaint();
eGraph.repaint();
output.repaint();
}
void checkThermalizing () {
if (thermalizing && mcs == thermSteps) {
thermalizing = false;
mGraph.clearValues();
eGraph.clearValues();
initializeSums();
}
}
public void step () {
for (int step = 0; step < stepsToSkip; step++) {
checkThermalizing();
Metropolis();
updateGraphData();
}
checkThermalizing();
Metropolis();
updatePlots();
}
public void reset () {
L = LReader.readInt();
antiFerromagnetic = ferroReader.readBoolean();
H = HReader.readDouble();
T = TReader.readDouble();
thermSteps = thermStepsReader.readInt();
hotStart = startReader.readBoolean();
stepsToSkip = skipReader.readInt();
stepsToGraph = graphReader.readInt();
initial();
mGraph.vMin = mMinReader.readDouble();
mGraph.vMax = mMaxReader.readDouble();
eGraph.vMin = eMinReader.readDouble();
eGraph.vMax = eMaxReader.readDouble();
mGraph.clearValues();
eGraph.clearValues();
updatePlots();
}
public static void main (String[] args) {
Ising ising = new Ising();
ising.frame("Metropolis Simulation of 2-d Ising Model", 430, 700);
}
}
Questions or comments:
phygons@acsu.buffalo.edu