next up previous contents
Next: A DMC Algorithm Up: Diffusion Quantum Monte Carlo Previous: Diffusion Quantum Monte Carlo

The Method


The basis of DMC is to write the Schrödinger equation in imaginary time, taking


tex2html_wrap_inline6333 can be expanded as:




The tex2html_wrap_inline6335 's and tex2html_wrap_inline5849 '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.(gif),


noting that tex2html_wrap_inline6339 is the imaginary-time evolution operator. Furthermore, expressing tex2html_wrap_inline6333 in the form of Eq.(gif) implies that


Letting tex2html_wrap_inline6343 be tex2html_wrap_inline6345 and hence tex2html_wrap_inline6347 , the evolution time for the system implies that, as long as tex2html_wrap_inline6345 is not orthogonal to the ground-state, tex2html_wrap_inline6351 , then


where tex2html_wrap_inline6353 is the ground-state energy. This can be seen by remembering that all other states have higher energies, tex2html_wrap_inline5849 , than the groundstate energy tex2html_wrap_inline6353 , and will therefore decay away faster. In the tex2html_wrap_inline6359 representation, Eq.(gif) becomes:




We now introduce an arbitrary energy offset term tex2html_wrap_inline6195 , such that the imaginary-time Schrödinger equation, Eq.(gif), is recast as:


Then if tex2html_wrap_inline6195 is adjusted to be the true ground-state energy, tex2html_wrap_inline6353 , 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.(gif), gives:


This is just a 3N-dimensional diffusion equation, with tex2html_wrap_inline6369 playing the role of the density of diffusing particles. The tex2html_wrap_inline6371 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.(gif) 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 tex2html_wrap_inline6373 , can diverge to tex2html_wrap_inline6375 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 gif).

We will follow the scheme of Ref.[29] for the introduction of importance sampling. The first step is to introduce a guiding function, tex2html_wrap_inline6377 . We now define a new distribution tex2html_wrap_inline6379 which, if tex2html_wrap_inline6369 satisfies the Schrödinger equation, is also a solution of the Schrödinger equation. Substituting tex2html_wrap_inline6383 into Eq.(gif) yields


Where tex2html_wrap_inline6385 can now be interpreted as a ``quantum force''. We follow [30] in defining this term as


and tex2html_wrap_inline6387 , 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 tex2html_wrap_inline6383 . The branching term is proportional to the ``excess local energy'', tex2html_wrap_inline6391 , which with a good choice of tex2html_wrap_inline6393 need not become singular when tex2html_wrap_inline6373 does. To control branching we need to choose tex2html_wrap_inline6393 such that tex2html_wrap_inline6399 is everywhere as smooth as possible, i.e. we want as little variance as possible in tex2html_wrap_inline6401 . Methods for optimising trial wavefunctions with exactly this property are described in detail in chapter gif. In general, the trial wavefunction, Eq.(gif), 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, tex2html_wrap_inline6393 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 tex2html_wrap_inline6393 to ensure that f is always of the same sign (see section gif). The accuracy of the position of the nodes in tex2html_wrap_inline6393 therefore determines how good the estimate of the ground-state energy, tex2html_wrap_inline5997 , is. This can be seen by considering the fact that at long (imaginary) times the distribution tex2html_wrap_inline6383 approaches tex2html_wrap_inline6415 , up to the constraint (imposed by the fixed-node approximation) that tex2html_wrap_inline6417 must vanish at the nodes of tex2html_wrap_inline6377 . This implies that the long-time limit is the true fermionic ground-state if and only if the nodes of tex2html_wrap_inline6377 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.(gif) can be written in integral form, in doing this we follow the procedure of Ref. [30]:


where tex2html_wrap_inline6423 is the Green's function for the case tex2html_wrap_inline6425 . The energy shift tex2html_wrap_inline6427 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 tex2html_wrap_inline6431 .

The three terms on the left-hand side of Eq.(gif) describe respectively, diffusion, drift and growth/decay. An approximate Green's function can be formed, with an error of tex2html_wrap_inline6433 for small tex2html_wrap_inline6435 , 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 gif). 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, tex2html_wrap_inline5881 ) from the distribution f, using


next up previous contents
Next: A DMC Algorithm Up: Diffusion Quantum Monte Carlo Previous: Diffusion Quantum Monte Carlo

Andrew Williamson
Tue Nov 19 17:11:34 GMT 1996