Writing a simple VMC code in Python

[edit] Variational Monte Carlo

Variational Monte Carlo is a method to evaluate the integral

To accomplish this, we use markov chain monte carlo. This works in the following manner:

1. Choose a new trial location for the particles
2. Evaluate the ratio squared of the wavefunctions for the new and old trial locations
3. If the ratio squared is greater then some random number, then keep the new trial location, otherwise keep the old trial location.

Typically the integration part of Variational Monte Carlo (what we are starting our focus on) is coupled with an optimization technique for the parameters. For experts, you can look here to see how that entire structure goes together: VMC Object Diagram

[edit] Trial wave functions for a hydrogen molecule

We will be writing a simple variational monte carlo code to calculate the bond length of a Hydrogen molecule.

For sake of simplicity, this lab will focus on the hydrogen molecule. Since it has only two electrons (one spin up and one spin down), we need not be concerned with computing determinants. We will be optimizing two very simple trial wave functions, and . Let our molecule be centered at the origin with protons at .

The first trial function is a simple bond-centered Gaussian with width controlled by the parameter ,

where is to be optimized variationally.

The second trial function add electron correlation and the correct cusp condition through the use of Jastrow factors. Let us define

where

where , are choses to satisfy the electron-electron and electron-proton cusp conditions, respectively. The parameter gives additional variational freedom through the relations,

The data structure that stores the coordinates of the ions (electrons) will be a 2d numpy array and called coords (ions). The first element will be particles and the second element will be coordinates. So to access the second particle you will call coords[1] and to access the y element you would call coords[1,1] (note that in python, all arrays start counting from 0)

[edit] Working procedure

In your VMC directory, you should find VMCHelp.py, which has will have useful functions for you throughout the tutorial) and the VMC_template.py (which is the file you will mainly work in). There will also be a number of test files (linked immediately after their relevant section name) that you will run throughout the tutorial. These files are all a few (<10) lines and are tests you can run to verify that pieces of your variational monte carlo code is working. Throughout this tutorial, there will be a number of code snippets like this:

Here is a code snippet

with a filename at top to indicate which file you would look in to modify or write this piece of code.

[edit] Building the wavefunction object

Test File: WFTest.py

We will start by building an object for our wavefunction. The paradigm of object oriented programming is a very powerful one. It will allow us to couple the computation of the wave-function (and the local energy) with the data (parameters, ion locations) that goes along with it. This will allow us to simplify the optimization and VMC runs. Our object will look like the following:

def __init__(self):
    self.I = zeros((2,3),float)
    self.SetBondLength(1.4)
    self.SetParams([0.32])
    self.name = 'H2_WF1'
 
  def SetParams(self,params):
    self.alpha = params[0]
 
  def SetIons(self,ions):
    self.ions=ions.copy()
 
  def SetBondLength(self,BondLength):
    self.BondLength=BondLength
    self.I[0] = [-0.5*length, 0.0, 0.0]
    self.I[1] = [+0.5*length, 0.0, 0.0]
 
  def WaveFunction(self,R):
    # Write code to evaluate your wavefunction here
    # Recall that self.alpha stores the WF parameter
    return 0.0 # you will return something other then 0
 
  def LocalEnergy(self,R):
    KE=-0.5*LaplacianPsiOverPsi(R,self.WaveFunction)
    V=0.0 #change this to actually calculate V
    return V+KE

You should note a couple things. The class has a number of internal variables (params, ions). These internal variables are able to be accessed by prepending self to them (as in self.params). The class also has a number of functions. All functions that live inside a python class must take the parameter self as the first parameter. Add a line to SetParams that stores the bond length as an internal variable in the class (remember you would have to access it using self.BondLength). Another is that in python we must make copies of objects we want to clone because the equals operator, when operating on two objects (left=right), simply makes the left object point towards the right. Afterwards changes to the left object will result in changes to the right, they are the same object.

Once we have finished writing the class, we then proceed to initialize it as follows:

H2=H2Class()
H2.SetParams([0.5])
H2.SetBondLength(1.4)

and we can evaluate the wavefunction by calling H2.WaveFunction(R) where R is our coordinate.

[edit] Evaluating the wavefunction

Now, that we have the scaffolding in place for our class, let us actually write the function to evaluate the wavefunction. The definition for the wavefunction is of the form:

VMC_template.py

def Wavefunction(self,coords):
   alpha=self.params[0]
   #you will fill this in

and it should return the value of the wavefunction (remember that you will eventually have to square this in order to evaluate the probability of making a move.) The parameter "ions" is acessible by calling self.ions and is a numpy array that specifies the location of the ions. The params is simply a list that contains the parameters of the function.

We now need to specify the form of the wavefunction that we will use. We will start by having a product of guassian's for each electron where the gaussian's are centered on the bond between them. (i.e. the wave-function are of the form where C is the center of the bond between the hydrogen atom and i is the product of the wavefunctions.)

To get you started, you might want to look at this example of a wavefunction that has guassians centered on the ions: PythonWaveFunctionGaussianIons

Useful things to know include the fact that

1. numpy.dot calculates dot products and
2. math.sqrt does square roots

Write your function to evaluate this quantity.

We will check your wavefunction by evaluating with a bond length of 1.4 and alpha=0.5 is 0.496585 through the lines

WFTest.py

R = numpy.array([[1.0, 0.5, 0.3], [-0.2, 0.1, -0.1]])
energyVal=H2.WaveFunction(R)
val = H2.WaveFunction(R)
corrval = 0.499574
print 'My val      = %9.6f' % (val)
print 'Correct val = %9.6f' % (corrval)

As you can see in WFTest.py, we set up the "coordinate", R by calling numpy.array() on a list of lists representing the particle positions. This produces a 2x3 array, which is what the rest of the code expects for coordinates.

Call ./WFTest.py to run this and verify your wavefunction is working! (If you don't get this value, stop here and fix the problem. Building on a broken base will make things much harder to debug in the future).

[edit] Building a monte-carlo code

Test File: VMCTest.py

Now that we can evaluate the wavefunction let's write a VMC function

VMC_template.py

def VMC(WFClass,numSteps):
   #pseudoCode:
   #Looping over steps:
      #choose new cooridinates
      #evaluate the wave function on these new cooridinates
      #if (newWaveFunction^2/oldWaveFunction^2>random_number):
         #accept the move making sure coords and oldWaveFunction has the up to date information
      #else:
         #reject the move making sure coords and oldWaveFunction has the up to date information

which will be called in the following way:

VMCTest.py

#first set up our class
H2=H2Class()
H2.SetParams([0.5])
H2.SetBondLength(1.4)
#then call the VMC
numSteps=100
VMC(H2,numSteps)

Write your VMC function. Useful things to know include the fact that

1. random.random() will return a random number uniformly between 0 and 1
2. to square a number do a**2 (the ^ will not work)

The body of the function currently has some pseudocode to get you going.

Congratulations! You now have a (pretty useless) quantum monte carlo code.

[edit] Calculating the energy

Test File: EnergyTest.py

Now, we wish to calculate the energy. The energy is not only an important observable in it's own right, but being able to calculate it is the key ingredient to the successful usage of variational monte carlo. We know that the exact ground state wave function has the minimum energy over all possible wave functions. Consequently, by comparing the energy of different wavefunctions, you are able to establish which one is a better representation of the ground state of your system.

To accomplish this, we need to evaluate The term is relatively easy to evaluate as it just reduces to

(for a system of electrons (e) and ions(I)). While in a real code we would want to use an analytic expression for the calculating , we will save you time by giving you a routine which computes it with finite-differences: LaplacianPsiOverPsi(coords,WaveFunction). This will also allow you to quickly play with wave-function variants without recalculating the laplacian for each variant. If you want your code to go faster, you can write your own laplacian function, but it is a good idea to use a numerical derivative as a check.

Using this information we will add a function to our object of the form

VMC_template.py

 
def LocalEnergy(self,coords):
   KE=-0.5*LaplacianPsiOverPsi(coords,self.WaveFunction)
   V=...

to our PathClass.

You should write this function now. Remember, to access the ion locations you have to use self.ions and to access the coordinates use coords.

To verify that your function works correctly, we use the following chunk of code:

EnergyTest.py

H2=H2Class()
H2.SetParams([0.5])
H2.SetBondLength(1.4)
R=numpy.zeros((2,3),float)
R[0]=[1.0,0.3,0.2]
R[1]=[2.0,-0.2,0.1]
print H2.LocalEnergy(R)

You should get an answer of -1.819 (maybe?? 0.932??).

Now add some code to your VMC function to append the local energy to the EnergyList every ten steps.

We test this VMC run by having the lines

EnergyTest.py

energyList=VMC(H2,100000)
print CalcStatistics.stats(energyList)
pylab.plot(energyList)

where the second line is a call to the statistics package you've written and the third call is to pylab to plot the energy. Run EnergyTest.py. You should see a plot like:

You should notice a couple things.

1. It is quite spiky
2. The data has some equilibration
3. your answer should be around -0.86

Now that we can successfully integrate a wavefunction, we will start down our path of optimizing wavefunction: VMC Hydrogen Optimizing

Retrieved from "http://www.qmcwiki.org/index.php/Writing_a_simple_VMC_code_in_Python"