QMC_H2

  N = 100

N = 10000

 

Run I/O:

----jGRASP exec: java qmcH2.MonteCarlo 100

Integrating...
r1 in Stabilising: -0.8105120797046247,0.06259845647383633,0.5586589974492433
x1 mit r1.getComponent(0): -0.8105120797046247
pot. Energy: -0.9274765038394915 Var: -0.012578247210878135
Energy (Virial Theorem): -0.46373825191974577

----jGRASP: operation complete.

----jGRASP exec: java qmcH2.MonteCarlo 10000

Integrating...
r1 in Stabilising: 0.8317117155546266,-2.298805229125598,2.0149717993222036
x1 mit r1.getComponent(0): 0.8317117155546266
pot. Energy: -1.1783991953195483 Var: -2.3591177085201525E-4
Energy (Virial Theorem): -0.5891995976597741
 

MonteCarlo.java

package qmcH2;

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 H2-molecule 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
   private DoubleVector r1;  // Vector of electron1
   private DoubleVector r2;
   private DoubleVector nucR1;  // Vector of nucleus1
   private DoubleVector nucR2;
   private double R;

   PlotFrame frame = new PlotFrame("x", "y", "QMC H2-Molecule");

   /**
   * 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;
     R = 1.4;  // distance between the two nuclei in atomic units; experimental value: 1.40
     r1=new Double3Vector(Math.random(),Math.random(),Math.random());
     r2=new Double3Vector(-Math.random(),-Math.random(),-Math.random());
     nucR1=new Double3Vector(R/2,0.0,0.0); // nucleus1 is set to coordinates x=R/2, y=0, z=0
     nucR2=new Double3Vector(-R/2,0.0,0.0)// nucleus2 is set to coordinates x=-R/2, y=0, z=0
     //energy=new double[N];
     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.setXYColumnNames(1,"x2","y2");
     frame.setRowNumberVisible(true);
   }

   /**
   * Compute the local potential energy.
   */

   private void compute()
   {
     DoubleVector tmpr1,tmpr2;
     double prob,tmpprob;
     // Stabilising
     for(int i=0;i<1000;i++)
     {
       tmpr1=r1.mapComponents(this);
       tmpr2=r2.mapComponents(this);
       tmpprob=trialWF(tmpr1, tmpr2);
       tmpprob*=tmpprob;
       prob=trialWF(r1, r2);
       prob*=prob;
       if(tmpprob/prob>Math.random())
       {
         r1=tmpr1;
         r2=tmpr2;
       }
     }

     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);
       tmpprob=trialWF(tmpr1, tmpr2);
       tmpprob*=tmpprob;
       prob=trialWF(r1, r2);
       prob*=prob;
       if(tmpprob/prob>Math.random())
       {
         r1=tmpr1;
         r2=tmpr2;
       }
       //energy[i]=localEnergy(r1,r2);
       potE[i]=localpotE(r1, r2);

       frame.append(0, r1.getComponent(0), r1.getComponent(1));
       frame.append(1, r2.getComponent(0), r2.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, DoubleVector r2)
   {
     double modR1=r1.norm();
     double modR2=r2.norm();
     double modR12=r1.subtract(r2).norm();
     //return Math.exp(-1*modR1)*Math.exp(-1*modR2)*Math.exp(modR12/2); // wavefunction including el-el repulsion
     return Math.exp(-1*modR1)*Math.exp(-1*modR2);
   }

   private double localpotE(DoubleVector r1, DoubleVector r2)
   {
     DoubleVector r12=r2.subtract(r1);  // interelectronic distance
     DoubleVector nucR12=nucR2.subtract(nucR1)// distance between the two nuclei
     DoubleVector r1R1=r1.subtract(nucR1)// distance between el1 and nuc1
     DoubleVector r2R2=r2.subtract(nucR2)// distance between el2 and nuc2
     DoubleVector r1R2=r1.subtract(nucR2)// distance between el1 and nuc2
     DoubleVector r2R1=r2.subtract(nucR1)// distance between el2 and nuc1

     double termr1R1=-1/(r1R1.norm());
     double termr2R2=-1/(r2R2.norm());
     double termr1R2=-1/(r1R2.norm());
     double termr2R1=-1/(r2R1.norm());
     return termr1R1 + termr2R2 + + termr1R2 + termr2R1 + 1/r12.norm() + 1/nucR12.norm();
   }


   /**
   * Update electron co-ordinates.
   */

   public double map(double x)
   {
     stepW = 0.5;    
     return x+(2*Math.random()-1.0)*stepW;
   }

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