The basis of DMC is to write the Schrödinger equation in imaginary time, taking
can be expanded as:
where
The 's and 's are the energy eigenvectors and eigenvalues respectively of the time-independent Schrödinger equation.
We can write down the formal solution of the imaginary-time Schrödinger equation, Eq.(),
noting that is the imaginary-time evolution operator. Furthermore, expressing in the form of Eq.() implies that
Letting be and hence , the evolution time for the system implies that, as long as is not orthogonal to the ground-state, , then
where is the ground-state energy. This can be seen by remembering that all other states have higher energies, , than the groundstate energy , and will therefore decay away faster. In the representation, Eq.() becomes:
where
We now introduce an arbitrary energy offset term , such that the imaginary-time Schrödinger equation, Eq.(), is recast as:
Then if is adjusted to be the true ground-state energy, , the asymptotic solution is a steady-state solution.
We now use the fact that the Schrödinger equation in imaginary time looks like the diffusion equation. Explicitly writing out the Hamiltonian in Eq.(), gives:
This is just a 3N-dimensional diffusion equation, with playing the role of the density of diffusing particles. The term is a rate term, and describes the branching (or creation/annihilation) processes. The entire equation can be simulated by a combination of a diffusion and branching processes, in which the number of diffusing particles increases or decreases at a given point proportional to the density of diffusers and the potential energy at that point in configuration space.
It turns out that solving Eq.() this way is a very inefficient way to simulate the Schrödinger equation on a computer. This is because the branching rate, which is proportional to , can diverge to for systems of particles interacting via the Coulomb interaction. This leads to large fluctuations in the number of diffusing particles which leads to a large variance in the estimate of the energy. These fluctuations can be dramatically reduced by the introduction of importance sampling [28] in a similar way to the implementation in the VMC algorithm (see section ).
We will follow the scheme of Ref.[29] for the introduction of importance sampling. The first step is to introduce a guiding function, . We now define a new distribution which, if satisfies the Schrödinger equation, is also a solution of the Schrödinger equation. Substituting into Eq.() yields
Where can now be interpreted as a ``quantum force''. We follow [30] in defining this term as
and , the branching term as
which is defined in terms of the local energy of the guiding function
We now have a drift-diffusion equation for . The branching term is proportional to the ``excess local energy'', , which with a good choice of need not become singular when does. To control branching we need to choose such that is everywhere as smooth as possible, i.e. we want as little variance as possible in . Methods for optimising trial wavefunctions with exactly this property are described in detail in chapter . In general, the trial wavefunction, Eq.(), used as the input to a VMC calculation, makes a suitable choice of guiding wavefunction.
As well as reducing the fluctuations in the number of diffusing particles, also has another important role for fermionic systems. It determines the position of the nodes of the final wavefunction, due to the necessity of using the fixed-node approximation where the nodal structure of the exact groundstate wavefunction is assumed to be the same as the nodal structure of to ensure that f is always of the same sign (see section ). The accuracy of the position of the nodes in therefore determines how good the estimate of the ground-state energy, , is. This can be seen by considering the fact that at long (imaginary) times the distribution approaches , up to the constraint (imposed by the fixed-node approximation) that must vanish at the nodes of . This implies that the long-time limit is the true fermionic ground-state if and only if the nodes of correspond to the exact nodes of the ground-state wavefunction. The fixed-node energy is an upper bound to the exact fermionic energy[31]. We simulate the equation within only a small number of nodal regions. Each walker moves within one nodal region and rejects all moves that attempt to cross the nodal surface into another nodal region. This contradicts the requirement that we stipulated earlier that the walk must be ergodic. The tiling theorem [32, 33], however, states that the nodal regions of the true ground-state eigenfunction of a system of identical fermions are all related by permutation symmetry. Furthermore the nodal regions of the determinant of LDA eigenfunctions are also related by permutation symmetry because the wavefunction is the ground-state of a Hamiltonian, although it is not the many-body Hamiltonian of the system we are interested in. This means that we can simulate the equation within one nodal region and still be guaranteed to obtain the ``best'' variational result.
Eq.() can be written in integral form, in doing this we follow the procedure of Ref. [30]:
where is the Green's function for the case . The energy shift plays the role of an arbitrary time-dependent renormalisation, chosen in such a way that the probability distribution f remains finite and non-vanishing in the limit .
The three terms on the left-hand side of Eq.() describe respectively, diffusion, drift and growth/decay. An approximate Green's function can be formed, with an error of for small , by the product of Green's functions for diffusion, drift and growth/decay[34]:
In order to deal with the nodes in the fermion ground-state we must use the fixed-node approximation, as stated earlier (this and the various schemes to try to improve upon it are discussed in section ). What this actually entails is that if a move is such that a walker would cross the node then it is immediately rejected. In other words, the node acts as an infinite potential barrier. It should be noted that the use of importance sampling does not introduce any extra approximations beyond the fixed node approximation already mentioned. With the exception of the fixed node approximation it is possible to calculate exact energies (and expectation values of any other operator that commutes with the Hamiltonian, ) from the distribution f, using