A Variational Monte Carlo Approach to Atomic Structure

Research: Science and Education edited by

Advanced Chemistry Classroom and Laboratory

Joseph J. BelBruno Dartmouth College Hanover, NH 03755

A Variational Monte Carlo Approach to Atomic Structure

W

Stephen L. Davis Department of Chemistry and Biochemistry, George Mason University, Fairfax, VA 22030; [email protected]

Chemistry undergraduates are usually introduced to quantum mechanics as part of a physical chemistry course taken during the junior or senior year. The treatment in such a course typically begins with the exact solution of the Schrödinger equation for a few simple, model systems, such as the infinite square well. Wave functions and energy levels are then presented, with or without proof, for a few more exactly solvable systems, such as the simple harmonic oscillator and the hydrogenic atom. The more interesting topics in multi-electron or multi-center systems often receive only qualitative treatment. Early on, students learn how to perform a few simple calculations with wave functions, such as expectation values of position in the simple model systems. However there is a large gap between this sort of calculation and a quantitative treatment of atoms with more than one electron. This of course is due to the difficulty of evaluating the multi-dimensional integrals that appear in the expectation values required for a variational treatment. Even if analytic formulas for these integrals are made available, they are too daunting and too ad hoc to be very satisfying. “Black-box”, ab initio, electronic structure software can be unsatisfying as well. One possible way out of this dilemma is to use a relatively simple variational Monte Carlo method (1) (VMC) to compute expectation values of the Hamiltonian for atoms containing a few electrons. The literature reveals much interest in quantum Monte Carlo methods, including VMC (2–5), due to its ability to exploit functional forms for the trial wave functions that explicitly include dependence on the inter-electronic coordinates, thus yielding a more efficient treatment of electron correlation. This paper’s interest in VMC, however, is pedagogical. VMC will be used here not to obtain highly accurate, correlated energies via sophisticated trial functions, but to permit undergraduates to perform quantitative calculations on atoms with more than one electron, while retaining some understanding of how the numerical results emerge. As described below, such calculations can be done without the need for computer programming; a current spreadsheet program, such as Microsoft Excel, suffices. The variational Monte Carlo method depends upon generating a random sample of electron positions distributed according to a given probability function involving ψ*ψ, where ψ is some trial function. Applying the method to more than a small number of particles requires very efficient generation of these random samples. In the literature this is done using sophisticated algorithms such as Metropolis sampling (6). However, since the topic of this paper is a method suitable for undergraduate, physical chemistry students, we use instead a simpler and much less efficient, “brute-force” approach to generating the nonuniform, random position samples. This approach is suitable for the few-electron systems www.JCE.DivCHED.org

•

treated here, but of course would be inadequate for larger systems or for higher statistical accuracy. An earlier paper in this Journal (7) described a diffusion Monte Carlo approach. Although preferable in many ways, the diffusion Monte Carlo approach may not be as easily understood by the student as the VMC approach used here. The sections below first describe the VMC method and then discuss the method used for generating nonuniform random numbers. This is followed by a series of computational examples on the electronic energies of the first four atoms. Variational Monte Carlo Method In the variational method, the energy of the atomic state is calculated as the expectation value of the Hamiltonian in some trial function. For an atomic system this may be written as E

=

∑ Tˆ i

i

+

∑ Vi

+

i

∑ Vi j

i< j

(1)

where the sums run over all N electrons in the system, and where Tˆi, Vi, and Vij represent, respectively, the operator for the kinetic energy of electron i, the potential energy of attraction between the nucleus and electron i, and the energy of repulsion between electrons i and j. In spherical coordinates, the kinetic energy operator can be written

1 d2 2 d Tˆi = − + ri dri 2 dri 2

+

ᐉ(ᐉ + 1) 2ri 2

(2)

where l is the angular momentum quantum number and ri is the radial position of the electron. Equation 2 assumes that Tˆ operates on an eigenfunction of the angular momentum with eigenvalue l. This is the case for all of the trial functions for the atomic systems treated in this paper. Note that this expression, as well as all others in this paper, uses so-called atomic units (8a). This is the system in which and the electron charge (e) and mass (me ) take values of unity. The atomic unit of energy is the hartree (Eh) and the unit of length is the bohr (a0). In these units the potential energy Vi is given by

Vi =

−Z ri

(3)

where Z is the number of protons, and the electron–electron repulsion potential by Vij =

1 rij

(4)

where rij is the distance between electrons i and j.

Vol. 84 No. 4 April 2007

•

Journal of Chemical Education

711

Research: Science and Education

Evaluation of these expectation values normally involves integration over the coordinates of either one electron (for ˆ i and Vi ) or two electrons (for Vij ). Though these integrals T can be done analytically for simple trial functions such as those used here, the complexity of the resulting expressions tends to make the student lose sight of the physical significance behind the mathematics. This paper makes use of a Monte Carlo method for evaluating the necessary expectation values. This avoids much of the mathematical complexity and also helps the student see more clearly the probabilistic meaning of the wave function. Using VMC, the calculation of the energy proceeds according to the following steps: 1. Choose a trial function Ψ(r1, r2, ... , rN) of the N position vectors ri of the electrons. In this paper the trial function will be taken as a product (or antisymmetrized product) of one-electron hydrogenic orbitals, but any trial function is possible. The trial function will usually involve one or more parameters that can be varied to minimize the resulting energy expectation value. 2. Choose a set of random position vectors for the electrons r1, r2, ... , rN that follow the probability distribution |Ψ(r1, r2, ... , rN)|2. There are standard methods for choosing random values distributed according to a given distribution function, which are described in the next section. ˆ and V 3. Evaluate all the “local energy” contributions T i

i

when electron i is at random position ri and Vij when electrons i and j are at random positions ri and rj. The ˆ Ψ)/Ψ, evaluated local kinetic energy is defined by (T i at ri. The local potential energies are simply Vi (ri ) and Vij (rij ). 4. Repeat steps 2 and 3 a sufficiently large number of times and average the resulting local energy contribuˆ , V , and V . These averages are approximations tions T i i ij to the corresponding expectation values and may be added together according to eq 1 to give E . These Monte Carlo estimates of E approach the analytical value as the number of random positions in the sample becomes larger and larger. However for any finite sample size, there is statistical error in E . This is of course in addition to any error arising from a nonexact trial function.

According to the variation theorem, the expectation value E is an upper bound to the true energy. Therefore it is possible to improve the energy by varying any parameters in the trial function so as to minimize the resulting E . Nonuniform Random Numbers Many software applications and almost all programming languages allow one to generate random (or pseudorandom) numbers uniformly distributed on an interval. For instance, Microsoft Excel includes the worksheet function RAND(), which returns an evenly distributed random number greater than or equal to zero and less than one. Similarly the Visual Basic language used for Excel macros includes the function RND.

712


•

Applying the VMC method to atomic structure problems requires generating a sequence of random values for the electron positions that follow a nonuniform probability distribution. In general the probability distribution of each electron is dependent on the distribution of all the others, rather than following an “orbital approximation” in which each electron is described by its own wave function. Though certainly not required by the Monte Carlo method, all the examples in this paper will use an orbital approximation for simplicity. Therefore the spherical coordinates of electron i are assumed distributed according to the function 2

P (ri , θi , φi ) = ψi (ri ) Yᐉm (θ , φ ) ri2 sin θi 2

= ψi (ri ) Θᐉm (cos θi ) ri2 sinθi

(5)

where ψ(r) is the radial part of the wave function, Ylm is a spherical harmonic, and Θlm an associated Legendre polynomial. Here we have assumed that the one-electron trial functions are eigenfunctions of the angular momentum operators L2 and Lz, as appropriate for atomic systems. Since the resulting probability distribution is independent of the angle φ, this coordinate is actually uniformly distributed on the interval 0–2π. Thus the RAND() function in Excel is sufficient for generating this coordinate. However, generating random values of the nonuniform coordinates r and θ requires another approach. There are two standard methods for generating nonuniform random numbers: the transformation method and the rejection method (9). Both are briefly described in this section. The transformation method uses a mapping of uniform random numbers p onto the cumulative probability distribution. For instance, suppose that Θlm(cos θ) is the wave function describing the θ-dependence of a certain electron and that p is a uniform random number in the range 0 to 1. Then θ

Θᐉm2 (cos θ ′ ) sin θ ′ dθ ′

p =

(6)

0

relates the uniform random number p to a corresponding nonuniform random θ. Inverting this relation gives properly distributed values of θ in terms of uniform random numbers p. As an example, for an s-electron with l m 0 and Θ00 (1/2) this gives, after integrating eq 6 and solving for θ, θ = cos −1 (1 − 2 p )

(7)

Generating an ordinary, uniform random number p thus yields a corresponding nonuniform random number θ. Obviously the transformation method is most useful when the probability distribution can be integrated analytically and the result inverted, as is the case for the angular functions used in this paper. An alternate method, applicable to any distribution, is the so-called rejection method. Suppose that we wish to generate random values of the radial coordinate r, distributed according to the trial function φ(r). This may be done as follows:


•

www.JCE.DivCHED.org


Hydrogen Atom Ground State

1. Choose a random value r0, uniformly distributed over the range 0 ≤ r0 ≤ rmax. This may be done using an ordinary random number generator.

Although the exact energy levels and wave functions of the hydrogen atom are well known, it is useful to illustrate the VMC method by calculating the ground state energy of hydrogen using the known 1s function as a “trial” function. Since the hydrogen atom can have no electron–electron reˆ It is easily pulsion, its energy is given by E V T. shown from eq 2 that the local kinetic energy is

2. Choose a random value p, uniformly distributed over the range 0 ≤ p ≤ Pmax, where Pmax is the maximum value of the probability function P (r) | φ(r)|2r 2. 3. If p ≤ P (r0 ) accept the value r0. Otherwise reject it. Successive values of r0 chosen according to this “acceptance criterion” will produce a sample distributed according to P (r).

One complication in applying this method to r, representing the distance between the electron and the atomic nucleus, is the choice of rmax. The domain of the variable r of course is actually 0 ≤ r ≤ ∞. Since it is impossible to choose random numbers uniformly distributed over an infinite range, one must treat the domain as 0 ≤ r ≤ rmax, where rmax is a value large enough that the energy expectation value E is essentially unaffected by the finite upper limit. One advantage of the rejection method is that the trial functions need not be normalized. This is because the normalization constant scales both P(r) and Pmax in the same way, leaving the acceptance criterion in step 3 unaffected. Note also that it is not necessary to determine the value of Pmax analytically. It is often simpler to equate Pmax to the largest value of |ψ(r)|2r 2 on the finite domain of random r-values actually generated in step 1 of the rejection method. As an illustration of the rejection method, suppose we wish to generate random values of r for the 1s electron in hydrogen, assuming the (un-normalized) trial function ψ1s(r ) = e − ς r

(8)

This gives the radial probability distribution

P1s (r ) = ψ1s (r )

2

= r 2 e −2 ς r

(9)

ζ is a parameter that can be used to minimize E . It is easily shown that Pmax 1/(e2 ζ2). Choosing rmax 5/ζ, guarantees that P(r) has fallen to less than 1% of its maximum value by the rmax cutoff. Now we generate two random numbers r0 and p0, which are uniformly distributed within 0 ≤ r0 ≤ rmax and 0 ≤ p0 ≤ Pmax respectively. The distance r0 is then accepted as part of the desired random sample if it satisfies the criterion P1s (r0 ) = r02 e −2 ς r0 ≥ p 0

•

ς ς2 − r 2

(11)

ˆ and V , which involve only The expectation values T the radial coordinate, are now obtained, not by explicit integration, but by Monte Carlo averaging over a random sample of r-values generated by the rejection method. For each value of r in the sample, the local kinetic (eq 11) and potential (eq 3) energies are evaluated and then the overall average is taken. This process is easily implemented within a simple, sixcolumn spreadsheet, illustrated here with the Microsoft Excel formulas in Table 1. The first four columns generate the required random value of r, using the rejection method. Column A generates a uniform random value r0. Column B then generates the corresponding values of P(r) calculated from eq 9. Column C generates a uniform random value p for use in the acceptance criterion of eq 10. Column D will produce the numerical values of accepted r0, with the rejected values appearing as blanks. Column E evaluates the local kinetic energy in eq 11 for the specific random r0 in column A. Similarly, column F evaluates the local potential energy using eq 3. The symbol ZETA in column E represents the parameter z in the trial function, while the symbol Z in the column F represents the value of the actual nuclear charge Z. ZETA and Z should be replaced with the desired numerical value or with the address of a cell containing the value. After replicating these formulas through as many rows as desired, the expectation values may be calculated using the Excel AVERAGE function. T would be the average of all the values in column E, while V would be the average of all the values in column F. Note that because the Excel AVERAGE function simply Table 1. Excel Formulas for Calculating Expectation Valuesa T and V for the H Ground Stateb Column

(10)

Otherwise it is rejected. Repeated choices of random pairs (r0, p0) followed by application of the acceptance criterion in eq 10 will generate a set of r-values that approaches the distribution P(r) as the sample size becomes larger and larger. In the examples described in the rest of this paper, the rejection method is always used to generate random values of radial positions r. The transformation method is usually used to generate random values of the spherical coordinate θ. The angle φ, since it is uniformly distributed, can be easily generated from the Excel formula 2*PI()*RAND().

www.JCE.DivCHED.org

Tˆ ψ1s (r ) = ψ1s (r )

Formula 5*RAND( )/ZETA $A1^2*EXP(2*ZETA*$A1) RAND( )/(EXP(2)*ZETA^2) IF(C1>$B1,”,$A1) IF($C1>$B1,”,ZETA/$A1ZETA^2/2) IF($C1>$B1,”,Z/$A1)

A B C D E F

Values for T and V are given in columns E and F, respectively. 1s hydrogenic trial function was used. The symbol ZETA represents the parameter from eq 8; the symbol Z represents the actual atomic number. a

bA


•


713


ignores the blank rows, only the accepted values of the position contribute. Table 2 shows the data resulting from 10 recalculations of a 10,000-row version of the spreadsheet specified by Table 1, with ζ 1, Z 1, and rmax 5 a0. Each recalculation gives roughly a 37% acceptance rate, yielding a total sample size of about 37,000 for the aggregate ten recalculations. The resulting values of V and T are within about 0.5% of the known analytical values. Furthermore the statistical errors in V and T cancel out, giving the exact value E 0.500 Eh. This is because the exact hydrogenic wave function was used, with ζ Z, and because the r 1 expectation values in T and V cancel out exactly. This will not be true in general. The optimal value of ζ could also be determined “experimentally” by minimizing E . The relatively low acceptance rate, which gets worse as the number of electrons increases, could be improved by storing only accepted position values in the rows of the spreadsheet. This could be done, for instance, using Excel macros. The Supplemental MaterialW includes a version of the spreadsheet implementing such a macro. Because most students are unlikely to be proficient in writing Excel macros, the examples below do not employ macro methods. Instead sufficient accuracy is obtained by recalculating a spreadsheet multiple times and averaging the results, as in Table 2. In this way, with a little patience, students can generate sample sizes as large as needed.

the 2s state ψ2s(r ) = (c − ς r ) e

ςr 2

(12)

This has the same form as the exact hydrogenic wave function, but replaces the nuclear charge Z with a variable parameter, ζ. The constant c is determined by the requirement that the wave function be orthogonal to the 1s function in eq 8. For the exact hydrogenic wave function, ζ 1 and c 2. However, here we use the more general form so that the parameter ζ can be determined from the variational theorem. This flexibility will be especially useful for later calculations on atoms that contain both 1s and 2s electrons. In such a situation, it is useful to allow the ζ parameters to take different values for the 1s and 2s functions, corresponding to different effective nuclear charges. If we designate the 1s parameter in eq 8 by ζ1 and the 2s parameter in eq 12 by ζ2, it can be shown that orthogonality of the two functions requires c =

6ς 2 2ς1 + ς 2

(13)

Applying the kinetic energy operator to ψ2s (l 0) gives a local kinetic energy of ς 23 r − (8 + c ) ς 22 + Tˆ ψ2 s = 8 (c − ς 2 r ) ψ 2s

Hydrogen Excited States Approximate energy levels for the excited states of hydrogen can also be obtained using VMC, if one ensures that the trial functions are orthogonal to lower states of the same symmetry. For instance, one might use as a trial function for

−

4( 2 + c ) r

(14)

Monte Carlo computation of E now requires averaging the local kinetic energy of eq 14 and the local potential energy of eq 3 over a random sample of r-values distributed according

Table 2. Recalculated Values of the Ground State Energy for Ha

Table 3. Recalculations of the 2s Energy for H

Trial

Number of Hits

T / E h

V / E h

Trial

Number of Hits

T / E h

V / E h

1

3592

0.5077

1.0077

1

4366

0.1331

0.2582

2

3675

0.5044

1.0044

2

4360

0.1174

0.2424

3

3681

0.5406

1.0406

3

4247

0.1339

0.2590

4

3706

0.4982

0.9982

4

4269

0.1144

0.2393

5

3676

0.4917

0.9917

5

4343

0.1219

0.2469

6

3732

0.4677

0.9677

6

4336

0.1218

0.2468

7

3684

0.4890

0.9890

7

4298

0.1248

0.2498

8

3718

0.5110

1.0110

8

4333

0.1247

0.2497

9

3700

0.5265

1.0265

9

4405

0.1260

0.2510

10

3630

0.4888

0.9888

10

4394

0.1260

0.2511

Average

3679

0.5026

1.0026

Average

4335

0.1244

0.2494

Std. Dev.

42

0.0207

0.0207

Std. Dev.

51

0.0061

0.0061

0.0026

0.0026

Error

0.0006

0.0006

0.51

0.26

Error Error (%)

Using the 1s trial function with 1. All energies are in atomic units (hartrees). a

714


•

Error (%)

0.50

0.20

Using the 2s hydrogenic trial function with 1. All energies are in atomic units (hartrees). a


•

www.JCE.DivCHED.org


to |ψ2s|2r 2. This random sample is generated using the rejection method, just as described above for the hydrogen ground state. Similarly, using a hydrogen 2p trial function −

ψ 2 p (r ) = r e

ςr 2

(15)

the local kinetic energy is given by

E

Tˆ ψ2 p ς2 ς = − + ψ2 p 8 r

(16)

Note that eq 16 includes the centrifugal contribution (l 1) from eq 2, which has the effect of canceling out the r2 term. Averaging the local kinetic energy and the potential energy over the sample of r-values distributed according to |ψ2p|2r 2 gives the Monte Carlo estimate of the 2p energy. Table 3 shows the averaged results from 10 recalculations of a 10,000-row spreadsheet using the 2s trial function, with ζ taken as its optimal value of unity. The value of rmax was somewhat arbitrarily chosen as 12/ζ, while Pmax was taken as the actual maximum value of |ψ2s(r)|2r 2 on the set of rvalues. Since the acceptance rate is about 43% for this value of rmax, this represents a random sample of about 43,000 rvalues. The resulting E 0.125 Eh agrees with the exact 2s energy, again because the exact trial function was used. It would be easy to verify that ζ 1 is indeed the optimal value, by using this spreadsheet to minimize E with respect to ζ. Similar results are obtained using the 2p trial function in eq 15. However, for the 2p function one can no longer use eq 7, which assumed an isotropic orbital. Instead, the probability distribution P(θ) in eq 6 becomes, using Θ10(cos θ) in place of Θ00(cos θ),

P (θ ) =

3 2

with ψ1s given by eq 8 and where r1 and r2 are radial coordinates of the two electrons and ζ is a variational parameter corresponding to the effective nuclear charge experienced by the electrons. The one-electron expectation values T and V are computed just as described previously for the 1s state of hydrogen, but with Z 2. The new feature in helium is the electron–electron repulsion (eq 4), so that E becomes

cos 2 θ sin θ

(20)

It can be shown that r12 =

r12 + r22 − 2r1 r2 × cos θ1 cos θ2 − sin θ1 sin θ2 cos (φ1 − φ2 )

(21)

where (r1, θ1, φ1) are the spherical coordinates of electron 1, and similarly for electron 2. Since r12 depends on the angular positions as well as the radial distances, the random electron positions must now include specification of all three spherical coordinates for each electron. Thus r121 must be calculated as follows: 1. Use the rejection method to choose random values of r1 and r2 that satisfy the radial probability distribution in eq 9. 2. Choose random values of θ1 and θ2 using the transformation method in eq 7. 3. Choose uniform random values of φ1 and φ2. 4. Evaluate the local electron–electron repulsion energy, r121, for these coordinates, using eq 21. 5. Repeat steps 1–4 for a large sample of random coordinates and average to obtain V12

This could be implemented in a spreadsheet with formulas as shown in Table 4. Columns A and D generate random

(17)

Inverting eq 6 then gives, in place of eq 7

Table 4. Excel Formulas for Calculating Expectation Valuesa T , V , and V 1 2 for the He 1s Ground Stateb

1

θ = cos −1 (1 − 2 p ) 3

(18)

Helium Ground State Energy

Formula

Column A B C D E F G H I J K

where p is a uniform random number 0 ≤ p ≤ 1. Alternatively, a properly distributed sample of θ values for the pelectron could be chosen using the rejection method. This would have the advantage of greater conceptual simplicity for the student. The disadvantage is that the percentage of accepted electron positions would be much lower, with correspondingly greater statistical errors in E . One complication in applying VMC as described above is the need to derive the local kinetic energy anew for each different trial function, as in eq 14, for instance. This could be avoided by computing the required derivatives numerically within the spreadsheet, allowing the use of the same Excel formulas regardless of the trial function.

5*RAND( )/ZETA ACOS(12*RAND( )) 2*PI( )*RAND( ) 5*RAND( )/ZETA ACOS(12*RAND( )) 2*PI( )*RAND( ) $A1^2*EXP(2*ZETA*$A1)*$G1^2*EXP(2*ZETA*$G1) PMAX*RAND( ) IF($H1>$G1,”,ZETA/$A1ZETA^2/2) IF($H1>$G1,”,Z/$A1) IF($H1>$G1,”,1/SQRT($A1^2$G1^22* A1*G1*(COS($E1)*COS($K1)SIN($E1)* SIN($K1)*COS($F1$L1))))

Values for T , V , and V1 2 are given in columns I, J, and K, respectively. b An orbital trial function (eq 19) was used. The symbol ZETA represents the parameter for the 1s radial function (eq 8); the symbol Z represents the actual atomic number. The symbol PMAX represents the actual maximum of all of the values in column G. a

The VMC method can also be used to obtain an energy estimate for the helium 1s2 ground state, taking as a trial function the un-normalized orbital product ψ (r1, r2 ) = ψ1s (r1 ) ψ1s (r2 ) www.JCE.DivCHED.org

= 2 T1s + 2 V1s + V1s :1s

(19) •


•


715


ates the local electron–electron repulsion energy for the accepted coordinates, using eq 21. The expectation values T , V , and V12 are obtained by averaging columns I, J, and K over the entire random sample of the six coordinates. Minimization of E then gives the optimal value of the effective nuclear charge parameter ζ. Calculations were performed using a 10,000-row spreadsheet constructed according to Table 4. The value of rmax was taken to be 5/ζ, just as for the hydrogen ground state. The acceptance rate is only 13%, much lower than for the oneelectron case. This is because acceptance depends on two random variables, r1 and r2, each meeting their criterion at the same time for the two-electron averages. Figure 1 shows the resulting E values as a function of ζ. The points were fitted to a quadratic in ζ. Due to the rather small sample of ∼1300 accepted points, the statistical errors are quite pronounced. Nevertheless, the quadratic fit gives an optimized ζ of 1.70, in rather good agreement with the analytical value of 27/16. Table 5 shows the results from 10 recalculations of the 10,000-row spreadsheet at the fixed, optimal value of ζ . Thus the final averages in Table 5 are based on a total sample of about 13,000 random positions of the two electrons. The resulting ground state energy E 2.834 Eh is within 0.5% of the exact, analytical value (8b). For the repulsion between two spherically-symmetric charge distributions, as in the 1s2 case treated here, the usual computational approach employs an expansion in spherical harmonics to replace r121 inside the expectation value with r>1, where r> is defined as the greater of r1 and r2 (8c). This method can be implemented by modifying Column K in Table 4. In addition, the random angles generated in Columns B,

Figure 1. E versus ζ for the helium 1s2 ground state, from ∼1300 random r1, r2 pairs distributed according to the trial function of eq 19. All values are given in atomic units.

values of r1 and r2, respectively. Columns B and E generate properly distributed values of θ1 and θ2. Columns C and F generate uniformly distributed values of φ1 and φ2. Column G contains the value of |ψ(r1) |2 |ψ(r2)|2 r12 r22, for use in the rejection test. Column H is the uniform random number p0 used to accept the pair r1, r2 according to the criterion that p0 < |ψ(r1) |2 |ψ(r2)|2 r12 r22. Note that the angular variables do not enter into the acceptance test at all, since properly distributed values of the angles are already guaranteed by use of the transformation method. Columns I and J evaluate the local kinetic and potential energies for the accepted positions, exactly as in the case of a one-electron atom. Column K evalu-

Table 5. Recalculated Values of the Ground State Energy Contributions for Hea Trial

Number of Hits

T / E h

V / E h

V 1 2 / E h

E / E h

1

1359

1.5439

3.5174

1.0798

2.8672

2

1335

1.2704

3.1931

1.0781

2.7673

3

1339

1.2668

3.1888

1.0475

2.7965

4

1317

1.5080

3.4748

1.0730

2.8606

5

1314

1.4021

3.3493

1.0494

2.8449

6

1407

1.3213

3.2535

1.0661

2.7982

7

1393

1.3353

3.2701

1.0307

2.8387

8

1349

1.3346

3.2692

1.0566

2.8126

9

1313

1.5875

3.5691

1.0592

2.9040

10

1305

1.4673

3.4265

1.0663

2.8522

Average

1343

1.4037

3.3512

1.0607

2.8342

Std. Dev.

35

0.116

0.1375

0.0153

0.0405

0.0201

0.0238

0.0060

0.0135

0.70

0.60

0.50

Error Error (%) a

1.40

Using the 2s orbital product trial function in eq 19, with 1.68. All energies are given in atomic units (hartrees).

716


•


•

www.JCE.DivCHED.org


C, E, and F are no longer needed and can be omitted. This approach would be preferred in actual calculations due to its greater efficiency. However there are some pedagogical advantages to the simpler approach used in Table 4. Helium Excited States The VMC method correctly describes the shielding between the two 1s electrons in the helium ground state. The same approach can be used to reveal the difference in shielding between 2s and 2p. This section describes simple VMC calculations of the energies and optimal screening parameters for the 1s2s and 1s2p excited states of helium. The antisymmetry requirement on the wave functions is ignored in this section but will be treated in the next one. Here we use a simplistic orbital product for these excited states and their corresponding radial probability distributions

Figure 2. E versus ζ for the 1s2s and 1s2p excited states of helium, from ∼60,000 random r1, r2 pairs. All quantities are given in hartrees.

2

P1s 2s (r1, r2 ) = ψ1s (r1 ) ψ2 s (r2 ) r12 r22 P1s2p (r1, r2 ) = ψ1s (r1 ) ψ 2 p (r2 )

2 2 2 r1 r2

(22)

where the one-electron orbitals are defined in eqs 8, 12, and 15. Eq 20 is replaced by E

= T1s + T2 s + V1s + V2s + V1s : 2 s

(23)

and similarly for 1s2p. Monte Carlo calculations for 1s2s and 1s2p were implemented by adapting the spreadsheet used for the 1s2 state. Columns were added for T2 (eqs 14 or 16) and V2 (eq 3). In addition, random values of θ2 now have to be generated using either eq 7 or eq 18. Different values of rmax are used for the 1s and either or both of the 2s and 2p electrons, 5/ζ for the former and 12/ζ for the latter. Note also that distinct values of the orbital parameters are used (ζ1 for the 1s and ζ2 for either or both of the 2s and 2p electrons), allowing for different effective nuclear charges. As in the case of the helium ground state, the acceptance rate is only about 15%. The Monte Carlo results were based on 20 averaged recalculations of a 20,000-row spreadsheet, giving a total sample size of about 60,000 random r1, r2 pairs. Energy expectation values were obtained in this way with ζ1 fixed at 2.00 (for the 1s electron) and ζ2 varied between 1.00 and 1.50 (for either or both of the 2s and 2p electrons). Figure 2 shows the resulting E versus ζ2 for both the 1s2s and 1s2p states. The curves show a quadratic fit to the actual Monte Carlo energies. Although some statistical deviation from the smooth curve is evident even for the 60,000 sample size, the scatter is noticeably less than in Figure 1, which used a much smaller sample size. Figure 2 convincingly demonstrates the difference in 2s and 2p energies caused by the presence of the inner 1s electron. The minimized energies occur at ζ2 1.25 for the 2s electron and at ζ2 1.0 for the 2p electron, corresponding to E values of 2.162 Eh for the 1s2s state and 2.126 Eh for the 1s2p state. For comparison, the analytical expressions show a minimum energy of E 2.155 Eh for 1s2s at ζ2s 1.20 and E 2.128

www.JCE.DivCHED.org

•

Eh for 1s2p at ζ2s 1.05, in fair agreement with the VMC results. In the absence of the 1s electron, the exact ζ values would be 2.0 for both the 2s and 2p states, corresponding to the actual value of Z. The fact that optimized ζ values for the outer electron are considerably less than 2.0 in both the 1s2s and 1s2p states shows the effect of the shielding by the inner 1s electron. The fact that optimization yields ζ2s > ζ2p shows that the 2s electron is shielded less effectively than the 2p electron. The higher effective nuclear charge in the 2s case is equivalent to a lower energy. Thus, the simple VMC spreadsheet calculation gives a quantitative illustration of the usual argument for the dependence of electron energies on the orbital angular momentum quantum number. Helium Singlet and Triplet States According to the Pauli principle, wave functions must be antisymmetric with respect to exchange of identical particles. In the previous section this antisymmetry requirement was ignored; the goal was simply to understand the difference in shielding between 2s and 2p electrons. However, the VMC approach gives a simple way to calculate singlet–triplet energy differences. It is only necessary to choose properly antisymmetrized trial functions in place of the orbital products used in the previous section. It is well known that configurations such as helium 1s2s give rise to triplet and singlet terms that differ in energy. In the 1s2s 3S term, the spatial part of the wave function is antisymmetric, while the spin part is symmetric. In the 1S term this is reversed, with the spatial part being symmetric. Therefore we take as trial functions for the 1s2s 3S and 1S terms

ψ triplet = ψ1s (r1 ) ψ 2 s(r2 ) − ψ 2s (r1) ψ1s(r2 ) ψ singlet = ψ1s (r1 ) ψ2s(r2 ) + ψ 2s (r1 ) ψ1s(r2 )

(25)

The one-electron contributions T and V are unchanged by the antisymmetry requirement; their calculation


•


717

Research: Science and Education Table 6. Excel Formulas for Calculating the Expectation Valuea V 1 2 for the He 1s2s 3S Excited Termb Formula

Column

12*RAND( )/ZETA1 ACOS(12*RAND( )) 2*PI( )*RAND( ) 12*RAND( )/ZETA2 ACOS(12*RAND( )) 2*PI( )*RAND( ) (EXP(ZETA1*$A1)*EXP(ZETA2*$G1/2)* (CZETA2*$G1)*EXP(ZETA1*$G1)*EXP(ZETA2* $A1/2)*(CZETA2*$A1)^2*$A1^2*$G1^2 PMAX*RAND( ) IF($H1>$G1,”,1/SQRT($A1^2$G1^22* $A1*$G1*(COS($E1)*COS($K1)SIN($E1)* SIN($K1)*COS($F1$L1))))

A B C D E F G

H I

a The value for V1 2 is given in column I. V1 2 is obtained from eq 21 and properly distributed random angle variables. bAn antisymmetrized trial function was used. The symbols ZETA1 and ZETA2 are, respectively, the values of for the 1s and 2s functions. The symbol C is the value given by eq 13. The symbol PMAX represents the maximum value of the probability distribution in column G. Note that the corresponding 1s2s 1S energy may be calculated simply by changing the sign between the orbital products in column G.

proceeds just as before. Table 6 shows the spreadsheet formulas needed to calculate V12 from the antisymmetric (triplet) trial function. The corresponding singlet trial function is handled similarly. Note that the larger value rmax 12/ζ2 must now used for both r1 and r2, since either electron may occupy the 2s orbital. The acceptance rate in the rejection method drops to about 7%, due to the use of the larger rmax for the 1s electron. Energy expectation values E were obtained by averaging 20 recalculations of a 20,000-row spreadsheet, yielding a total sample size of about 30,000 r1, r2 pairs. The parameter ζ1 in the 1s trial function was held constant at ζ1 2.0, physically equivalent to taking Zeff Z for the inner 1s electrons. The parameter ζ2 for the 2s trial function was varied between 1.0 (complete shielding) and 1.6. The results are shown in Figure 3. The variation of E was fit to a quadratic in ζ2 for both the 3S and 1S. This fit was then used to determine the optimal ζ2 and minimum energy for each term. For the 3S term, this gives ζ2 1.37 and E 2.172 Eh; for the 1S term, ζ2 1.17 and E 2.155 Eh. Table 7. Energy Contributions for the 2s Electron in the He 1s2s 3S and 1S Termsa Triplet Values, in E h

a

Singlet Values, in E h

2

1.37

1.17

T2

0.152

0.097

V2

0.594

0.501

V1 2

0.271

0.247

All energies are given in atomic units (hartrees).

718


•

Figure 3. E vs ζ for the 1s2s 1S and 3S excited terms of helium, from ∼30,000 random r1, r2 pairs All quantities are in atomic units.

This is in good agreement with the corresponding analytical values: E 2.172 Eh for 1s2s 3S at ζ2 1.39; and E 2.155 Eh for 1s2s 1S at ζ2 1.17. For comparison, the experimental singlet–triplet energy difference is 0.029 Eh, considerably larger than the 0.017 Eh value obtained with this simple trial function. Significantly, even these simple calculations reproduce the correct energy ordering for the antisymmetric and symmetric versions of the trial function. They also convincingly demonstrate the origin of the singlet–triplet energy difference in the antisymmetry requirement. The lower energy of the triplet state can be understood by examining the numbers from the spreadsheet. Table 7 shows the various energy contributions for the 2s electron. Usually the lower energy of the triplet is attributed to higher electron repulsion in the spatially symmetric state. However the electron repulsion is necessarily lower for the triplet only if the same value of ζ2 is used for both states. Figure 3 shows that the optimized ζ2 is greater for the triplet state than for the singlet. This can be understood as representing a contraction of the triplet 2s orbital (greater effective nuclear charge) relative to the singlet state due to the greater repulsion in the singlet. Table 7 shows that the electron repulsion contribution V1s:2s is actually greater in the triplet than in the singlet when the optimized ζ2 values are used. But the increased ζ2 in the triplet also increases the (negative) potential energy of attraction between the 2s electron and the nucleus. This more than compensates for the increase in V1s:2s. Lithium Ground and First Excited States As a further example we consider the energies of the ground and first excited states of the lithium atom. The expectation value for the ground state energy may be written

E

= 2 T1s + 2 V1s + T2s + V 2 s + V1s:1s + 2 V1s:2 s

(26)

with a similar expression for the 1s22p state. The trial functions are taken as simple orbital products for the three electrons, where the one-electron orbitals are defined in eqs 8, 12, and 15. The one-electron orbitals depend on the param-


•

www.JCE.DivCHED.org


Figure 4. E versus ζ2 for the lithium 1s22s and 1s22p states, with ζ1 2.70. All quantities are given in atomic units.

Figure 5. E versus ζ2 for Be 1s22s2 and Be+ 1s22s2 states, with ζ1 2.70. All quantities are given in atomic units.

eters ζ1 and ζ2, which may be understood as the effective nuclear charges for the 1s and 2s or 2p electrons. Even though lithium is a three-electron system, the probability distribution needed to calculate the expectation values will only involve, at most, two electrons at a time. Thus, each of the expectation values on the right side of eq 26 has already appeared in the earlier calculations of this paper. For instance, the V1s:1s term is handled just as in the helium 1s2 ground state, while the V1s:2s term is handled just as in the helium 1s2s excited state. Random pairs (r1, r2) distributed according to |ψ1s(r1) ψ1s(r2)|2 are chosen to compute V1s:1s or according to |ψ1s(r1) ψ2s(r2)|2 to compute V1s:2s . Random values of the angular coordinates of each electron are also chosen in order to evaluate the inter-electron distance (eq 21). Fifty recalculations of a 20,000-row spreadsheet were averaged to obtain energy expectation values. The acceptance rate was about 13% for both of the two-electron distribution functions, yielding a total sample size of 130,000. The acceptance rates and sample sizes were of course much larger for the one-electron distributions. The 1s effective nuclear charge ζ1 was held fixed at a value of 2.7, reflecting the shielding expected of each 1s electron by the other. The 2s and 2p effective nuclear charge ζ2 was varied to minimize the energy. Different values of rmax were used for the 1s (5/ζ1) and 2s (14/ζ2) electrons. The resulting variation of E with ζ2 is shown in Figure 4 for both states. For each state the points were fitted to a quadratic, yielding an optimum ζ2 value of 1.41 for 1s22s and 1.02 for the 1s22p state. This shows the expected differential screening of 2s and 2p electrons. The higher ζ2 value is due to the higher probability that the 2s electron can penetrate the 1s2 core. The corresponding minimum energies E are 7.403 Eh for the 1s22s and 7.344 Eh for the 1s22p state. The excitation energy for the 1s22s → 1s22p transition from this simple VMC calculation is found to be 0.059 Eh, compared with the experimental value of 0.068 Eh.

third ionization energies of beryllium. The energy for Be may be written as E

+ V1s:1s + 4 V1s:2 s + V2 s :2s

www.JCE.DivCHED.org

•

(27)

Only the final term in eq 27 is new to the Be atom, and it may be calculated using methods already described. Furthermore, eq 26 gives the energy for ground state of the Be ion, which is computed just as for the Li ground state. Similarly, eq 20 gives the energy for the Be2 ground state, computed just as for the He ground state. In each case one must use the proper nuclear charge (Z 4) to evaluate the potential energy. Of course the change in Z will also change the energy minimization and the optimal values for the effective nuclear charges ζ relative to the isoelectronic states of Li and He. Table 8 shows optimized ζ and E values for Be, Be, and Be2, obtained using the VMC method described in earlier sections. The sample size was the same as for the Li calculations in the previous section. For the Be2 species ζ1 was optimized; there is no ζ2 parameter. For simplicity, ζ1 was not optimized for the Be and Be states; instead the optimal value of ζ1 from the Be2 calculation was used. The ionization energy for Be2 assumed the exact energy for one-electron Be3. Table 8 also shows the corresponding ionization energies for each species along with the experimental values (10). Figure 5 shows the variation of E with ζ2 for the Be and Be ground state. These results show the increase in the Table 8. Calculated Ground State Energiesa and Calculated and Experimental Ionization Energies for Be, Be, and Be2 Ionization Energies/ kJ mol1

Exp. Ionization Energies/ kJ mol1

1 / E h

2 / E h

E / E h

2

Be

3.70

2.52

14.503

735

899

Be2

3.70

2.18

14.233

1575

1757

Be2

3.70

—

13.594

14,670

14,848

Species

Ionization Energies of Beryllium Ionization energies are often used to illustrate the periodic properties of the elements. As a final example of the simple VMC approach, we calculate the first, second, and

= 2 T1s + 2 V1s + 2 T2s + 2 V2s

a

Ionization energies are in kJ/mol; all others are in atomic units (hartrees).


•


719


ionization energy after the first 2s electron is removed. This is due to the decrease in the effective nuclear charge for the 2s electron from 2.52 in Be to 2.18 in Be and the consequent expansion of the 2s cloud. The first, second, and third ionization energies calculated from these simple orbital product trial functions differ from experiment by 18%, 10%, and 1%, respectively. Of course the merit of the approach is not its accuracy, but the fact that it allows realistic calculations to be done in a simple way, using tools and techniques readily accessible to undergraduates.

this VMC approach is that it reinforces the idea of the wave function as a probability distribution. Seeing the actual generation of a random sample of properly distributed electron positions helps make concrete the concept of the “electron cloud”.

Summary

Literature Cited

This paper has demonstrated the practicality and usefulness of variational Monte Carlo calculations at the undergraduate level. Using only a commonly available spreadsheet program such as Microsoft Excel, it is possible for students to obtain numerical energies for states of simple, many-electron atoms. In order to keep the method as transparent as possible, a simple rejection method is used to generate the required random positions. In spite of the low efficiency of such a method, it easily succeeds in quantitatively illustrating a number of key concepts in atomic structure theory. These concepts include electron shielding, effective nuclear charge, l-dependence of the orbital energies, singlet–triplet energy splitting, and ionization energy trends. Though not presented in this paper, this spreadsheet approach can also be successfully used to obtain quantitative estimates of the L–S term energies in the ground state of carbon and the bond energy and equilibrium bond distance in H2. Usually these concepts are introduced merely qualitatively, but VMC calculations allow students to see them emerge from numerical results they can obtain for themselves using only familiar tools. Another pedagogical advantage of

1. Hammond, B. L.; Lester, W. A., Jr.; Reynolds, P. J. Monte Carlo Methods in Ab Initio Quantum Chemistry; World Scientific Publishing Co.: Singapore, 1994. 2. Galvez, F. J.; Buendia, E.; Sarsa, A. Chem. Phys. Lett. 2003, 378, 330–336. 3. Galvez, F. J.; Buendia, E.; Sarsa, A. J. Chem. Phys. 2001, 115, 1166–1171. 4. Alexander, S. A.; Coldwell, R. L. Int. J. Quantum Chem. 1997, 63, 1001–1022. 5. Alexander, S. A.; Coldwell, R. L. J. Chem. Phys. 1995, 103, 2572–2575. 6. Metropolis, N.; et. al. J. Chem. Phys. 1953, 21, 1087. 7. Cuthbert, H. L.; Rothstein, S. M. J. Chem. Educ. 1999, 76, 1378–1379. 8. Lowe, J. P. Quantum Chemistry, 2nd ed.; Academic Press: San Diego, CA, 1993; (a) pp 106–107; (b) p 191; (c) pp 598– 602. 9. Newman, M. E. J.; Barkema, G. T. Monte Carlo Methods in Statistical Physics; Clarendon Press: Oxford, UK, 1999; pp 396–404. 10. Emsley, J. The Elements; Clarendon Press: Oxford, UK, 1989.

720


•

W

Supplemental Material

The various spreadsheets used to obtain the results in this paper are available in this issue of JCE Online.


•

www.JCE.DivCHED.org

A Variational Monte Carlo Approach to Atomic Structure

Recommend Documents