A general program for computing matrix elements in atomic structure with nonorthogonal orbitals

Computer Physics Communications ELSEVIER

Computer Physics Communications 98 (1996) 235-254

A general program for computing matrix elements in atomic structure with nonorthogonal orbitals Oleg Zatsarinny Institute of Electron Physics, Uzhgorod 294016, Ukraine

Received 10 August 1995; revised 29 February 1996

Abstract To solve some problems concerning the atomic structure and transition probabilities, one needs to evaluate the matrix elements of various operators with respect to the configuration wave functions formed from the nonorthogonal one-particle orbitals. The program presented performs the angular integration necessary for expressing such matrix elements as weighted sums of radial integrals. Any amount of nonorthogonality between the orbitals may be presented leading to overlap integrals for the matrix elements. The calculations follow the method based on the representation of configuration wave functions through the Slater determinants. The program, therefore, can also be used for obtaining the corresponding vector-coupling coefficients. A distinguishing feature of the present procedure is the preliminary generation of vector-coupling coefficients for the individual subshells, which reduces considerably the amount of calculations. The operators included in the present version of the program are the one- and two-particle electrostatic interaction along with the one-particle tensor operator and simple overlap. Matrix elements of these operators provide the basis for studying the formation or decay of the inner-shell vacancies. Keywords: Matrix elements; Nonorthogonality; LS-coupling; Vector-coupling coefficients; Complex atoms; Wave function; Bound state; States integrals; Multipole transitions

PROGRAM SUMMARY Title of program: ZAP_ NO Catalogue identifier: ADDV Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland Licensing provisions: none

Programming language used: FORTRAN 90 (MS FORTRAN 5.1) High speed memory required to execute with typical data: 229 Kbytes Peripherals used: terminal, disk No. of bytes in distributed program including test data, etc.: 253368 Distribution format: ASCII

Computers: 386~486-based PCs Operating systems under which the program has been tested: MS-DOS

Keywords: matrix elements, nonorthogonality, LS-coupling, vector-coupling coefficients, complex atoms, wave function, bound state, Slate integrals, multipole transitions

0010-4655/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PI1 S 0 0 1 0 - 4 6 5 5 ( 9 6 ) 0 0 0 7 9 - 3

236

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

Nature of physical problem In many atomic processes involving inner atomic shells, the relaxation of eleclyon orbitals plays a important role. The relaxation can be most efficiently taken into consideration by using the nonorthogonal orbitals for the initial and final state. This requires the evaluating of the matrix elements of various operators with respect to the nonorthogonal one-particle orbitals. Initially, these matrix elements can be expressed as weighted sums of relevant radial integrals, possibly multiplied by overlap integrals. The program computes all the arising coefficients of the radial integrals and the corresponding overlap factors.

Method of solution At first, the configuration wave functions are expanded over the Slater determinants [1], using the combination of the vector-coupling and fractional parentage methods. Then the coefficients of the radial integrals and their overlap factors are obtained from integration over all spin and angular coordinates for the separate Slater functions which is a much simpler task. Integration over the radial coordinates is defined either as unity, zero, or as an expression for radial integrals. Due to the additional task of f'mding the determinant expansion, the method used is more laborious than those based on the Racah techniques [2,3] and is widely used now for calculating the matrix elements with orthogonal orbitals. On the other hand, the present technique admits a simple extension to the case of nonorthogonal orbitals by the most general way. Besides, considerable reduction of calculations has been achieved by using the tables of vector-coupling coefficients for the individual subshells.

Restrictions on the complexity of the problem Any number of s, p, or d electrons are allowed in a shell, but no more than two electrons or two holes in any shell of higher orbital angular momentum. In principle, the program admits the simple extension to the case of any f-subshell, but this requires the use of a very large auxiliary file (about 2 Mbytes) for the corresponding vector-coupling coefficients. LS-coupling is used, but its extension to jj-coupling is straightforward. The following parameters have also been adopted in the present program: the maximum number of electrons is 80, the maximum number of different subshells in one configuration is 20, and up to 20 configurations may be considered simultaneously in each run. These parameters may be augmented by changing dimension statements, if there is enough memory in the computer.

Typical running time The typical running time is 0.01 to 10 seconds on a 486-based PC for a single matrix element, and depends crucially on the complexity of the configurations involved, i.e. primarily on the size of the determinantal expansion and the number of electrons.

References [1] E.U. Condon and G.H. Shortley, The Theory of Atomic Spectra (Cambridge Univ. Press, Cambridge, 1935). [2] U. Fano, Phys. Rev. A 140 (1965) 67. [3] A. Hibbert and C. Froese Fischer, Comput. Phys. Commun. 64 (1991) 417.

LONG W R I T E - U P

1. Introduction

The calculation of many atomic and nuclear properties depends on the evaluation of matrix elements of various one- and two-particle tensor operators. The evaluation of the matrix elements requires integration over all the radial, angular, and spin coordinates. This can be done in stages by first integrating over all angular and spin coordinates, leaving the radial integrals in an unevaluated form. Several general programs are now available for performing such angular integration. Among the various techniques proposed, the method developed by Fano [1], based on the Racah algebra, turned out to be best suited for the systematic angular calculation including the complex configurations with several open shells. The matrix elements in this method are expressed as a product of two recoupling coefficients, spin and orbital, coefficients of fractional parentage, and one- or two-particle reduced matrix elements. Hibbert [2], following closely the scheme described by Fano [1] and using a general recoupling coefficient package developed by Burke [3], has written a program for the evaluation of the electrostatic interaction between arbitrary LS-coupled configurations. A similar program has been developed by Robb [4] to evaluate the reduced matrix elements of one-particle tensor operators, and by Klotz [5] to evaluate the matrix elements of the spin-orbit interaction. This method has then been generalized by Glass [6] to the case of two irreducible double tensorial sets of operators, of any order, and is used by Glass and Hibbert [7] for computing the angular integrals of the Breit-Pauli Hamiltonian. The above method has also been extended [8,9] to the case of jj-coupling to be used in relativistic calculations. Notwithstanding the great advance of the above technique, its applications have one restriction. All the above programs assume that both configurations in the matrix element under consideration are formed from the same orthogonal set of one-particle orbitals. The general formalism does not require orthogonality but its imposition results in a considerable simplification of the calculations. In many cases, this leads

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

to no loss of generality. However, at studying the number of atomic processes, the nonorthogonality arises naturally from the relaxation of electron orbitals, which is a consequence of unavoidable change in the screening of the nucleus by the surrounding electrons whenever there is a change of the electron configuration. First of all, this concerns the processes involving inner-shell vacancies. A vivid example of the use of the nonorthogonal orbitals is the creation probability of inner-shell holes during the radioactive decay of an atomic nucleus, or during the high-energy scattering. In the sudden shake-up approximation, this probability can be directly calculated through the overlap integral between the initial and final states [10]. It has also long been recognized that the consideration of nonorthogonality leads to a more correct structure for interaction matrix elements and has great importance for the accurate determination of radiative decay [11] or autoionization [12]. A comprehensive review of the application of nonorthogonal orbitals to atoms has been recently given by Kupliauskiene [13]. The nonorthogonality may be necessary not only to calculate the transition matrix elements but for the calculation of the separate atomic states as well. One important example has been discussed by Hibbert et al. [14] and lies in the following. In a configuration interaction or multi-configurational Hartree-Fock (MCHF) study of correlation effects, the most important additions to the Hartree-Fock scheme arise through one- and two-electron changes in the electron assignment of the initial configuration. The success of these calculations depends on the rate of convergence of the many-configuration expansion to the acceptable accuracy. It was found that convergence is faster if the orbitals associated with different pair correlations are not necessarily required to be orthogonal. Just for the above type of calculations, Hibbert and Froese Fischer [15] have developed a general program for computing angular integrals of the nonrelativistic Hamiltonian with nonorthogonal orbitals. The program is based on an extension of the Fano scheme [1], which formed the basis of the earlier WEIGHTS program [2] for the orthogonal orbitals. But it has the following serious restrictions (when we consider the matrix element ( ~ I H 1~')): (i) there are at most two subshells in containing electrons whose orbitals are not or-

237

thogonal to orbitals of ~ ' ; (ii) if all the spectator electrons with nonorthogonal orbitals have the same /-value, then there are at most two such electrons in each of • and ~ ' . For the purpose of the MCHF calculations, these restrictions do not cause a serious limitation in the flexibility of the wave functions. However, in calculations of transition matrix elements involving electrons from inner shells, the imposed restrictions limit strongly the wave functions used. Another approach to the problem of nonorthogonality is one where the wave functions are transformed to a new radial basis so that the evaluation of operators between states can proceed as though the radial bases were orthonormal. This imposes certain requirements on the wave function. A few years ago, an efficient method was found [16] which solves this problem in the case of states closed under deexcitation: for some preferred orbital ordering, if one or more orbitals in a determinant are replaced by lower-numbered orbitals, the results must be in expansion space. The above technique has been recently generalized [17] to the atomic, symmetry adapted case requiring the treatment of degenerate shells nl s, with arbitrary occupation numbers 0 < N < 41 + 2. The above approach is most efficient for the extensive multiconfiguration calculations based on a restricted active space configuration selection scheme (see [17,18] for details). From our point of view, the further extension of the Racah technique, used in the above-cited programs, leads to a rather complicated scheme for the angular calculations in the case of nonorthogonal orbitals. It is much more convenient to return to the methods which have been used in the beginning of the atomic calculations [19] and use, in particular, the determinantal representation of the atomic wave function. From a computational point of view, the latter technique is more cumbersome, because it requires the additional computation of determinant expansion coefficients and takes into consideration the magnetic quantum numbers. On the other hand, this technique admits a simple extension to the case of nonorthogonal orbitals in the most general way. The latter provokes further interest to use this straightforward (even though computational tedious) way to solve the problem of nonorthogonality by expanding a given atomic state function into a linear

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

238

combination of Slater determinants. A good example in this respect is the recent multiconfiguration calculations of the intercombination transition in Be-like ions [20] in the framework of the Dirac-Fock approximation. In this work, the matrix elements of the radiative multipole operator have been calculated including electron relaxation by using the determinantal expansion for the initial and final states obtained in the framework of a state-specific optimization procedure. The authors also mention that a program module for the complete expansion of an atomic state function in Slater determinants has been developed as an extension to the GRASP code [21], but, to our best knowledge, no further details have been published yet. The main problem in this technique is the effective way for obtaining the determinantal representation of an arbitrary atomic wave function. However, the task of obtaining the determinant expansion coefficients may be considerably simplified by using just the same Racah technique. This approach has been realized in the present program and is described in the following sections. The resulting formulae are very simple and convenient for computer programming.

where I L'l -- L31
2. Outline of the method

2.1. Matrix elements between the determinantal wave functions

In the configuration model of a many-electron system, a nonrelativistic wave function for an atomic state is expanded in terms of a configuration-state function, q~(3,LS), where y refers to the configuration and coupling of shells, and L, S are the total orbital momentum and the total spin. The configuration is defined by the q pairs of principal quantum numbers, n i l i . . , nqlq, where q is the number of subshells. Each open subshell is defined additionally by its term a i L i S ~. The total orbital momentum and total spin are obtained by coupling the momenta of all occupied subshells according to the given scheme, completely specified by the intermediate angular quantum numbers, L'i, ~.,

(--(((r,L2) (...

q~(TLSMLMs) = ( N ! )

. . . )s,

i z / 2 E i c i lu~i . . . uN),

(2) where the individual determinant l u~ " - - u N) has the form

= (ttlUl) [u I "'" u N) (tulul)

... ...

(tl luN) "

(tslus)

(3)

Here t = (r, s) stands for spin and position of the individual electrons, and u stands for the one-electron quantum numbers, I u) = [ nlmtz)). At this stage we have to introduce the magnetic quantum numbers in the description of states, in contrast to the abovementioned programs based on the Racah technique. This shows the prevalence of the Racah procedure for the construction of matrix elements of quantummechanical operators between many-particle states. On the other hand, matrix elements between determinantal wave functions can be easily generalized to the case of nonorthogonal orbitals.

Let us consider the matrix elements associated with two Slater determinants,

U = ( N l ) - 1 / 2 l ul, u 2 . . . . . u s ) and

V= (N!)-'/elui,

Pz . . . . . ~u),

which are built up from two basic set of spin-orbitals u l, u 2. . . . . u N and ~'l, ~'2. . . . . uN. For the sake of completeness, we will not impose any orthogonality condition on the set u~ and vj, and further assume that they have mutual overlap integrals ( u i l v j ) , which may differ from zero for i 4: j. To obtain the general expression for the matrix elements, it is convenient to express the determinantal wave function in the form

C, L3)C2L,)L'3 • - )L,

((( s,s )s',sOs' sOs'

-

N

(1)

U= (N!)-l/2E(-1)Pl--Iui(Pi), p

i~l

(4)

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

where P is an operator which permits the coordinates of the electrons, p is the parity of the permutation, and the sum is over all N! possible permutations. Using the above representation, one can show (see, for example, Ref. [22] or Ref. [23], ch. 6) that for any operator F, which is symmetric with respect to the electron interchange, the matrix element between the Slater functions U and V has the form

239

spherical harmonic and may be expressed through the Clebsch-Gordan coefficients as

ck( lm, l'nt) m' {(21 + 1)(21' + 1)}1/2 =(--1)

2k+l

X C(/01'0; kO) * C( 1 - ml'm'; km' - m). (9)

N

(UIFIV)=

E(-1)" f l-lu;(i) Fvei(i p =

) d'r i.

(3) Operator F

(5)

=

Y',i
1 (UI E - - I v ) i < j rij

Consider now the special cases for the operator F, involved in the present program.

N

=

E p < o - , p ' < oJ

(1) Operator F = 1.

1

X ( u , ( 1 ) u~(2) l ~ l

N

(UIV) =

~,,(1) ~ , ( 2 ) )

E(-1)PI-I (u, lup,) i

i=1

XO( p, tr, p', or'),

= Duv = det{duv (ij)}.

(6)

Here D vv is the determinant of the matrix formed from all overlap integrals between the one-electron orbitals u i and vj. If there is no risk for confusion, we will further omit the indices U and V. (2) Operator F = Ziu__,f~, where f~ is the one-electron operator with multipole k, acting on the electron i, ,~lf/' "=

=

~

(--l)(P+d)(u, lf'IPp,)

× D( p,

where D( p, tr; p', o-') is the second order cofactor of the determinant Dvv and corresponds to the matrix of overlap integrals with the orbitals p, odeleted from the initial set of orbitals and with orbitals p', o-' deleted from the final one; Pp,~,, permutes the two states nlmp~ labeled p' and o-'; the two-electron matrix element is expressed through the Slater integrals R k by 1

(up(l) u~(2) l -~-j215,(1) v,,,(2)) = Ex

p,p'= 1

(10)

(

k

p'),

(7)

where D( p, p') denotes the cofactor of determinant

Duv and corresponds to the matrix of overlap integrals with orbitals up and vp, deleted from the initial and final sets, respectively. For the k-pole transition operator r k included in the present program, the one-electron integral in (7) is reduced to (up I rkl vo,) = ck(lomp, l,,mp,)(PplrkIPp,)a(/xp,/~¢)

(8)

where (Ppl r k IPp, ) denotes the corresponding radial integral, c k is the matrix element of the renormalized

Xa(p,p, P'd) a(P'o" P'o-'),

(11)

and the angular integral Xk is po-,

= ck( Ipmp, Ip,, mp, )c'( l~m~, l~,ma, ) X 6 ( m p + m ~ , my+me,, ).

(12)

We see that the above formulae differ from the case of orthogonal orbitals only by the determinantal factor Dvv or its cofactors. In comparison to the existing codes for angular integrals based on the Racah technique, the sum in (7) or (10) is taken over all electrons but not over shells. It increases consider-

240

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

ably the amount of operations required. On the other hand, the present formulae are much simpler and convenient for programming. Only the ClebshGordan procedure is required. But the main advantage is that the above 'nlml,-representation' procedure allows the direct and rather simple extension to the case of nonorthogonal orbitals to be made in the most general way. The main problem here is the effective way for obtaining the determinantal representation of an arbitrary atomic wave function, including the configurations with several open shells. One of the ways for solving this problem is proposed in the present program.

2.2. Vector-coupling coefficients for the determinantal representation The expansion coefficients c i in (2) may be calculated by direct diagonalization of the operators L 2 and S 2 or by sequential coupling of the one-electron states, using the Clebsch-Gordan coefficients. The first method has been programmed by Nussbaumer [24]. The main shortcoming of this method is that it does not lead to a fully defined phase of the total wave function which is necessary in the manyconfiguration calculations. Besides, this method fails to give the vector-coupling coefficients for a given coupling scheme (1) with the indication of subshell and intermediate terms. It provides, in some arbitrary order, all eigenfunctions for the given configurations, but this information is often redundant and not convenient for applications. The effective application of this method, therefore, is restricted only by configurations with one open shell at most. The second method, the sequential coupling of one-electron states, is not applicable directly to the equivalent electrons due to the Pauli principle. The sequential coupling, therefore, may be used only with respect to the subshell terms, whereas the subshell eigenfunctions must be calculated beforehand. The latter can be obtained by using the fractional parentage method or again, by direct diagonalization of the operators L 2 and S 2. Due to the above-mentioned reasons, the fraction parentage method is seen to be more workable. This method is able to give the expansion coefficients for an individual spectroscopic term aLS, including the seniority number a , and it does not require any diagonalization proce-

dure. What is more important, the fractional parentage method allows one to introduce the self-consistent phase of eigenfunctions for different subshells and terms. Besides, as the number of different subshells is restricted, the considerable reduction of relevant calculations may be achieved by creating the tables of the subshells expansion coefficients. The latter is the main modification proposed in the present program.

2.3. Vector coupling coefficients for the individual subshells Let [INaLS) denotes the antisymmetrized wave function for the l N subshell with term aLS. It may be expanded in terms of the fractional parentage, [ INaLS) = ~,~'L'S' [ Iu- Ia'ES'INLS) X (l N- Ia'E~ILS]}INaLS),

(13)

where the subscript N specifies that the electron N has been removed from the antisymmetrized part of the wave function and II N- Ia'ES'INLS') is a coupled but not antisymmetrized wave function with respect to the electron N,

[I N- Ia'ES'INLS ) = ~_, Y'. C(EM'Llm; LML) C(S'M'ssiz; SMs) g'~ m gb, xll N- Ia'ES'M'LM's) I I,).

(14)

The function [ 1N- Za'ES') may be further expanded in terms of the Slater determinants l u'j • • • u s_ l),

I I N- 'oe'L'S'M'LM's) 1

-

~_,

C(u'; a'ES'M'LM's)

I / ( N - 1)! u(M't,M's)

x lu',

...

U'N_I > '

(15)

where u' specifies the corresponding set of one-electron states and u(M'LM's) denotes the Slater states with the same resulting M~ and M~. Taking together the relations (13)-(15) and considering the sums of M~, M~, m, and /, being replaced by the single summation over all possible states of the 1N electron, we obtain the final expression of the I lNaLS) wave

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

function in terms of the Slater determinants for the l/V- 1 subshell,

l l/VaLSMLMs ) -

E

+ 2-, + 2 ÷ .....

EE

¢( N - I) ! u(M~,M,s) u~ a'ES'

x (I/v- 'a'ES'ILSl}IUaLS) C( EM'LIm; LML) XC( S'M'sslx; SMs) XC(u", a'ES'g'Lg's) lu'~ ... u'/v_, I u~v). (16) On the other hand, the direct determinant expansion for the I INaLS) wave function has the form

I I/VaLSMLM s) 1

-

y"

I'-~!

C(u;aLSMLMs)

U(ML,Ms)

×lu, "'" us),

(17)

and by expanding the involved determinants over last row, we may extract the electron l/v,

I I/VctLSMLMS) 1

-

~.,

Y'.(--1)N-iC(u; aLSMLMs)

(I/-~ ! u(ML,Ms) ui

X l u I " " ui_lui+ 1 "'" u/v) ui).

(18)

Comparing now the expressions (16) and (18), we obtain the following relation between the vector-coupling coefficients for l/v and I u- 1 subshells,

C( u; aLsmLms) =V~

involve a f a c t o r - 1 to the vector coupling coefficients. In the present program, the following order of one-electron states nlmtm s was adopted: 0 - , 0+, - 1 - , + 1 - , + 1 + , - 2 - , - 2

1

Y'. (lU-la'ES'ILSl}l/vaLS) ot'L'S'

×C(CMFm; LMD

SMs)

X C(u'; a'L'S'M'LM's) ×(--1)N-i~(uI

u', . . . u'u),

"'" Ui_lUi+l "'" UNUi;

(19)

where the ~-function denotes that between the oneelectron states u and u' there is unique correspondence in the indicated order. In general, when we use the determinantal representation, we must adhere to some standard order of the one-electron orbitals, otherwise eventual reordering to standard order might

241

+, (20)

where an integer stands for the quantum number m and the superscript specifies the quantum number 1 /x = + 7. This order is convenient to use, because it has the same structure for the subshells with different l and, therefore, it is sufficient to indicate only the corresponding position in the above list to specify the quantum numbers m and /z. Expression (19) gives the rule for obtaining the vector-coupling coefficients for the subshell l N when the coefficients for the subshell l N- 1 are known. Vector-coupling coefficients may now be obtained by repeated application of (19), starting with N = 2. Note that in doing so we obtain the self-consistent phases for all terms and subshells, which also reflects the sign convention accepted for the coefficients of fractional parentage. This procedure was carried out for all p/v and d/v subshells, and the corresponding coefficients were stored in the form of two-dimension arrays C(i, j), where i specifies the Slater state and j specifies the subshell term. Each element C(i, j) is represented as a square root from the ratio of two integer numbers. The full relevant information, namely the list of all possible determinants, the list of subshell terms and the arrays C(i, j), are stored in file DET_ ARC and used in the present program as external information. In the present version of the program, file DET_ARC contains the data concerning only the shell p3, d 3, d 4, d 5, d 6, d 7, whereas the expansion coefficients for the more simple cases of the shells, containing not more than two electrons or two holes, are obtained by the program. It allows us to remove the restriction imposed on the possible electron orbital momentum in such shells. In the case of oneelectron shells, the expansion coefficients are trivial and equal to 1. For the two-electron shells, they can be obtained through the product of two ClebschGordon coefficients,

C( ulu2, LSMLMs) = (2)l/2C( llmll2m2; LMt.) N C ( S l / x l s 2/x2; SMs),

(21)

242

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

provided that//1 ~ u2- For shells with two vacancies, the expansion coefficients were found to be expressed through the same coefficients, ( - 1 ) M~C(ulu 2, LSMIMs), but here the symbols u I and u 2 denote the holes. In the case of the one-hole shells, the expansion coefficients are ( - 1 ) ME+I, while for the closed shells they are - 1. The above sing factors result from our choice of the standard order (20) for the one-electron states.

turn, depends on the chosen quantum numbers M L and M s . For the scalar operators, such as simple overlap (6) or inter-electron interaction (10), the matrix elements do not depend on the M L and M s. In this case, M L = L and M s = S are used, providing a minimum number of the relevant Slater states. For the tensor one-particle operators (8), our task is to obtain the reduced matrix element, which is defined through the Wigner-Eckhart theorem as

( yLSML M s I O~k) I y'L'S'M'L M's)

2.4. Final remarks

= 6(S, S') 6 ( M s, Ms) C(EM'Lkq; LML) With the availability of the determinantal representations for individual subshells, the task of finding the vector coupling coefficients for the configuration wave function is reduced to the sequential coupling of the subshell terms according to the given coupling scheme (1), and the corresponding coefficients have the form

× .

(23)

In this case, to determine the reduced matrix element, we use the q = 0 component of the matrix element with the maximum possible M s and M L compatible with the condition that the multiple factor at the right side of (23) is not equal to zero.

C( aLSML Ms; u) 3. Program structure and operation

q

=

E

I-I C(aiLiSiMLMs,; ui)

Mi,M~,ui i= 1 q-I

× H C(Ei-IM'L,_,LiML,; EiM'L,) i

× C(S~_ iM's,_SiMs~; S'iM's~).

(22)

Here the sum is taken over all possible magnetic quantum numbers for the subshell and intermediate terms, provided that T.ML, = M L and S, Ms, = M s, and over all associated subshell determinants. The full determinant [u) is a straightforward superposition of the subshell ones, l u) = l u I • -- Uq), careful following a given order of subshells and one-electron states. From the last equation we see that to obtaine the vector-coupling coefficients one has to use only the Clebsch-Gordan procedure. The present algorithm for matrix element calculations consists of two independent steps: the determination of the determinantal representation for two input configuration wave functions with certain quantum numbers M L and M s, according to Eq. (22), and the calculation of matrix elements between all obtained Slater states by the formulae (6)-(12). The amount of calculations crucially depends on the number of corresponding Slater states, which, in

The structure of the program is depicted in Figs. 1 and 2. The central subprogram is the subroutine ZNO_ coef, constructed in such a way that it can be used as an additional independent block for angular calculations in other, more general, programs. A description of the main program ZAP_ NO follows. All other subprogram modules are subsequently described in the order in which they are called.

Program ZAP_ NO This is a test driving program for running the central subroutine ZNO_ coef. All other routines provide only the input/output of the required information. Execution begins by attempting to open FORTRAN units 5 and 6 for input/output. Then the

Fig. 1. Top level block diagram for the ZAP_ NO program.

O. Zatsarinny/ Computer Physics Communications 98 (1996) 235-254

program reads the information strings and the full list of input configurations (subroutines R_core, R_ conf, Add_ core, and Pri_ conf). This latter information is stored in the archive array IC by subroutine Red_conf. All data are tested to the extent possible. If an error occurs at any stage, a message is issued to FORTRAN unit 6 and the execution is terminated. Then the program reads the additional orthogonality conditions for the radial orbitals and parameters, which specify the type of operators under consideration. The last input parameter ipr_d indicates that the output of the vector-coupling coefficients is also of interest. If any operator is not specified to be considered,then the program works in the mode of generating the vector-coupling coefficients. It is the additional possibility of the present program. Further, the principal loop of program ZAP_ NO over all input configurations is executed, and the matrix elements between two configurations are treated subsequently by the subroutine ZNO_coef. The resulting vector-coupling coefficients and the structure of matrix elements are printed to FORTRAN unit 6 by subroutines Pri_ dete and Pri_ coef, respectively. Note also that angular momenta are either integer or half-integer quantities. So to be able to use the FORTRAN INTEGER data type, all angular momenta and their projections are used internally by ZAP_NO in the representation ( 2 / + 1).

!

3.1. I n p u t / o u t p u t

243

routines

S u b r o u t i n e R _ core

This routine reads the number of common shells, NZ, and, if NZ > 0, reads the core configuration, which is common for all input states. This information is stored in C O M M O N / C O R E / .

Subroutine R _ c o n f

This routine reads all data which define the configuration and coupling scheme for one state. The description of the configuration contains the number of shells, N, and the nl q values for each shell. The quantum numbers of the shells coupling scheme are read into the array LS(k, i). k = 1, 2, 3 corresponds to the seniority, 2 L + 1 and 2 S + 1 quantum numbers, respectively. The first N columns (1 < i < N) in LS correspond to the quantum numbers of the occupied shells. These numbers are followed by ( N - 1) further numbers, which give the intermediate terms resulted from the coupling of the occupied shells (here seniority = 0). If the 'seniority' of the final total terms is zero, then the successive coupling of shells are assumed: the first intermediate quantum numbers are obtained by coupling the first and the second occupied shells, then the second intermediate coupling is obtained by coupling the third occupied shell with the first intermediate term, and so on. If the 'seniority' of the final total term is equal to 2 or 3, then subroutine R_conf reads additionally the

zNo_coEF

!

|

Fig. 2. Block diagram for ZNO_ COEF. Boxes in dotted face represent the two main independent steps of calculations: determining the determinantal expansion and computing the angular coefficients between two determinants. Double boxes indicates the exchange with the disk file.

244

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

array JJ(k, i), which defines the coupling of the shells. The value of each elements in the array JJ gives the position in the LS array, and each row in JJ describes one elementary coupling: LS(JJ(1, i ) ) + LS(JJ(2, i ) ) ~ LS(JJ(3, i)). If the 'seniority' of the final total terms is equal to 1 or 3, then R_ conf reads additionally the array KN(i), which indicates the set number of radial orbitals for each subshell i. These quantities are used to treat the nonorthogonality conditions for the radial wave functions in the further calculations. Subroutine Add_ conf The data read by R_ conf are only supplemented to the core information, and if the fixed core is used, R_ conf also calls the subroutine Add_ conf to obtain the full description of the state. The final data are stored in C O M M O N / I N P / . Subroutine Pri_ conf Prints the spectroscopic notation of one from the two configuration placed in C O M M O N / I N P / . Subroutine AL Provides some spectroscopic symbols needed for printing the configurations and angular coefficients.

tween the orbitals i and j for two input configurations, placed in COMMON / I N P / . It is assumed that only the orbitals of one set are mutually orthogonal. Subroutine Pri _ dete Prints the determinant expansion for one configuration from / I N P / with given M L and M s. In addition to the real coefficients, their integer representation is also output. Subroutine Pri _ coef Prints the resulting angular coefficients with the corresponding overlap factors from the COMMON / O U T N / . The integer representation of angular coefficients is also output. For the individual determinantal multiplier, the following abbreviation is used: nil,

n21

n'll,

d21

•

•

•

• ..

nql n'ql

( n l l ln'll)

(n211n'2l)

-..

(nqlln~l)

(n211n'll>

(n2lln'21)

• ..

(n21ln~l)



•"

(nqlln~l)

•

.

.

.

( n q l l n'll>

(24) Subroutine Red_ conf Stores in or extracts one input configuration from the archive IC(mc,261) to COMMON / I N P / . The parameter mc defines the maximum of configurations, which can be considered simultaneously. The full description of one state requires 261 (INTEGER * 1) locations, and in the present programs the parameter mc was adopted to 20. The user may change it to any required value, based on the available memory. Subroutine R _ iort Reads the pointer JORT, and, if JORT > 0, reads the additional nonorthogonality conditions imposed on radial orbitals and places the information in COMMON / Z O R T / (see description o f / Z O R T / for details). Subroutine Pre_ iort Based on the set number KN, this routine prepares the orthogonality conditions IORT(i, j) be-

Subroutine NUM(Z, K1, K2) For given real number Z, this routine finds such integers K1 and K2 that Z = ( K 1 / K 2 ) 1/2. 3.2. Subprogram ZNO_ coef This is the central routine which performs the calculations of the angular integrals and corresponding overlap factors for the matrix elements between two input configurations in the case of nonorthogonal orbitals. The description of two input configurations is given by C O M M O N / I N P / , the orthogonality constrains imposed on radial wave functions are placed in COMMON / Z O R T / , the type of operators under consideration is given by the formal parameters list, and the resulting angular coefficients with the corresponding overlap factors are placed in COMMON / O U T N / . Thus, this routine can be used as an independent angular block in other, more general, CI or MCHF programs.

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

The structure of the subprogram ZNO_coef is depicted in Fig. 2, and with some auxiliary routines contains two main independent blocks, indicated by dotted frames: the block of calculation of the determinantal expansion (at the left of Fig. 2) and the block of determination of angular coefficients between two Slater determinants (at the right of Fig. 2). Function I_ diort

This routine controls the orthogonality of two input configurations states according to the number of electrons, the triangle rule for the total spin and orbital momentum, and parity conservation. Subroutine Red _fact

According to the Wigner-Eccart theorem, this routine gives the normalization factor which connects the reduced matrix element with the matrix element for the maximum value of M L and M s. If this factor is equal to zero for maximum M L, the other value of M L is selected. Subroutine D e t _ arc

This routine serves to store in or extract the vector-coupling coefficients and corresponding determinants from the archive arrays CA(mdet) and NBA(320,mdet). The parameter mdet defines the maximum of determinants, the description of which can be stored in memory. In the present program, the parameter mdet was adopted to 40. In the cases with several open shells, the number of relevant determinants may be very large and varies considerably from case to case. For their storage, therefore, it is better to use the disk file. If the number of determinants exceeds the parameter mdet, then to store these additional determinants the scratch file DET_ ARC is created. This is a binary file with direct access. The record length is equal to 328. Such file structure allows us to extract quickly the required determinant. As the operation with the memory is more effective, the user can increase parameter mdet to the maximum possible value based on the available memory and the complexity of the configurations under consideration. Subroutine ZNO_ LS

This routine prepares the auxiliary internal arrays LS and N on the basis of the determinant description

245

NB (see COMMON / I N P N / ) . LS(i) specifies the angular symmetry of the ith orbital, N(i) gives the corresponding principal quantum number. The angular symmetry of the one-electron wave functions is stored in the LS array in the packed form l X 1000 + (m t + 50) X 100 + (2m~ + 1). It offers savings in memory and allows us to compare quickly the full angular symmetry of orbitals. Besides, the LS and N arrays differ from NB by the ordering of the orbitals. In LS and N, the orbitals are ordered according to their angular symmetry, which allows us to compare directly the determinants and define most effectively the interacting electrons and the corresponding determinant factors. Note that the parity of the above reordering must be included in the resulting angular coefficients. Subroutine Press_ c

Sequential consideration of the matrix elements between the individual Slater functions gives rise to the same radial integrals. After a full summation over determinant expansions, many angular coefficients become equal to zero. Subroutine Press_c excludes such coefficients from the output list placed in C O M M O N / O U T N / . Subroutine Order_ c

This routine sorts the resulting angular coefficients according to the type of radial integrals and overlap factors, for the sake of conveniency of their further consideration. 3.3. Block f o r calculating the determinant expansion coefficients

This block can be used as a fully independent one for determining the vector-coupling coefficients for one configuration state. Subroutine D e t _ expn

According to Eq. (22), this routine calculates the determinant expansion coefficients for one configuration using the expansion coefficients for the individual subshells provided by the subroutine C_ shell. The full description of the configuration, with additional indication of the quantum numbers M L and M s , is accepted through the formal parameters list. This subroutine returns only the total number of

246

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

resulting determinants, whereas the description of individual determinants and corresponding expansion coefficients is stored in the memory or in the file DET_ ARC by means of the subroutine Det_ arc (see above). Function C_shell (l, q, IT, ID, M L, M s, IDET) This function is equal to the expansion coefficient of the determinant ID for the subshell I q with term IT. The parameters ID and IT are only the corresponding pointers, relative to the adopted common list of shell determinants and terms. The function C_shell also returns the corresponding quantum numbers M L and M s and the description of the determinant as a list IDET of one-electron orbitals (more precisely, their positions in the common list of one-electron orbitals (20)). For the subshell p3, d 3, d 4, d 5, d 6, d 7, the required information is read from file Det_ exp.dat, where the subshell expansion coefficients obtained according to the procedure of Section 2.3 are stored. For other more simple cases, ii, 12, 121+4, 12/+3, and 121+4, the data are obtained by the program (see the last paragraph of Section 2.3). Thus, for many configurations which are widely met in practice and consist of the closed shell and shells with no more than two electrons or two vacancies, the present program may be used without the data file Det_ exp.dat. Function K_term(1, q, a, L, S) This function returns the position of a given term a L S according to the adopted total list of terms, which can be formed from the 1q subshell. The list of the subshell terms is stored in some internal arrays of the subroutine Term, and its ordering must coincide with the term ordering used to store the subshell expansion coefficients in file Det_exp.dat (see the subroutine C_ shell above). Function Iglq(l, q) This function gives the total number of determinants for the given I q subshell. This value is simply equal to the total statistical weight of the subshell (41 + 2 ) ! / q ! / ( 4 l + 2 - q)!. Subroutine Clebsh This function evaluate the Clebsch-Gordon coefficients by using their expression through the 3jsymbols.

Function Z_ 3 j This function calculates the 3j-symbols by means of an ingenious method not requiring the direct use of factorials. The 3j-symbol is expanded as a sum of integer-product ratios. Then, the corresponding numerators and denominators are compared and all equal integer components are cancelled. The resulting values are proved to be rather small. The above procedure allows one to consider the large input orbital momenta (the procedure has been tested for the orbital momenta up to 1000) without any troubles arising due to extremely large values of the factorials, through which the 3j-symbols are generally expressed. Subroutines M L _ id, MS_ id These routines return the m t and m s values for the one-electron orbital with the pointer ID (according to the adopted common list (20)). 3.4. Block for the calculation o f the matrix elements between the determinantal wave functions

This block is depicted at the right of Fig. 2 by a dotted frame. It contains a set of routines one for each operator under consideration: ZNO_ ms treats a simple overlap operator, ZNO_ mi treats a one-electron scalar operator, Z N O md treats a one-electron tensor operator, and ZNO_ mee treats a two-electron operator for electrostatic interaction. The first three operators may, in principal, be treated by one routine, but their separate consideration gives rise to a much more effective procedure. All these routines have the same structure and we consider, as an example, only one, Z N O mee. Subroutine ZNO_ mee According to Eqs. (10)-(12), this routine computes the angular coefficients for electrostatic interelectron interaction between two given determinants, which are now defined by the arrays LS 1 and LS2 in BLANK COMMON (see the subroutine ZNO_ LS). The calculations are repeated in the four-fold sum over all one-electron states. The sequence of operations and callings in these loops is as follows. First of all, the l, m, and /x quantum numbers of the interacting electrons are determined by using the subroutine ZNO_mls. Then a set of orthogonality

O. Zatsarinny/ ComputerPhysics Communications98 (1996)235-254

247

restrictions are verified: the restrictions imposed on the interacting electrons according to the &functions in Eqs. (11)-(12); the restriction due to the angular symmetry of the spectator electrons, which must be the same for two considered determinants; and the restriction due to the additional orthogonality of the radial wave functions, imposed by the input data (by using the subroutine I_ iort). If all orthogonality conditions mentioned above are satisfied, then the determinantal overlap factor is determined by the subroutine Det_ fact. Then the angular coefficients for all possible multipoles k are calculated according to Eq. (12). The resulting coefficients and overlap factors are stored in C O M M O N / O U T N / b y the subroutine Add_ coef.

the description of the overlap factors requires different memory and may be the same for different angular coefficients, it appeared more appropriate to describe these factors in their own arrays. Thus, the output information in C O M M O N / O U T N / c o n t a i n s the following information: the list of angular coefficients and the corresponding radial integrals in the arrays AK and KAR, the description of the structure of the corresponding overlap factors in the array ND, and the list of the principal quantum numbers for the involved one-electron states in the array NN (for further details, see the description of COMMON / O U T N / ) . The above structure of the output information allows the resulting matrix elements to be stored in the most memory-saving way.

Subroutine ZNO_ mls For a given orbital i, this routine extracts the l, m t, and m s values from LS(i), where these three quantum numbers are packed into a single integer.

Subroutine Det_ simple If the input data imply the partial orthogonality of the radial wave functions, then the overlap factor can be further simplified. Subroutine D e t simple performs this procedure. As a result, some determinantal components D i are reduced to unity, or their dimensions are decreased.

Function I_ ort According to the information in COMMON / Z O R T / , this routine analyzes the determinantai overlap factor with respect to the additional orthogonality constrains on the radial wave functions. I_ ort = 0 if this factor is zero. Function Zcklm By calling the function Clebsh, this routine evaluates the value of c k coefficients according to Eq. (9). Subroutine Det_ fact Determines the overlap factors by analyzing the spectator electrons. The overall overlap factor D( p, tr; p', or') in Eq. (10) is the determinant which is formed from all overlap (monopole) integrals and from which the rows p, tr and columns p', tr' were deleted. Due to the angular symmetry of the one-electron states involved, this determinant has the block structure. The subroutine Det fact reduces D ( p , or; p', tr') to the product of much smaller determinants, FI i Di(l, q, k, n, n'), where each component D i is characterized by the following quantities: the orbital momentum l and the principal quantum numbers n, n' of the involved one-electron states, and the dimension q and the extent k of the corresponding determinant. In the output listing, such determinants have the two-strings notation (24). As

Function INOINL Locates position of the nl orbital in the description of the input configuration. This routine is used by I ort and D e t simple to obtain the correspondence between the overlap integral ( n l [ n ' l ) and the elements of array IORT(i, j) which defines the imposed orthogonality conditions. Subroutine Add_ coeff As already indicated, the sequential calculation of the matrix elements between the individual Slater states gives rise to the same radial integrals. Subroutine Add_ coef controls such cases and adds to or modifies the resulting angular coefficients in the output list in C O M M O N / O U T N / , as the situation requires.

4. Description of COMMON BLOCKS ( 1 ) / C O R E / N Z , NN(20), L(20), IQ(20), LS(3, 40), KN(20), JJ(3, 20) Contains the common part of input configurations: NZ - - number of shells.

248

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

NN(I), L(I), I = 1, NZ - - (n, l) values for each shell. IQ(I), I = 1, NZ - - number of electrons in each shell. LS(K, 1), K = 1, 3; I = 1, NZ + N Z - 1 - - angular momentum quantum numbers. K = 1 is seniority; K = 2 is ( 2 L + l ) ; K---3 i s ( 2 S + l ) . I = I , NZ is the set of quantum numbers for each shell, the remaining ( N Z - 1) values define the intermediate quantum numbers resulting from the coupling of the shells. KN(I), I = 1, NZ - - optional parameters, which indicate the set numbers of radial orbitals and are used only to define the nonorthogonality conditions. JJ(K, I), K = 1, 3; I = 1, NZ - 1 - - optional parameters which define the coupling of shells in the case when this coupling is not sequential. The value of each element in JJ gives the position in the LS array and each row in JJ describes one elementary coupling: L(JJ(1, I)) + L(JJ(2, I) L(JJ(3, I)). (2) / I N P / N1, NNI(20), Ll(20), IQl(20), LSI(3, 40), KNI(20), JJl(3, 20), N2, NN2(20), L2(20), IQ2(20), LS2(3, 40), KN2(20), JJ2(3, 20) Contains the full description of two configuration states q~t and ~2- The corresponding quantities have the same meaning as in the description of the common core (see above). This information is the only input data (with the orthogonality conditions in / Z O R T / ) needed for calculating the single matrix elements (~1 IO1~). ( 3 ) / Z O R T / I O R T ( 2 0 , 20), JORT Contains the orthogonality conditions between the radial orbitals of two input configuration wave functions ~ l and q~2, described i n / I N P / : JORT - - indicates the orthogonality mode: JORT = - 1 - - all radial orbitals are considered to be orthogonal; JORT = 0 - - all radial orbitals are considered to be nonorthogonal; JORT = 1 - - orthogonality conditions are defined by the array IORT; JORT = 2 - - orthogonality conditions are derived from the arrays KN1 and KN2 in / I N P / . It is assumed

IORT(I, J), I = number of IORT(I, J ) = integral (Pi[ any value.

that only the orbitals of one set are mutually orthogonal. 1, NI; J = 1, N2 - - NI(N2) is the shells for q~l (~2), respectively. 0 (or 1) indicates that the overlap P~) = 0 (or 1), otherwise (Pi [ Pj) is

(4) / I N P N / NE, NV, C1, C2, NBI(4, 80), NB2(4, 80), ILl, IL2, IS1, IS2, MS, ML, CN, LS3(3, 80), LS4(3, 80), N3(80), N4(80), KZ1, KZ2 Lists the quantum numbers and parameters for full description of two Slater determinants: NE - - number of electrons (the determinant dimension). C1, C2 - - expansion coefficients for these Slater functions. NB(K, I), K = 1, 4; I = 1, NE - - describes a Slater determinant as a set of (n, l, m t, m s) values for each electron, m s being represented as a 2m s + 1 integer ( = 0, 2). IL, IS, ML, MS total orbital momentum, total spin, and the corresponding magnetic quantum numbers. CN - - the ratio of the matrix element under consideration to the reduced one for the given ML, MS. NC1, NC2 - - number of Slater functions in the expansion of qb I or ~2, respectively. LS3(I), LS4(I), N3(I), N4(I), I = 1, NE - - the internal working variables. The LS arrays describe the angular part of the one-electron wave-function in the parked form [l × 10000 + m t + 50) × 100 + (2m s + 1)], and N3, N4 contains the corresponding n values. These arrays are ordered according to the increase of LS, and then that of N. KZ1, KZ2 - - the number of permutations needed for transformation from the initial order of the one-electron orbitals given in NB to the order of the orbitals in LS and N. -

-

(5) / O U T N / KIR(8, mc), RK(mc), ND(4, md), NN(2, ran), KRK, m_ id, m in Contains the output angular coefficients and the corresponding overlap factors which fully define the single matrix element (q0j IO 1~2) for ~ l and q~2 given i n / I N P / . The parameter mc limits the maximum number of coefficients, md = mc × 10, and mn = md x 5. KRK number of output radial integrals; -

-

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

KIR(K, I), K = 1,8; I = 1, KRK - - describes the Ith radial integral: KIR(1, I) is its multipole; K(2-5, I) are the (n × 100+ l) values for the involved radial orbitals; KIR(6, I) defines the type of integral; KIR(7, I) is the number of overlap determinants; KIR(8, I) locates the description of the overlap factor in array ND. RK(I), I = 1, KRK - - values of the corresponding angular integrals. ND(K, J), K = 1, 4; J = KIR(8, I) + 1, KIR(8, I) + KIR(7, I) - - the list of the overlap determinants for the Ith angular coefficient: ND(1, J) is the /-value for the involved one-electron orbitals; ND(2, J) is the dimension of the determinant; ND(3, J) is its extent; ND(4, J) locates the description of one-electron principal quantum numbers in list NN. NN(K, N), K = 1, 2; N = ND(4, J) + 1, ND(4, J) + ND(2 J) - - list of the principal quantum numbers for the orbital, involved in the overlap determinant J.

249

NT - - number of shell terms; MS(i), ML(i) - - magnetic quantum numbers for the ith determinant; IDET(i) - - description of the Slater determinant as the list of pointers for the involved one-electron orbitals, according to the common list (20); NP(j), IP(ij) - - give the integer representation of the expansion coefficient, corresponding to the ith determinant and jth terms, in the form

(IP( ij) /NP( j) ) I/2. To save memory, the latter arrays are written as INTEGER*2 variables. The file Det_exp.dat is opened only in the case when the input configurations contain the above-indicated open shells. The size of this file is 33126 bytes.

File Det_ arc This is the scratch internal file which is used for recording the large determinant expansions. It is connected with unit 1 and opened only in the case when the number of determinants exceeds the parameter mdet (see the subroutine Det_ arc).

6. Input data (6) B L A N K C O M M O N - - contains only the inter-

nal working variables.

5. Files used

File Det_ exp.dat This file contains the vector coupling coefficients and the list of the corresponding determinants for the p- and d-subshells. The expansion coefficients have been obtained according to the procedure described in Section 2.3. The use of these tables reduces considerably the amount of calculations and considerably simplifies the final algorithm. File Det_ exp.dat is a binary, sequential file, consisting of 6 records, one for each following subshells: p3, d 3, d 4, d 5, d 6, and d 7. Each record has the following structure: IQ, KD, NT, (MS(i), ML(i), IDET(i), i = 1, KD), (NP(j), (IP(ij), i = 1, KD), j = 1, NT), where IQ - - number of electrons in the subshell; KD - - number of all possible determinants;

Except for the first three information strings, which may be ignored, all the data required by the program are integer in form and the only format statement used is FORMAT(3012). The variables to be read in for one run are as follows: TEXT, repeat 3 times NZ NN(i), L(i), 1Q(i), i = 1, NZ

optional, i f N Z > 0

LS(k, i), k= I,3, i= I , N Z + N Z - 1 KN(i), i = 1, NZ JJ(k, i), k = I, 3 NC N l , ( N N I ( i ) , L l ( i ) , I Q l ( i ) , i = 1, N l ) L S l ( k , i). k = 1,3, i = 1,NI + N I - 1 K N I ( i ) , i = 1, NI JJl(k, i), k = 1,3, repeat i = I , N I - 1 JORT, I I, J 1 IORT(i, j ) , j = 1, J l , repeat i = 1,11 MS, MI, ME, MD, KK MET IPR_ D

repeat NC times

TEXT is an information string (72 characters long) for the identification of the case.

250

O. Zatsarinny/ ComputerPhysics Communications98 (1996) 235-254

NZ is the number of common shells, NC is the number of configurations to be considered during one run. The arrays NN, L, IQ, LS, KN, and JJ contain the full description of input configurations and corresponding coupling schemes (see, for details, the description of common blocks / C O R E / and / I N P / ) . Note that the description of the common shells is optional (only for NZ > 0). The input of the additional specification for the radial orbitals, KN, and specific coupling scheme, JJ, is also optional and is driven by the total 'seniority' number LS(1, NZ + NZ - 1). JORT indicates the orthogonality mode. Array IORT(i, j) contains the orthogonality constrains on the radial wave functions (see, for details, the description of common b l o c k / Z O R T / ) . The input of IORT is optional: in the case of full nonorthogonality of the radial wave functions or for the orthogonal orbitals, this array is not required, as well as in the case of the additional specification of the orbitals, which is given in array KN (see the description of the subroutine Pre_ iort). MS, MI, ME, MD, and KK specify the operator to be considered and its order. If all these parameters are equal to zero, then the program generates only the vector-coupling coefficients for the input configurations. MET drives the calculation of the matrix elements between the separate configurations. MET = 0 indicates that the matrix elements are considered between all input states, MET = 1 - - only diagonal matrix elements, MET = 2 - - only nondiagonal matrix elements, MET = 3 - - only transitions from the first states. IPR_D > 0 causes the additional print of the vector-coupling coefficients for each input configuration.

toionizations and Auger decays) the present program is intended.

7. T e s t r u n

References

Example 1. This example presents the dipole transition ls22s22p63s~ ls2s22p63sep, involving the deep shells. Such matrix elements arise, for example, in the relaxed Hartree-Fock calculations of the sodium ls photoionization [25], and present the clear example of full nonorthogonality mode. Namely, for such deep shell transitions (including also the au-

Example 2. In this example we consider four of the configurations analyzed by Hibbert and Froese Fischer [14] in the test run of their program. These configurations are (1) lS2pl(IP)3dl 2p, (2) ls2pl(ap)3dl 2p, (3) lS3PE(Ip)3d2 2p, (4) ls3P2(3p)3d2 2p. The additional subscripts here indicate the set numbers of orbitals. We have computed all the angular coefficients arising in the calculations of the electrostatic interaction matrix elements. This is an example of the partial nonorthogonality mode, when the nonorthogonality conditions are derived from the description of the configuration, on the basis of the set numbers of orbituls. This mode has been specially introduced to cover the requirements of the MCHF calculations. The obtained coefficients agree entirely with the values given by Hibbert and Froese Fischer [ 14]. Example 3. This example includes one matrix element between two arbitrary complex configurations which contain a set of open shells: (2p3(2p)3d9 (2D) 3p IV 12p 5 (2p)3d7 (]D) 3p). To do the calculations in this case, the program has to open and read the data file Det_arc and the auxiliary file Arc_det. Thus, this test run allows user to verify the correct execution of the program in this mode and check the correct structure of the data file Det_ exp.dat. Besides, here the parameter IPR_ D = 1 and the corresponding determinantal expansions are also printed. The times needed for running the above examples on the PC AT 4 8 6 / D X 2 / 6 6 are given in the Test Run Output.

[1] U. Fano, Phys. Rev. 140 (1965) A67. [2] A. Hibbert,Comput.Phys. Commun. 1 (1970) 359; 2 (1971) 80; 7 (1974) 318; 8 (1974) 329. [3] P.G. Burke, Comput. Phys. Commun. 1 (1970) 241. [4] W.D. Robb, Comput. Phys. Commun. 6 (1973) 132. [5] W.D. Klotz, Comput. Phys. Commun.9 (1975) 102. [6] R. Glass, Comput. Phys. Commun.(197) [7] R. Glass and A. Hibbert,Comput. Phys. Commun. 11 (1976) 125.

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

251

[18] J. Olsen, B.O. Roos, P. Jorgensen and H.J. Aa. Jensen, J.

[8] [9] [10] [11]

I.P. Grant, Comput. Phys. Commun. 5 (1973) 273. J.J. Chang, Comput. Phys. Commun. 7 (1974) 225. T. Aberg, Nucl. Instr. Meth. B 87 (1994) 5. Z. Kuplianskis, Izv. Akad. Nauk. SSSR Ser. Fiz. 41 (1977) 2626 [in Russian]. [12] C.A. Nicolaides and G. Aspomallis, J. Phys. B 16 (1983) L251. [13] A. Kupliauskiene, Lithunian J. Phys. 35 (1995) 113. [14] A. Hibbert, C. Froese Fischer and M.R. Godefroid, Comput. Phys. Commun. 51 (1988) 285. [15] A. Hibbert and C. Froese Fischer, Comput. Phys. Commun. 64 (1991) 417. [16] P.-,~. Malmqvist, Int. J. Quantum Chem. 30 (1986) 479. [17] J. Olsen, M.R. Godefroid, P. JSnsson, P.-,g,. Malmqvist and C. Froese Fischer, Phys. Rev. E 52 (1995) 4499.

Chem. Phys. 89 (1988) 2185. [19] E.U. Condon and G.H. Shortley, The Theory of Atomic Spectra (Cambridge University Press, London, 1935). [20] S. Fritzche and I.P. Grant, Physica Scripta 50 (1994) 473. [21] K.G. Dyall, I.P. Grant, C.T. Johnson, F.A. Parpia and E.P Plummer, Comput. Phys. Commun. 55 (1989) 425. [22] P.-O. LSwdin, Phys. Rev. 97 (1955) 1474. [23] H.A. Bethe, Intermediate Q.anmm Mechanics (Benjamin, New York, Amsterdam, 1964). [24] H. Nussbaumer, Comput. Phys. Commun. 1 (1969) 191. [25] B.I. Craig and F.P. Larkins, J. Phys. B 18 (1985) 3713.

Test Run Output

Test 1 dipole operator, full non-orthogonality ....................................................................... List of states: 1

2

2 2 6 1 n 1 q= is / 2s / 2p / 3s / Terms= IS0 IS0 IS0 2SI IS IS 2S 1 2 6 1 1 n 1 q= IS / 2s / 2p / 3s / 9p / Terms= 2SI IS0 1S0 2SI 2PI 2S 2S IS 2P

< state 1 II 0per. II state 2> ...................................................................... KRK= i0

full non-orthogonality

-.707107

[

i/

2]

< is Ir^ll 2p>

-.707107

[

I/

2]

< is Ir'll 2p>

.707107

[

i/

2]

< Is Ir~i[ 9p>

.707107

[

i/

2]

< IS Ir^ll

9p>

.707107

[

i/

21 < 2S Jr^if

2p>

.707107

[

i/

2]

< 2S

-.707107

[

i/

2]

< 2s Ir^ll 9p>

Ir^iI 2p>

-.707107

[

I/

2]

< 2s

-1.414214

[

2/

I]

< 3S Ir^iI 2p>

Ir^ll

1.414214

[

2/

i]

< 3s Ir^ll 9p>

overall time of calculations =

9p>

is, is, 2S, IS, IS, is, 2S, IS, is, is, is, Is, IS, IS, is, is, IS, IS, IS, IS,

2sl^l 2S[ 3SI^I 2Sl 281"1 2sl 3sIAl 2Sl 2Sl^l

281 3SI^I 2Sl 2s]^l 2SJ

3Sl^l

2Sl 2sl^I 2S] 281^I 2sJ

2S, 2S, IS, 2S, 2S, 2S, Is, 2S, IS, 2S, is, 2S, IS, 2S, IS, 2S, Is, 2S, IS, 2S,

3S 3S 2S 3S 3S 3S 2S 3S 3S 3S 2S 3S 3S 3S 2S 3S 2S 3S 2S 3S

^12p2P]^51 ^i 2Di^5 2DI "i 2DI^6 2Dr ^i 2D[^6 2DI

1

1 1 1 n 1 q= is / 2p / 3d / Terms= 2SI 2Pl 2DI IP 2P 1

1

1

2pJ'19p 2p ^i 9p

^1 2p2Pl^Sl I9p2pl^1 ^i ^i ^I ~I ^i

2D]A5 2DI 2Di^6 2DI 2Di^6 2DI 2DI^5 2DI 2Di^6 2Di

.06 sec

Test 2 (comparison with Hibbert and Froese Fischer, CPC 64 (1991) 417 ) ........................................................................ List of states:

I

2p "i 9p

1 2p ~i I 9p

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

252 2

n 1 q= is / 2p / 3d / Terms= 2SI 2PI 2DI 3P 2P 1 1 1 n 1 q= 2s / 3p / 3d / Terms= 2SI 2P1 2DI IP 2P

3

1 1 1 n 1 q= 2s / 3p / 3d / Terms= 2SI 2PI 2DI 3P 2P

4

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

......................................................................

KRK=

3

.173205 -.057735 -.222692

partial non-orthogonality [ 3/ i00] [ i/ 300] [ 243/ 4900]

R2( is, 3d; 3d, is) RI( 2p, 3d; 3d, 2p) R3( 2p, 3d; 3d, 2p)

...................................................................... < state 1 II 0per. II state 3> ...................................................................... KRK=

5

partial non-orthogonality

1.000000

[

1/

i]

1.000000

[

1/

I]

I( is, 2s) I 2Pl ^I 1 3dl'l l 3~I i 3p Ro( is, 2p; 2s, 3p] 1 3di'1

.333333

[

I/

9]

RI( ls, 2p; 3p, 2 s ) I

1.000000

[

1/

1]

so, is

-.lOOOOO

[

1/

IOO]

2o

3d 3di^1

3d l 2pI'i 3p

R2( IS, 3d; 3d, 2s,

2p ^I 3p

< state 1 I I Oper. }I state 4> ...................................................................... KRK=

1

.173205

partial non-or thogonality [

3/

i00]

R2( is, 3d; 3d, 2s)

, 2p,^l

I 13p overall

time of calculations

=

0.33 sec

Test 3 example for the testing file Det_exp.dat ........................................................................ List of states: 1

3 9 n 1 q= 2p / 3d /

O. Zatsarinny/ Computer Physics Communications 98 (1996) 235-254 Terms= 2P1

253

2DI 3P

5 . 7 2 n 1 q= 2p / 3d / Terms= 2Pl 2D3 3P .................................................................. 1

3 9 n 1 q= 2p / 3d / Terms= 2Pl 2D1 3P determinant expansion : N_det=

.54772

nlms= Coef= nlms= Coef=

.22361

3d2PI+0-I sqrt( 3d-2-2P 0sqrt(

6

ML=

3

MS=

3

3d-2+2P 0+ 3d2p-l+2- 3d3d2+0-1 3d 0+ 1 3d-l- I 3d-l+ 1 3d i- 1 3/ I0) 3d-2+2P 0+ 3d2Pi_+[ 3d3d2+0+1 3d-l-I 3d-l+l 3d i- I 3d i+ 1 i/ 20)

38730 o 1 ~ 3~-2-2p o+ 30-2+2p-12~3~~;1 3~30°:I 3~ 0+, 3~-~+, 3~ 1,3~ 1+, Coef= sqrt( 3/ 20) nlms= 3d-2-2P 0+ 3d-2+2p'l+ 3d2P2-I-I 3d3d2+0-1 3d 0+ 1 3d-l+ 1 3d i- I 3d i+ 1 Coef= sqrt( 3/ 20)

.38730

5.772 oloso 2~1-30 i+ 3~2+2~-I"3~2p~:13~3~°:I 3~ 0+, 301-, 3~-i+, 3~ 1, Coef= sqrt( 3/ 10) nlms= 2p-l+3d_2_3d-2+2P i- 3d2Pi_+l 3d3d2+0+1 3d-l-I 3d-l+l 3d 1- I 3d I+ I Coef= sqrt( 1/ 20) .................................................................. .22361

2

5 7 n 1 q= 2p / 3d / Terms= 2P1 2D3 3P determinant expansion:

.23905

.05976 .17928 .05976 .29277 .29277

nlms= 2p P-+1 3d Coef= sqrt( nlms= 2p 0+I 3d I-, Coef= sqrt( nlms= 2p 0+ 1 3d I-, Coef= sqrt( nlms= 2p 0+ 3d-2Coe£= sqrt( nlms= 2p 0+ 1 3d I+ I Coef= sqrt( nlms= 2p 0+ I 3d i+ I

2p-l-3d I÷ 2/ 35) 2p-l3d i+ 1/280) 2p-13d 1+ 9/280) 2p-l3d-2+ 1/280) 2p-13d-2+ 3/ 35) 2p-l3d-2-

N_det =

19

ML=

3

MS=

3

3d-2+2p'l+3d2P2+i-1 2p i+ 1 3d 0- I 3d 0+ I 3d-l- 1 3d-2-2p-l+ 2P3d2+i-1 2p i+ I 3d 0- I 3d 0+ 1 3d-l+ I 2p-l+3d_2+ 3d2P21[1 2p i+ 1 3d 0- I 3d 0+ 1 3d-l+ 1

2p-1+3d 2- 2P3~2÷I1 2p 1+I 31 0-I 31 0+I 31 1+I 2p-l+3d 2- 3d2P2+i-I 2p i+ 1 3d 0- I 3d-l- I 3d-l+ 1 3d2p-l÷2- 3d2P2+i-1 2p i+ 1 3d 0+ I 3d-l- I 3d-l+ 1

O. Zatsarinny / Computer Physics Communications 98 (1996) 235-254

254

Coef= sqrt( nlms= 2p O+ 3d-2Coef= sqrt( nlms= 2p O3d-l+ Coef= sqrt( nlms= 2p O3d i+ Coef= sqrt( nlms= 2p O3d lCoef= sqrt( nlms= 2p O3d lCoef= sqrt( nlms= 2p 03d-2Coef= sqrt( nlms= 2p 03d-2Coef= sqrt( nlms= 2p 03d-l+ Coef= sqrt( nlms= 2p O3d i+ Coef= sqrt( nlms= 2p O3d lCoef= sqrt( nlms= 2p O3d i+ Coef= sqrt( nlms= 2p 03d-2Coef= sqrt( nlms= 2p 03d i+ Coef= sqrt(

.17928 .16903 .16903 -.10351 .10351 -.10351 .10351 .42258 -.08452 -.25355 .33806 .41404 .25355

3/ 35) 2p-l3d-2+ 9/280) 2p O+ 3d-2+ 1/ 35) 2p O+ 3d-21/ 35) 2p O+ I 3d i+ 3/280) 2p O+ I 3d i+ 3/280) 2p O+ 3d-2+ 3/280) 2p O+ 3d-2+ 3/280) 3d

2p-l+ 3d 2-

2p I- I 2p i+ I 3d-l+ 1 3d I- I 3d I+ I 3d 2+

2p-l+

2p

3d 2-

3d 2+

2p-i+

2p

3d-2+

3d 2+

2p-l+ 3d-2-

2p l- I 2p i+ 1 3d 0+ 1 3d-l- I 3d-l+ I 3d 2+

2p-1+

2p

2p-l+ 3d 2-

2p I- I 2p i+ I 3d O+ 1 3d-l- I 3d i+ I 3d 2+

2p-l÷

2p

3d 2-

3d 2+

2p-l3d i+

2p-l+ 3d 2+

2p i+ 1 3d 0- I 3d O+ 1 3d-l- I

2p-l-

2p-1÷

2p l÷l 3d

3d 2-

3d 2+

2p-l3d 2-

2p-l+ 3d 2+

2p i+[ 3d O- I 3d O+ 1 3d-l+ I

2p-1-

2p-l+

2p i÷I

3d 2-

3d 2+

2p-l3d 2-

2p-1+ 3d 2+

2p i+ I 3d O+ 1 3d i- I 3d i+ I

2p-l3d 2-

2p-l+ 3d 2+

2p i+ 1 3d-l- I 3d-l+ I 3d I- I

1-1 2p l÷t 3d o-i 3d 0÷i 1-1 2p l+i

I 2p

1-1

3d O-I 3d O÷i 3d

1+I 3dO+l

3d-l-I

2

-.039678 [

3d 0-I 3d 0+I 3d-i+I

full non-orthogonality

27/17150]

R3( 3d, 3d; 2p, 2p)

overall time of calculations =

10.92 sec

~ ^3 2p

i-I

0-I 3d 0÷I 3d-1-I

< state 1 II Oper. II state 2> ...................................................................... KRK=

1-I

3d-1+I

2p 1+I 3d O+I 3d-1+I 3d

5/ 28) 2p O+ 3d-2+ 1/140) 2p O+ 3d-2+ 9/140) 2p O+ 3d-24/ 35) 2p O+ 3d-2+ 6/ 35) 2p O+ 3d-2+ 9/140)

3d-I-I

3d 3d