A Program Library for Computing Pure Spin–Angular Coefﬁcients for One-and Two-Particle Operators in Relativistic Atomic Theory

: A program library for computing pure spin-angular coefﬁcients for any one-and scalar two-particle operators is presented. The method used is based on the combination of the second quantization and quasi-spin techniques with the angular momentum theory and the method of irreducible tensorial sets. A relativistic approach is assumed. This program library is integrated in the General Relativistic Atomic Structure Package but it can be implemented in other program packages, too.


Introduction
The improved accuracy of modern experiments challenges theorists to match or exceed experimental precision. Models of many-electron atoms and ions require both relativistic and correlation effects to be taken into account; this can be performed, for example, by using various versions of perturbation theory (PT) [1][2][3], the configuration interaction method (CI) [4,5], a combination of the many-body perturbation theory with the configuration interaction method [6,7], the multiconfiguration Hartree-Fock method (MCHF) [5,8] or the multiconfiguration Dirac-Hartree-Fock (MCDHF) method [5,9].
The MCDHF method [10] is probably the most efficient and consistent way to account simultaneously for correlation and relativistic effects in complex energy spectra and in other properties of many-electron atoms with open subshells. However, there are practical and theoretical difficulties related with a large set of configuration state functions (CSF) which need to be taken into account for obtaining accurate results. One of them is related with the integrations over spin-angular variables, which constitute a considerable part of the computation, especially when atoms with many open shells are treated [10]. Therefore, the most efficient approach of spin-angular integration is welcome for implementation in computer codes.
The method [11,12] used in the present program library is based on the combination [13] of the second quantization [14] and quasi-spin techniques [15,16] with the angular momentum theory [17,18] and with the generalized graphical method [19] and the method of irreducible tensorial sets [15,20,21], which has many advantages compared with traditional approaches. The background of this theory [11,12] is presented in Section 2. The structure of the program library, its documentation with a description of main routines, and the library itself are presented in Section 3. This library is integrated in the GRASP-2018 package [22][23][24] for calculating energy structure [5] and such atomic properties as hyperfine structures, transition parameters, and isotope shift [23], but it can be implemented in other program packages, too.

Recoupling Coefficients and Second Quantization
Racah algebra [25][26][27][28] (at the coefficients of fractional parentage level) based on the Fano approach [29] was implemented in atomic structure codes in the seventies, but the calculation of recoupling coefficients with NJSYM [30] remained the bottleneck. The evaluation of matrix elements of complex operators for electronic configurations involving many open subshells was rather time consuming. The performances for calculating recoupling coefficients have been sensitively increased with NJGRAF [31] using the graphical approach to first transform them into structureless graphs from which zero-valued angular momenta are taken out and minimal loops to generate the optimal expression as a sum over products of 6j-coefficients are searched for, but these improvements were not enough. When using Slater determinants, any matrix element of a physical operator can be calculated from the second quantization form of the one-body F = ∑ N i f (i) and two-body G = 1 2 ∑ i =j g(ij) operators [14] where a µ is the electron creation, a † σ is electron annihilation operators, and µ ≡ (n i i ) j i m i , η ≡ (n j j ) j j m j , σ ≡ (n i i ) j i m i , and ≡ (n j j ) j m j . This second quantization formalism [14] was adapted for dealing with symmetry-adapted configuration state functions instead of Slater determinants, leading to a more efficient approach for spin-angular integrations [11,12]. This method uses the coupled tensorial form of the various operators [11], allowing a generalized graphical method [19] based on quasi-spin and on the reduced coefficients of fractional parentage [32,33] and is implemented in the nonrelativistic [34] and relativistic codes [22,35]. This approach is realized in the program library presented in the paper and is described below in more detail.

Quasi-Spin Notation for Antisymmetric Subshell States
In the relativistic theory, each electron shell nl (apart from the ns subshells) is of course known to split into two subshells with j = ± 1/2 = ± , which then affects the representation of the relativistic configuration state function. Instead of the antisymmetric (LS-coupled) shell states |n w αLSJ , as often used in nonrelativistic theory, one then has to deal with the (antisymmetrized product) functions |(n ) j 1 w 1 j 2 w 2 α 1 J 1 α 2 J 2 J . Here, the quantum number still refers to the parity of the configuration state but no longer to the orbital angular momenta of the equivalent electrons as in the nonrelativistic theory.
A (relativistic) subshell state of w equivalent electrons with quasi-spin Q and total angular momentum J is written as [16] |(n ) j w ανJ = |(n ) j w αQJ = (n ) jαQJ; M Q where α refers to all additional quantum numbers that are needed for a unique classification of all subshell states. For any subshell (n j) ≡ (nκ); therefore, the quasi-spin momentum Q helps to encode the seniority quantum number ν by Q = 2j+1 2 − ν /2 , while its z-component characterizes the electron occupation N of the subshell state by M Q = w − 2j+1 2 /2 . The set of two quantum numbers J, Q of a subshell with j ≤ 7/2 and additional α ≡ Nr for j ≥ 9/2 defines the subshell term αQJ, which identifies the state |(n ) j w αQJ constructed with one subshell [33].

Quasi-Spin Formalism and Second Quantization
Although crucial for the most recent code developments, we do not go into much detail in spin-angular algebra. For further explanations, we refer to [11,13,32]. We would like, however, to stress the relation between the fractional parentage coefficients and the completely reduced matrix elements of the creation/annihilation operators appearing in (1) and (2). In the quasi-spin formalism, the operators of second quantization are the components of an irreducible tensor of rank q = 1/2, in a quasi-spin space as well a (q j) In (3), a (j) m j is an electron creation operator, and in the meantime, the tensorã (j) m j is defined as where a † (j) −m j is the electron annihilation operator. Such double tensors a (q j) m q m j are basic elements in modern atomic spectroscopy. Using the Wigner-Eckart theorem in quasi-spin space, we obtain reduced CFPs in quasi-spin space [32] (n ) j w α QJ a where the coefficients Q 1 2 M Q m q | Q 1 2 Q M Q are known as Clebsh-Gordan coefficients and [Q] means (2Q + 1).
The relation between CFPs and the reduced CFPs is [33] So, applying the quasi-spin method, we can use the reduced matrix elements j α QJ a (q j) j α Q J of a (q j) tensor operator, which are independent of the occupation number of the subshell instead of the usual fractional parentage coefficients j w−1 α Q J |} j w αQJ . Thus, an amount of their numerical values is much smaller in comparison with that of the CFP.

Matrix Elements for One-Particle Operator between Complex Configurations
The most simple operator in atomic theory is the one-particle scalar operator in jjcoupling. Let us start the analysis with it.
So, the matrix elements of a one-particle scalar operator F (0) between configuration state functions with any number of open subshells can be expressed as a sum over oneelectron contributions where In jj-coupling, all states are defined. γ α J| and γ β J are, respectively, bra and ket functions with any number of open subshells, n i i ) j i f (0) n j j ) j j is the one-electron reduced matrix element of the operator denote the respective sets of active subshell total angular momenta. Some selection rules for the matrix element of one-particle scalar operators in jj-coupling come from (8). They are presented in Table 1 as the first group of selection rules. Table 1. Selection rules for one-particle scalar operator in jj-coupling.

The Matrix Element Diagonal
Off-Diagonal The second part:

THE FOURTH GROUP OF SELECTION RULES (COMING FROM EFFECTIVE INTERACTION STRENGTH (16))
For example, the matrix element of the Dirac operator H D has (17): The triangular delta δ(J 12 ... a , j a , J 12 ... a ) is nonzero only when the following conditions are fulfilled .. a and J 12 ... a + j a + J 12 ... a is an integer.
The recoupling matrix R j i , j j , Λ bra , Λ ket in (8) is particularly simple. It is either a product of delta functions of the set of intermediate angular momenta 1 (see Equation (18) [11]) or a combination of delta functions and 6j-coefficients (see Equation (22) [11]) when i = j. For example, in the last case (i = j) when i = 1 and j = 2, the recoupling matrix is proportional to In the more complex case, when γ consists of at least 3 subshells (i, j ≥ 3 and j = u and a = min(i, j) and b = max(i, j)), the recoupling matrix is proportional to There are more special cases in matrix element calculation. All these cases can be found in [11]. Some other selection rules come from the recoupling matrix R j i , j j , Λ bra , Λ ket . They are presented in Table 1 as the second group of selection rules.
The operator a (q j) m q in (8) is the second quantization operator in the quasi-spin space of rank q = 1/2 (for more details, see Section 2.3). By applying the Wigner-Eckart theorem in j α j Q j J j we obtain RCFP, for which we can use the tables of RCFP. The submatrix element of the simplest compound tensor operator of type where j αQJ W (kq k j ) j α Q J denotes the reduced matrix element of the tensor operator W (kq k j ) (nj, nj) = a (q j) × a (q j) (k q k j ) in quasi-spin space. In terms of the fully reduced coefficients of fractional parentage j αQJ a (qj) j α Q J , we find This construction has the advantage that the completely reduced matrix elements on the right-hand side of (6) and (12) are independent of the occupation number of the subshell. The last selection rules come from the calculation of the submatrix element of the operator of the second quantization or its combinations. They are presented in Table 1 as the third group of selection rules.
The phase factor ∆ in (8) arises from the reordering needed to match the recoupled creation and annihilation operators in the bra and ket vectors. We have when n i κ i = n j κ j ; otherwise, where w r is the occupation number of subshell r.
The general expression (8) can be used for any scalar one-particle physical operator. It only remains to define the one-electron interaction matrix element (the effective interaction strength) in (8). The only operator required in this implementation is the matrix element of the Dirac operator, a tensor operator of rank zero, Therefore, the matrix element of Dirac operator (7) can be expressed through spin-angular coefficients t αβ ab and radial integrals I(a, b) [5,10] where a ≡ n i i j i and b ≡ n j j j j . The value of pure spin-angular coefficients for one-particle scalar operators is obtained with the help of (8), keeping in mind that the effective interaction strength is equal to one A similar expression as in (7) and (8) has the matrix elements of the one-particle nonscalar operator F (k) between configuration state functions with any number of open subshells where where Γ in R j i , j j , Λ bra , Λ ket , Γ, k refers to the array of all shell terms and intermediate quantum numbers of the bra and ket functions. These recoupling matrices have analytical expressions in terms of just 6j-and 9j-coefficients (see Equations (30) and (34) [36]). Some selection rules for the matrix element of a one-particle nonscalar operator in jj-coupling come from (21). They are presented in Table 2 as the first group of selection rules. Other selection rules come from the recoupling matrix and the tensorial part of the matrix element. The expression (21) is fairly general and covers all cases of the one-particle nonscalar operators in relativistic atomic theory. It only remains to define the value of k rank and the one-electron interaction matrix element (the effective interaction strength) For example, this can be taken from [37] or [38] for relativistic radiative transitions or for hyperfine interaction. Table 2. Selection rules for one-particle nonscalar operator with rank k in jj-coupling.

The Matrix Element Diagonal
Off-Diagonal

THE FIRST GROUP OF SELECTION RULES (COMING FROM (21))
The second part:

Matrix Elements for Scalar Two-Particle Operator between Complex Configurations
According to the approach in [11], a general expression of the submatrix element for any two-particle operator between functions with any number of open subshells can be written as follows: with γ α J G (k j k j 0) n i i j i , n j j j j , n i i j i , n j j j j γ β J = ∑ κ 12 (−1) ∆ Θ n i i j i , n j j j j , n i i j i , n j j j j , Ξ T j i , j j , j i , j j , Λ bra , Λ ket , Ξ, Γ R j i , j j , j i , j j , Λ bra , Λ ket , Γ , where Γ refers to the array of coupling parameters connecting the recoupling matrix R to the submatrix element T, and Ξ refers to the whole array of parameters that connect the amplitude Θ to the submatrix element T. Some selection rules for the matrix element of a two-particle scalar operator in jj-coupling come from (24). They are presented in Table 3 as the first group of selection rules.
To calculate the spin-angular part of a submatrix element of this type, one has to obtain: 1.

4.
Θ n i i j i , n j j j j , n i i j i , n j j j j , Ξ .
Some important points to note are the following: 1. The recoupling matrices R j i , j j , j i , j j , Λ bra , Λ ket , Γ in our approach are much simpler than in other known approaches. We obtained their analytical expressions in terms of just 6j-and 9j-coefficients. That is why we chose a special form of operator in the second quantization, where second quantization operators acting upon the same subshell are tensorially coupled together. Some other selection rules come from this recoupling matrix. They belong to the second group of selection rules (see Table 3).
2. The tensorial part of a two-particle operator is expressed in terms of (products of) operators of the type We denote their submatrix elements by T j i , j j , j i , j j , Λ bra , Λ ket , Ξ, Γ . The parameter Γ represents the whole array of parameters connecting the recoupling matrix It is worth noting that each of the tensorial quantities (25)-(29) act upon one and the same subshell. So, all the advantages of tensor algebra and the quasi-spin formalism may be efficiently exploited in the process of their calculation. We obtain the submatrix elements of operator (25) by using straightforwardly the Wigner-Eckart theorem in quasi-spin space (5).

The Diagonal Matrix Element
operator acts on one subshell (n i i ) j w i i operator acts on two subshells (n i i ) j w i i (n i j ) j w j j for the radial integral R k (aa, aa) 1 ( [5], (89) and (90)) for the radial integral R k (ab, ab) 1 ( [5], (89) and (90)) The second part: where a = min(i, j) and b = max(i, j) additional triangular delta from 6j-coefficients; it depends on the case [11]

THE THIRD GROUP OF SELECTION RULES (COMING FROM
For example, the matrix element of Coulomb operator H C has: δ(l i , k, l i ) and k = even δ(l i , k, l i ) δ(l j , k, l j ) and k = even

THE FIRST GROUP OF SELECTION RULES (COMING FROM (24))
The second part: additional triangular delta from 6j-coefficients; it depends on the case [11]

THE THIRD GROUP OF SELECTION RULES (COMING FROM
THE FOURTH GROUP OF SELECTION RULES (COMING FROM Θ n i i j i , n j j j j , n i i j i , n j j j j , Ξ ) For example, the matrix element of Coulomb operator H C has: δ(j i , k, j j ) δ(l i , k, l j ) and l i + k + l j = even As it is seen, by using this approach, the calculation of spin-angular parts of matrix elements between functions with u open subshells ends up in obtaining the submatrix elements of tensors (25) and (26) within one subshell of equivalent electrons. As these completely reduced (reduced in the quasi-spin, orbital, and spin spaces) submatrix elements do not depend on the occupation number of the subshell, the tables of them are reduced considerably in comparison with the tables of ordinary matrix elements and the tables of fractional parentage coefficients (CFP) [33]. That is why the expressions obtained are very useful in practical calculations. This lets us exploit all the advantages of Racah algebra [11]. The last selection rules (see Table 3) come from the calculation of the submatrix element of the operator of the second quantization or its combinations T j i , j j , j i , j j , Λ bra , Λ ket , Ξ, Γ . They belong to the third group of selection rules.
The amplitude Θ n i i j i , n j j j j , n i i j i , n j j j j , Ξ is proportional to the two-electron submatrix element (the effective interaction strength) of a two-particle operator (2) To obtain the expression of a specific physical operator, analogous to expression (23), the tensorial structure of the operator and the two-electron matrix elements (31) must be known. For example, for a Coulomb operator, the amplitude Θ is proportional to the radial integral where a ≡ n i i j i , b ≡ n j j j j , c ≡ n i i j i , and d ≡ n j j j j . Therefore, the matrix element of the Coulomb operator (23) can be expressed through spin-angular coefficients v αβ abcd; k and radial integrals R k (ab, cd) [5,10] We do not present details on obtaining phase factors ∆ and Θ n i i j i , n j j j j , n i i j i , n j j j j , Ξ , since no essential generalizations may be made here; these are possible only after a particular operator is chosen (for more details see [11,39]). The selection rules, which are dependent on the particular operator, come from Θ n i i j i , n j j j j , n i i j i , n j j j j , Ξ . For example, for the Coulomb operator, they are listed in the fourth group of selection rules in Table 3. For other two-particle physical operators, they are others.

The Method Implementation in Software Packages
The spin-angular integration method [11,12] was implemented in the following software packages: • A general-purpose relativistic atomic structure program (GRASP) [22,35]. The current library (the library presented in this paper) is implemented in GRASP packages. The library version written in FORTRAN 77 programming language is installed in [35] and written in FORTRAN 95 and is installed in [22]. • An MCHF atomic-structure package (ATSP) [34]. The library is written in FORTRAN 77 programming language. It has a similar structure to that described in this paper, but it uses an LS-coupling scheme. Some additional information is published in [40,41]. • A program for relativistic configuration interaction calculations (RELCI) [42]. The library is written in FORTRAN 90/95. It has a similar structure to that described in this paper. Some additional information is published in [43][44][45]. • Jena atomic calculator (JAC) [46]. The library is written in the JULIA programming language.
• A flexible atomic code (FAC) [47]. This program uses an adapted version of the library [42][43][44][45]. • THE RACAH program presents Maple procedures for the coupling of angular momenta [48,49]. The program is written in the MAPLE programming language. This implementation is not suitable for large-scale calculations. It serves for the manipulation of reduced matrix elements and some simple expressions of spin-angular integrations from this theory in both LSand jj-couplings.

•
The HFS program presents Maple procedures as an environment for hyperfine structure parametrization [50]. The program is written in the MAPLE programming language. This implementation is realized in LS coupling and is not suitable for ab initio large-scale calculations.

The Spin-Angular Coefficients for Some Simple Cases and Average Energy of a Configuration
In general, a matrix element of a Dirac-Coulomb Hamiltonian can be expressed through spin-angular coefficients and radial integrals (18) and (33) where a ≡ n i i j i , b ≡ n j j j j , c ≡ n i i j i , and d ≡ n j j j j . The one-body interactions give the spin-angular coefficients t αβ ab and the I(a, b) integrals defined by ( [5], (88)), and the two-body Coulomb interactions give the spin-angular coefficients v αβ abcd; k and the relativistic radial integrals R k (ab, cd). The present program library serves for the calculation of all these spinangular coefficients t αβ ab and v αβ abcd; k for configuration state function with any open subshells using the expressions (7) and (23), respectively. As it is seen from Section 2.4, the analytical expressions for these coefficients are complicated. However, for one subshell case, they have simple expressions [16,[51][52][53]. The spin-angular coefficient t αβ aa is equal to the occupation number for the subshells where α denotes the index of the CSF of bra function and β denotes the index of the CSF of ket function on the left side of (35) and additional quantum numbers of subshell on the right side. Meanwhile, the spin-angular coefficients v αβ aaaa; k can be expressed over the unit tensors [52,53]. However, for k = 0, the spin-angular coefficients v αβ aaaa; 0 depend only on occupation number w i for the subshell with any j and have only the diagonal matrix element v αβ for k > 0 and w i = 2j i + 1 depend on occupation number w i and j i quantum number and for k > 0 and The expression (34) can be rewritten for a diagonal matrix element in respect of configurations with two open subshells as [52,53] where the radial integrals are (39) has the same value as in (35). The spin-angular coefficients v αβ abab; k with w i = 2j i + 1 depend on occupation numbers for subshells with any j and have only the diagonal matrix element and v αβ abba; k with w i = 2j i + 1 have the following expression: where a ≡ n i i j i and b ≡ n j j j j or a ≡ n j j j j and b ≡ n i i j i . The level energies are defined in the relativistic theory by the Coulomb and Breit interactions between and within open subshells. The Breit energy is merely a correction to the electrostatic energy, and therefore is not considered in the average energy of a configuration having a specific electron distribution in the subshells (the one-configuration jj-average energy). Therefore, taking into account the above-presented expressions, the average energy E of a configuration [16,53] is the following: where a ≡ n i i j i and b ≡ n j j j j . The (42) and the orthogonality condition of P nκ and Q nκ lead to relativistic Dirac-Fock equations for finding the one-electron orbitals P nκ and Q nκ [10].

Structure of the Library
A brief description of the algorithm of the library is presented in this section. The program library is divided into three routine groups, METWO, REC, and SQ, according to the peculiarity of the expressions (8), (21), and (24). The purpose of METWO is to calculate the spin-angular coefficients of the matrix elements of any one-and scalar twoparticle operators. It calls REC and SQ along with the subroutines CLRX, CXK, DRACAH, SPEAK, SNRC,TALK, and some modules with the definition of variables and arrays from the GRASP-2018 package [22], which this program library is intended to supplement. The REC routines group calculates the required recoupling coefficients, and SQ is a routines group of standard reduced matrix elements of second-quantized operators. METWO computes the spin-angular coefficients for the Coulomb and Breit interactions. It can be easily extended to calculate the spin-angular coefficients for the second-order effective operator in perturbation theory or for any three-particle operator. All routines use the Fano-Racah phase convention, while GRASP-2018 [22] uses that of Condon-Shortley. The final results are transformed to the Condon-Shortley convention before being output by SPEAK and TALK. It is assumed that all of the first group of angular momentum selection rules (see Tables 1-3) are checked already before calling this library.
The program library uses reduced matrix elements (6) of the operator of second quantization a (q j) instead of coefficients of the fractional parentage. All reduced matrix elements of a (q j) tensor operator and all reduced matrix elements of the tensor operator [ a (q j) × a (q j) ] (k 1 k 2 ) [33] required by METWO and SQ are stored in memory. The routine RMEAJJ extracts the first one from it, and the routine RWJJ extracts the second one.

The METWO Routines Group
This group contains the following routines: ONEPARTICLEJJ1, ONEPARTICLEJJ1, ONESCALAR1, ONESCALAR2, EL1, EL2, EL3, EL4, and EL5, each of which calculates particular spin-angular coefficients of the matrix elements of one-or two-particle operators. The arrays and variables from module m_C of GRASP-2018 [22] library LIBMOD must be set before applying METWO.

The Subroutine ONEPARTICLEJJ1
The purpose of the subroutine ONEPARTICLEJJ1 is to calculate the spin-angular coefficients of the diagonal matrix elements for the one-particle operator in the basis of configuration state functions (21). Its structure is presented in Figure 1. The subroutines RECOP00 and RECOP1 first of all check the selection rules for the recoupling coefficient. If it is nonzero, then the routine RECOP1 computes it (see Section 3.2.7 for more detail), but before it, the routine WJ computes matrix element (91) (see the red frame W in Figure 1, and a detailed description of it is in Section 3.3.3). The subroutine PERKO2 is the interface between GRASP-2018 [22] and the SQ routines group. The subroutine has the following arguments: 1.
NS is a number of peel subshells from the module m_C.

2.
KA is a rank k of the operator (see (21)).

3.
JJA and JJB are the numbers of configuration state functions for the matrix element to be evaluated.

4.
JA is the index in the array JLIST of the orbital on which the creation operator acts.

5.
JB is the index in the array JLIST of the orbital on which the annihilation operator acts. 6.
COEFF is the value of the spin-angular part of the matrix element.

The Subroutine ONEPARTICLEJJ2
The purpose of the subroutine ONEPARTICLEJJ2 is to calculate the spin-angular coefficients of the off-diagonal matrix elements for the one-particle operator on the basis of configuration state functions (21). Its structure is presented in Figure 2. The subroutines RECOP00 and RECOP2 first of all check the selection rules for the recoupling coefficient. If it is nonzero, then the routine RECOP2 computes it (see Section 3.2.8 for more details), but before it, the routine C0T5S generates a Clebsch-Gordan coefficient depending on the occupation number of the subshell, and the routine RMEAJJ computes the reduced matrix element (89) (see Section 3.3.1). The subroutine PERKO2 is the interface between GRASP-2018 [22] and the SQ routines group. The subroutine has the following arguments: 1.
NS is a number of peel subshells from the module m_C.

2.
KA is a rank k of the operator (see (21)).

3.
JA is the index in the array JLIST of the orbital on which the creation operator acts.

4.
JB is the index in the array JLIST of the orbital on which the annihilation operator acts.

5.
COEFF is the value of the spin-angular part of the matrix element.

The Subroutines ONESCALAR1 and ONESCALAR2
The subroutines ONESCALAR1 and ONESCALAR2 have a similar structure as ONEPARTICLEJJ1 and ONEPARTICLEJJ2, respectively. They calculate the spin-angular coefficients of the matrix elements for the one-particle scalar operator in the basis of configuration state functions (8).

The Subroutine EL1
The purpose of the subroutine EL1 is to calculate the spin-angular coefficients of the diagonal matrix elements for the two-particle operator on the basis of configuration state functions (24). Its structure is presented in Figure 3. The subroutines RECO and RECO2 first of all check the selection rules for the recoupling coefficient. If it is nonzero, then the routine RECO2 computes it (see Section 3.2.4 for more details). The subroutine PERKO2 is the interface between GRASP-2018 [22] and the SQ routines group. The subroutines ITREXG and IXJTIK organize the calculation. The subroutine SIXJ calculates 6j-coefficient. The subroutines SPEAK and TALK prepare the output of spin-angular coefficients for Coulomb and Breit interactions, respectively.  Figure 3. Structure of the subroutine EL1. The box "GRASP" routines correspond to the routines ITRIG, SNRC, CXK, SPEAK, and TALK from the GRASP92 package [54]. The blue frame WW corresponds to the routine WW1 (computing the matrix element (94)) and shows its structure. Additionally, the red frame W shows the structure of the routine WJ (computing the matrix element (91)).
WJ1 calls the subroutines in the red box W (in red frame) in Figure 3; a detailed description appears in Section 3.3.3.
The subroutine EL1 has the following arguments:

1.
JJA and JJB are the numbers of configuration state functions for the matrix element to be evaluated.

2.
JA and JB locate the position of the two interacting orbitals in the array JLIST [22] from the module m_C of GRASP-2018 [22] library LIBMOD. These parameters are used for the combinations (43)-(47).

3.
IIRE must be set equal to 0 when the matrix element is diagonal, 1 when the matrix element is off-diagonal with respect to configuration state functions.

The Subroutine EL2
This calculates the spin-angular coefficients of the off-diagonal matrix elements of the two-particle operator (24) in the case when the bra and ket configurations have a pair of interacting shells whose occupation numbers differ by two a (j 1 ) (n 1 l 1 ) a (j 1 ) (n 1 l 1 )ã (j 2 ) (n 2 l 2 )ã (j 2 ) (n 2 l 2 ).
The structure of this subroutine is the same as that of EL1, omitting the subroutine WW1. The subroutine has the following arguments: 1.
JJA and JJB are the numbers of configuration state functions for the matrix element to be evaluated.

2.
JA is the index in the array JLIST of the orbital on which the two creation operators act.

3.
JB is the index in the array JLIST of the orbital on which the two annihilation operators act.

The Subroutine EL3
EL3 calculates the spin-angular coefficients of the off-diagonal matrix elements of the two-particle operator (24). The subshell occupation numbers of the bra configuration differ from those of the ket configuration by at most one. There are two interacting subshells on each side. The structure of this subroutine is presented in Figure 4.    The subroutine has a similar list of arguments as the routine EL2, but instead of the two arguments JA and JB it has four arguments, JA, JB, JC, and JD, which point to the interacting subshells of the matrix element of interest in the JLIST array.
The structure of EL31 is very similar to that of EL1. The recoupling coefficient is the same as in (44)- (47), and the only additional subroutines are the C0T5S, which generates a Clebsch-Gordan coefficient depending on the occupation number of the second subshell, and the RMEAJJ from the routine GG1222, which provides the reduced matrix element of the a (q j) tensor operator. The subroutine AWP1 handles the more difficult calculation for the spinangular part of the first shell. The structure of this subroutine is defined by the AW box (in red frame) in Figure 5; more details are given in Section 3.3.5.  Figure 5. Structure of the subroutine EL31. The box "GRASP" routines corresponds to the routines SNRC, CXK, SPEAK, and TALK from the GRASP92 package [54]. The red frame AW shows the structure of the routine AWP1 (computing the matrix element (92)) to which the subroutine EL31 explicitly is calling.
The organization of EL32 is shown in Figure 6. The subroutine WAP1 organizes the more complicated calculation for the second subshell in this case. The structure of this subroutine is defined by the WA box (in red frame) in Figure 6; more details are given in Section 3.  Figure 6. Structure of the subroutine EL32. The box "GRASP" routines correspond to the routines SNRC, CXK, SPEAK, and TALK from the GRASP92 package [54]. The red frame WA shows the structure of the routine WAP1 (computing the matrix element (93)) to which the subroutine EL32 explicitly is calling.
The routine RECO3 calculates the recoupling coefficient; more details are given in Section 3.2.5. The subroutine JFAZE determines the phase factor arising from operator permutations. The routine EILE reorders the subshells JA, JB, JC and places the ordered pointers in JAA, JBB, and JCC. The routine GG1233 organizes the calculation of the spinangular coefficients. C0T5S determines two Clebsch-Gordan coefficients: one dependent on the occupation of the first subshell, the other on the occupation of the second subshell. RMEAJJ provides two reduced matrix elements of the a (q j) tensor operator, and WJ1 determines the spin-angular factors for the third shell (see Section 3.3.3).

The Subroutine EL4
The subroutine EL4 organizes the calculation of the spin-angular coefficients of offdiagonal matrix elements of the two-particle operator (24) involving three interacting subshells. The operator combinations are a (j 3 ) (n 3 l 3 ) a (j 3 ) (n 3 l 3 )ã (j 1 ) (n 1 l 1 )ã (j 2 ) (n 2 l 2 ); (57) and It uses the subroutine EL41, which has the same structure as EL33 (Figure 7). The subroutine EL4 has the arguments JJA and JJB, which define the numbers of configuration state functions for which the matrix element is to be evaluated, four arguments JA, JB, JC, and JD, which locate the interacting subshells in the array JLIST, and the argument ICOLBREI, which determines the calculation of spin-angular coefficients.

The Subroutine EL5
The EL5 calculates spin-angular coefficients of the off-diagonal matrix elements of the two-particle operator (24) when there are four interacting subshells. The structure of this routine is presented in Figure 8. The subroutine EL5 has the arguments JJA and JJB, which define the numbers of configuration state functions for which the matrix element is to be evaluated, four arguments JA, JB, JC, and JD, which locate the interacting subshells in the array JLIST, and the argument ICOLBREI, which determines calculation of spin-angular coefficients.

The REC Routines Group
This contains the following subroutines RECOONESCALAR, RECOP00, RECO, RECO2, REC3, RECO4, RECOP1, and RECOP2, all of which are concerned with the calculation of recoupling coefficients of one-or two-particle operators. The first checks that all of the second group of angular momentum selection rules coming from the recoupling coefficients R j i , j j , Λ bra , and Λ ket or R j i , j j , j i , j j , Λ bra , Λ ket , and Γ (see Tables 1-3) are satisfied and, if so, computes the recoupling coefficients for one, two, three or four subshells, respectively. The arrays and variables from the module m_C of GRASP2018 [22] library LIBMOD must be defined before calling REC.

The Subroutine RECOONESCALAR
This routine checks the first part of the second group of angular momentum selection rules coming from the recoupling coefficients R j i , j j , Λ bra , Λ ket , and Γ (see Table 1) for the one-particle scalar operator.
The subroutine RECOONESCALAR has the following arguments: 1.
NS is the number of peel subshells from the module m_C for the diagonal case and NS = -1 for the off-diagonal matrix element.

2.
JA1 identifies the first subshell, on which the creation operator a (j) or annihilation tensorã (j) acts in the array JLIST.

3.
JA2 identifies the second subshell, on which the creation operator a (j) or annihilation tensorã (j) acts in the array JLIST. The subshells must be numbered so that the arguments JA1 and JA2 are in an increasing order.

4.
KA is the parameter which determines the number of subshells coupled by the interaction, taking the values KA = 0 if one subshell, KA = 1 if two subshells. 5.
The subroutine returns the value of IAT, which is 0 if the selection rules are not satisfied and the recoupling coefficient is zero, 1 if the recoupling coefficient is to be calculated.

The Subroutine RECOP00
where A (k) and B (k) are simple or composite tensor operators of rank k. A (k) acts only on the first active shell and B (k) on the second active shell in the order in which they are coupled in the configuration. The structure of the routine is presented in Figure 10. In most cases, the recoupling coefficients factorize in three parts which are generated by the subroutines DIAGA1, DIAGA2, and DIAGA3. When the argument IRE = 0, all triads of the recoupling coefficient are checked using the routine IXJTIK and perform the actual calculation when IRE = 1. The necessary 6j-coefficients are generated using SIXJ. . The green frame A1 shows the structure of the routine DIAGA1, the red frame A2 shows the structure of the routine DIAGA2, and the blue frame A3 shows the structure of the routine DIAGA3 to which the subroutine RECO2 explicitly is calling.
The subroutine RECO2 has the following arguments: 1.
JA1 identifies the first subshell on which the operator A (k) acts in the array JLIST.

2.
JA2 identifies the second subshell on which the operator B (k) acts in the array JLIST.

3.
KA is the intermediate rank k.

4.
IRE takes the input value 0 if only the coupling triads are to be checked 1 if the recoupling coefficient is to be calculated.

5.
When IRE = 0, the subroutine returns the value of IAT, which is 0 if the selection rules are not satisfied and the recoupling coefficient is zero, 1 if the recoupling coefficient is to be calculated. 6.
REC is the value of the recoupling coefficient computed when IRE = 1.

The Subroutine REC3
This subroutine checks the second part of the second group of angular momentum selection rules coming from R j i , j j , j i , j j , Λ bra , Λ ket , and Γ or calculates the recoupling coefficients for the two-particle scalar operator acting in three different subshells, the tensorial structure of which can be represented as .
As in (85), A (k 1 ) , B (k 2 ) , and C (k) are simple or composite tensor operators which act on subshells i, j, and m, respectively.
The structure of the routine is presented in Figure 11. The recoupling coefficients can be factorized into several simple parts, which are generated by the subroutines DIAGA1, DIAGA2, DIAGA3, and DIAGA4. REC3 checks all triads of the recoupling coefficient using IXJTIK if IRE = 0 and performs the calculation if IRE = 1. The subroutine has the following arguments: 1. JA1, JA2, and JA3, which point to the orbitals i, j, and m in the array JLIST.

3.
IRE takes the input value 0 if only the coupling triads are to be checked 1 if the recoupling coefficient is to be calculated.

4.
When IRE = 0, the subroutine returns the value of IAT, which is 0 if the selection rules are not satisfied and the recoupling coefficient is zero, 1 if the recoupling coefficient is to be calculated.

5.
REC is the value of the recoupling coefficient computed when IRE = 1.

The Subroutine RECO4
This subroutine checks the second part of the second group of angular momentum selection rules coming from R j i , j j , j i , j j , Λ bra , Λ ket , and Γ or calculates the recoupling coefficients for the two-particle scalar operator acting in four different subshells, the tensorial structure of which can be represented as where A (k 1 ) , B (k 2 ) , C (k 3 ) , and D (k 4 ) may be simple or composite tensor operators of the orders indicated, corresponding to the structure of (87). The subshells must be ordered so that A (k 1 ) operates on the first and D (k 4 ) on the last in order. The structure of the subroutine RECO4 is the same as that of RECO3 ( Figure 11). The subroutine has the following arguments: 1. JA1, JA2, JA3, and JA4, which point to the orbitals 1, 2, 3, and 4 in the array JLIST.

3.
IRE takes the input value 0 if only the coupling triads are to be checked 1 if the recoupling coefficient is to be calculated.

4.
When IRE = 0, the subroutine returns the value of IAT, which is 0 if the selection rules are not satisfied and the recoupling coefficient is zero, 1 if the recoupling coefficient is to be calculated.

5.
REC is the value of the recoupling coefficient computed when IRE = 1.

The Subroutine RECOP1
This subroutine checks the second part of the second group of angular momentum selection rules coming from R j i , j j , Λ bra , and Λ ket (see Table 2)) or calculates the recoupling coefficients for the one-particle nonscalar operator acting in one subshell. The structure of the subroutine RECOP1 is similar to that of RECO2 (Figure 10), except that DIAGA5 is used instead of DIAGA2.
The subroutine RECOP1 has the following arguments: 1.
NS is the number of peel subshells from the module m_C.

2.
JA1 identifies the subshell on which the operator a acts in the array JLIST.

3.
KA is the intermediate rank k.

4.
IRE takes the input value 0 if only the coupling triads are to be checked 1 if the recoupling coefficient is to be calculated.

5.
When IRE = 0, the subroutine returns the value of IAT, which is 0 if the selection rules are not satisfied and the recoupling coefficient is zero, 1 if the recoupling coefficient is to be calculated. 6.
RECC is the value of the recoupling coefficient computed when IRE = 1.

The Subroutine RECOP2
This subroutine checks the second part of the second group of angular momentum selection rules coming from R j i , j j , Λ bra , and Λ ket (see Table 2) or calculates the recoupling coefficients for the one-particle nonscalar operator acting in two different subshells, the tensorial structure of which can be represented as where A (k) and B (k) are simply the creation operator a (j) or annihilation tensorã (j) . A (k) acts only on the first active subshell and B (k) on the second active subshell in the order in which they are coupled in the configuration. The structure of the subroutine RECOP2 is similar to that of RECO3 (Figure 11), except that DIAGA5 is used instead of DIAGA2.
The subroutine has the following arguments: 1.
NS is the number of peel subshells from the module m_C.

2.
JA1 and JA2, which point to the orbitals 1 and 2 in the array JLIST.

4.
IRE takes the input value 0 if only the coupling triads are to be checked 1 if the recoupling coefficient is to be calculated.

5.
When IRE = 0, the subroutine returns the value of IAT, which is 0 if the selection rules are not satisfied and the recoupling coefficient is zero, 1 if the recoupling coefficient is to be calculated. 6.
RECC is the value of the recoupling coefficient computed when IRE = 1.

The SQ Routines Group
The Section SQ (standard quantities) is a collection of utilities used by the routines groups METWO and REC. Most of them are independent and may be used in other programs. The routines check the third group of angular momentum selection rules (see Tables 1-3) by mainly calling the routines ITJJ, IXJTIK, C0T5S, C1E0SM, C1E1SM, and CLE0SM and, if so, calculates the matrix or reduced matrix elements of standard quantities.
Most of the subroutines use arrays from the module CONS_C of GRASP2018 [22] library LIBMOD. Single subshell data needed for calculation are stored in the two arrays I and B. The former consists of: Table 4. Allowed couplings |(n ) j w ανJ of states for j = 1/2 − 9/2, which are in the array MT from the module MTJJ_C. The subshell quasi-spin angular momentum, seniority of the coupling, and the subshell angular momentum are denoted by Nr, Q, ν, and J, respectively (see [33]).

No.
Nr

The Subroutine RMEAJJ
This subroutine determines the value of the reduced matrix elements of operator (25) j α QJ a (q j) j α Q J .
The subroutine uses the table of reduced matrix elements of the a (q j) tensor operator for j = 1/2 − 9/2, which are stored in memory. The subroutine has the following arguments: 1.
LL is the quantum number j multiplied by two.

2.
IT is the state number of the bra function (see column No. in Table 4).

3.
LQ is the quasi-spin Q for the bra function multiplied by two.

4.
J is the total angular momentum J for the bra function multiplied by two.

5.
ITS is the state number of the ket function. 6.
LQS is the quasi-spin Q for the ket function multiplied by two. 7.
J1S is the total angular momentum J for the ket function multiplied by two. 8.
COEF is the value of the reduced matrix element (89) which is returned by the subroutine.

The Subroutine RWJJ
The subroutine determines the value of the reduced matrix element of operator (26) The subroutine uses the table of reduced matrix elements of the tensor operator a (q j) × a (q j) (k 1 k 2 ) stored in memory for j = 1/2 − 7/2. The subroutine has the following arguments:

1.
J is the quantum number j multiplied by two.

2.
J1 is the state number of the bra function.

3.
J2 is the state number of the ket function.

6.
COEF is the value of the reduced matrix element (90) which is returned by the subroutine.

The Subroutine WJ1
Before determining the value of the matrix element of operator (26) this subroutine checks the main part of the third group of angular momentum selection rules coming from T j i , j j , j i , j j , Λ bra , Λ ket , Ξ, Γ for j = 1/2 − 37/2 (see Table 3). The subroutines calls W1JJG for the calculation of the matrix element (91) in cases j = 9/2 − 37/2. The subroutine finds the Clebsch-Gordan coefficient, which gives the dependence on the subshell occupation number in cases j = 1/2 − 7/2. If the tensor product (91) consists of either two electron creation operators or two annihilation operators, then C1E1SM is called. Otherwise, C1E0SM is called. The subroutine RWJJ finds the reduced matrix elements of the operator a (q j) × a (q j) (k 1 k 2 ) . The structure of the subroutine WJ1 is given in the red box labeled W in Figure 3.
The subroutine has the following arguments:

8.
AW is the value of the reduced matrix element (92) which is returned by the subroutine.

The Subroutine AWP1JJG
The subroutine determines the value of the matrix elements (92) for j = 9/2 − 37/2. The subroutine has the following arguments: 3.3.9. The Subroutine WW1 Before determining the value of matrix elements of operator (29) (n ) j w α QJ a the subroutine checks the main part of the third group of angular momentum selection rules coming from T j i , j j , j i , j j , Λ bra , Λ ket , Ξ, Γ for j = 1/2 − 37/2 (see Table 3). The subroutine WW1 uses the subroutines ITJJ, RUMT, and IZAS1 for calculation of this sort of matrix element. The subroutine WJ1 calculates the first and the second parts of the operator to be calculated. The structure of the subroutine WW1 is presented in Figure 3 (blue block WW).
The subroutine has the following arguments: 1.
IK is the array I for the bra function.

2.
BK is the array B for the bra function.

3.
ID is the array I for the ket function.

4.
BD is the array B for the ket function.
WW is the value of the reduced matrix element (94) which is returned by the subroutine.

Description of the New Modules for Arrays Used in the Program Library
The principal new developed modules for arrays used in the program library are listed in the Table 5. They are located in the GRASP-2018 [22] library LIBMOD. Table 5. Arrays used in this program library for computing pure spin-angular coefficients for oneand two-particle operators in the relativistic atomic theory.

Name Dimension Function mtjj_C
The arrays for the atomic states of |(n ) j w ανJ with any occupation of subshells MT 63 the arrays for j = 1/2 − 9/2 subshells (see Table 4) mtjj2_C The arrays for the atomic states of |(n ) j w ανJ in case w = 1, 2 MT9 6 The arrays for j = 9/2 subshell MT11 189 The arrays for j = 11/2 − 37/2 subshells trk_C The data of the orbitals for interacting subshells BD1 3 The array B (see Section 3.3) for the first subshell of the ket function BD2 3 The array B for the second subshell of the ket function BK1 3 The array B for the first subshell of the bra function BK2 3 The array B for the second subshell of the bra function ID1 7 The array I (see Section 3.3) for the first subshell of the ket function ID2 7 The array I for the second subshell of the ket function IK1 7 The array I for the first subshell of the bra function IK2 7 The array I for the second subshell of the bra function BD3 3 The array B (see Section 3.3) for the third subshell of the ket function BD4 3 The array B for the fourth subshell of the ket function BK3 3 The array B for the third subshell of the bra function BK4 3 The array B for the fourth subshell of the bra function ID3 7 The array I (see Section 3.3) for the third subshell of the ket function ID4 7 The array I for the fourth subshell of the ket function Let us first consider simple cases with a small number of CSFs (3SD, 3SDT, and 4SD with J = 0, 1, 2, 3, 4, 5, 6,7,8,9). Although GRASP(NEW) [22,35] with the present spinangular library generates the full set of "pure" coefficients for both one-and two-particle operators, the calculations run 1.4-2.3 times faster than in equivalent calculations with GRASP92 [54] because of the small number of computational operations. Table 6 compares the performance of the GRASP92 and GRASP(NEW) codes for the larger-scale 4SD and 5SD problems. The number of CSFs is presented in column two, and column three lists the number of nontrivial t αβ ab coefficients; this number is the same for both calculations. In addition, the GRASP(NEW) code calculates pure spin-angular coefficients, which are fully sufficient to find the v αβ abcd; k of any physical operator. Therefore, the number of pure spin-angular coefficients calculated by GRASP92 and GRASP(NEW) is different. Columns four and five of Table 6 list the number of v αβ abcd; k coefficients from GRASP92 [54] and GRASP92(NEW) [22,35] calculations, respectively. Column "with a different number of v αβ abcd; k " presents the speed-up of GRASP(NEW) in case it calculates the different number of v αβ abcd; k (see column four) as GRASP92 (see column five). The fact that GRASP(NEW) calculates approximately twice the number of spin-angular coefficients as GRASP92 therefore increases the efficiency per coefficient by a factor of two. The actual speed-up of the GRASP(NEW) program is listed in the column "actual". From the results presented in Table 6, we conclude that the new program is not much faster for simple cases but shows better performance for more complicated cases.

Limitations of the Program Library
All orbitals in a wave function expansion are assumed to be orthonormal; meanwhile, the nonorthogonal orbitals are not supported by the program library.
Configuration state functions with any distribution of electrons in subshells with j ≤ 9/2 are allowed in this library. For all subshells with j ≥ 11/2 (i.e., h−, h+, i−, i+, . . . electrons), the maximum number of equivalent electrons is restricted to two. This permits the user to take into account the single, double, triple, and quadruple excitations from open d and f shells for the systematic MCDHF or RCI studies of heavy and superheavy elements (Z > 95). Other limitations can come from the package in which the library is installed.

Conclusions
The approach to matrix element evaluation presented in the paper, and realized in the program library, is based on the combination [13] of the angular momentum theory [17], as described in [18], on the concept of irreducible tensorial sets [14,15,20]), on a generalized graphical approach [19], on the second quantization in coupled tensorial form [15], on the quasi-spin approach [16], and on the use of reduced coefficients of fractional parentage [33]. It introduces a number of new features in comparison with traditional approaches:

1.
A number of theoretical methods known in atomic physics facilitate the treatment of spin-angular parts of matrix elements, among them are the theory of angular momentum, its graphical representation, the quasi-spin, and the second quantization in its coupled tensorial form. However, while treating the matrix elements of physical operators in general, including the off-diagonal ones, with respect to configurations, the above-mentioned methods are usually applied either only partly and/or inefficiently. An idea to combine all of these methods in order to optimize the treatment of the general form of matrix elements of physical operators in atomic spectroscopy is presented and carried out in this work. It allows us to investigate even the most complex cases of atoms and ions efficiently in relativistic approaches.

2.
The general tensorial expressions of one-and two-particle operators, presented in this work, allow one to exploit all the advantages of tensorial algebra. In particular, this is not only a reformulation of spin-angular calculations in terms of standard quantities but also the prior determination from symmetry properties in which matrix elements are equal to zero without performing further explicit calculations (the first and third grout of selection rules (see Tables 1-3)).

3.
The tensorial forms of a one-and two-particle operator (see (8) and (24)), allow one to obtain simple expressions for the recoupling matrices. Hence, the computer code based on this approach would immediately use the analytical formulas for the recoupling matrices R j i , j j , j i , j j , Λ bra , Λ ket , Γ . Among the rest, this feature saves computing time because (i) complex calculations lead finally to simple analytical expressions [11,12], and (ii) a number of momenta triads (triangular conditions (the second group of selection rules (see Tables 1-3)) can be checked before the explicit calculation of a recoupling matrix leading to a zero value. These triangular conditions may be determined not only for the terms of subshells that the operators of the second quantization act upon but also for the rest of the subshells and resulting terms.

4.
The tensorial form of any operator presented as products of tensors of the second quantizaqtion a , and a (q j) allows one to exploit all the advantages of a new version of Racah algebra based on quasi-spin formalism. So, the application of the Wigner-Eckart theorem in quasispin space for the submatrix element of the operator of the second quantization or its combinations provides an opportunity to use the tables of reduced coefficients of fractional parentage and tables of the other standard quantities, which do not depend on the occupation number of a subshell of equivalent electrons. Thus, the volume of tables of standard quantities is reduced considerably in comparison with the analogous tables of submatrix elements of T (k) (for jj-coupling) [16,56] and the tables of coefficients of fractional parentage. These tables cover all the electronic configurations needed in practice. Therefore, the process of selecting the standard quantities from the tables becomes simpler. It also allows to determine the third group of selection rules (see Tables 1-3). 5.
In this approach, which is both diagonal and off-diagonal with respect to configurations, matrix elements are considered in a uniform way and are expressed in terms of the same quantities, namely, reduced coefficients of fractional parentage or reduced submatrix elements of standard tensors, which are independent of the number of electrons in a subshell. The difference is only in the values of the projections of the quasi-spin momenta of separate subshells. The complete numerical tables of these quantities allow practical studies of any atom or ion in the periodical table.
Funding: No funding for this research was granted.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, [G.G.], upon reasonable request.