The Method

The Schrödinger equation for the He atom in atomic units is,

$$\left( - \frac{1}{2} \nabla_{r_1}^{2} - \frac{1}{2} \nabla_{... ..., \bi{r}_{2} ) = E \Psi ( \bi{r}_{1}, \bi{r}_{2} ) \; \; \; ,$$ (4.35)

where $\bi{r}_{1}$ and $\bi{r}_{2}$ are the position vectors of the two electrons, $r_{1} = \vert \bi{r}_{1} \vert$, $r_{2} = \vert \bi{r}_{2} \vert$, and ${r}_{12} = \vert \bi{r}_{1} - \bi{r}_{2} \vert$. Energies are in units of the Hartree energy (1 Hartree = 2 Rydbergs) and distances are in units of the Bohr radius. The ground state spatial wavefunction is symmetric under exchange of the two electrons (the required antisymmetry is taken care of by the spin part of the wavefunction, which we can forget about otherwise).

The expression for the energy expectation value of a particular trial wavefunction, $\Phi_T( \bi{r}_1 , \bi{r}_2 )$, is,

$$E_{T} = { \int\ldots\int\Phi_T^* \hat{H} \Phi_T \d^{3}r_{1... ...int\ldots\int\Phi_T^* \Phi_T \d^{3}r_{1} \d^{3}r_{2} } \quad .$$ (4.36)

In the variational Monte Carlo method, this equation is rewritten in the form,

$$E_{T} = \int \ldots \int \left( \frac{1}{\Phi_T} \hat{H} \... ...) f\left(\bi{r}_1, \bi{r}_2\right) \d^3 r_1 \d^3 r_2 \quad ,$$ (4.37)

where

$$f( \bi{r}_1, \bi{r}_2 ) = { \Phi_T^*( \bi{r}_1, \bi{r}_2 ) ... ...2) \over \int \ldots \int \Phi_T^* \Phi_T \d^3 r_1 \d^3 r_2 }$$ (4.38)

is interpreted as a probability density which is sampled using the Metropolis algorithm. Note that the Metropolis algorithm only needs to know ratios of the probability density at different points, and so the normalisation integral,

$$\int \ldots \int \Phi_T^* \Phi_T \d^3 r_1 \d^3 r_2 \quad ,$$ (4.39)

always cancels out and does not need to be evaluated.

The mean of the values of the ``local energy'',

$${1\over\Phi_T} \hat{H} \Phi_T \quad ,$$ (4.40)

at the various points along the Monte Carlo random walk then gives an estimate of the energy expectation value. By the variational principle, the exact energy expectation value is always greater than or equal to the true ground state energy; but the Monte Carlo estimate has statistical errors and may lie below the true ground state energy if these are large enough. Anyway, the better the trial wavefunction, the closer to the true ground state energy the variational estimate should be.

In this project you are given a possible trial wavefunction,

$$\Phi = e^{- 2 r_1} e^{- 2 r_2} e^{r_{12}/2} \quad .$$ (4.41)

Use the variational Monte Carlo technique to calculate variational estimates of the true ground state energy.

Before you start programming, you will need the analytic expressions for the local energy. This involves some nasty algebra, but the answer is,

$${1\over\Phi} \hat{H} \Phi = - \frac{17}{4} - { \bi{r}_1 \c... ...12}} + { \bi{r}_2 \cdot {\bf r}_{12} \over r_2 r_{12}}\quad.$$ (4.42)

where $\bi{r}_1$ and $\bi{r}_2$ are the coordinates of the 2 atoms relative to the nucleus and $\bi{r}_{12} = \bi{r}_2 - \bi{r}_1$.

The Monte Carlo moves can be made by generating random numbers (use a library routine to do this) and adding them to the electron coordinates. I suggest that you update all six electron position coordinates $( x_{1}, y_{1}, z_{1}, x_{2}, y_{2}, z_{2})$ each move, and so you will need six random numbers each time. The accepted lore is that the Metropolis algorithm is most efficient when the step size is chosen to keep the acceptance probability close to 0.5. However, the method should work in principle no matter what the step size and you should try a few different step sizes to confirm that this is indeed the case. The starting positions of the two electrons can be chosen randomly, but remember that the Metropolis algorithm only samples the probability distribution exactly in the limit as the number of moves tends to infinity. You will therefore have to throw away the results from the moves near the beginning of the run and only start accumulating the values of the local energy once things have settled down. You should experiment to find out how many moves you need to throw away.

The statistical errors in Monte Carlo calculations decrease like $1/\sqrt{N}$, where $N$ is the total number of moves after the initial equilibration period. The errors therefore improve only slowly as the length of the run is increased. You will not be able (and should not attempt) to attain great accuracy. However, you should think hard about the magnitude of the statistical errors involved. Calculating the variance of the values in the list of energies accumulated during the random walk is easy and you should certainly do it.