08 z=z*exp(-W)

09 s=sum(z)

10 E=-log(s)

11 y=z/s

12 return [y,E]

Comments

Line 03: the ¬rst two and last two points are zero and are kept zero.

Line 06: di¬usion step: each point becomes the average of its neighbors; zeros added at ends.

Line 07: second and before-last points set to zero.

Line 08: evolution due to potential.

Line 10: Energy (—”„ / ) needed to keep y normalized. This converges to the ground state energy.

Line 12: Use y = SRstep[0] as input for next step. Monitor E = step[1] every 100 steps for

convergence. Last y is ground state wave function.

Using the values for HF given in Table 1.1, using 1000 points (0, 3b) (b

is the bond length), and starting with a Gaussian centered around b with

σ = 0.01 nm, the energy value settles after several thousand iterations to

24.495 kJ/mol. The resulting wave function is plotted in Fig. 4.2. It very

closely resembles the Gaussian expected in the harmonic approximation.

The algorithm is only of second-order accuracy and the discretization error

is proportional to (”x)2 , amounting to some 3 J/mol for 1000 points.

4.3 The few-particle problem 89

The ¬rst excited state can be computed when the wave function is kept

orthogonal to the ground state. Denoting the ground state by ψ0 , we add a

proportion of ψ0 to ψ such that

ψψ0 dx = 0, (4.23)

which can be accomplished by adding a line between line 8 and 9 in the

program SRstep:

08a z=z-sum(z*y0)*y0

where y0 is the converged and normalized ground state wave function. Fig-

ure 4.2 includes the ¬rst excited state, with energy 71.866 kJ/mol. This

is accurate within 3 J/mol. The vibration wave numbers then are 4139.8

(harmonic oscillator) and 3959.9 (Morse function). Higher excited states

can be computed as well, as long as the wave function is kept orthogonal to

all lower state wave functions.

4.3.5 Di¬usional quantum Monte Carlo methods

If we look carefully at (4.12), we see that the ¬rst term on the r.h.s. is a

di¬usion term, as in Fick™s equation for the time dependence of the concen-

tration c of di¬using particles:

‚c

= D∇2 c. (4.24)

‚t

This ¬ts with what this term does after discretization (4.22): the function

splits in two halves located at the neighboring points. This is what would

happen with a probability distribution of particles after each particle has

made a random step over one grid element, either to the left or to the right.

This equivalence suggests a di¬erent way to solve the SE, which has been

pioneered by Anderson (1975, 1976).6 Suppose we wish to obtain the ground

state of a system of n particles. Consider an ensemble of a large number

of replicas of the system, each with its own con¬guration of the n particles.

The members of the ensemble are called psi-particles (psips) or walkers.

Then evolve each member of the ensemble as follows:

(i) Displace each of the coordinates of the particles in a “time” ”„ with

a random displacement with variance

(”x)2 = 2D”„ = ”„. (4.25)

m

6 See reviews by Anderson (1995) and by Foulkes et al. (2001), the latter with applications to

solids.

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

o

This can be done by sampling from a Gaussian distribution with that

variance, or by displacing the coordinate by ± ”„ /m.

(ii) Duplicate or annihilate the walker, according to the probability

P = exp[’(V ’ E)”„ / ] (4.26)

in order to satisfy the second term (source term) on the r.h.s. of

(4.12). This, again, can be done in several ways. One way is to let

the walker survive with a probability P if P < 1, and let it survive

but create a new walker with probability P ’ 1 if P ≥ 1. This is

accomplished (Foulkes et al., 1995) by de¬ning the new number of

walkers Mnew as

Mnew = integer(P + ·), (4.27)

where · is a uniform random number between 0 and 1. A higher-

order accuracy is obtained when the potential energy for the cre-

ation/annihilation (in (4.26)) is taken as the average before and after

the di¬usion step:

V = 1 [V („ ) + V („ + ”„ )]. (4.28)

2

Another way is to enhance or reduce a weight per walker, and then

applying a scheme like: duplication when the weight reaches 2 or

annihilation with 50% probability when the weight falls below 0.5;

the surviving and new walkers will then start with weight 1. Applying

weights only without creation/annihilation scheme does not work, as

this produces uneven distributions with dominating walkers.

(iii) The energy E determines the net gain or loss of the number of walkers

and should be adjusted to keep that number stationary.

The advantage of such a stochastic scheme above relaxation on a grid is

that it can more easily be expanded to several dimensions. For example, a

four-particle system in 3D space (a hydrogen molecule) involves 12 degrees of

freedom, reducible to six internal degrees of freedom by splitting-o¬ overall

translation and rotation. A grid with only 100 divisions per dimension

would already involve the impossible number of 1012 grid points, while an

ensemble of 105 to 106 walkers can do the job. The method remains exact,

and includes all correlations between particles. The result has a statistical

error that reduces with the inverse square root of the number of steps.

Di¬usional Monte Carlo methods cannot straightforwardly handle wave

functions with nodes, which have positive as well as negative regions. Nodes

occur not only for excited states, but also for ground states of systems

4.3 The few-particle problem 91

containing more than two identical fermions, or even with two fermions with

the same spin. Only two fermions with opposite spin can occupy the same

all-positive ground state wave function. The di¬using particles represent

either positive or negative wave functions, and when a positive walker crosses

a nodal surface, it would in principle cancel out with the negative wave

function on the other side. One scheme that avoids these di¬culties is to

keep the nodal surface ¬xed and annihilate the walkers that di¬use against

that surface; however, by ¬xing the nodal surface the method is no longer

exact. Schemes that use exact cancelation between positive and negative

walkers, instead of ¬xed nodal surfaces, have been successfully applied by

Anderson (1995).

The implementation of di¬usional Monte Carlo methods is much more

e¬cient when an approximate analytical wave function is used to guide

the random walkers into important parts of con¬guration space. This im-

portance sampling was introduced by Grimm and Storer (1971) and pur-

sued by Ceperley and Alder (1980). Instead of sampling ψ, the function

f (r, „ ) = ψ(r, „ )ψT (r) is sampled by the walkers, where ψT is a suitable

trial function. The trial function should be close to the real solution and

preferably have the same node structure as the exact solution; in that case

the function f is everywhere positive. Noting that the exact energy E is the

ˆ

eigenvalue of H for the solution ψ at in¬nite time

ˆ

Hψ = Eψ (4.29)

ˆ

and that H is Hermitian:

ˆ ˆ

ψ HψT dr = ψT Hψ dr = E ψT ψ dr, (4.30)

we can write the energy as the expectation of the local energy for the trial

wave function (which can be evaluated for every member of the generated

ensemble) over an ensemble with weight f = ψψT :

ˆ

f HψT dr

ˆ ˆ

ψ HψT dr HψT