Home
|
Course Outline
|
Lectures
|
Homework
|
Files
The potential in this simulation is V(x) = 0.5 x2 + b x3. The red dots along the x axis in both graphs represent the positions of the ensemble of random walkers. The top graph in this applet contains a plot of the potential, and the current value of the energy E, with the scales of the x and V(x) axes fixed. The bottom graph shows the wave function: the scales of the x and y axes are adjusted as the applet runs to show all of the random walkers and the wave function.
Here is the code which generates this applet:
// QMWalk.java
import comphys.graphics.*;
import comphys.Easy;
public class QMWalk extends Animation {
// define potential function
double b = 0; // anharmonic coefficient
double V (double x) {
return 0.5 * x * x + b * x * x * x;
}
double[] x; // walker positions
double[] psi; // wave function
double E; // current estimate for energy
int N0 = 50; // desired number of walkers
double width = 1; // width of initial region
int N; // number of walkers
double ds = 0.1; // step length
double dt; // time step
double Vref; // reference potential
double Vave; // mean potential
double Esum; // energy sum
int mcs; // Monte Carlo steps
int nequil = 100; // thermalization steps
double binsize; // binsize for psi array
void initial () {
if (x == null)
x = new double[2001];
for (int i = 0; i < 2001; i++)
x[i] = 0;
if (psi == null)
psi = new double[1001];
for (int i = 0; i < 1001; i++)
psi[i] = 0;
dt = ds * ds; // time step
mcs = 0;
N = N0; // initial # of walkers = desired #
// choose initial positions of walkers at random
Vref = 0;
for (int i = 1; i <= N; i++) {
x[i] = (2 * Math.random() - 1) * width;
Vref += V(x[i]);
}
Vref /= N;
Vave = Vref;
binsize = 2 * ds; // binsize for psi array
E = Esum = 0;
data();
}
void walk () {
int Nin = N; // # of walkers at beginning of trial
double Vsum = 0;
for (int i = Nin; i > 0; i--) {
if (Math.random() < 0.5)
x[i] += ds;
else
x[i] -= ds;
double potential = V(x[i]); // potential at x
double dV = potential - Vref;
if (dV < 0) { // check to add walker
if (Math.random() < - dV * dt) {
N++;
x[N] = x[i]; // new walker
// factor of 2 since two walkers at x
Vsum += 2 * potential;
} else {
Vsum += potential; // only do old walker
}
} else {
if (Math.random() < dV * dt) {
x[i] = x[N];
--N;
} else {
Vsum += potential;
}
}
}
Vave = Vsum / N; // mean potential
Vref = Vave - (N - N0) / (N0 * dt); // new reference energy
++mcs;
}
void zeroData () {
Esum = 0;
for (int i = 0; i < 1001; i++)
psi[i] = 0;
}
void data () {
// accumulate data
Esum += Vave;
// bin walkers
for (int i = 1; i <= N; i++) {
int bin = (int) Math.round(x[i] / binsize);
bin += psi.length / 2;
psi[bin] += 1;
}
int n = mcs;
if (n > nequil)
n -= nequil;
if (n > 0)
E = Esum / n;
}
double Vmin = 0; // minimum energy for plotting
double Vmax = 2; // maximum energy for plotting
double xmin = -3; // minimum x for plotting
double xmax = 3; // maximum x for plotting
class Graph extends Plot {
Graph () {
setSize(400, 250);
}
public void paint () {
clear();
setColor("black");
drawAxes(xmin, xmax, Vmin, Vmax);
int points = 100;
double delta = (xmax - xmin) / points;
double xOld = xmin;
double yOld = V(xmin);
setColor("green");
for (int i = 0; i < points; i++) {
double x = xOld + delta;
double y = V(x);
plotLine(xOld, yOld, x, y);
xOld = x;
yOld = y;
}
setColor("blue");
plotLine(xmin, E, xmax, E);
setColor("red");
double rx = 0.01 * (xmax - xmin);
double ry = 0.01 * (Vmax - Vmin);
for (int i = 1; i <= N; i++)
floodCircle(outer.x[i] - rx, outer.x[i] + rx, E - ry, E + ry);
setColor("blue");
double y = Vmax + 0.02 *(Vmax - Vmin);
plotString("Energy E", xmin, y);
setColor("green");
plotStringLeft("V(x)", xmax, y);
plotString("V = " + Easy.format(Vmax, 4), 0, y);
y = -0.08 * (Vmax - Vmin);
setColor("black");
plotStringCenter("x", 0, y);
plotStringCenter(Easy.format(xmin, 3), xmin, y);
plotStringCenter(Easy.format(xmax, 3), xmax, y);
}
}
class Output extends Plot {
Output () {
setSize(400, 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);
plotString("E = " + Easy.format(E, 4), 30, 1);
plotString("V_ref = " + Easy.format(Vref, 4), 55, 1);
plotString("N = " + N, 80, 1);
}
}
class PsiGraph extends Plot {
PsiGraph () {
setSize(400, 250);
}
public void paint () {
clear();
int i0 = psi.length / 2;
int i = i0;
double pmax = psi[i] * psi[i];
double sum = pmax;
double P = pmax;
double Vmin = 0;
double Vmax = 0;
do {
++i;
P = psi[i] * psi[i] + psi[2 * i0 - i] * psi[2 * i0 - i];
if (P > pmax)
pmax = P;
double x = binsize * (i - i0);
double v = V(x);
if (v < Vmin)
Vmin = v;
if (v > Vmax)
Vmax = v;
v = V(-x);
if (v < Vmin)
Vmin = v;
if (v > Vmax)
Vmax = v;
} while (P > 0);
int imax = i - i0;
double norm = Math.sqrt(sum * binsize);
double ymax = Math.sqrt(pmax) / norm;
setColor("black");
drawAxes(-imax, imax, 0, ymax);
setColor("magenta");
for (i = -imax; i <= imax; i++) {
plotLine(i - 1, psi[i0 + i - 1] / norm, i, psi[i0 + i] / norm);
}
setColor("green");
double xOld = binsize * (- imax - 1);
double yOld = (V(xOld) - Vmin) / (Vmax - Vmin) * ymax;
for (i = -imax; i <= imax; i++) {
xOld += binsize;
double y = (V(xOld) - Vmin) / (Vmax - Vmin) * ymax;
plotLine(i - 1, yOld, i, y);
yOld = y;
}
if (Vmin < 0) {
yOld = - Vmin / (Vmax - Vmin) * ymax;
plotLine(-imax, yOld, imax, yOld);
}
setColor("red");
double rx = 0.01 * 2 * imax;
double ry = 0.01 * ymax;
for (i = 1; i <= N; i++) {
xOld = outer.x[i] / binsize;
floodCircle(xOld - rx, xOld + rx, -ry, +ry);
}
double y = (E - Vmin) / (Vmax - Vmin) * ymax;
setColor("blue");
plotLine(-imax, y, imax, y);
y = 1.02 * ymax;
setColor("magenta");
plotString("psi(x)", -imax, y);
setColor("green");
plotStringLeft("V(x)", imax, y);
plotString("V = " + Easy.format(Vmax, 4), 0, y);
y = -0.08 * ymax;
setColor("black");
plotStringCenter("x", 0, y);
double x = imax * binsize;
plotStringCenter(Easy.format(-x, 3), -imax, y);
plotStringCenter(Easy.format(x, 3), imax, y);
}
}
Graph graph;
Output output;
PsiGraph psiGraph;
Reader bReader, N0Reader, widthReader, dsReader, nequilReader, skipReader;
int skip;
public void init () {
initial();
add(graph = new Graph());
add(output = new Output());
add(psiGraph = new PsiGraph());
add(bReader = new Reader("b = ", b, 4));
add(N0Reader = new Reader("N_0 = ", N0));
add(widthReader = new Reader("width = ", width, 4));
add(dsReader = new Reader("ds = ", ds, 4));
add(nequilReader = new Reader("Therm steps", nequil));
add(skipReader = new Reader("Skip steps", skip));
addControlPanel();
}
public void step () {
for (int s = 0; s <= skip; s++) {
walk();
if (mcs == nequil)
zeroData();
else
data();
}
graph.repaint();
output.repaint();
psiGraph.repaint();
}
public void reset () {
b = bReader.readDouble();
N0 = N0Reader.readInt();
width = widthReader.readDouble();
ds = dsReader.readDouble();
nequil = nequilReader.readInt();
skip = skipReader.readInt();
if (skip < 0)
skip = 0;
initial();
graph.repaint();
output.repaint();
psiGraph.repaint();
}
QMWalk outer; // hack to get around problem on unix
public QMWalk () {
super(); // create comphys.graphics.Animation
outer = this; // allows access to x[i] in inner classes
}
public static void main (String[] args) {
QMWalk qmWalk = new QMWalk();
qmWalk.frame("Random Walk Quantum Monte Carlo", 430, 730);
}
}
Questions or comments:
phygons@acsu.buffalo.edu