# Algorithms

## Monte Carlo Part

### Metropolis Algorithm

A generic Metropolis Monte Carlo algorithm draws samples from a target distribution $P(x)$ by following a random walk. The algorithm can be described in succinct steps:
- Initialise a starting point for the random walk $x_0$
- For each step in the chain, propose a change to the configuration $x_i \rightarrow x_{i+1}$
- Calculate the acceptance probability of this transition
$$
a = min\left(1 , \frac{P(x_{i+1})}{P(x)} \right)
$$
- Draw from a uniform distribution $u = U(0,1)$ and:
-- if $u < a$ accept the transition
-- if $u > a$ reject the transition, set $x_{i+1} = x_i$
- Calculate observables (position, energy, ...)

The output of this algorithm is a chain of points $({x_0 , x_1 , ... , x_n})$ exploring the target distribution $P(x)$ , and lists of observable quantities from the chain.

#### Acceptance Rate

The acceptance rate (how often transitions to new points are accepted in the chain) for this algorithm must be tuned carefully. A low acceptance rate will not explore the target space, whereas a high acceptance rate will ignore the target space entirely. This rate can be effectively tuned by the size of the proposed change $x_i \rightarrow x_{i+1}$ , and we aim for an acceptance rate around the theoretically optimal value of approximately 30%.

#### Thermalisation and Resampling

Initial points in MCMC walk may not have thermalised around the correct value - it is entirely possible that the initial configuration is far from the average. Removal of the first points in the chain will give more accurate results.

Subsequent steps in the Metropolis MCMC are correlated and therefore not independent. It is required to resample the points in the chain into independent points. This can be done by calculating the correlation time $t_{corr}$ within the chain, and then thinning or resampling the data so that statistics/further analyses are performed on a new set of points, which would be separated by $t_{corr}$ in the original chain.

To find $t_{corr}$ we compute the autocorrelation function for an observable $Z_{i}$ measured at the $i^{th}$ point in the chain, with mean $\mu$ for varying lag length $k$

$$
A(k) = \sigma^{-2} (Z_{i} - \mu)(Z_{i+k} - \mu)
$$

where $\sigma$ is the standard deviation of the observable distribution, and then search for the value of $k=t_{corr}$ such that $A(t_{corr}) \approx 0$.

Once $t_{corr}$ is found, the data is thinned so that one in every $t_{corr}$ values is 'kept' / observables are only calculated for those values.

An initial run for each problem will be performed to find the correct number of points to burn to ensure thermalisation, and the correlation length, and these will be passed as constant values to the full simulation.

## Quantum Part

### Quantum Harmonic Oscillator

The first system to be implemented within the code will be a quantum harmonic oscillator. This system has a simple analytical solution for both the wavefunction and the energy, and so will be useful as an early test against analytic results to check if the VMC algorithm is working properly.

The analytic solution of this problem (in atomic units) takes the form:
$$
\psi_{0}(x)=e^{-\frac{x^{2}}{2}}
$$
with Hamiltonian:
$$
\hat{H}=-\frac{1}{2}\frac{d^{2}}{dx^{2}}+\frac{1}{2}x^{2}
$$
and ground state energy:
$$
E_{0}=\frac{1}{2}
$$

A trial wavefunction can therefore be suggested:
$$
\psi^{\alpha}_{T}(x)=e^{-\alpha x^{2}}
$$

The local energy $E^{\alpha}_{L}(x)$ is evaluated to be
$$
E^{\alpha}_{L}(x)=\frac{\hat{H}\psi^{\alpha}_{T}(x)}{\psi^{\alpha}_{T}(x)}=\alpha+x^{2}\left ( \frac{1}{2}-2\alpha^{2} \right )
$$
and the probability distribution to be used in the Metropolis algorithm is
$$
\frac{\psi^{\alpha}_{T}(x')}{\psi^{\alpha}_{T}(x)}=e^{-2\alpha(x^{'2}-x^{2})}
$$

By iterating over walkers and updating positions with the specified probability, the energy of the trial wavefunction can be evaluated. The ground state wavefunction can then be found by varying $\alpha$ until a minimum energy is found. The code should evaluate $\alpha$ to be equal to 0.5, in agreement with the analytic solution of the problem.

### $H_{2}^{+}$

This problem requires that a single electron at position $\mathbf{r}$ be interacting with two nuclei, $A$ and $B$, frozen at positions $\mathbf{R}_{A}$ and $\mathbf{R}_{B}$ respectively. The Hamiltonian of this system is as follows:
$$
\hat{H}=-\frac{\hbar^{2}}{2m}\nabla^{2} - \frac{e^{2}}{4\pi\epsilon_{0}\left | \mathbf{r}-\mathbf{R}_{A} \right |} - \frac{e^{2}}{4\pi\epsilon_{0}\left | \mathbf{r}-\mathbf{R}_{B} \right |}
$$
The simplest way of forming a trial wavefunction for this system is to take a linear combination of two hydrogenic $1s$ orbitals:
$$
\psi^{c}_{T}(\mathbf{r})=c\phi_{1s}^{A}(\mathbf{r}) + \sqrt{1-c^{2}}\phi_{1s}^{B}(\mathbf{r})
$$
where $\phi_{1s}^{A}(\mathbf{r})$ and $\phi_{1s}^{B}(\mathbf{r})$ take the form
$$
\phi_{1s}^{A}(\mathbf{r})=\frac{1}{\sqrt{\pi}a_{0}^{3/2}}e^{-\frac{\left |\mathbf{r-\mathbf{R}_{A}} \right |}{a_{0}}} \ \ \ \ \ \ a_{0}=\frac{\hbar^{2}}{2m}
$$

If this trial wavefunction proves to be too inflexible for the purposes of the variational quantum Monte Carlo code, further hydrogenic wavefunctions of increasing principal quantum numbers could be added to the linear combination, for example:
$$
\psi^{\mathbf{c}}_{T}(\mathbf{r})=c_{1}\phi_{1s}^{A} + \sqrt{1-c_{1}^{2}}\phi_{1s}^{B} + c_{2}\phi_{2s}^{A} + \sqrt{1-c_{2}^{2}}\phi_{2s}^{B}
$$
where there are now two coefficients to be optimised, and the $2s$ hydrogenic wavefunction takes the form
$$
\psi^{A}_{2s} = \frac{1}{4\sqrt{2\pi}a_{0}^{3/2}}\left [ 2 - \frac{\left |\mathbf{r-\mathbf{R}_{A}} \right |}{a_{0}} \right ] e^{-\frac{\left |\mathbf{r-\mathbf{R}_{A}} \right |}{2a_{0}}}
$$

From this, formulae for the local energy $E_{L}^{c}(\mathbf{r})$ and the acceptance probability $a$ can be calculated and implemented directly as functions to be called by the Monte Carlo algorithm. If differentiability of the hydrogenic wavefunctions becomes a problem, it may be better to approximate the exponential by a Taylor expansion in a polynomial basis. This will be investigated as the code is implemented.

### $H_{2}$

This problem requires that two electron at positions $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ interact with two nuclei, $A$ and $B$, frozen at positions $\mathbf{R}_{A}$ and $\mathbf{R}_{B}$ respectively. The Hamiltonian of this system is as follows:
$$
\hat{H}=-\frac{\hbar^{2}}{2m}\nabla^{2} + \frac{e^{2}}{4\pi\epsilon_{0}} \left ( -\frac{1}{\left | \mathbf{r}_{1}-\mathbf{R}_{A} \right |} - \frac{1}{\left | \mathbf{r}_{1}-\mathbf{R}_{B} \right |} - \frac{1}{\left | \mathbf{r}_{2}-\mathbf{R}_{A} \right |} - \frac{1}{\left | \mathbf{r}_{2}-\mathbf{R}_{B} \right |} + \frac{1}{\left | \mathbf{r}_{1}-\mathbf{r}_{2} \right |} \right )
$$

In order to satisfy antisymmetry under exchange, the simplest trial wavefunction for this problem is
$$
\psi^{c}_{T}(\mathbf{r}_{1}, \mathbf{r}_{2})=c\phi_{1s}^{A}(\mathbf{r}_{1})\phi_{1s}^{B}(\mathbf{r}_{2}) + \sqrt{1-c^{2}}\phi_{1s}^{A}(\mathbf{r}_{2})\phi_{1s}^{B}(\mathbf{r}_{1})
$$
Similarly to the $H_{2}^{+}$ problem, this can be extended to hydrogenic wavefunctions of higher prinicpal quantum numbers if desired:
$$
\psi^{\mathbf{c}}_{T}(\mathbf{r}_{1}, \mathbf{r}_{2})=c_{1}\phi_{1s}^{A}(\mathbf{r}_{1})\phi_{1s}^{B}(\mathbf{r}_{2}) + \sqrt{1-c_{1}^{2}}\phi_{1s}^{A}(\mathbf{r}_{2})\phi_{1s}^{B}(\mathbf{r}_{1}) + c_{2}\phi_{2s}^{A}(\mathbf{r}_{1})\phi_{2s}^{B}(\mathbf{r}_{2}) + \sqrt{1-c_{2}^{2}}\phi_{2s}^{A}(\mathbf{r}_{2})\phi_{2s}^{B}(\mathbf{r}_{1}))
$$

Again, $E_{L}^{c}(\mathbf{r})$ and $a$ can be evaluated from this point to be called directly from the Monte Carlo algorithm.