3.3 Variational Monte Carlo

Variational Monte Carlo (VMC) is based on a direct application of Monte Carlo integration to explicitly correlated many-body wavefunctions. The variational principle of quantum mechanics, derived in the following section, states that the energy of a trial wavefunction will be greater than or equal to the energy of the exact wavefunction. Optimised forms for many-body wavefunctions enable the accurate determination of expectation values.

3.3.1 The variational principle

The variational principle of quantum mechanics may be derived by expanding a normalized trial wavefunction, $\psi_{T}$, in terms of the exact normalized eigenstates of the Hamiltonian,

$$\psi_T=\sum_{i=0}^{\infty} c_{i} \psi_i \;\;\;,$$ (3.11)

where the expansion coefficients, $c_{i}$, are normalised

$$\sum_{i=0}^{\infty} \vert c_{i}\vert^2=1 \;\;\;.$$ (3.12)

The expectation of the many-body Hamiltonian, $\hat{H}$, may be evaluated

$$<\psi_T\vert\hat{H}\vert\psi_T > = <\sum_i c_i\psi_i\vert\hat{H}\vert\sum_j c_j\psi_j >$$

$$= \sum_i \sum_j c_{i}^{*} c_{j} <\psi_i\vert\hat{H}\vert\psi_j >$$

$$= \sum_i \vert c_{i}\vert^{2} \epsilon_{i} \;\;\;,$$ (3.13)

where

$$\epsilon_{i} = <\psi_i\vert\hat{H}\vert\psi_i > \;\;\;.$$ (3.14)

The expectation value of a trial wavefunction with the Hamiltonian must therefore be greater than or equal to the true ground state energy.

Variational calculations depend crucially on the form of trial wavefunction used. By selecting trial wavefunctions on physically motivated grounds, accurate wavefunctions may be obtained. Commonly, wavefunctions obtained from a Hartree-Fock or similar calculations are used and additional parameters added, building-in additional physics such as known limits and derivatives of the many-body wavefunction. The additional variational freedom is then exploited to further optimize the wavefunction. These techniques are briefly described in the remainder of this section, and in more detail in chapters 4 and 5.

3.3.2 Monte Carlo integration

Before the expectation value of a trial wavefunction with the many-body Hamiltonian may be computed, the integral must be transformed into a form suitable for Monte Carlo integration.

Trial wavefunctions, $\psi_T$, are dependent on the set of $N$ electron positions, ${\bf R}=\{{\bf r}_1,{\bf r}_2, \ldots {\bf r}_N\}$. The expectation value is given by

$$E=\frac{\int \psi_{T}^{*}\hat{H}\psi_{T} d{\bf R}} {\int \psi_{T}^{*}\psi_{T}d{\bf R}} \;\;\;,$$ (3.15)

which may be rewritten in an importance sampled form in terms of the probability density $\vert\psi_{T}\vert^{2}$:

$$E=\frac{\int \vert\psi_{T}\vert^{2}\frac{\hat{H}\psi_{T}}{\psi_{T}} d{\bf R}} {\int \vert\psi_{T}\vert^{2}d{\bf R}} \;\;\;.$$ (3.16)

The Metropolis algorithm samples configurations - sets of electron positions ${\bf R}$ - from the probability distribution $\vert\psi_{T}\vert^{2}$, and the variational energy is obtained by averaging the ``local energy'' $E_{L}$ over the set of configurations $\{{\bf R}\}$:

$$E_{L}({\bf R})=\frac{\hat{H}\psi_{T}({\bf R})}{\psi_{T}({\bf R})}$$ (3.17)

and

$$E=\frac{1}{N} \sum E_{L}({\bf R}_i) \;\;\;.$$ (3.18)

3.3.3 The local energy

The local energy, $E_{L}({\bf R})$, equation 3.17 is one of the central quantities in QMC methods. It occurs in both the variational and diffusion Monte Carlo algorithms and its properties are exploited to optimise trial wavefunctions.

The local energy has the useful property that for an exact eigenstate of the Hamiltonian, the local energy is constant. For a general trial wavefunction the local energy is not constant and the variance of the local energy is a measure of how well the trial wavefunction approximates an eigenstate. The spatially averaged variance of the local energy is therefore a quantity suitable for optimisation, and methods exploiting this observation are presented in chapter 5.

Determination of the local energy is one of the most computationally costly operations performed in QMC calculations. Application of the Hamiltonian to the trial wavefunction requires computation of the second derivatives of the wavefunction and the calculation of the electron-electron and electron-ion potentials. Efficient methods for the evaluation of $E_{L}({\bf R})$ are given in chapter 4.

3.3.4 Trial wavefunctions

The choice of trial wavefunction is critical in VMC calculations. All observables are evaluated with respect to the probability distribution $\vert\Psi_T({\bf R})\vert^2$. The trial wavefunction, $\Psi_T({\bf R})$, must well approximate an exact eigenstate for all ${\bf R}$ in order that accurate results are obtained. Improved trial wavefunctions also improve the importance sampling, reducing the cost of obtaining a certain statistical accuracy.

Quantum Monte Carlo methods are able to exploit trial wavefunctions of arbitrary forms. Any wavefunction that is physical and for which the value, gradient and laplacian of the wavefunction may be efficiently computed can be used.

The power of Quantum Monte Carlo methods lies in the flexibility of the form of the trial wavefunction. In early studies of bosonic $^{\rm 4}$He by McMillan [21] the wavefunction was taken to be a Jastrow or two-body correlation function, [22]

$$\psi=\exp \left[ \sum_{i<j}^{N} -u(r_{ij}) \right] \;\;\;.$$ (3.19)

The function $u$ was chosen to miminise the energy of the system under consideration, by choosing $u$ to increase the probability of particles being at a distance that minimises their interaction energy. Variations on this idea have been successfully applied to fermionic systems by multiplying a determinantal wavefunction by a two-body or higher body correlation functions. [23,24] Well chosen correlation functions include correlation effects more efficiently than CI-based approaches.

It is important that the trial wavefunction satisfies as many known properties of the exact wavefunction as possible. A determinantal wavefunction is correctly anti-symmetric with respect to the exchange of any two electrons. An additional local set of constraints which may be readily imposed are for electron-electron and electron-nucleus coalescence. These constraints are the ``cusp conditions''[25], and are a constraint on the derivatives of the wavefunction. For particle-particle coalescence, it may be shown that [25]

$$\left.\frac{d\Psi}{dr}\right\vert _{r=0} = \zeta \Psi\vert _{r=0} \;\;\;,$$ (3.20)

where $\Psi$ is the many-body wavefunction and $r$ is either an electron-electron or electron-ion separation. $\zeta$ is $-Z$ for coalescence of an electron with a nucleus of charge $Z$, and $1/2$ or $1/4$ for spin parallel and anti-parallel electron-electron coalescence respectively.

A commonly used and simple form of Jastrow factor suitable for solids is

$$u(r) = \frac{A}{r}(1-e^{-{r\over F}}) \;\;\;,$$ (3.21)

where $F$ is parameterised in terms of $A$ and chosen to satisfy the electron-electron cusp conditions,

$$\left.\frac{du}{dr}\right\vert _{r=0} = \left\{ \begin{array... ... & \mbox{for anti-parallel spins} \end{array} \right. \;\;\;,$$ (3.22)

and the complete trial wavefunction is then taken to be

$$\Psi({\bf R})=D({\bf R}) \exp \left[ \sum_{i<j}^{N} -u(r_{ij})\right] \;\;\;,$$ (3.23)

where $D$ is a determinant of Hartree-Fock or DFT orbitals. The variational freedom in the value of $A$ is used to optimise the trial wavefunction.

The development of better approximations to exact many-body wavefunctions, and practical methods for obtaining them will remain an important area of research due to the improved levels of accuracy and efficiency that result. The Jastrow factors used in this thesis contain of the order of 10-200 parameters and details of parameter selection and optimisation are given in the next chapter.

3.3.5 The VMC algorithm

The VMC algorithm consists of two distinct phases. In the first a walker consisting of an initially random set of electron positions is propagated according to the Metropolis algorithm, in order to equilibrate it and begin sampling $\vert\Psi\vert^2$. In the second phase, the walker continues to be moved, but energies and other observables are also accumulated for later averaging and statistical analysis.

To simplify the notation, $\Psi({\bf r}^\prime)$ is used to indicate the wavefunction with a single electron moved to point ${\bf r}^\prime$ and all other electrons unmoved:

1. Equilibration phase:
   1. Generate initial configuration using random positions for the electrons
   2. For every electron in the configuration:
      1. Propose a move from ${\bf r}$ to ${\bf r}^\prime$
      2. Compute $w=\vert\Psi({\bf r}^\prime)/\Psi({\bf r})\vert^2$
      3. Accept or reject move according to Metropolis probability $\mathrm{min}(1,w)$
   3. Repeat configuration moves until equilibrated
2. Accumulation phase:
   1. For every electron in the configuration:
      1. Propose a move from ${\bf r}$ to ${\bf r}^\prime$
      2. Compute $\vert\Psi({\bf r}^\prime)/\Psi({\bf r})\vert^2$
      3. Accumulate the contribution to the local energy, and other observables, at ${\bf r}^\prime$ and at ${\bf r}$, weighted by the Metropolis acceptance and rejection probabilities respectively. [26]
      4. Accept or reject move according to Metropolis acceptance probability
   2. Repeat configuration moves until sufficient data accumulated

In this algorithm, the electrons are moved individually and not as a whole configuration. This improves the efficiency of the algorithm in larger systems, where configuration moves require increasingly small steps to maintain the acceptance ratio. Observables are also accumulated on a per electron basis, weighted by the acceptance and rejection probabilities:

$$<O>=\frac{1}{M}\sum_{i=1}^{M} \left[p_{i}O_{i}({\bf r}^\prime)+(1-p_{i})O_{i}({\bf r})\right] \;\;\;,$$ (3.24)

where $p_{i}$ indicates the Metropolis acceptance probability of electron $i$ and $O_{i}({\bf r})$ its contribution to the value of the observable, and there are $M$ moves in the simulation. This formula, which can be applied to any observable, has improved statistics compared with averages computed with reference to the unmoved (or moved) positions alone. [26]