This is a simple one dimensional VMC program that I use for didactic purposes.
It is written in standard FORTRAN 77 and implements the plain Metropolis algorithm to estimate the energy of a trial function for a user defined one dimensional quantum system.

You need

var.f               the main VMC program
var.h               include file
random.f            random number generator
seed.dat            random seed
file1               walkers file
system.f            define the quantum system you are studying
Makefile            Makefile
psi.inp             parameter of trial wave function
var.inp             main input file for var.
 

The code performs a simple VMC simulation on a user defined one dimensional system using a user defined trial wave function.
var.f  is the main program, it performs the simulation using the parameters read from the input file.

Example of input file

10     blocks
200    steps
500    walkers
2.5    side. size of the box
1      ifill: 0 = random 1 = from file1

As usual in Monte Carlo techniques, the simulation is divided into blocks, each block subdivided into steps.
In a single step all the walkers (configurations of the system) are updated using the Metropolis scheme, moving them within a box of a given side.

After compiling the program using the Makefile, you can run it with var < var.inp
The system under study is defined in the file system.f.
You must define the potential energy (in the present file is a simple harmonic oscillator)

      double precision function pot(x)

This function should return the value of the user defined potential at a given position.

You also need to write a subroutine

      subroutine ftrial(x,psi,psi1,psi2)

which, for a given x, returns the value of the trial wave function and its first and second derivatives.
Currently, the wave function coded in the program is a gaussian:  Exp[a x²]. Note that you do not need to normalize the wave function, since it is inessential for the calculation of the energy in a Monte Carlo program.
The program estimates the variational energy <H>, the square of the Hamiltonian <H²>, and some other
useful quantity.

The wave function, presently a gaussian, can have an optimizable parameter which is saved in file psi.inp
Of course when you study a new system, you do not know what is the best value of the trial parameter, so we need to optimize it in some way. There is a second program that takes care of the optimization called opt
 

Optimization

you need the files (in addition to the previouses):

opt.f
brent.f
mnbrak.f
opt.inp
 

opt.f is the main optimization program, and the optimization is started by the command  opt < opt.inp

A typical input file is

500             nwalks, number of walkers
1               iweight 1 = weighted estimates, 0 unweighted
1               iopt: 0 = energy 1 = sigma 2 = mu
0.5             eref: reference energy
0.00001         tol: tolerance for minimization

The most important parameter is iopt: if 0 the optimization process will try to minimize the energy, but beware that numerically this is a very unstable process. It is better to try to optimize the variance of the local energy, with iopt=1    S(H) = <H²> - <H>²
which is always bounded and positive. If we set iopt to 2  the code will try to minimize the variance of the local energy with respect to eref that is <(H-eref)²>. The remaining parameter, iweight determines if the optimization should be done by reweighting the walkers or not. This is a technicality and in most cases it does not matter. If you are confused, you should probabily read my review on VMC optimization.

Things to do

• Generalize the program for a particle of a generic mass m.
• The current version allows only one optimizable parameter into the wave function.
• Implement the Langevin algorithm instead of the plain Metropolis.
• Collect statistics on other quantities besides the energy.
• Collect the density of the walkers by dividing the real axes into bin. This is a way to reconstruct a numerical version of the wave function.

Things to try

• The Metropolis algorithm generates walkers that are statistically correlated. Compute the autocorrelation function for the energy as a function of the number of steps. Check how the number of walkers influence the exponential decay of the autocorrelation function. It is important that the number of steps in a block be much larger than the number of steps is takes the autocorrelation function to reach zero.
• For the Harmonic oscillator, plot the variational energy as a function of the variational parameter. Check the result against the analytical calculation
• Plot the Sigma of the Hamiltonian S(H) = <H²> - <H>² against the variational parameter
• Change the trial wave function. For example, use Exp[-a x⁴]
• Try to study the first excited state. Check that the symmetry of you trial wave function is the correct one.
• Change potential. You might try an anharmonic potential a x² + b x⁴ a>0, b>0
• If you allow, for the potential a x² + b x⁴ , a to be negative, you have a double well potential. How do you change your trial wave function to best accomodate the change?
• Implement the particle in a box potential. This can be tricky. How do you keep the walkers confined between the walls? do you need to hack the potential?
• Implement other systems. Needs a suggestion? a van del Waals potential, the Morse potential, etc ....