A critical review of programs in atomic structure for bound states

COMPUTER PHYSICS COMMUNICATIONS 5 (1973) 147-152. NORTH-HOLLAND PUBLISHING COMPANY

A CRITICAL REVIEW OF PROGRAMS IN ATOMIC STRUCTURE FOR BOUND STATES Charlotte FROESE FISCHER * Department of Applied Analysis & Computer Science, University of Waterloo, Waterloo, Ontario, Canada

Received 7 November 1972

Programs in the Computer Physics Communications library for the calculation of atomic structure properties of bound states are summarized. Some of the difficulties encountered in inplementing many of these programs are mentioned and suggestionsare made for improvement of both the program and the published write-up.

1. Programs for atomic structure calculations The Computer Physics Communications Program Library now includes a fair number of programs for atomic structure calculations for bound states. These programs divide, quite naturally, into three classes: (i) programs for evaluating angular momentum integrals; (ii) programs for determining and manipulating radial wave functions; (iii) programs for special problems which may indude both (i) and (ii). 1.1. Programs for angular momentum integrals In atomic structure calculations, the total wavefunction of the system is approximated by antisymmetric function i,Ii constructed from products of oneelectron orbitals of the form [1, 2]

~r(0, ~)X(a; .L)

~(r, 0, 0, a) = r1P (ni; r)

in the non-relativistic approximation, or [31

/

ir’Q(nl;r)X_k,m

0(r) Q(nl;r) Xk

m

/ *

\ 1 / J

This research was supported by Grant No. A-3136 from the National Council of Canada.

in the relativistic approximation, where Xkm angular momentum eigenfunction Xkm

=

~,

iS

the

C(l~I; rn—a a) Y~ (0, ~)X(u; ~).

Here C(l~/;rn—a a) is a Clebsch-Gordon coefficient directly related to the 3 / symbol of Wigner, (0, ~)a spherical harmonic, X(u; -~-)a spin function; also k = —(j+.~)awith 1j —~aand a = ±1. Hence the evaluation of matrix elements such as (i~iIH~ i/i’) for the energy, or (i/i R (1)1 iii’) for a transition probability, require the evaluation of integrals involving spherical harmonics or, more generally, eigenfunctions of angular momentum operators. Two basically quite different approaches are possible. In the non-relativistic approximation, the total wave functions are eigenfunctions of the angular momentum operatorsL 2, ML, S2 and MS. Nussbaumer’s program [4] is based essentially on Condon and Shortley’s method for generating such eigenfunctions. The state is represented as a linear combination of determinants and the vector coupling coefficients are computed either by direct diagonalization of the operatorsL2 and S2 or by applying step-up and step-down operators. The advantage of this ap-

~r

—

once the linear combination of determinants is known, matrix elements can be evaluated quite simply by using the rules for determinants. The disadvantage is that the number of determinants may become quite large and in the evaluation of matrix elements, cancellation could occur which would reduce the accuracy. The seriousness of such cancellation has proach is that

148

C. Froese Fischer, Review of programs in atomic structure for bound states

not been considered. Of less importance, from a computational aspect, is the fact that the determinantal approach does not lead to any general formulas. To some extent, the shortcomings of the determinantal method have been overcome by Racah’s [5] powerful tensor operator methods for evaluating matrix elements. In the so-called “Racah algebra”, angular momenta are coupled, and recoupled, using 3 /, 6 /, and 9 / symbols (n / symbols is general). Furthermore, the concept of “coefficients of fractional parentage” is introduced so that a single electron may be uncoupled from a group of equivalent electrons. General rules for evaluating matrix elements have been developed and are summarized in Shore and Menzel [2] for example. The evaluation of the dcctrostatic interaction between configurations has been analyzed by Fano [6]. Rohrlich [7] has performed a similar analysis for transition probabilities. Wybourne [8] has derived formulas for many atomic properties with particular reference to configuration with f electrons. Several basic programs are now available for performing computations using Racah algebra techniques. Subroutines for 3 —/~6 j, and 9 —/ symbols have been published by Tamura [91 with an improvement by Wills [10]. These subroutines store tables of log (n!) so that calculations can be performed with large angular momentum quantum numbers. Burke’s subroutine [111 NJSYM calculates a general recoupling coefficient for an arbitrary number of integer or halfinteger angular momenta. Allison [121has written a pair of subroutines for fractional parentage coefficients for equivalent p shell and d shell electrons; a similar program for equivalent f shell electrons is under development. Finally, Hibbert [13] has integrated many of these subroutines into a program for evaluating the electrostatic interactions between configuration. A similar program has been developed by Robb [14] to evaluate the reduced matrix elements of one-particle tensor operators. Programs for other atomic properties like Hibbert’s and Robb’s program would form a valuable addition to a library of atomic structure programs. The programs mentioned above all assume an LS coupling scheme, but other coupling schemes are also receiving same attention. Recently, Grant [15] has published a program for fractional parentage coeffi—

—

—

—

cients in the! / coupling scheme. These coefficients are useful in relativistic atomic structure calculations irvolving open shells. —

1.2. Programs fbr radial wavefuncrions A variety of approximations are commonly used in atomic structure calculations for generating radial functions. One of the simplest, from a computational point of view, is the Hartree—Fock----Slater approximation used by Herman and Skillman [16] in their well-known program published in 1963. They use a linear r variable in their calculations with more than 440 points in the range of the functions. This program has recently been rewritten by Desclaux [17] to use the p = log r variable so that the number of points can be reduced to at most 220. The program by Klapisch [18] computes radial wavefunctions by the parametric potential method as well as many other related quantities such as Slater integrals, spin—orbit coupling, mixing of configurations, hyperfine structure and transition probabilities. But information about the angular integrals is part of the input data. A general Hartree—Fock program [19] is also available which can perform calculations for ground states and excited states with several open shells. It can also perform a variety of multi-configuration calculations though the configurations are restricted to those for which the interactions can be expressed as Fk, G k or R k integrals. This includes all two-electron substitutions and one-electron substitutions of the type nl—*n’l’, l’:~l. The program includes a large number of subroutines for evaluating many of the commonly occurring radial integrals; for example, the Slater integralsFk, Gk and R k; integrals Mk, Nk, and Vk which occur in the evaluation of spin—orbit interaction; the transition integral ~ ‘J-~i-P r dr ,)

ni ~ ~

ni

mean value (r’~)integrals; and the one-electron integral =

r

~

r~ + 2Z — 1(1 + ~\P I ~~dr2 r r2 /

nl’

nl

(r) dr

Finally, Liberman et al. [20] have published a pro-

C. Froese Fischer, Review of programs in atomic structure for bound states

gram which performs self-consistent field calculations

for atoms and ions using the Dirac operator in the oneelectron equations. A number of alternatives are permitted for the exchange potential. These include (I) the Slater type exchange, (2) Coulomb tails, (3) free-electron Hartree---Fock exchange, and (4) a modified Kohn—Sham potential. A p = log r variable is again used for the integration of the differential equations and the output includes the total energy, electron binding energies, radial wavefunctions, charge density, and potential function. Subroutines for evaluating atomic properties other than those associated with the energy are not included. The programs developed by Lewis [21] differ in that they deal with hydrogenic functions which are the zero order approximations in a perturbation theory expansion of atomic properties in terms of liz. Lewis has developed a set of routines for evaluating R k integrals using series expansion techniques (as distinct from numerical integration), which give results of high accuracy; iirthe test case the exact and-computed values agree to nine significant digits. Another program evaluates the single particle substitution contribution to the second-order term in the Z-expansion of the total energy. 1.3. Programs for special classes ofproblems The programs summarized so far can be applied to many configurations of most atoms. Though still quite general, the program by Beck and Zare [22] is restricted to atoms having a closed core and two valence dcctrons. Radial functions are obtained using the Hartree— Fock—Slater equations which are solved using a niodification of the Herman and Skillman program. The total wavefunction is expressed as a linear combination of single configuration functions which are eigenstates of the operatorsS2,L2,J2,J~and parity. The atomic hamiltonian is approximated by means of the low-Z Pauli approximation, though the program allows for the exclusion of any relativistic operators. Upon application of the variational principle, a secular equation is obtained whose solution then defines the n~ixing of the eigenstates and hence the total energy. In effect, this program combines the evaluation of matrix

149

cialized. It computes total wavefunction for the helium iso-electronic sequence using a configuration interaction type expansion with sturmian functions for the radial functions.

2. A critical review of programs and their write-up This summarizes some of the programs in the iibrary for performing atomic structure calculations. For the purpose of this review, several of these programs were implemented on the IBM 360/75 at the University of Waterloo. All of these programs were written in FORTRAN. The rest of this review will consist of some general recommendations to future authors based on the experience gained in implementing these programs, attempting to run and modify them. 2.1. Compilation of the programs and portability The first step in a computational task is to obtain a program which will compile correctly under the opcrating system that is being used. Such a program is said to be syntactically correct for that system. As was expected, not even programs written for other System 360 computers would compile without change. The most common fault was undefined variables which can result in underfiow/overfiow problems in our MVT operating system. For this reason, programs were run using WATFIV [24] which checks for undefined variables at execution time. In each case,

correct results were obtained by setting these variables to zero. It seems programmers assume that the operating system will clear memory when the pro. gram is loaded which may not be the case. With the compilers in use today, it is virtually impossible to write completely portable programs for all machines. By restricting the language to ASA Fortran [25], syntactically portable programs can be written which will compile with at most minor changes on most machines. However, it will not guarantee coin. pilation without error. For example A(I) B I)) = I 10) ‘

‘

elements with a Herman and Skiliman type program

valid list for inpul/oiitput in ASA Fortran. but will not compile in I B1\1/.~()Fort ran IV. which cx-

and an energy diagonalization process. The program by Knox [23J is somewhat more spe-

(A (I) B(t) i

is

il

peetS .

.

=

C.

150

Froese Fischer, Review of programs in atomic structure for bound states

One of the greatest shortcomings of ASA Fortran is

gram with the test data supplied with the program from

the extremely restrictive form of the DATA statement. Many programs in atomic structure rely on basic ta-

the library. All of the above problems were encountered but

bles of data such as the average energy of interaction, for example. Generally, these are quite lengthy multidimensional arrays. ASA Fortran requires that every element be listed as in

even so, they were mainly minor inconveniences. In each case a running program, which succesfully processed the test data, was obtained in a relatively short time by an assisstant who knew nothing about the program itself. A more difficult question than syntactic portabili-

DIMENSION ARRAY (2,2) DATA ARRAY (1,1), ARRAY (2,1), ARRAY (1,2), ARRAY (2,2) /1.0, 2.0, 3.0, 4.0/ Many compilers have a convenient way of implementing the initialization of arrays either by implied DO’s or by initializing the whole array in a pre-defined or-

der. It is a simple matter to convert from a program which uses the implied DO, to IBM FORTRAN which initializes the whole array in an order where the first subscript varies most rapidly and the last least rapidly and vice versa, as long as the implied DO’s are consistent with the IBM convention. With the implied DO the above DATA statement becomes DATA((ARRAY(I,J),I = 1,2)/1.0, 2.0, 3.0, 4.0/ whereas the IBM Fortran version is DATA ARRAY

/1.0,2.0,3.0,4.0/

When the implied DO’s are in the wrong order, conversion can become a major task since the order of the da-

ty is semantic portability whichon requires that programs execute in an identical manner all machines. The problems here are related primarily to the internal word length of the machine — the different ranges for integers, different precision and different ranges for floating point arithmetic.

The differing precision among computers is a factor not usually considered when developing a program for a particular computer. For example, the basic data arrays which define the average energy of a configuration consist of rational numbers such as 1/14 = 0.071 428 571 428 57• .To what precision should such data arrays be initialized? This question is difficult to answer without a particular computer in mind. On the other hand, the question does not arise at all if initialization is done at execution time by means of executable statements such as B(I,J)= 1.0/14.0 The constant will automatically be represented to ma-

ta itself must be changed. Special care must then be

chine accuracy. This fact should be kept in mind when

taken to ensure that errors are not introduced. Many programs did initialize “by rows” rather than “by columns” so that data cards had to be repunched. My own opinion is that ASA Fortran is too restrictive, that implied DO’s should be accepted for programs in the library, but that they should adhere to the IBM Fortran order and initialize the whole array. By the same token, IBM Fortran programs should indicate, by means of comment cards, what implied DO’s should be used to initialize the data. One area in which it is important to keep within ASA Fortran is in input/output statements. For example, to convert the statement

data arrays are initialized.

READ (IREAD, 1) ((A(I,J), J= 1, K(I)), I = 1, 10) 1 FORMAT (10F8.3)

2.2. Description of the input data Users refer to a library program for one of two reasons: either to implement a program which calculates a desired result (such as the coefficient of an angular momentum integral) or to adapt the program to solve another related problem. The user in the first category probably hopes he never has to become involved in the details of the program; however, he certainly must know how to prepare the data for the program, what the meaning of the variables are that form part of the data (or calling sequence), and what the limitations and restrictions of the program may be, Some programwriteups imbedded a description of the data requirements in

the general description of the program and the test da-

to ASA Fortran requires that the input data itself be

ta. This made it more difficult than necessary to locate

be repunched. The user then cannot even test his pro-

these data requirements. A card-by-card type descrip-

C. Froese Fischer, Review ofprograms in atomic structure for bound states tion (where the “card” may in fact extend over several physical cards) should be available to the user either in the published write-up or at the beginning of the listing of the program. The latter alternative has the added flexibility in that when an adaptation is published which changes the input format, an adaptation for the description of the input data can be submitted at the same time, One feature which distinguishes programs in atomic structure from, say, mathematical subroutines, is that the nature of the problem is more complex. It is more difficult to state precisely, in a manner readily understood by prospective users, exactly what a program does, what its limitations are, and what the input data requirements are. As a result, there is much greater possibility for error on the part of the user. Because of this I would recommend that programs check the input data for consistency wherever possible, inform the user of these errors, and terminate what could otherwise be a lengthy calculation with incorrect data. Several programs did have a large number of error comments incorporated into the program but error checking was by no means universal, 2.3. Numerical accuracy of the program

.

Checking the numerical accuracy of a large complex program can be a formidable task, and, if the errors are due mainly to round-off, may be quite machine dependent. In many instances though, exact values are known and it would be reassuring to know that for a particular machine a program can reproduce these results to a certam accuracy. Certainly, a program which solves a radial equation should check out its accuracy on the hydrogen equation varying both the principal and angular quantum numbers. One program, written for an S360 where floating point arithmetic has slightly more than six significant decimal digits, gave an energy parameter for the Is function which was accurate to only four significant digits. On the other hand, the energy parameter for the 7s function was also accurate to four significants digits which is somewhat surprising since generally the higher eigenfunctions cannot be determined as accurately as the lower ones. In a letter to the editor, Grant [26] has already emphasized the need for checking. This point may seem obvious but cannot be stressed too strongly. All special properties of the problem not directly incorporated into the method of solution should be exploited in de-

15 1

signing tests — symmetry, sum rules, orthogonality of solutions, and so on. 2 4 Efficiency -

The correctness of a program is, of course, of prime importance, but for library programs which may be used many times, efficiency is also desirable. Several measures of efficiency may be applied to a program but storage requirements and execution time are probably the two most obvious ones. Often there is a trade-offbetween storage and execution time; given more storage it generally is possible to write a program for which execution is faster. Other desirable programming practices (such as readability, modular design, etc.) also compete with efficient execution. However, once a program has been written it is extremely helpful to know where a program spends most of its time. By concentrating on those sections, it may then be possible to improve the efficiency with a relatively small amount of additional effort. For example, Gentleman and Douglas [27] have written a general Fortran timing routine for the S360 which keeps track of the number of times a subroutine is called and the amount of time spent in each subroutine. With the help of this routine, it was discovered that the MCHF routine [19] may spend 2/3 of the execution time computing functions which enter primarily into the calculation of the exchange function. It was then a fairly simple matter to check the statements in the DO-loops of the subroutines involved and re-arrange the most important arithmetic statement to reduce the number of multiplications. The resulting decrease in execution time was significant. As the same time, it has been helpful to know just how time-consumingan exchange calculation may be.

3. Conclusion Like many authors who have submitted programs to the CPC program library, I myself have been developing programs for many years without relying significantly on the work of others. The experience gained in implementing and testing library programs has been benefical to me in that I now have an extremely useful library program to supplement my own, which I probably would never have written myself. At the same time, some broad, general goals for a library of programs in the atomic structure field has become apparent. Firstly, the fact that a library program does not corn-

152

C. Froese Fischer, Review of programs in atomic structure for bound states

pile without change is not too important provided the changes are minor syntactical changes and do not require a re-arrangement of data and provided the program is of high quality. Program developers must ensure that the advantages of the program far outweigh the inconvenciences encountered in implementing the program for a particular system. That is, the program should be accurate, efficient, and well documented. Even errors in the program can be tolerated when the author is willing to cooperate in the debugging, Secondly, authors should not assume that users will be particularly knowledgeable about the problem; the program should be tolerant of user mistakes and detect as many errors as possible giving appropriate, helpful error messages. (One program actually detected an error and did nothing but give a variable, a somewhat larger value than usual!) Finally, significant benefits could accrue if programmers would begin to think in terms of a software systern for atomic structure calculations, In such a system, programs would make use of other programs in the ]i brary as the output of one program becomes part of the input for another. The Hibbert and Robb programs are examples of programs which make extensive use of other programs in the library. As already mentioned, similar programs for the evaluation of angular integrals for other atomic properties would be useful additions to such a software system. A modification of the Hib bert program to punch (or store on disk) part of the input data for the multi-configuration Hartree—Fock program has already been implemented and will be submitted as an adaptation. The concept of a software system may also have psychological benefits. We all know what we expect from a manufacturer’s software — accuracy, reliability, efficiency, relevant error messages, and clear documentation. If we apply the same standards to the programs in the library which we expect from other software systems, then the quality and usefulness of the library will improve.

Acknowledgements I gratefully acknowledge the programming assistance given to me by Mr. Alan 0 ‘Dacre who implement-

ed many of the programs mentioned in the paper on the IBM System 360/75.

References Condon and G.H. Shortley, The theory of atomic spectra (Cambridge Univ. Press, London, 1935). 121 B.W. Shore and D.H. Menzel, Principles of atomic spectra (J. Wiley, New York, 1968). [3] l.P. Grant, Proc, Roy. A 262 (1961)555. [1] E.U.

[4] H. Nussbaumer, Computer Phys. Commun. 1(1969)191. [5] G. Racah, PhysRev. 61(1942)186; 62 (1940) 438; 63 (1943) 367. [61 Phys. Rev. 140 J.(1965)A67. [7] U. F. Fano, Rohrlich, Astrophys. 129 (1959)44 L [8] B.G. Wybourne, Spectroscopic properties of rare earths (Willey — Interscience, New York, 1965). [9] T. Tamura, Computer Phys. Commun 1 (1970) 337. [101 J.G. Wills, Computer Phys. Commun. 2 (1971) 381. [11] P.G. Burke, Computer Phys. Commun. 1(1970)241. [121 D.C.S. Allison, Computer Phys. Commun. 1(1969)15. [13] A, Hibbert, Computer Phys. Commun. 1(1970)359; 2 (1971) 180. [14] W.D. Robb, Computer Phys. Commun, to be published,

[151 I.P. Grant, Computer Phys.

Commun. 5 (1973) 279. F. Herman and S. Skillman, Atomic structure calculations (Prentice-Hall, Englewood Cliffs, 1963). J.P. Desclaux, Computer Phys. Commun. 1(1969)216. M. Klapisch, Computer Phys. Commun. 2 (197 1)239. C. Froese Fischer, Computer Phys. Commun. 1(1969) 151;4 (1972) 107. [20] D.A. Liberman, D.T. Cromer and J.T. Waber, Computer Phys. Commun. 2 (1971) 107. [211 M.N. Lewis, Computer Phys. Commun. 1(1970)325, 265. [221 D.R. Beck and R.N. Zare, Computer Phys. Commun. 1(1969)113. [23] H.O. Knox, Computer Phys. Commun. 1 (1969) 167. [241P. Cress, P. Dirksen and J.W. Graham FORTRAN IV with WATFOR and WATFIV (Prentice-Hall, Englewood Cliffs, 1970). [16] [17] [18] [19]

[25] A Programming Language for Information Processing on Automatic Data Processing Systems, Commun. ACM 7 (1964)591. [26] I.P. Grant, Computer Phys. Commun. 2 (1971)383. [27] W.M. Gentleman and L.B. Douglas, FORSTAT — A FORTRAN Timing Program, Technical Report CSTR ioio, Department of Applied Analysis and Computer Science, University of Waterloo, Waterloo, Canada (1971).