python program 4.4 walkers

Three functions to be used for simple QMD of the helium atom

01 def initiate(N):

02 # create array of N walkers (helium atom)

03 # returns [walkers,weights]

04 walkers=zeros((N,3,3), dtype=float)

05 sigma=0.5

06 walkers[:,1,:]=sigma*randn(N,3)

07 walkers[:,2,:]=sigma*randn(N,3)

7 See Hetherington (1984), Sorella (1998) and Assaraf et al. (2000) for a discussion of noise and

bias related to stochastic recon¬guration.

4.3 The few-particle problem 95

08 r1=walkers[:,1]-walkers[:,0]; r2=walkers[:,2]-walkers[:,0]

09 r1sq=sum(r1**2,1); r2sq=sum(r2**2,1)

10 weights=exp(-2.*(sqrt(r1sq)+sqrt(r2sq))+0.5*(r1sq+r2sq)/sigma**2)

11 ratio=N/sum(weights)

12 weights=ratio*weights

13 return [walkers,weights]

14

15 def step(walkers,weights,delt,M):

16 # move walkers one time step delt; M=nuclear/electron mass

17 N=len(walkers)

18 r1=walkers[:,1]-walkers[:,0]; r2=walkers[:,2]-walkers[:,0]

19 r1norm=sqrt(sum(r1**2,1)); r2norm=sqrt(sum(r2**2,1))

20 for i in range(3):

21 r1[:,i]=r1[:,i]/r1norm; r2[:,i]=r2[:,i]/r2norm

22 alphadelt=2.*M/(M+1.)*delt

23 d1=-alphadelt*r1; d2=-alphadelt*r2

24 d0=-(d1+d2)/M

25 sd0=sqrt(delt/M); sd1=sqrt(delt)

26 walkers[:,0,:]=walkers[:,0,:]+d0+sd0*randn(N,3)

27 walkers[:,1,:]=walkers[:,1,:]+d1+sd1*randn(N,3)

28 walkers[:,2,:]=walkers[:,2,:]+d2+sd2*randn(N,3)

29 # adjust weights one time step

30 V=1./sqrt(sum((walkers[:,1]-walkers[:,2])**2,1))

31 weights=weights*exp(-V*delt)

32 ratio=N/sum(weights)

33 E=log(ratio)/delt

34 weights=ratio*weights

35 return [walkers,weights,E]

36

37 def renew(walkers,weights,Nnew):

38 # select Nnew new walkers with unit weight

39 wtacc=cumsum(weights)

40 s=wtacc[-1]

41 index=[]

42 for i in range(Nnew):

43 u=s*rand()

44 arg=argmax(where(greater(wtacc,u),0.,wtacc)) + 1

45 index=index+[arg]

46 wa=take(walkers,index)

47 wt=ones((Nnew))

48 return [wa,wt]

Comments

The coordinates of the walkers form an array walkers[n, i, j], where n numbers the walkers, i

numbers the particles (0 = nucleus, 1,2 = electrons) in each walker and j = 0, 1, 2 indicates the

x, y, z coordinates of a particle.

Function initiate creates a number of N walkers: lines 06 and 07 assign a normal distribution

to the electron coordinates (randn generates an array of normally distributed numbers); line 10

adjusts the weights to make the distribution exponential, and lines 11 and 12 normalize the total

weight to N .

The function step moves the particles in lines 26“28 by a simultaneous di¬usion and drift dis-

placement through sampling a normal distribution with prescribed mean and standard deviation.

Then the weights are adjusted according to the computed excess energy 1/r12 of each walker in

lines 30“31, and restored in lines 32“34 to the original total weight, which yields the “energy” E.

The function renew reconstructs a new set of walkers each with unit weight, representing the

same distribution as the old set. It ¬rst constructs the cumulative sum of the original weights

(in line 39); then, for each new walker, a random number, uniformly distributed over the total

weight, is selected, and the index of the original walker in whose range the random number falls,

is determined. All indices are appended in one list index (line 45) which is used in line 46 to copy

the walkers that correspond to the elements of that list.

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

o

4.3.7 Green™s function Monte Carlo methods

The di¬usional and drift step of a walker can be considered as random

samples from the Green™s function of the corresponding di¬erential equa-

tion. The Green™s function G(r, „ ; r 0 , 0) of a linear, homogeneous di¬eren-

tial equation is the solution ψ(r, „ ) when the boundary or initial condition

is given by a delta-function ψ(r, 0) = δ(r ’ r 0 ); the general solution is

the integral over the product of the Green™s function and the full boundary

function:

ψ(r, „ ) = G(r, „ ; r 0 , 0)ψ(r 0 , 0) dr 0 . (4.42)

For the di¬usion equation the Green™s function is a Gaussian.

There is an alternative, iterative way to solve the time-independent Schr¨-

o

dinger equation by Monte Carlo moves, by iterating ψn according to

2m 2m

’∇2 ’ E ψn+1 = V ψn . (4.43)

2 2

The function ψn is sampled, again, by walkers, who make a step sampled

from the Greens function of the di¬erential operator on the left-hand-side,

which in this case involves a modi¬ed Bessel function of the second kind. The

iterations converge to the “exact” solution. We shall not further pursue these

Green™s function Monte Carlo methods, which were originally described by

Kalos (1962), and refer the reader to the literature.8

4.3.8 Some applications

Quantum Monte Carlo methods have been used to solve several few-particle

and some many-particle problems.9 Particular attention has been paid to the

full potential energy surface of H3 in the Born“Oppenheimer approximation:

a nine-dimensional problem (Wu et al., 1999). In such calculations involving

a huge number of nuclear con¬gurations, one can take advantage of the fact

that the ¬nal distribution of walkers for one nuclear con¬guration is an

e¬cient starting distribution for a di¬erent but nearby con¬guration. All

electron correlation is accounted for, and an accuracy of better than 50

J/mol is obtained. This accuracy is of the same order as the relativistic

correction to the energy, as calculated for H2 (Wolniewicz, 1993). However,

the adiabatic error due to the Born“Oppenheimer approximation is of the

order of 1 kJ/mol and thus not all negligible.

8 See also Anderson (1995); three earlier papers, Ceperley and Kalos (1979), Schmidt and Kalos

(1984) and Schmidt and Ceperley (1992), give a comprehensive review of quantum Monte Carlo

methods. Schmidt (1987) gives a tutorial of the Green™s function Monte Carlo method.

9 See Anderson (1976) and the references quoted in Anderson (1995).

4.4 The Born“Oppenheimer approximation 97

There is room for future application of QMC methods for large systems

(Foulkes et al., 2001; Grossman and Mitas, 2005). Systems with up to a

thousand electrons can already be treated. Trial wave functions can be

obtained from density functional calculations (see below) and the QMC

computation can be carried out on the ¬‚y to provide a force ¬eld for nuclear

dynamics (Grossman and Mitas, 2005). Because QMC is in principle exact,

it provides an ab initio approach to integrated dynamics involving nuclei and

electrons, similar to but more exact than the “ab initio” molecular dynamics

of Car and Parrinello (see Section 6.3.1). But even QMC is not exact as it

is limited by the Born“Oppenheimer approximation, and because the nodal

structure of the trial wave function is imposed on the wave function. The

latter inconsistency may result in an incomplete recovery of the electron

correlation energy, estimated by Foulkes et al. (2001) as some 5% of the total

correlation energy. For elements beyond the second row of the periodic table,

the replacement of core electrons by pseudopotentials becomes desirable for

reasons of e¬ciency, which introduces further inaccuracies.

Finally, we note that a QMC program package “Zori” has been developed

by a Berkeley group of scientists, which is available in the public domain

(Aspuru-Guzik et al., 2005).10

4.4 The Born“Oppenheimer approximation

The Born“Oppenheimer (B“O) approximation is an expansion of the be-

havior of a system of nuclei and electrons in powers of a quantity equal

to the electron mass m divided by the (average) nuclear mass M . Born