VMC Hydrogen Optimizing
From PIMC++
Return to: VMC Tutorial
Contents |
[edit] Optimizing the wavefunction: A First Attempt
Test File: OptimizeTest.py
Now that we can evaluate the energy, we should use this information to allow us to choose more effective wavefunctions. We will approach this in two different ways. We will examine both quantitatively optimizing paramaters in the wave functions as well as qualitatively different forms of wavefunctions.
We will start by writing some code that allows us to numerically optimize the paramaters of our two gaussians.
We know that we can change the paramaters of our wavefunction by doing H2.SetParams([alpha])
Let's use this fact to write a loop that samples over values of alpha, runs VMC on them , computes the appropriate statistics and keeps the average energy and average error in a list (i.e. have a function like
def Optimize(H2):
optimizeList=[]
loop of alpha's:
Calculate the average energy and error for this alpha
optimizeList.append(alpha,energy,error)
return (optimizeList,bestAlpha)
Helpful hint: A useful way to loop in python is to call for alpha in numpy.arange(0,4,0.1) which will loop over numbers between 0 and 4 with a step size of 0.1
We then call our optimizer like
#set up wave function as always (optimizeData,bestAlpha)=Optimize(H2) print "The optimum alpha is ",bestAlpha pylab.errorbars(optimizeList[:,0],optimizeList[:,1],optimizeList[:,2])
to test it. (Hint: You should be near an optimum value for alpha of around 0.5) From your plot, you should be able to tell if you've selected enough VMC steps to give you small enough error bars that mean your optimum value of alpha is meaningful.
This optimization is a little slow. Ideally we would like to find the optimum value of alpha somewhat faster. One approach is to calculate a couple points and then fit them to a curve. In our test file, we show you how this is done by doing it after we return the optimizeList.
Although this isn't a necessary step for the rest of the tutorial, if you find you are moving quickly, it is a nice exercise to modify your Optimize function to use this technique.
[edit] Calculating the bond length and the vibrational energy: Attempt 1
Test File: TestBondLength.py
We now want to calculate the bond length of our system. In order to calculate the vibrational energy and the average bond length, we will evaluate the energy at a variety of different bond lengths with this fixed optimized value of alpha (in a moment we will do something better). Write a function
def BondLength(): bondList=[] #write stuff here bondList.append([bondLength,energy,error]) return (optimum_bondLength,bondList)
that returns the optimum bond length and the list of alpha's and energy.
We will now call it using
... (opBond,bondList)=BondLength() print "The optimum bondlength is ",opBond pylab.errorbar(bondList[,:0],bondList[:,1],bondList[:,2])
which should give us our answer and plot the values.
Let's now improve on this. Ideally we want the value of alpha to be optimized for each independent bond length. You should be able to couple your two functions (Optimize, BondLength) If you haven't sped up your alpha optimization up yet, this might turn out to be a little slow. Another option is to replace your calls to VMC with a call to fast.VMC which is a version of VMC we've written that does exactly what yours does (but is written in C++ and so is an order of magnitude faster). See here for a discussion of python and speed.
You may notice that the minimum bond length is reasonably good. On the other hand, the binding energy is -1.0 and so you will find that there actually isn't any binding
[edit] Adding a jastrow
Test File: WFJastrowTest.py
We will now make a change to our wavefunction by adding a jastrow.
We will start working with H2JastrowClass. Since the way in which the local energy is calculated doesn't change (since we are using finite differences), make a copy of your current LocalEnergy function from your H2Class
Now we will write the wavefunction for our class.
The appropriate wave function is of the form:
where 
and 
.
Let's change our WaveFunction in our H2JastrowClass to reflect this change. Our list of params will now contain 
To verify our wavefunction is working, we will do the following checks
-
-
- the VMC integral is approximately -1.14 with an acceptance ratio of 0.38.
by calling
H2=H2JastrowClass()
H2.SetParams([0.5,0.3])
H2.SetBondLength(1.4)
coords=numpy.zeros((2,3),float)
coords[0,:]=[1.0,0.5,0.3]
coords[1,:]=[-0.2,0.1,-0.1]
energyVal=H2.WaveFunction(coords)
print "Your wavefunction is ",energyVal
print "Your local energy is ",H2.LocalEnergy(coords)
energyList=VMC(H2,10000)
pylab.plot(energyList)
pylab.savefig("JastrowEnergy.png")
pylab.show()
This will also make a trace plot of the energy. You will notice that the spikes that were so prevalant when we were using the wavefunction without the jastrow are no longer manifested here. This is because the jastrow satisfies the cusp condition.
[edit] Bond Length and Optimization: Take 2
Test File: OptimizeJastrowTest.py
You could again optimize these two parameter in a similar method as earlier. This will be fairly slow (we will do it anyway). Although some of the reason for this is that we are using python, larger systems will run into problems (speed and otherwise) independent of the language. If you have time, you can start learning how to better deal with VMC optimization. Otherwise, the next step will be to modify your Optimize function to successfully optimize the parameters alpha and beta. Copy it over and call it
OptimizeJastrow. Modify it to optimize both parameters: alpha and beta.
(because python is slow, you will probably need to switch over to fast.VMC) Instead of returning a listing you should now have it return three 2-dimensional numpy arrays. They should respectively be "a grid of alpha values", "a grid of beta values" and "a grid of Energy values." (i.e. your function should contain the arrays
def OptimizeJastrow(H2):
bestAlpha=0.0
bestBeta=0.0
N=10 # amount of grid points you want on one side
alphaMat=numpy.zeros((N,N),float)
betaMat=numpy.zeros((N,N),float)
EMat=numpy.zeros((N,N),float)
...
return (alphaMat,betaMat,EMat,(bestAlpha,bestBeta))
OptimizeJastrowTest.py will return a 2d color plot of the values of these parameters. (Useful note: Set N to be small (say 3) and the number of VMC steps you do to be small (say 100) the first time through so you make sure you don't have any bugs or typos in your code. Since python only catches bugs at runtime, you don't want to do an entire 20 minute 2d optimization and then discover you have a typo.
Once you've produced a 2d color plot, you could this new optimization method into your BondLengthJastrow() function. (Of course, any steps you can take to speed up your optimization before this point would be helpful).
Copy and modify your function BondLength() to be BondLengthJastrow()
Again, at each step solve for the optimum parameters alpha and beta and then calculate the bond length at these parameters.
We haven't supplied a test function for this. Write your own test file that plots the energy versus bond length. You should find that the system actually bonds this time.
[edit] Extension 1: Qualitative Change of Wave Functions
As the base for our wave function, we've been assuming that there are gaussians on the bond orbital. Certainly, as the bond gets sufficiently large, it makes more sense to locate the gaussians on the bonds and not on the center of the orbitals. Earlier in this tutorial, we showed an example of a wave function where that was the case. You could modify your code so that it calculates both wavefunctions and then chooses the one with the minimum energy. One could also imagine having a parameter that slowly interpolates between these two.
