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