QMC_Li
100 Steps
|
10000 steps
|
Run I/O
----jGRASP exec: java QMC_Li.MonteCarlo 100
Integrating...
r1 in Stabilising: -0.38014904947160455,-0.9210851025964849,0.2729655174999872
x1 mit r1.getComponent(0): -0.38014904947160455
pot. Energy: -16.541845941968024 Var: -2.9650144775619016
Energy: -8.270922970984012
----jGRASP: operation complete.
----jGRASP exec: java QMC_Li.MonteCarlo 10000
Integrating...
r1 in Stabilising: -0.12759974449962552,0.26781908269656074,-0.8764658835492126
x1 mit r1.getComponent(0): -0.12759974449962552
pot. Energy: -17.141635374902123 Var: -0.041293387061769694
Energy: -8.570817687451061
MonteCarlo.java
package QMC_Li;
import java.io.*;
import JSci.maths.*;
import JSci.maths.statistics.*;
import JSci.io.*;
import java.awt.Color;
import org.opensourcephysics.frames.PlotFrame;
import org.opensourcephysics.display.*;
/**
* Monte Carlo calculation of Helium ground state energy.
* @author Mark Hale
* @version 1.0
*/
/**
* Monte Carlo calculation of Lithium ground state energy.
* @author Guenter Haucke
* @version 1.0
*/
public class MonteCarlo implements Mapping
{
private int N;
private double potE[];
private double stepW;
// walker step width in atomic units
(bohr)
private DoubleVector r1;
private DoubleVector r2;
private DoubleVector r3;
PlotFrame frame = new PlotFrame("x", "y", "QMC
Lithium");
/**
* Instantiate class.
*/
public static void main(String arg[])
{
if(arg.length==1)
{
int n=Integer.valueOf(arg[0]).intValue();
new MonteCarlo(n);
}
else
{
System.out.println("Need to specify command line arguments:");
System.out.println("<number of iterations>");
}
}
/**
* Constructor.
* @param n number of iterations
*/
public MonteCarlo(int n)
{
N=n;
r1=new Double3Vector(Math.random(),Math.random(),Math.random());
r2=new Double3Vector(-Math.random(),-Math.random(),-Math.random());
r3=new
Double3Vector(-Math.random(),-Math.random(),-Math.random());
potE = new double[N];
frame.setConnected(0, true);
frame.setMarkerShape(0, Dataset.NO_MARKER);
frame.setConnected(1, true);
frame.setMarkerShape(1, Dataset.NO_MARKER);
frame.setConnected(2, true);
frame.setMarkerShape(2, Dataset.NO_MARKER);
frame.setVisible(true);
// frame.repaint();
compute();
saveResults();
frame.setDefaultCloseOperation(javax.swing.JFrame.EXIT_ON_CLOSE);
frame.setXPointsLinked(true); // default is true
frame.setXYColumnNames(0,"x1","y1");
frame.setXYColumnNames(1,"x2","y2");
frame.setXYColumnNames(2,"x3","y3");
frame.setRowNumberVisible(true);
}
/**
* Compute the ground state energy.
*/
private void compute()
{
DoubleVector tmpr1,tmpr2,
tmpr3;
double prob,tmpprob;
// Stabilising
for(int i=0;i<1000;i++)
{
tmpr1=r1.mapComponents(this);
tmpr2=r2.mapComponents(this);
tmpr3=r3.mapComponents(this);
tmpprob=trialWF(tmpr1,tmpr2,tmpr3);
tmpprob*=tmpprob;
prob=trialWF(r1,r2,r3);
prob*=prob;
if(tmpprob/prob>Math.random())
{
r1=tmpr1;
r2=tmpr2;
r3=tmpr3;
}
}
System.out.println("Integrating...");
System.out.println("r1 in Stabilising: " + r1);
System.out.println("x1 mit r1.getComponent(0): " +
r1.getComponent(0));
for(int i=0;i<N;i++)
{
tmpr1=r1.mapComponents(this);
tmpr2=r2.mapComponents(this);
tmpr3=r3.mapComponents(this);
tmpprob=trialWF(tmpr1,tmpr2,tmpr3);
tmpprob*=tmpprob;
prob=trialWF(r1,r2,r3);
prob*=prob;
if(tmpprob/prob>Math.random())
{
r1=tmpr1;
r2=tmpr2;
r3=tmpr3;
}
//energy[i]=localEnergy(r1,r2);
potE[i]=localpotE(r1,r2,r3);
frame.append(0, r1.getComponent(0),
r1.getComponent(1));
frame.append(1, r2.getComponent(0),
r2.getComponent(1));
frame.append(2, r3.getComponent(0), r3.getComponent(1));
frame.setVisible(true);
frame.repaint();
}
}
/**
* Trial wavefunction.
* @param r1 position vector of electron 1
* @param r2 position vector of electron 2
* @param r2 position vector of electron 3
*/
private double trialWF(DoubleVector
r1,DoubleVector r2,DoubleVector r3)
{
double modR1=r1.norm();
double modR2=r2.norm();
double
modR3=r3.norm();
double modR12=r1.subtract(r2).norm();
double
modR13=r1.subtract(r3).norm();
double
modR23=r2.subtract(r3).norm();
//return Math.exp(-3*modR1)*Math.exp(-3*modR2)*(1-modR3/2)*
Math.exp(-3*modR3/2)*Math.exp(modR12/3)*Math.exp(modR13/3)*Math.exp(modR23/3);
return Math.exp(-3*modR1)*Math.exp(-3*modR2)*(1-modR3/2)*Math.exp(-3*modR3/2);
// without e-e interaction terms
}
/**
* potential energy calculation.
* @param r1 position vector of electron 1
* @param r2 position vector of electron 2
* @param r3 position vector
of electron 3
*/
private double localpotE(DoubleVector r1,DoubleVector r2,DoubleVector
r3)
{
DoubleVector r12=r2.subtract(r1);
DoubleVector r13=r3.subtract(r1);
DoubleVector r23=r3.subtract(r2);
double termR1=-3/(r1.norm());
double termR2=-3/(r2.norm());
double termR3=-3/(r3.norm());
return 1/(r12.norm())+ 1/(r13.norm())+ 1/(r23.norm()) +
termR1 + termR2 + termR3;
}
/**
* Update electron co-ordinates.
*/
public double map(double x)
{
stepW = 0.5;
return x+(2*Math.random()-1)*stepW;
}
/**
* Log results to disk.
*/
private void saveResults()
{
File file=new File("results.dat");
char buf[]=null;
//NormalDistribution norm=new NormalDistribution(energy);
NormalDistribution normpot = new NormalDistribution(potE);
double data[][]=new double[1][3];
//double mean=norm.getMean();
//double var=norm.getVariance();
double meanpot = normpot.getMean();
double varpot = normpot.getVariance();
System.out.println("pot. Energy: " + meanpot + " Var: " + varpot);
System.out.println("Energy: " + meanpot/2);
data[0][0]=N;
data[0][1]=meanpot/2;
data[0][2]=varpot;
// Read in existing data
try {
FileReader in=new FileReader(file);
buf=new char[(int)file.length()];
in.read(buf);
in.close();
} catch(Exception e) {
System.out.println("No previous data - new file.");
}
// Save all to file
try {
TextWriter out=new TextWriter(file,',');
if(buf!=null)
out.write(buf);
out.write(data);
out.close();
} catch(Exception e) {
System.out.println("Failed to save.");
}
}
}