# What are the types of Quantum Monte Carlo?

Similar to:

I would like to ask this time, if people can summarize the types of QMC in up to 3 paragraphs:

## Stochastic Series Expansion (SSE) Monte Carlo

Theory: SSE is a finite-temperature, discrete-time technique that works well for quantum spin problems (e.g. Heisenberg model) and other lattice Hamiltonians in any number of dimensions. The method works by expanding the partition function in a Taylor series
$$\tag{1} Z = \mathrm{Tr}[ \rho] = \mathrm{Tr}[e^{-\beta H}] = \mathrm{Tr} \sum \limits_{\alpha_0} \sum \limits_{n=0}^\infty \left\langle \alpha_0 \left| \frac{(-\beta H)^n}{n!} \right| \alpha_0 \right\rangle$$ and then inserting a complete set of states in some basis between each term in the Taylor expansion
$$\tag{2} Z = \sum \limits_n \frac{(-\beta)^n}{n!} \sum \limits_{\alpha_0} \sum \limits_{\alpha_1} ... \sum \limits_{\alpha_{n-1}} \langle \alpha_0 | H | \alpha_1 \rangle \langle \alpha_1 | H | \alpha_2 \rangle ... \langle \alpha_{n-1} | H | \alpha_0 \rangle$$

The resulting matrix elements $$\langle \alpha_{i} | H | \alpha_{i+1} \rangle$$ are usually simple to evaluate.

Algorithm: The goal of the algorithm is to sample the sum in the previous equation. This is sum in an extremely high-dimensional space, perfect for importance sampling Monte Carlo. In general, the program stores a MC configuration as a starting state $$\alpha_0$$ and the list of local operators that act upon that state (the operator string). There are many different methods of updates, but they generally consist of two steps:

1. A diagonal update that adds and removes diagonal matrix elements, sampling the order ($$n$$) of the Taylor expansion.
2. An off-diagonal update which leaves $$n$$ fixed and samples the configurations $$\{\alpha\}$$. These are commonly called operator loop updates and there are many different algorithms for different situations.

Use case: SSE is efficient for lattice Hamiltonians at finite temperature (although $$T=0$$ can be reached as a limit). To measure any quantity that is diagonal in the working basis is each, since the measurements can be directly computed from the SSE configurations, e.g. $$\tag{3} \langle S^z \rangle = \mathrm{Tr} \left[ S^z e^{-\beta H} \right]$$ If an observable is off-diagonal in the simulation basis, but appears in the Hamiltonian, it can sometimes be calculated using some clever formula, e.g. $$\tag{4} \langle E \rangle = -\frac{\langle n \rangle }{\beta}$$ $$\tag{5} \langle C_v \rangle = -\frac{\langle n(n-1) \rangle }{\beta }$$

Computing arbitrary off-diagonal observables often requires complicated secondary sampling procedures.

References The main source for learning about SSE is this big review article by Sandvik: arXiv:1101.3281. There is also an excellent recorded lecture by Roger Melko from the 2010 Boulder Summer School (the link to the videos is right below the title).

## DMC (Diffusion Monte Carlo)

Theory. Consider the SchrÃ¶dinger equation in imaginary time $$\tau=it$$: $$-\hbar\frac{\partial\psi(x,\tau)}{\partial\tau}=\hat{H}\psi(x,\tau).$$ For a time-independent Hamiltonian $$\hat{H}$$, the $$\tau$$-dependence can be solved in a way analogous to the usual time dependence to obtain: $$\psi(x,\tau)=\sum_nc_n(0)e^{-E_n\tau/\hbar}\psi_n(x),$$ where $$\hat{H}\psi_n(x)=E_n\psi_n(x)$$. The function $$\psi(x,\tau)$$ at imaginary time $$\tau$$ is a sum over an exponentially decaying superposition of energy eigenstates with the exponential decay rate proportional to $$E_n/\hbar$$. This means that in the limit of large $$\tau$$: $$\psi(x,\tau\gg1)\simeq c_0(0)e^{-E_0\tau}\psi_0(x).$$ In this limit, the ground state $$n=0$$ is "projected out" of the initial state, because the corresponding exponential decay is the slowest one. Therefore, by evolving the system in imaginary time we can obtain the ground state of the Hamiltonian $$\psi_0(x)$$ as the long imaginary time limit.

Algorithms. So why is the method called "diffusion" Monte Carlo? The kinetic energy term of the Hamiltonian together with the imaginary time dependence is mathematically a diffusion equation, which is simulated using stochastic methods to evolve a collection of "walkers" or samples of the wave function. The potential term is then treated as a "branching" term, in which walkers are created or annihilated. Will the algorithm converge to the ground state? If the initial wave function has some overlap with the ground state, then the ground state will be projected out. Otherwise, the lowest energy state with non-zero overlap with the initial wave function will be projected out.

Reference. An excellent review paper of the method was published some time ago in Reviews of Modern Physics.

## FN-DMC (Fixed-node diffusion Monte Carlo)

Theory. See my answer about DMC. The only addition for FN-DMC is that the ground state of an arbitrary Hamiltonian will not be antisymmetrized, and therefore DMC will not converge to the fermionic ground state of interest in electronic systems. To force the system to project out the fermionic ground state, then the nodes of the wave function are fixed during the simulation to those of a fermionic wave function.

Algorithms. How are these nodes obtained? A typical approach is to first solve the problem using a different method (typically DFT or Hartree-Fock, but sometimes with post-SCF methods like CISD) and then fixing the nodal surface to that predicted by these methods. Although fixing the nodal surface introduces an approximation to the method, the overall methodology is still one of the most accurate to solve the electronic structure problem in periodic systems.

References. An excellent review paper of the method was published some time ago in Reviews of Modern Physics, with an emphasis on the fixed-node formulation and applications to solids. Codes implementing this method include CASINO, QMCPACK, and QWalk.