Home
|
Course Outline
|
Lectures
|
Homework
|
Files
The potential in this simulation is the simple harmonic oscillator V(x) = 0.5 x2.
Here is the code which generates this applet:
// PathIntegralMC.java
import comphys.graphics.*;
import comphys.Easy;
public class PathIntegralMC extends Animation {
// Use path integral Monte Carlo to find ground state in a potential
double V (double x) { // potential energy function
return 0.5 * x * x; // harmonic oscillator
}
double dVdx (double x) { // derivative of V(x) used in virial theorem
return x; // harmonic oscillator
}
int N = 10; // number of "atoms"
double tau = 15; // imaginary time = inverse temperature
double deltaTau = tau / N; // imaginary time step
double T = 1 / tau; // temperature
double[] x; // displacements from equilibrium of "atoms"
double E0; // ground state energy
double xmin = -4; // minimum value of x to bin
double xmax = 4; // maximum value of x to bin
double dx = 0.1; // binsize for probability array
int nBins; // number of bins
double[] P; // probability array
double[] phi0; // ground state wave function
int nequil = 2000; // number of thermalization steps
int mcs; // Monte Carlo step number
double delta = 1; // Metropolis step size
void initial () {
deltaTau = tau / N;
T = 1 / tau;
if (x == null || x.length < N)
x = new double[N];
nBins = (int) Math.round((xmax - xmin) / dx);
if (P == null || P.length < nBins) {
P = new double[nBins];
phi0 = new double[nBins];
}
zeroProbability();
// set random displacements
for (int i = 0; i < N; i++) {
x[i] = (2 * Math.random() - 1) * delta;
updateProbability(x[i]);
}
E0 = energySum() / N;
computeVArray(); // for plotting
computePhi0();
mcs = 0;
}
double energySum () {
double sum = 0;
for (int bin = 0; bin < nBins; bin++) {
double x = xmin + (bin + 0.5) * dx;
// use virial theorem
double kineticEnergy = 0.5 * x * dVdx(x);
sum += P[bin] * (kineticEnergy + V(x));
}
return sum;
}
void updateProbability (double x) {
int bin = (int) Math.round((x - xmin) / (xmax - xmin) * nBins);
if (bin < 0)
bin = 0;
if (bin >= nBins)
bin = nBins - 1;
P[bin] += 1;
}
void zeroProbability () {
for (int n = 0; n < nBins; n++)
P[n] = 0;
}
double phi0Max; // used for plotting phi0
void computePhi0 () {
phi0Max = 0;
for (int i = 0; i < nBins; i++) {
phi0[i] = Math.sqrt(P[i] / dx);
if (phi0[i] > phi0Max)
phi0Max = phi0[i];
}
}
void MonteCarloStep () { // one Monte Carlo step per atom
for (int i = 0; i < N; i++) {
// choose an atom j at random
int j = (int) (N * Math.random());
// trial displacement
double xTrial = x[j] + (2 * Math.random() - 1) * delta;
// periodic boundary conditions
int jMinus1 = j - 1;
if (jMinus1 == -1)
jMinus1 = N - 1;
int jPlus1 = j + 1;
if (jPlus1 == N)
jPlus1 = 0;
// change in energy
double deltaE =
(x[jPlus1] - xTrial) * (x[jPlus1] - xTrial) +
(xTrial - x[jMinus1]) * (xTrial - x[jMinus1]) -
(x[jPlus1] - x[j]) * (x[jPlus1] - x[j]) -
(x[jMinus1] - x[j]) * (x[jMinus1] - x[j]);
deltaE /= 2 * deltaTau * deltaTau;
deltaE += V(xTrial) - V(x[j]);
if (deltaE < 0 || Math.random() <= Math.exp(-deltaTau * deltaE))
x[j] = xTrial;
updateProbability(x[j]);
}
++mcs;
// compute energy
double norm = mcs; // number of MC steps * number of atoms
if (mcs > nequil)
norm -= nequil;
norm *= N;
if (norm != 0)
E0 = energySum() / norm;
computePhi0();
}
int nPlot = 100; // number of points for plotting V(x)
double Vmin, Vmax; // used for plotting
double[] VArray; // store values for plotting
void computeVArray () {
if (VArray == null || VArray.length < nPlot)
VArray = new double[nPlot];
Vmin = 1e10;
Vmax = -1e10;
double dx = (xmax - xmin) / (nPlot - 1);
for (int i = 0; i < nPlot; i++) {
double x = xmin + i * dx;
VArray[i] = V(x);
if (VArray[i] < Vmin)
Vmin = VArray[i];
if (VArray[i] > Vmax)
Vmax = VArray[i];
}
}
class Graph extends Plot {
Graph () {
setSize(500, 400);
}
public void paint () {
clear();
setColor("black");
drawAxes(xmin, xmax, 0, N);
plotStringCenter(Easy.format(xmin, 3), xmin, -0.07 * N);
plotStringCenter("x", 0, -0.07 * N);
plotStringCenter(Easy.format(xmax, 3), xmax, -0.07 * N);
// draw potential
double xOld = xmin;
double yScale = N / (Vmax - Vmin);
double yOld = (VArray[0] - Vmin) * yScale;
double dx = (xmax - xmin) / (nPlot - 1);
setColor("green");
for (int i = 1; i < nPlot; i++) {
double x = xOld + dx;
double y = (VArray[i] - Vmin) * yScale;
plotLine(xOld, yOld, x, y);
xOld = x;
yOld = y;
}
plotString("V(x) = 0.5 x^2", xmin, 1.04 * N);
// draw path
xOld = outer.x[0];
yOld = 0;
setColor("cyan");
for (int i = 1; i < N; i++) {
double x = outer.x[i];
plotLine(xOld, yOld, x, yOld += 1);
xOld = x;
}
plotLine(xOld, yOld, outer.x[0], N);
// draw atoms
double rx = 0.01 * (xmax - xmin);
double ry = 0.01 * N * 5.0 / 4.0;
setColor("red");
for (int i = 0; i < N; i++) {
floodCircle(outer.x[i] - rx, outer.x[i] + rx, i - ry, i + ry);
}
floodCircle(outer.x[0] - rx, outer.x[0] + rx, N - ry, N + ry);
if (mcs == 0)
return;
// draw energy
setColor("orange");
double y0 = (E0 - Vmin) * yScale;
plotLine(xmin, y0, xmax, y0);
// draw probability
double Pmax = 0;
for (int i = 0; i < nBins; i++)
if (Pmax < P[i])
Pmax = P[i];
yScale = 0.8 / Pmax * N;
xOld = xmin;
dx = (xmax - xmin) / nBins;
yOld = P[0] * yScale;
setColor("magenta");
for (int i = 1; i < nBins; i++) {
double x = xOld + dx;
double y = P[i] * yScale;
plotLine(xOld, yOld, x, yOld);
plotLine(x, yOld, x, y);
xOld = x;
yOld = y;
}
plotStringCenter("Probability", 0, 1.04 * N);
// draw wave function
xOld = xmin;
yScale = N / (3 * phi0Max);
yOld = y0 + phi0[0] * yScale;
setColor("blue");
for (int i = 1; i < nBins; i++) {
double x = xOld + dx;
double y = y0 + phi0[i] * yScale;
plotLine(xOld, yOld, x, y);
xOld = x;
yOld = y;
}
plotStringLeft("Wave function", xmax, 1.04 * N);
}
}
class Output extends Plot {
Output () {
setSize(500, 30);
}
public void paint () {
clear();
setWindow(0, 100, 0, 3);
setColor("blue");
boxArea(0, 100, 0, 3);
setColor("white");
if (mcs > nequil)
plotString("MC step: " + (mcs - nequil), 5, 1);
else
plotString("Thermalizing " + mcs, 5, 1);
plotStringCenter("E_0 = " + Easy.format(E0, 5), 50, 1);
plotStringLeft("T = " + Easy.format(1 / tau, 3), 95, 1);
}
}
Graph graph;
Output output;
Reader NReader, tauReader, xminReader, xmaxReader, dxReader;
Reader deltaReader, nequilReader, skipReader;
int skip;
public void init () {
initial();
add(graph = new Graph());
add(output = new Output());
add(NReader = new Reader("N = ", N));
add(tauReader = new Reader("tau = ", tau, 4));
add(deltaReader = new Reader("delta = ", delta, 4));
add(xminReader = new Reader("xmin = ", xmin, 4));
add(xmaxReader = new Reader("xmax = ", xmax, 4));
add(dxReader = new Reader("dx = ", dx, 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++) {
MonteCarloStep();
if (mcs == nequil)
zeroProbability();
}
graph.repaint();
output.repaint();
}
public void reset () {
N = NReader.readInt();
tau = tauReader.readDouble();
delta = deltaReader.readDouble();
xmin = xminReader.readDouble();
xmax = xmaxReader.readDouble();
dx = dxReader.readDouble();
nequil = nequilReader.readInt();
skip = skipReader.readInt();
initial();
graph.repaint();
output.repaint();
}
PathIntegralMC outer; // provide access to inner classes
public PathIntegralMC () {
super();
outer = this;
}
public static void main (String[] args) {
PathIntegralMC pimc = new PathIntegralMC();
pimc.frame("Path Integral Monte Carlo", 530, 650);
}
}
Questions or comments:
phygons@acsu.buffalo.edu