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.");
}
}
}