E= = = . (4.31)

ψT

ψT ψ dr f dr

f

This has the advantage that the energy follows from an ensemble average

instead of from a rate of change of the number of walkers.

The time-dependent equation for f follows, after some manipulation, from

the de¬nition of f and the Schr¨dinger equation in imaginary time (4.12:

o

ˆ

‚f 1 HψT

[∇2 f ’ 2∇ · (f ∇ ln ψT )] ’ ’E

= f. (4.32)

‚„ 2m ψT

92 Quantum chemistry: solving the time-independent Schr¨dinger equation

o

Note that in the multiparticle case the mass-containing term must be sum-

med over all particles. The ¬rst term on the r.h.s. of this equation is a

di¬usion term; the second term is a drift term: it is (minus) the divergence

of a ¬‚ux density f u with a drift velocity

∇ ln ψT .

u= (4.33)

m

The term ’ ln ψT acts as a guiding potential, steering the walker in the

direction of large values of ψT . The third term replaces the strongly varying

potential V (r) by the much weaker varying local energy for the trial function.

QMC algorithms of this type su¬er from a time-step error and extrapo-

lation to zero time step is needed for a full evaluation. The accuracy can

be considerably improved, allowing larger time steps, by inserting an accep-

tance/rejection step after the di¬usion-drift step has been made (Reynolds

et al., 1982; Umrigar et al., 1993). The procedure without the reactive step

2

should lead to sampling of ψT , which can be made exact (for any time step)

by accepting/rejecting or reweighting moves such as to maintain detailed

2

balance under the known distribution ψT .

4.3.6 A practical example

Let us, for the sake of clarity, work out the programming steps for a realistic

case: the helium atom. We shall use atomic units (see page xvii) for which

the electron mass, elementary charge, and are all unity. The helium atom

consists of three particles: a nucleus with mass M (equal to 6544 for He-4),

charge +2 and coordinates R, and two electrons with mass 1, charge ’1 and

coordinates r 1 , r 2 . The Hamiltonian, with nine degrees of freedom, reads

12

ˆ

H=’ ∇R ’ 1 ∇2 ’ 1 ∇2 ’ ’

2 2 1

+ r12 , (4.34)

21 22 |r 1 ’R| |r 2 ’R|

2M

where

r12 = |r 1 ’ r 2 |. (4.35)

In the Born“Oppenheimer approximation we can eliminate R as a variable

by assuming M = ∞ and R ≡ 0, but there is no pressing need to do so.

Neither is there a need to separate center-of-mass and rotational coordinates

and reduce the number of internal degrees of freedom to the possible min-

imum of three, and rewrite the Hamiltonian. We can simply use all nine

coordinates; the e¬ect of reduced masses using internal coordinates is im-

plied through the ¬rst term concerning the di¬usion of the heavy nucleus.

The center-of-mass motion (corresponding to a free particle) will now also

4.3 The few-particle problem 93

be included, and will consist of a simple random walk, leading to an indef-

initely expanding Gaussian distribution. All relevant data will concern the

relative distribution of the particles.

The attractive Coulomb terms in the Hamiltonian cause troublesome be-

havior when electrons move close to the nucleus. The time step should be

taken very small in such con¬gurations and it becomes virtually impossi-

ble to gather su¬cient statistics. However, the use of importance sampling

with a trial function that is the product of the single-electron ground-state

solutions for the two electrons eliminates these attractive Coulombic terms

altogether. Take as trial function

ψT = exp[’±|r 1 ’ R| ’ ±|r 2 ’ R|]. (4.36)

Choosing

2M

±= , (4.37)

M +1

we ¬nd that

ˆ

HψT 4M 1

=’ + . (4.38)

ψT M + 1 r12

We note that this trial function is a very rough approximation to the real

wave function. For realistic applications it is necessary to use much bet-

ter trial functions, e.g., obtained from variational Monte Carlo or density

functional theory (see Section 4.7).

The time evolution of f according to (4.32) is solved by a collection of

walkers, each consisting of a nucleus and two electrons, that:

1

(i) di¬use with a di¬usion constant 1/2M , resp. 2;

(ii) drift with a drift velocity (see (4.33))

r1 ’ R r2 ’ R

±

u0 = + (4.39)

|r 1 ’ R| |r 2 ’ R|

M

for the nucleus, and

ri ’ R

ui = ’± (4.40)

|r i ’ R|

for the two electrons i = 1, 2;

(iii) create/annihilate according to

‚f 4M 1

’

= E+ f. (4.41)

‚„ M + 1 r12

94 Quantum chemistry: solving the time-independent Schr¨dinger equation

o

If the creation/annihilation step is implemented at each time by a stochastic

process according to (4.27), additional noise is introduced into the process.

It is better to assign a weight to each walker and readjust the weights ev-

ery step, while some form of population control is carried out at regular

intervals.7 The latter may involve duplication of heavy and annihilation of

light walkers (according to Grassberger, 2002: if the weight exceeds a given

upper threshold, then duplicate, giving each copy half the original weight;

if the weight is less than a given lower threshold, draw a random number

· between 0 and 1, annihilate if · < 1 , but keep with the double weight

2

if · ≥ 1 ), or a complete random reassignment of walkers chosen from the

2

weighted distribution of the existing walkers.

In the following Python program a number of functions are de¬ned to

realize the initial creation of walkers, the drift-di¬usion step with readjust-

ments of weights, and the population control. It is left to the reader to

employ these functions in a simple program that computes the ground state

energy of the helium-4 atom. There are two di¬erent ways to compute the

energy in excess of ’4M/(M + 1): ¬rst by monitoring the factor by which

the weights must be readjusted to keep their sum constant (E), and second

the average of 1/r12 over the ensemble of walkers (V ). When the time step

is small enough, both energies tend to be equal; their di¬erence is a good

indicator for the suitability of the time step. One may choose 1000 walkers,

a time step of 0.002 and make 1000 steps with weight adjustment before

the walkers are renewed. The excess energy above the value of ’3.999 455

hartree for ’4M/(M + 1) should be recovered by E or V ; this value should

equal +1.095 731, given the exact energy of the helium atom of ’2.903 724

(Anderson et al., 1993). With this simple approach one may reach this

value within 0.01 hartree, much better than the Hartree“Fock limit (see

page 101) which is 0.04 hartree too high due to lack of correlation (Clementi

and Roetti, 1974).