Computational Physics

12.3 Diffusion Monte Carlo 391
Next we analyse the helium atom using the diffusion Monte Carlo method. This
is less successful. The reason is that writing the time-evolution operator as a product
of a kinetic and potential energy evolution operator

e−

τ (K+V−ET)=e−

τKe−

τ (V−ET)+O(

τ^2 ) (12.57)

is not justified when the potential diverges, as is the case with the Coulomb potential
atr=0. Formally, this equation is still true, but the prefactor of theO(
τ^2 )term
diverges. However, even if the potential does not diverge but varies strongly, the
statistical efficiency of the simulation is low. This is due to the fact that if a walker
moves to a very favourable region, it will branch into many copies. But these are
all the same, and together they form a rather biased sample of the distribution in
that region. It requires some time before they have diffused and branched in order
to form a representative ensemble. Frequent occurrence of such strong branching
events will degrade the efficiency considerably. Quite generally one can say that
the efficiency increases with the flatness of the potential.
There exist, in principle, two ways to solve the divergent potential problem.
The first one consists of finding a better alternative to the simple approximation
to the time-evolution operator than in(12.57). Such approximations have been
devised and we shall consider these in the context of path-integral Monte Carlo
(seeSection 12.4). The common procedure, however, is to use aguide function,
which transforms the original Schrödinger equation into a new one with a flatter
potential, just as in the case of the Fokker–Planck variational Monte Carlo method.
This method will be described in the next section.

12.3.3 Guide function for diffusion Monte Carlo

As we have just seen, the diffusion Monte Carlo method causes problems if the
potential is unbounded, and this is the case in almost every many-particle system.
Sampling some other function instead of the ground state wave functionψmight
cure this problem.
A suitable function isρ(R,τ)=ψ(R,τ)T(R)whereT(R)is some trial func-
tion which models the exact wave function in a reasonable way. It turns out thatρ
satisfies a Fokker–Planck type of equation:

∂ρ(R,τ)

∂τ

=

1

2

∇R[∇R−F(R)]ρ(R,τ)−[EL(R)−ET]ρ(R,τ). (12.58)

Here, the ‘force’F(R)is again given as 2∇RT(R)/T(R). This form differs
from (12.49) because (12.58) is not a ‘pure’ Fokker–Planck equation: it contains a
‘potential term’EL(R)−ET. The ‘local energy’EL(R)is given as usual by

EL(R)=

HT(R)

T(R)

=

−∇^2 T(R)/ 2 +V(R)T

T(R)

. (12.59)