QMC H-Atom
Run main
Run I/O
----jGRASP exec: java qmcH.MonteCarlo 10000
Integrating...
r1 in Stabilising:
-0.08848529339583566,0.059584609091649754,-0.032991675658692876
x1 mit r1.getComponent(0): -0.08848529339583566
pot. Energy: -0.9757972282367787 Var: -1.9808595321628946E-4
Energy (Virial Theorem): -0.48789861411838936 (exakte Energie: -0.5 a.u.)
MonteCarlo.java
package qmcH;
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 H-atom ground state energy.
* @author Guenter Haucke
* @version 1.0
*/
public class MonteCarlo implements Mapping
{
private int N;
private double energy[], potE[];
private DoubleVector r1;
PlotFrame frame = new PlotFrame("x", "y", "QMC H-Atom");
/**
* 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());
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.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.setRowNumberVisible(true);
}
/**
* Compute the ground state energy.
*/
private void compute()
{
DoubleVector tmpr1;
double prob,tmpprob;
// Stabilising
for(int i=0;i<1000;i++)
{
tmpr1=r1.mapComponents(this);
tmpprob=trialWF(tmpr1);
tmpprob*=tmpprob;
prob=trialWF(r1);
prob*=prob;
if(tmpprob/prob>Math.random())
{
r1=tmpr1;
}
}
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);
tmpprob=trialWF(tmpr1);
tmpprob*=tmpprob;
prob=trialWF(r1);
prob*=prob;
if(tmpprob/prob>Math.random())
{
r1=tmpr1;
}
potE[i]=localpotE(r1);
frame.append(0, r1.getComponent(0),
r1.getComponent(1));
frame.setVisible(true);
frame.repaint();
}
}
/**
* Trial wavefunction.
* @param r1 position vector of electron 1
* @param r2 position vector of electron 2
*/
private double trialWF(DoubleVector r1)
{
double modR1=r1.norm();
//double modR2=r2.norm();
//double modR12=r1.subtract(r2).norm();
return Math.exp(-1*modR1);
}
/**
* Local energy calculation.
* @param r1 position vector of electron 1
* @param r2 position vector of electron 2
*/
private double localpotE(DoubleVector r1)
{
//DoubleVector r12=r2.subtract(r1);
double termR1=-1/(r1.norm());
//double termR2=-2/(r2.norm());
return termR1;
}
/**
* Update electron co-ordinates.
*/
public double map(double x)
{
return x+(Math.random()-0.5);
}
/**
* Log results to disk.
*/
private void saveResults()
{
File file=new File("results.dat");
char buf[]=null;
NormalDistribution normpot = new
NormalDistribution(potE);
double data[][]=new double[1][3];
double meanpot = normpot.getMean();
double varpot = normpot.getVariance();
System.out.println("pot. Energy: "+ meanpot + " Var: "
+ varpot);
System.out.println("Energy (Virial Theorem): " +
meanpot/2);
data[0][0]=N;
data[0][1]=meanpot;
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.");
}
}
}