Scalable Codes for Precision Calculations of Properties of Complex Atomic Systems

Abstract

High precision atomic data are indispensable for studies of fundamental symmetries, tests of fundamental physics postulates, developments of atomic clocks, ultracold atom experiments, astrophysics, plasma science, and many other fields of research. We have developed a new parallel atomic structure code package that enables computations that were not previously possible due to system complexity. This code package also allows much quicker computations to be run with higher accuracy for simple systems. We explored different methods of load-balancing matrix element calculations for many-electron systems, which are very difficult due to the intrinsic nature of the computational methods used to calculate them. Furthermore, dynamic memory allocation and MPI parallelization have been implemented to optimize and accelerate the computations. We have achieved near-perfect linear scalability and efficiency with the number of processors used for calculation, paving the way towards the future where most open-shell systems will finally be able to be treated with good accuracy. We present several examples illustrating new capabilities of the newly developed codes, specifically correlating up to all 60 electrons in the highly charged Ir${}^{17+}$ ion and predicting certain properties of Fe${}^{16+}$.

1. Introduction

Studies of the fundamental symmetries with atoms and ions require knowledge of atomic properties as well as calculations to extract potential new physics from the experiments [1]. For example, the most accurate atomic calculation for a heavy atom was carried out to study Cs parity violation [2,3]. A precision (0.35% accuracy) experiment [4] had to be supplemented by theory computation of a comparable accuracy in order to extract the value of the weak charge and test the standard model of elementary particles. The theory accuracy is more than an order of magnitude worse in Yb, not sufficient for a similar analysis of the Yb experiments [5,6], motivating studies with multiple isotopes [6]. While Yb has two valence electrons, low-lying excitations from the closed shells cannot so far be accurately treated by any high-precision methods with significantly lower attainable precision. Theory accuracy is even worse for more compacted systems such as Dy [7] or Sm of interest to parity violation and other studies. Studies of CP-violation with atoms also require computations of the enhancement factors [8]. Studies of Lorentz symmetry violations need computations of the Lorentz invariance violating matrix elements [9]. Moreover, the development of experiments with new systems requires accurate knowledge of many atomic properties. For example, development of even more precise atomic clocks [10] that are used to test the Einstein equivalence principle, search for the variation of the fundamental constants and dark matter [1], is accelerated by strong theory-experimental collaborations.

In the present work, we improve on the capabilities of an ab initio atomic structure code package for calculating atomic properties of complex many-electron systems. The methods used here, including configuration interaction (CI) and the combination of CI with either many body perturbation theory (CI+MBPT) [11,12] or the linearized coupled cluster method (CI+all-order) [13], are very broadly applicable to any atom in the periodic table.

A common problem of the CI method is that the number of configurations for a many-electron atom grows exponentially on the number of correlated electrons, and for systems with more than 4–5 electrons, it is possible to include only a small fraction of them. A standard approach is to divide all electrons into core and valence electrons. In the framework of the CI method, only valence–valence interactions are treated explicitly, while core excitations are neglected. This essentially limits the accuracy of the method when it is applied to a heavy many-electron atom or ion. A way to circumvent this problem was suggested by Dzuba, Flambaum, and Kozlov in Ref. [11]. They used an idea that a CI method treats the interactions between valence electrons, while MBPT can be effectively used to account for core–core and core–valence correlations. Combining these two methods, it is possible to acquire benefits from both approaches and attain higher computational accuracy. A similar approach, combining CI and MBPT, was suggested by Savukov and Johnson [14]. Another program package based on the CI+MBPT method discussed in Ref. [11] was developed by Dzuba and co-authors (see Refs. [15,16] and the references therein). The version of the CI+MBPT code package applied by us was modified for public use and published in Computer Physics Communications in 2015 by Kozlov et al. [12]. Although not published yet, the inclusion of the all-order method which provides accurate solutions for a large number of properties of atoms and ions with up to 5–6 valence electrons has been completed and made fully compatible with this code package.

Numerous problems that were not tractable before have been solved with our newly developed parallel codes, including problems in astrophysics, metrology of highly charged ions, neutral atoms, and negative ions. Our focus will be on the computational developments that enable these new large-scale calculations. In particular, we implemented a message passing interface (MPI) to parallelize computationally expensive portions of the programs. This allows us to fully take advantage of modern high performance computing facilities, such as the University of Delaware high performance Caviness community cluster, where we developed and tested our programs. We use the new package of parallel programs to consider two cases of particular experimental significance: Ir${}^{17+}$, which was proposed for the development of novel atomic clocks, and Fe${}^{16+}$, which has lines essential for plasma diagnostics tools for astrophysics. We found great success in the parallelization efforts, as the new programs enabled precision calculations of these systems beyond what was desired.

The organization of the rest of this paper is as follows: in Section 2, a summary of atomic structure theory and the methods used in this work is presented. In Section 3, we describe main programs of the CI/CI+MBPT/CI+all-order code package and the computational developments done in this work. In Section 4, the application of the newly developed code package to highly charged Ir${}^{17+}$ and Fe${}^{16+}$ ions is discussed. Finally, we conclude and describe future prospects of this work in Section 5.

2. Theory and Methods

In this section, we describe the methods used in our programs and applications. For any many-electron system, we can divide all electrons into core and valence electrons. In this way, we can separate the electron–electron correlation problem into one describing the valence–valence correlations under the frozen-core approximation, and another describing the core–core and core–valence correlations. In the initial approximation, we start from the solution of the restricted Dirac–Hartree–Fock (HFD) equations in the central field approximation to construct one-electron orbitals for the core and valence electrons. Virtual orbitals can be constructed from B-splines or by other means to account for correlations. The valence–valence correlation problem is solved using the CI method, while core–core and core–valence correlations are included using either MBPT or the all-order method. In either case (CI+MBPT/CI+all-order), we form an effective Hamiltonian in the valence CI space, then diagonalize the effective Hamiltonian using the CI method to find energies and wave functions for the low-lying states.

2.1. The CI Method

The CI method is a standard ab initio method for calculating atomic properties of many-electron systems. In the valence space, the CI wave function is constructed as a linear combination of all distinct states of a specified angular momentum J and parity

$$\psi =\sum _i{c}_i{\Phi }_i,$$

(1)

where the set $\left\{{\Phi }_i\right\}$ are Slater determinants generated by exciting electrons from some reference configurations to higher orbitals.

Varying the coefficients $c_i$ results in a generalized eigenvalue problem

$$\sum _j\langle {\Phi }_i\left|H\right|{\Phi }_j\rangle c_j=E{c}_i,$$

(2)

which can be written in matrix form and diagonalized to find the energies and wave functions of the low-lying states. The energy matrix of the CI method can be obtained as a projection of the exact Hamiltonian H onto the CI subspace $H^{\mathrm{CI}}$ [11]

$$H^{\mathrm{CI}}=E_{\mathrm{core}}+\sum _{i>N_{\mathrm{core}}}{h}_i^{\mathrm{CI}}+\sum _{j>i>N_{\mathrm{core}}}{V}_{ij},$$

(3)

where $E_{\mathrm{core}}$ is the energy of the frozen core, $N_{\mathrm{core}}$ is the number of core electrons, $h_i^{\mathrm{CI}}$ accounts for the kinetic energy of the valence electrons and their interaction with the central field, and $V_{ij}$ accounts for the valence–valence correlations.

Having obtained from CI the many-electron states $|JM\rangle$ and $|J^{\prime }{M}^{\prime }\rangle$ with the total angular momenta $J,J^{\prime }$ and their projections $M,M^{\prime }$, one can form density transition matrix in terms of the one-electron states $|nljm\rangle$ [12]

$$\widehat{\rho }={\rho }_{nljm,n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }}|nljm\rangle \langle n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }|,$$

(4)

where

$${\rho }_{nljm,n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }}=\langle J^{\prime }{M}^{\prime }|a_{n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }}^{†}{a}_{nljm}|JM\rangle .$$

(5)

Here, un-primed indices refer to the initial state and primed indices refer to the final state. The many-electron matrix element can then be written as

$$\langle J^{\prime }{M}^{\prime }|T_q^L|JM\rangle =\mathrm{Tr}{\rho }_{nljm,n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }}\langle n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }|T_q^L|nljm\rangle ,$$

(6)

where the trace sums over all quantum numbers $\left(nljm\right)$ and $\left(n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }\right)$, and $T_q^L$ is the spherical component of the tensor operator of rank L. Using the Wigner–Eckart theorem, one can reduce the many-electron matrix element to

$$\langle J^{\prime }\parallel T^L\parallel J\rangle =\mathrm{Tr}{\rho }_{nlj,n^{\prime }{l}^{\prime }{j}^{\prime }}^L\langle n^{\prime }{l}^{\prime }{j}^{\prime }\parallel T^L\parallel nlj\rangle ,$$

(7)

where  

$${\rho }_{nlj,n^{\prime }{l}^{\prime }{j}^{\prime }}^L={(-1)}^{J^{\prime }-M^{\prime }}{\left(\begin{array}{ccc}{J}^{\prime }& L& J\\ -M^{\prime }& q& M\end{array}\right)}^{-1}\sum _{m{m}^{\prime }}{(-1)}^{j^{\prime }-m^{\prime }}\left(\begin{array}{ccc}{j}^{\prime }& L& j\\ -m^{\prime }& q& m\end{array}\right){\rho }_{nljm,n^{\prime }{l}^{\prime }{j}^{\prime }{m}^{\prime }}.$$

(8)

We have developed new parallel programs based on these methods: conf realizes the CI method, which forms the CI Hamiltonian and uses Davidson’s algorithm of diagonalization [17] to find low-lying energies and wave functions; dtm calculates reduced matrix elements (7) of one-electron operators by forming the reduced density transition matrices (8). We discuss the challenges of developing the parallel programs in Section 3.

2.2. Selecting Important Configurations with Valence Perturbation Theory (CI+PT)

As the number of configurations contributing to the CI wave function grows exponentially with the number of valence electrons, efficient selection of the most important configurations from a set of configurations becomes the main challenge of accurate computations. To significantly reduce the number of configurations, we further developed a method suggested in Ref. [18] to predict important configurations based on a set of configurations with known weights. This method can be used to optimize the CI space by identifying the most important configurations from a list of CI configurations using perturbation theory (PT) [18].

All second-order corrections are taken into account and added to the energy calculated from CI to obtain the total energy, $E^{\mathrm{CI}}=E_0+E_1$, while first-order corrections to the wave functions are stored for use in subsequent CI calculations. This process of using CI on a small subspace, calculating corrections via PT, and reordering the list of configurations in descending weights, can be repeated several times to form the most optimal CI subspace. Once the energy differences between subsequent CI calculations are relatively small, it can be assumed that convergence has been met.

We have developed a new parallel program conf_pt that realizes the CI+PT method. The parallel version enables computations of extremely large problems, with tests running up to 400 million determinants.

2.3. Including Core Correlations with Other Methods (CI+MBPT/CI+All-Order Method)

The CI+MBPT [11,13,19] and CI+all-order [13] methods include core–core and core–valence correlations using the combination of CI with second-order MBPT and the linearized coupled cluster method, respectively. The many-electron Schrödinger equation can be written as

$$H_{\mathrm{eff}}\Psi =E\Psi ,$$

(9)

where the effective Hamiltonian has the form

$$H_{\mathrm{eff}}=H^{\mathrm{CI}}+\Sigma .$$

(10)

Here, $H^{\mathrm{CI}}$ is the CI Hamiltonian described by Equation (3), and $\Sigma ={\Sigma }_1+{\Sigma }_2$ is the core–valence correlation potential obtained from either MBPT or the all-order method, where ${\Sigma }_1$ and ${\Sigma }_2$ are the one- and two-electron parts of the core–valence correlation potential, respectively. Following the methods described in Section 2.1, Equation (9) is solved by diagonalizing the effective Hamiltonian using the CI method to obtain energies and wave functions of the low-lying states, then the wave functions can be used to calculate matrix elements of desired one-electron operators.

3. Computational Developments

A schematic of the main programs in the CI/CI+MBPT/CI+all-order code package is illustrated in Figure 1. First, the program hfd solves the restricted HFD equations to construct one-electron orbitals, and the program bass constructs the basis set required by all subsequent codes. If including core correlations, either the all-order block of codes (allcore-ci, valsd-ci, sdvw-ci and second-ci) or the second-order MBPT code (second-ci) are used to calculate the corrections to the Hamiltonian. The set of configurations are generated using the program add and then radial integrals are calculated with the program basc. The main program conf, which realizes the CI method as discussed in Section 2.1, is then ran to calculate wave functions and energies. The CI+PT code conf_pt can optionally be utilized here to optimize the set of configurations or improve the results of CI. Matrix elements of desired observables can then be calculated using the program dtm and polarizabilities using the program pol-ci.

Of the programs listed in the scheme in Figure 1, the CI program (conf) and the matrix elements program (dtm) have been parallelized with the message passing interface (MPI) [20,21]. The CI+PT program (conf_pt) is an auxiliary program that has also been parallelized with MPI. There are several key features of our recently developed parallel codes over previous serial (non-parallel) programs published in Ref. [12]. For all parallel programs, we focused on improving readability and usability, reducing the total memory footprint, implementing dynamic memory allocation where possible, and parallelizing the computationally expensive portions of code using MPI.

The newly developed parallel programs enable computations that were not possible before due to lack of memory or prohibitive computational times. Calculations for smaller and simpler systems can now be done with much more accuracy in a shorter amount of time, allowing a broader range of effects to be investigated in a short period of time. With the original serial version of the programs, calculations of properties of very small systems of 1–2 valence electrons could last up to a few days, while calculations for more complex systems with 3–6 valence electrons could last up to weeks. Some problems involving systems with more than six valence electrons were simply intractable with this code due to the large amount of memory required or extremely long computation times. The main objectives of the new developments were to optimize memory usage in the programs while speeding up computational times by implementing MPI.

3.1. Parallelization of Codes

Our main focus will be on the developments of the parallel CI, CI+PT, and matrix element programs. Each of these methods requires calculating matrix elements between a pair of determinants, which is very difficult to parallelize due to the nature of the problem. In the CI program, the CI Hamiltonian matrix, as well as the matrix of the operator ${\mathit{J}}^{\mathbf{2}}$, are constructed by calculating these matrix elements. Similarly, the CI and PT blocks of the Hamiltonian matrix in the CI+PT program, and the density or transition matrix in the matrix element program, are constructed in the same manner. Construction of each of these matrices is essentially a different variant of the same problem.

We will describe only the parallel implementation to the CI program, but the methods discussed here are applicable in general to the other parallel programs. The main difficulty here is the intrinsic unpredictability of the computational method, i.e., one cannot guess how many non-zero matrix elements there are without explicitly comparing all pairs of determinants. When forming the matrices and calculating the matrix elements, one has to compare electron occupancies between determinants for each configuration. The Slater–Condon rules [22,23] describe the number of operations that are done for each matrix element, which depends on the number of differences between the pair of determinants forming the matrix element.

The construction of the matrix follows a two-loop structure: an outer loop iterates over the total number of determinants and an inner loop iterates over the total number of configurations. Determinants belonging to each configuration are paired and the value of the matrix element is calculated if there are less than 3 differences between the pair of determinants. We separate the computation into two stages: (i) a comparison stage to find pairs of determinants with less than three differences between them and (ii) a calculation stage to calculate the value of the matrix element. After all determinants have been compared, one final loop is done through the total number of non-zero matrix elements to calculate their values. With this implementation, we have achieved near-perfect linear scalability and efficiency in our tests with up to 500 computing cores. It is also possible, and much simpler, to distribute the workload by dividing the outer loop by the number of cores. However, this results in very uneven load-balancing which is exacerbated for larger computations. This is not the case when dealing with smaller matrices such as in the case of the density and transition matrices.

In the serial implementation, non-zero matrix elements were written to disk due to lack of memory of the computers used when the codes were first developed. However, with modern computers, we are able to store the non-zero matrix elements in memory, which also has the advantage of reducing the computation time due to the removal of file input and output (I/O) to disk. The serial program wrote each matrix element to disk with 24 bytes: eight bytes for the counter, four bytes for each index, and eight bytes for the value of the non-zero matrix element. The parallel program removes the redundant eight bytes for the id, reducing the memory requirement of storing the Hamiltonian matrix by 33%.

Since the matrix elements of the CI Hamiltonian are distributed across cores, the Davidson procedure had to be modified to take advantage of the parallelization. The Davidson procedure diagonalizes the CI Hamiltonian matrix to find energies and wave functions of the low-lying states. In the present CI code, only the calculation of matrix-vector products has been parallelized due to it being the only computationally expensive portion of the algorithm. Since each core holds a portion of the Hamiltonian matrix in memory, the Davidson algorithm had to be modified so each core would only be responsible for computing products of their stored matrix elements with the wave function.

3.2. Speedup of Parallel Programs

The performance of the parallel conf program was tested using the Ir${}^{17+}$ ion. We describe Ir${}^{17+}$ as the ion with 30 valence electrons in the open $4f$, $4d$, and $4p$ shells with an [$8spdfg$] basis set. The basis set is designated by the highest principal quantum number for each partial wave included. For example, the [$8spdfg$] basis set includes $1-8s$, $2-8p$, $3-8d$, $4-8f$, and $5-8g$ orbitals. These test calculations included 24,895 relativistic even-parity configurations and 17,431,323 determinants. As reference, the total computation time of the original serial calculation was just short of two weeks.

In

Table 1. The runtime and speedup of the subroutines FormH and Dvdsn of the parallel CI program (conf) are presented for increasing number of cores N, relative to 50 cores (n = 50). The runtimes presented are given in seconds (s). These large test runs were done with Ir${}^{17+}$ with 24,895 relativistic configurations and $17.4\times {10}^6$ determinants. Note that the total times include all serial subroutines outside of FormH and Dvdsn.

	Runtime (s)			Speedup (from N = 50)
N	FormH	Dvdsn	Total	FormH	Dvdsn	Total
50	22,571	1343	24,218	1	1	1
100	12,843	1514	14,593	1.8	0.9	1.7
200	5800	957	6927	3.9	1.4	3.5
300	3810	678	4610	5.9	2.0	5.3
400	2913	596	3646	7.8	2.3	6.6
500	2292	535	2958	9.9	2.5	8.2

, we compare runtimes of the formation of the CI Hamiltonian matrix (FormH) and the Davidson iterative procedure (Dvdsn), as well as the total computational time, for increasing number of computing cores. Comparing a large number of cores to the base N = 50 case, we see that there is a near-perfect linear scalability up to 500 cores for the FormH subroutine, but very minimal performance gain for Dvdsn from 200 to 500 cores. In Figure 2, we plot the speedup of the main subroutines as a function of the total number of cores, with respect to the ideal speedup of 100% efficiency.

The implementation of the Davidson procedure in the parallel CI program does not scale well since only the calculation of matrix–vector products was parallelized, leaving a large majority of the algorithm serial. Since the Davidson procedure typically does not run as long as the formation of the Hamiltonian, the performance of the Davidson procedure was deemed sufficient for our problems, leaving remaining optimizations to a future project. Despite the inefficiency of Dvdsn, we can see its relatively small impact on the total speedup of the program in Figure 2.

The efficiency of the total speedup attained with the parallel CI program is about 70–80% with the number of cores. The speedup achieved with the parallel CI+PT program and the parallel matrix element program is very similar, with the parallel CI+PT program averaging around 75–85% efficiency and the parallel matrix element program averaging around 80–90% efficiency with the number of cores. These parallel programs have been routinely used in our group to calculate various atomic properties of many-electron systems.

3.3. Selection of Important Configurations

In addition to the parallel codes, we have developed algorithms and auxiliary codes for significantly reducing the number of configurations by efficiently identifying dominant configurations from a list of configurations. It is necessary to be able to select the most important configurations for the valence CI space, especially when the dimensionality of the CI problem becomes intractably large.

The first algorithm involves using the CI+PT method to predict important configurations based on a set of already selected configurations with known weights, as discussed in Section 2.2. Initially, a configuration list with no weights is constructed by allowing excitations from a few basic configurations. A small subspace is chosen from the top of the list to generate initial weights for each of the chosen configurations in the list using CI, then weights of all other configurations in the list are calculated using PT. We then reorder the initial configuration list by descending weights, then repeat this process until convergence is met, i.e., the CI space has been saturated as it has taken into account the most important configurations. The convergence is checked by comparing energies computed from consecutive CI and PT procedures.

The second algorithm explores the importance of the configurations based on the correlations between electrons that are present in the configurations. We find that contributions of electrons in orbitals of the same partial wave are very regular. For example, if our calculations of Fe${}^{16+}$ show that the contribution of the $1s^{2}2s2p^{5}6s^2$ configuration is negligible, then all similar configurations where the two last electrons have higher principal quantum number, e.g., $1s^{2}2s2p^{5}6s7s$, $1s^{2}2s2p^{5}6s8s$, $1s^{2}2s2p^{5}7s^2$, can be omitted. However, if the $1s^{2}2s2p^{6}6f$ configuration has a large weight, then other configurations with $7f$ and $8f$ electrons should be included.

4. Applications

4.1. Optical Clocks Based on Highly Charged Ions

Highly charged ions (HCI) such as Ir${}^{17+}$ are of particular interest for the development of novel atomic clocks due to its very high sensitivity to the variation of the fine structure constant $\alpha$ and related dark matter searches [24,25,26,27]. There are many advantages to creating optical clocks with HCI, such as enhanced sensitivity to the variation of $\alpha$, heavily suppressed systematic effects, and estimated potential clock uncertainty beyond the ${10}^{-18}$ limit [28,29,30]. Recent developments in quantum logic techniques for HCI spectroscopy have made rapid progress in the development of HCI possible [31].

In the case of Ir${}^{17+}$, theoretical calculations are particularly difficult due to atomic configurations with holes in the $4f$ shell leading to a very large uncertainty in prior theoretical results, with over 50% for the predicted clock transition energy [32,33]. The clock transitions were expected to be observed in recent experiments done in the last six years. However, no clock transitions or any $E1$ transitions have been found, even though the predicted transition rates were well within experimental capabilities. Our new results explain the lack of observations of the $E1$ transitions and provide a pathway towards detection of clock transitions based on highly charged ions. Here, we will only summarize the results of our application of the new parallel code to Ir${}^{17+}$. A detailed discussion of the computation can be found in Ref. [34].

We find that the best initial approximation is achieved by solving restricted DHF equations with partially filled shells, namely $[1s^2\cdots 4d^{10}]4f^{13}5s$. The hybrid CI+MBPT/CI+all-order approaches described in Section 2.3 can not be used with such an initial approximation. Therefore, we treat the inner shells with the CI method, leading to an exponential increase in the number of configurations. We found that, while the weights of most configurations are small, the number of important configurations are still very large.

We start with the most straightforward CI computation that includes single and double excitations from the $4f$ and $5s$ valence shells, similar to Ref [26]. The excitations are allowed to each of the basis set orbitals up to $\left[7spdfg\right]$. Next, we allow excitations of any of the 24 electrons from the $4d^{10}4f^{13}5s$ shells to the same basis set orbitals. We find drastic changes in the energies of the $E1$ transitions when excitations are allowed from the $4d$ shell. Due to such large contributions, we continued to include more and more electrons of the inner shells into the CI valence space, until all 60 electrons have been included. Both single and double excitations were allowed from the $4f,4d,4p,4s$ and $3d$ shells, while only single excitations are included for all other shells. We tested that the double excitation contribution is small for these inner shells, allowing us to omit them at the present level of accuracy. The results, obtained with different numbers of shells included in the CI valence space, are given in

Table 2. Energies of Ir${}^{17+}$ (in cm${}^{-1}$) obtained using CI with increasing number of open shells. Only configurations obtained by exciting $4f$ and $5s$ electrons are included in the “$5s4f$ only” column. Contributions from exciting electrons from the $4d$ shell are given in the column labeled “$4d$ contr.”, and contributions of all other shells are given separately in the next columns. The results with all 60 electrons correlated by the CI are listed in the column “All shells open”. Sum of all other corrections is given in the column labeled “Other”—see Ref. [34] for detailed explanation.

Configuration		$5\mathit{s}4\mathit{f}$ Only	$4\mathit{d}$	$4\mathit{p}$	$4\mathit{s}$	$3\mathit{d}$	$1\mathit{s}2\mathit{s}3\mathit{s}$	$3\mathit{p}$	$2\mathit{p}$	All Shells	Other	Final
			Contr.	Contr.	Contr.	Contr.	Contr.	Contr.	Contr.	Open
$4f^{13}5s$	${}^3{F}_4^o$	0	0	0	0	0	0	0	0	0	0	0
	${}^3{F}_3^o$	4714	31	15	14	8	−3	2	0	4781	−4	4777
	${}^3{F}_2^o$	25,170	−75	14	13	75	−2	25	−4	25,220	−34	25,186
$4f^{14}$	${}^1{S}_0$	9073	5797	−931	−1994	1097	−240	−54	9	12,757	−375	12,382
$4f^{12}5s^2$	${}^3{H}_6$	36,362	−8549	460	1848	−403	183	294	144	30,339	−56	30,283
	${}^3{F}_4$	46,303	−8680	−5	1858	−410	184	251	144	39,645	−81	39,564
	${}^3{H}_5$	59,883	−8638	454	1858	−326	183	324	143	53,882	−84	53,798
	${}^3{F}_2$	68,786	−8751	−188	1690	−384	253	191	64	61,662	−233	61,429

. Contributions from increased size of the basis set, triple excitations, full Breit, and QED, were found to be small at the present level of accuracy. We found that by far the largest effect comes from the inclusion of the inner electron shells into the CI. Three basis sets of increasing sizes $\left[7spdfg\right]$, $\left[8spdfg\right]$, and $\left[10spdfg\right]$ were used to test basis set convergence.

While allowing excitations from the $4d$ shell drastically reduced the energies of the $E1$ transitions, allowing excitations from the $4p$ shell drastically reduced several $E1$ matrix elements and rates to well below the detector ability. We have identified several other transitions at different wavelengths for the future $E1$ transition searches. As soon as any $E1$ transition is detected, we will be able to obtain a better prediction for the proposed clock transition wavelength.

The largest Ir${}^{17+}$ run with the latest version of the parallel CI program included 96,622 configurations and 58,224,918 determinants, while the largest number of determinants included in the old serial CI program was about 5,000,000. The largest run calculated and stored over 100 billion Hamiltonian matrix elements, with the whole run requiring a total of 2880 GB of memory, which was the maximum amount available to us at the time. Distributed across 80 cores, this run took two days and 18 h to complete. A larger number of cores were not used due to the large amount of memory required to store the basis set, the Hamiltonian matrix, and subsequent arrays for the Davidson procedure.

4.2. Calculation of the 3C/3D Line Intensity Ratio in Fe XVII

Some of the brightest lines of the spectra of many hot astrophysical objects arise from Fe${}^{16+}$ around 15 Å, namely the resonance line $3C$$({\left[{\left(2p^5\right)}_{1/2}3d_{3/2}\right]}_{J=1}\to {\left[2p^6\right]}_{J=0})$ and the intercombination line $3D$$({\left[{\left(2p^5\right)}_{3/2}3d_{5/2}\right]}_{J=1}\to {\left[2p^6\right]}_{J=0})$. They are crucial for plasma diagnostics of electron temperatures, elemental abundances, ionization conditions, velocity turbulences, and opacities. For the last four decades, there has been a persistent discrepancy between the observed intensity ratios and advanced plasma models, diminishing the utility of high-resolution X-ray observations. A recent experiment measured the most accurate $3C/3D$ oscillator strength to date, in an attempt to explain this puzzle [35].

We carried out a precision calculation with our newly developed MPI version of the CI program, allowing us to drastically increase the number of included configurations to over 230,000. This calculation correlated all 10 electrons, included full Breit and QED corrections, and predicted the transition rates with 1–2% accuracy. Our calculations ruled out incomplete inclusion of the electronic correlations in theoretical calculations as the potential explanation of the puzzle. A detailed study of the latest experiment and theoretical work done can be found in Ref. [35].

We started with all possible single and double excitations from the $2s^{2}2p^6$, $2s^{2}2p^{5}3p$ even and $2s^{2}2p^{5}3s$, $2s^{2}2p^{5}3d$, $2s2p^{6}3p$, $2s^{2}2p^{5}4d$, $2s^{2}2p^{5}5d$ odd configurations, correlating eight electrons with a short $\left[5spdf6g\right]$ basis set. Separate calculations were done to establish the effects of triple and quadruple excitations, as well as full correlation with the $1s^2$ shell. We found negligible corrections to both energies and matrix elements as illustrated by

Table 3. Contributions to the energies (in cm${}^{-1}$) of Fe${}^{16+}$ calculated with increased basis set sizes and number of configurations. Contributions from triple excitations, excitations from the $1s^2$ shell, and QED contributions are given in their respective columns. The last line gives the $3C-3D$ energy difference in eV.

Configuration		$\left[5\mathbf{spdf}6\mathit{g}\right]$	Triples	$1{\mathit{s}}^2$	$+\left[12\mathbf{spdfg}\right]$	$+\left[17\mathbf{dfg}\right]$	QED	Final	Diff. [36]
$2p^6$	${}^1{S}_0$	0	0	0	0	0	0	0
$2p^{5}3p$	${}^3{S}_1$	6,087,185	6	254	3876	772	67	6,092,159	0.02%
$2p^{5}3p$	${}^3{D}_2$	6,116,210	−21	24	2886	701	43	6,119,842	0.03%
$2p^{5}3p$	${}^3{D}_3$	6,129,041	−23	25	3015	711	94	6,132,864	0.03%
$2p^{5}3p$	${}^1{P}_1$	6,138,383	−11	41	2825	704	82	6,142,025	0.03%
$2p^{5}3s$	2	5,842,248	−10	108	3408	735	787	5,847,276	0.03%
$2p^{5}3s$	1	5,857,770	−10	70	3303	708	784	5,862,626	0.03%
$2p^{5}3s$	1	5,953,697	−10	74	3364	717	1042	5,958,883	0.03%
$2p^{5}3d$	${}^3{P}_1$	6,466,575	−11	16	2384	665	87	6,469,717	0.03%
$2p^{5}3d$	${}^3{P}_2$	6,481,385	−13	16	2250	658	86	6,484,383	0.03%
$2p^{5}3d$	${}^3{F}_4$	6,482,549	−12	27	1745	622	97	6,485,028	0.03%
$2p^{5}3d$	${}^3{F}_3$	6,488,573	−14	26	1740	607	84	6,491,016	0.03%
$2p^{5}3d$	${}^1{D}_2$	6,502,481	−17	21	1696	627	88	6,504,895	0.03%
$2p^{5}3d$	${}^3{D}_3$	6,511,163	−18	18	1762	604	87	6,513,617	0.02%
$2p^{5}3d$	${}^3{D}_1$	6,548,550	−16	−3	1747	616	134	6551029	0.02%
$2p^{5}3d$	${}^3{F}_2$	6,589,977	−16	22	1729	629	335	6,592,676	0.02%
$2p^{5}3d$	${}^3{D}_2$	6,596,316	−17	14	1947	641	334	6,599,235	0.03%
$2p^{5}3d$	${}^1{F}_3$	6,600,744	−17	19	1803	610	343	6,603,501	0.03%
$2p^{5}3d$	${}^1{P}_1$	6,656,872	−8	−52	1743	619	288	6,659,462	0.02%
$3C-3D$	(eV)	13.4302	0.0009	−0.0061	−0.0005	0.0004	0.0191	13.4440	0.15%

. The differences between our theoretical energies and experimental values were found to be less than those of the NIST database by 3000 cm${}^{-1}$. The energies from the revised analysis of the Fe${}^{16+}$ spectra given in Table 3 are estimated to be accurate to about 90 cm${}^{-1}$. The last line of Table 3 shows the difference of the $3C$ and $3D$ energies in eV, with the final value 13.44(5) eV.

4.3. Other Applications

The methods and codes developed in this work have also been utilized in various other calculations involving HCI, neutral atoms and negative ions. Recent works include exploring the development of optical clocks using highly charged Cf${}^{15+}$ and Cf${}^{17+}$ ions [37], predicting the lowest-lying odd-parity energy levels in neutral actinium [38], predicting quasibound states of negatively charged La${}^{-}$[39], and calculating $E2$ transition properties of Bi${}^{-}$ [40].

5. Conclusions and Further Developments

We have developed a new parallel atomic structure code package that has opened a lot of new possibilities for high-precision calculations of atomic properties of various complex many-electron systems. The efficiency and accuracy of our programs have been validated by solving many problems, involving astrophysics, in metrology of highly charged ions, and negative ions.

Although our parallel programs have displayed great success in efficiency and accuracy, there is still plenty of developmental work to be done. This includes optimization of current parallel algorithms, parallelization of remaining problematic serial codes (e.g., the Davidson procedure in the CI program), optimization of memory usage, completion of documentation, including check-pointing, making the codes user-friendly, and more.

The programs developed here are part of a larger project developing an online portal for high-precision atomic data and computation. The portal will feature a database of all high-precision data calculated in the last couple of decades, as well as wave functions of a large number of atoms and ions stored on a back-end server. Users will be able to request transition properties and polarizabilities for systems not in the database without having to download any codes or learn any inputs. Requested data will be generated automatically through code executions on a back-end server and then updated on the database. All codes will also be released to the public, optimized and user-friendly.