Quantum Monte Carlo (QMC)

Beim "klassischen" Diffusions-Quanten-Monte-Carlo-Verfahren wird die Potenzialhyperfläche mit einer Schar von "Random-Walkers" abgetastet. Die Wellenfunktion ergibt sich dann aus der Häufigkeit der Walker-Positionen.Ein Beispiel (H2-Molekül) hierfür findet man hier.
Im folgenden Programm wird eine Trial-Wellenfunktion mit den Elektronen als Walker abgetastet. Die potenzielle Energie wird als Mittelwert über die Walker-Positionen bestimmt. Daraus kann mit Hilfe des Virialtheorems die Gesamtenergie bestimmt werden. Das Verfahren hat den Vorteil, dass man die Bewegung der Elektronen anschaulich darstellen und dass man mit einem analytischen Ausdruck für die Wellenfunktion experimentieren kann.

Quantum-Monte-Carlo für das Heliumatom unter Verwendung der Bibliotheken JSci (A Science API for Java, http://sourceforge.net/projects/jsci/) und OSP (OpenSourcePhysics, www.opensourcephysics.org/) als Ausgangsbasis. JSci liefert 3DVektoren, OSP eine komfortable grafische Darstellung. Für jede Teilchen-Koordinate wird ein Walker verwendet (kein Entstehen und Vergehen von Walkern). Das Programm kann leicht auf andere Systeme erweitert werden und man kann damit experimentieren (Variation der Anzahl der Schritte, der Schrittweite, der Wellenfunktion). Beispiele: H, Li, H2.

Diese Seite kann auch als pdf-Datei heruntergeladen werden, Source-Code als jGRASP-Projekt.

Mit jGRASP:

mit n = 40 in Befehlszeile (Run Arguments)

   

Run main:

  Mit Views/Scale

 

Über Tools/DatasetTool kann man die Darstellung weiter variieren

   

Die Daten kann man sich auch als Tabelle ausgeben lassen:

   

Mit einer größeren Anzahl von Walker-Schritten wird die Elektronen-Verteilung deutlicher:

  N = 2000

N = 20000

 

Kommentierter Source-Code

Die blau-grau hervorgehobenen Passagen sind aus OSP eingefügt. Eigene Kommentare sind hell-grün hervorgehoben.

MonteCarlo.java

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
*/

public final class MonteCarlo implements Mapping
{
   private int N;  // walker steps
   private double energy[];   // 1D array of local energy
   private DoubleVector r1;   // radius vector of electron1
   private DoubleVector r2;
  
   PlotFrame frame = new PlotFrame("x", "y", "QMC Helium"); // (String xlabel, String ylabel, String frameTitle)

   /**
   * Instantiate class.
   */

   public static void main(String arg[])
   {
     if(arg.length==1)
     {
       int n=Integer.valueOf(arg[0]).intValue();  // for command line input
       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;  // walker steps
     r1=new Double3Vector(Math.random(),Math.random(),Math.random());  // components x1, y1, z1 with random numbers
                                                             // between 0 and 1 as initial walker position
     r2=new Double3Vector(-Math.random(),-Math.random(),-Math.random());  // same for electron2 but with opposite
                                                                          // sign
     energy=new double[N];

     frame.setConnected(0, true)// Points are connected by straight lines
     frame.setMarkerShape(0, Dataset.NO_MARKER); // Shapes are: NO_MARKER, CIRCLE, SQUARE, AREA, PIXEL, BAR, POST
                                             // weitere Methoden: setMarkerSize(int datasetIndex, int markerSize)
                                            
// setMarkerColor(int datasetIndex, Color color)
     frame.setConnected(1, true);
     frame.setMarkerShape(1, Dataset.NO_MARKER);
     frame.setVisible(true);

     compute();

     saveResults();

     frame.setDefaultCloseOperation(javax.swing.JFrame.EXIT_ON_CLOSE);
     frame.setXPointsLinked(true); // default is true  // X data for datasets > 0 will not be shown in a table view
     frame.setXYColumnNames(0,"x1","y1"); // Sets the column names when rendering this dataset in a JTable
     frame.setXYColumnNames(1,"x2","y2");
     frame.setRowNumberVisible(true);      // Table displays row number
   }


   /**
   * Compute the ground state energy.
   */

   private void compute()
   {
     DoubleVector tmpr1,tmpr2;
     double prob,tmpprob;
     // Stabilising
     for(int i=0;i<1000;i++)
     {
       tmpr1=r1.mapComponents(this); // mapComponents() applies a function on all vector components, here a new walker step
                                     // (-> map())
       tmpr2=r2.mapComponents(this);
       tmpprob=trialWF(tmpr1,tmpr2); // applies the trial wave function on the new coordinates of the two electrons
       tmpprob*=tmpprob;
       prob=trialWF(r1,r2); // applies the trial wave function on the old coordinates of the two electrons
       prob*=prob;
       if(tmpprob/prob>Math.random()) // now it depends on this condition whether the new coordinates will be accepted
       {
         r1=tmpr1;
         r2=tmpr2;
       }
     }

     System.out.println("Integrating...");

     System.out.println("r1 in Stabilising: " + r1); // for monitoring only
     System.out.println("x1 mit r1.getComponent(0): " + r1.getComponent(0));

     for(int i=0;i<N;i++)  // here the actual walk begins
     {
       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);

       //System.out.println("r1 = " + r1);   // for monitoring purposes
       //System.out.println("r2 = " + r2);
       //System.out.println("energy[ " + i + " ] = " + energy[i]);

      
       frame.append(0, r1.getComponent(0), r1.getComponent(1)); // append(int datasetIndex, double x, double y)
       frame.append(1, r2.getComponent(0), r2.getComponent(1));
       frame.setVisible(true)
       frame.repaint(); // Damit wird Grafik laufend dargestellt, was sehr rechenintensiv ist!     
     }
   }

   /**
   * 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();  // square root of (x^2 + y^2 +z^2)
     double modR2=r2.norm();
     double modR12=r1.subtract(r2).norm();  // sqrt of ( (x1^2-x2^2) + ...)
     return Math.exp(-2*modR1)*Math.exp(-2*modR2)*Math.exp(modR12/2);  // general formula of a two-electron system
   }

   /**
   * Local energy calculation.
   * @param r1 position vector of electron 1
   * @param r2 position vector of electron 2
   */

   private double localEnergy(DoubleVector r1,DoubleVector r2)
   {
     DoubleVector r12=r2.subtract(r1);
     double termR112=r1.scalarProduct(r12)/(r1.norm()*r12.norm());
     double termR212=r2.scalarProduct(r12)/(r2.norm()*r12.norm());
     return -17/4-termR112+termR212;  // expression for local energy
   }

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

   public double map(double x)
   {
     return x+(Math.random()-0.5);  // random walk steps between -0.5 and +0.5
   }

   /**
   * Not used, dummy implementation for Mapping interface.
   */

   public Complex map(Complex z)
   {
     return null;
   }

   /**
   * Log results to disk.
   */

   private void saveResults()
   {
     File file=new File("results.dat");
     char buf[]=null;
     NormalDistribution norm=new NormalDistribution(energy);
     double data[][]=new double[1][3];
     double mean=norm.getMean();
     double var=norm.getVariance();
     System.out.println("Energy: "+mean+" Var: "+var);

     data[0][0]=N;
     data[0][1]=mean;
     data[0][2]=var;
     // 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.");
     }
   }
}