David S. Alderdice and Ronald s. Watts University of New South Wales New South Wales, Australia 2033
I
I
Cclldati~nof Atomic Energy Levels A student exercise in theoretical chemistry
T h e past year or two has seen an increasing number of student exercises in various aspects of theoretical chemistry appear in THIS JOURNAL (1-3). This paper describes an exercise which serves as an example in introducing and coordinating three of the most important techniques required in theoretical chemistry; namely, the mathematical, numerical, and computer methods essential to most calculations using wavefunctions. At the same time, it is hoped that the exercise will yield some insight into the fundamental chemical concepts of atomic structure. Although the exercise is, in itself, simple, it is intended for advanced undergraduate students with.a working knowledge of basic quantum mechanics, Legendre polynomials (3), numerical analysis, and of Fortran or a similar scientific programing language. The exercise is intended to occupy two or three 3-hr sessions, spaced a week apart, to allow the computer program to be run. Mathematical
The aim of the exercise is to calculate the electronic energy levels of the isoelectronic sequence of the lithium atom, i.e., in the chemical notation Li, Be+, B2+, Ca+, etc., or in the physical notation Li I, Be 11,B 111, C IV, etc. In any atom comprising a single "series" electron outside a closed shell, the central field approximation may be used to give a Schrodinger equation of the form
VY
+ 12 E - V(r)l$ = 0
ing that the volume charge density p(r) at a point r is related to the electrostatic potential 6 by eqn. (4)
The charge density in atomic units (see Appendix) is directly proportional to the volume probability. Following Hartree (4), denote Y(n1;r)= TV(T)
and Y(nl;r)= N
where N is the nuclear charge, and q(n1) is the number of occupied wavefunctions in the closed shell given by n and I. Yo(nl;r) is a function whose form is found by expanding the right-hand side of eqn. (4) in terms of Legendre polynomials of argument
The generating function for Legendre polynomials is
If x
=
P(0) = P ( m ) = 0
(3)
and the normalization condition can be written JmP2(r)dr = 1
Now the central field potential V ( r ) experienced by the outer electron is due to both the nuclear charge and the Schrodinger distribution of electrons in the closed shell. The latter component can be computed by not-
1
1-2ooso-
If the angular part of the Laplacian is denoted by R, eqn. (1) may be written
and the requirement of stationary state solutions implies boundary conditions of the form
cos 0 then
-
(1)
Using the separation of variables technique with 2 = 0, 1 , 2, . . . r+ = P(r)Yz (@,.+) the radial function may be seen to satisfy an equation of the form
- q(n2) Ya(n2;r)
1;' + (f)"
therefore
where the effects of exchange have been omitted. I n order to calculate Yo,it is necessary to have a trial wavefunction P, for the closed shell electrons. The Hartree self-consistent field calculation employs an estimate for Pi, solves for the outer wavefunction, and calculates Yo, and thence PI until self-consistency is achieved. In this exercise, an accurate Hartree-Fock wavefunction for the inner electrons is used, and eqn. (2) is solved as an eigenvalue problem for the outer electron wavefunction P(r). A suitable analytic approximation to the Hartree-Fock function is given in (5). For example, for the lithium atom Volume 47, Number 2, February 1970
/
123
P&) = 2Za"re-"
But
(1 parameter)
or
P X ( R )= f ( R ) P ( R ) = C f ( R ) P d r ) = N?eCZr(l
+ kra)
( 2 parameter)
where Z, N, and k are variational parameters. Students are asked to derive eqns. (2) and (3) from the general Schrodinger equation, given the central field approximation and separation of variables. The formula in eqn. (4) for the electrostatic potential is then introduced and students are asked to use the generating function of eqn. (5) to derive the formula for Yo. The value of Yo for, say, the lithium atom, can then be evaluated analytically using one of the above wavefunctions (if tabulated or three parameter analytic wavefunctions are used, it would be preferable to perform this step by computer). The specific form of eqn. (2) for the lithium atom can then be written in closed form P"(v) = f ( r ) P ( r )
An approximate formula for the second differential of a function is
+ 0(h4)
PP(R) = Chaf(R)
Then, to a good approximation Y
=
[(l
+~)h~ff(R)I"~
(9)
The correct value of y to yield a solution in geometric progression is obtained by guessing a value (say y = O), substituting into eqn. (9) to obtain a better estimate, and repeating the process until y , + ~= y,. Once two values of P(r) have been obtained, the Numerov method is again used to integrate inwards to the matching point. The degree of mismatch is measured by the deviation of the quantity x
(6)
Numerical
Equation (6) is of the linear, homogeneous secondorder type, and the Numerov method (6) is applicable. Fourth-order differences in PXare neglected in this method, to yield a formula for P1 = P(rl) in terms of values of f(r) and P(r) at three successive values of r = r-1, TO,
from zero. If the solution does not match, the eigenvalue E (or c = 2E) must be made larger or smaller until x = 0. This corresponds to obtaining the correct number of nodes and the correct sign in the wavefunction for given values of the quantum numbers n and 1.
Ridley (7) has given a formula for the perturbation As to sin order to produce convergence to the matching
condition
x
=
0
Here, his the step interval in the independent variable r, i.e. The solution of the problem is complete only when the boundary conditions (eqn. (3)) have been satisfied as well as the radial eqn. (2). Hartree (4) h a shown that for small r the solution P(r) can be expanded in a power series
where A is an arbitrary constant (since the equation is homogeneous) and
where e is twice the energy of level of P(r); and vo, the potential (due to electrons) at the origin, -22. This approximation is valid, for the purpose of this exercise, up to r = 0.2. The solution can be started with h = 0.01 up to r = 0.2, and the Numerov method used up to some suitable value of r = ro called the "matching point" (say r = 8.0). The integrations have now to be started separately, from some large value of r to ensure that the solution is zero a t infinity. For large r, the solution P(r) will be approximately in geometric progression for successive values of r = R h, R, R - h, i.e.
+
where y is a small parameter and C a small (arbitrary) constant. Then
124
/
Journal o f Chemical Education
Thus, the numerical procedures to be followed are ( 1 ) Calculation of the function f(r) for all values of r in the range of integration. ( 2 ) Calculation of two initial vdues of P ( r ) far small values of r, from the formula. in eqn. (8). ( 3 ) Integration of eqn. (2) up to the matching point TO, using the initial values found above, and the Numerov method. (4) Calculation of two initial values of P ( r ) for large values of r by assuming the solution to be in geometric progression. (5) Integration of eqn. (2) down to the matching point TO, using the initial values found in step (4),and the Numerov method. (6) Calculation of the parameter y, from eqn. (10) to test the degree 01 mismatch. ( 7 ) Calculstion of A&) from the formula in eqn. ( 1 1 ) and correctine the value of e bv Ae. S t e.~ s(.1 .) to 16) . . are then repeated until .in = e..t, to the required accuracy.
-
Students are expected to understand the derivation of the Numerov formula in eqn. (7) and to be aware of the problems encountered in two-point boundary value problems. They should be able to derive the equations necessary for inward starting values, and to understand why normalization can be satisfied after integrating and matching. Computing
The procedures outlined above can easily be done on a computer. The programing is not complicated and can usually be surmounted by combining a series of simple Fortran subroutines available at most IBM computing establishments. For instance, the integrations required to calculate Ae and the normalization constant can be satisfactorily performed using Simpson's rule.
It is advisable to print out the second derivative of the solution wavefunction alongside the value of f(r) P(r),to ensure that the integration process is functioning correctly. If desired, the wavefunctions can be plotted as the matching process proceeds, and can then be checked for correct qualitative behavior such as number of nodes, sign a t infinity, etc. This is particularly easily done if the coding is in PL/l,as this work originally was. The student would usually be given a working program to run and would be expected to produce values of the energy eigenvalue for a particular atom, and given values of n and 1. These values of n and 1 can conveniently be done in one run, provided that the computer is fast enough. The values obtained by the integrations can then be compared with the experimental values of Moore (8). Example
For the lithium atom, a one-parameter radial function for one of the closed electrons is P t ( r ) = 2Zatgre-ZT
where Z = 2.960529
The function YOis then
Results o f Calculations for the Lithium Atom
n
Num. ber of Iter* 1 tions
-.
r,
X
r,
Ecx..
Ethssr.
Ps(r,) (em-') (cm-I)
A PL/1 program' has been written, for the IBM 360/50 computer, which integrates and matches the solutions for this equation. The matching condition x = 0 is tested using a simple numerical differential formula such as The values of the various parameters associated with the integrations for some states of the lithium atom are given in the table. The matching condition imposed here was 1x1 < 0.001. The agreement between experimental and theoretical energy eigenvalues is seen to be satisfactory for P and D states (1 = 1 and 2), but is not very good for the S state. This is because of the greater penetration of the core electrons by an s electron, and the consequent importance of exchange and polarization terms in the potential function; these were neglected in this calculation. Appendix Units: Atomic units are used throughout. This implies
The integration can be performed by parts and yields Yo(n1;r)= 1
- ~ - ' ~ ' (+l Z r )
and since
where me = electron rest mass, 6 = h/2* aa = 1st Bohr radius. Energy units are twice the Rydherg constant, and the general Schrbdinger equation is
q(n1) = 2
N = 3
(-'/z
Then Y(n1;r)= 3 - 2 ( 1 - ecmr(1 = 1 2ecam(1 Zr)
+
+
+ Zr))
Vs
+ U)$ = E+
Acknowledgments
R. S. W. wishes to acknowledge the financial assistance of the Colonial Sugar Refining Co. Ltd. Literature Cikd
The Schrodinger equation for the radial function is then
(1) 0 a r n s ~ 1 ; ~R. o . K.. J. C x ~ uEouc.. . 44, 286 (1967). ( 2 ) Pobbwow, G . F.,AND HOPPINOER, A. .I J. .. CHEX.EDUC.. 45, 528 (1968). (3) REBTE. L. F.. J . CHEM.EDUC.. 45. 151 (1969). 141 A ~ n m e e .D. R.. "Calculation of Atomic
Structure." John
Wilev R
Copies available upon request.
Volume 47, Number 2, February 1970
/
125
