Fock space coupled cluster study of the 11Π g state of the Li2 molecule

ABSTRACT Recently reported (M. Musiał, J. Chem. Phys., 136, 134111 (2012)) multireference coupled cluster method formulated in the (2,0) sector of the Fock space proved to be an efficient scheme in the treatment of potential energy curves (PECs) for alkali metal diatomics Me2. The (2,0) sector provides description of states obtained by an attachment of two electrons to the reference which for the Me2 is a doubly ionised Me+ 2 2 system. The latter has a very concrete advantage in the calculations of the PECs since it dissociates into closed shell fragments (Me+ 2 2 → Me++Me+); hence, the restricted Hartree–Fock reference can be used in the whole range of interatomic distances. In this work, an accurate PEC and vibrational energy levels are obtained for the state of the lithium dimer. The average deviation of the vibrational term values from experiment is 0.79 wavenumber and is significantly lower than that provided by other theoretical works. GRAPHICAL ABSTRACT


Introduction
In several recent papers [1][2][3][4][5][6][7][8], we studied an applicability of the coupled cluster (CC) theory [9][10][11][12][13] to the description of the dissociation process of the alkali metal diatomics. The correct reproduction of the potential energy curves (PECs) for the standard single bond formed by a pair of electrons is not a trivial task. The difficulty originates from the fact that in the vicinity of the equilibrium the molecule is correctly represented by a closed shell configuration easily obtained with the restricted Hartree-Fock (RHF) scheme. However, when the closed shell molecule dissociates into open shell fragments, the RHF leads to incorrect dissociation results and in such case the proper computational variant is the unrestricted HF (UHF). The latter, though, has some inconvenient features like poorly determined spin or CONTACT  certain convergence problems, in particular when correlated calculations are involved. The ideal and very welcome situation occurs when the closed shell system dissociates into closed shell fragments. Such a situation arises when we dissociate a double positive ion instead of the neutral molecule. Let us consider as an example the Li 2 molecule which dissociates into the open shell atoms each with an unpaired electron: However, the double positive ion of the Li 2 dissociates as and we obtain closed shell products, Li + , isoelectronic with the helium atom.
By focusing on the double cation rather than neutral molecule as a reference system, we may face a problem of inadequate form of the one-particle functions. A certain remedy for this is a utilisation of the large basis sets and exhaustive model spaces. In the current study, the size of the (2,0) model space is equal to 7396 which produces a large number of S amplitudes both in the (1,0) and (2,0) sectors, sufficient to ensure necessary relaxation effects in the correlated wave function.
In the current study, we focus on the accurate calculations of the lowest energy singlet g state, for which experimental data are precisely measured. The computations are based on the above mentioned strategy, i.e. in the first step, we calculate the RHF reference function for the Li +2 2 system, and in the next, the FS-MRCC(2,0) scheme is used to recover electronic states of the neutral Li 2 molecule. An important observation is that the standard FS calculations are plagued with the intruder state problem which prohibits using large basis sets. However, the FS-MRCC(2,0) version reported in [1] is based on IH formalism in which the iterative solution of the FS-MRCC equations is replaced with the diagonalisation of the IH matrix. Due to this, a relatively large basis set including 114 functions per atom could be used.
In the next sections, we present the synopsis of the theory and later on we discuss characteristics of the 1 1 g state, i.e. spectroscopic parameters, energies of the vibrational levels and inertial rotation constants. In the supplemental material, we include the set of the centrifugal distortion coefficients.

Synopsis of the theory
The CC [9][10][11][12][13] method is based on the exponential parameterisation of the wave function: where T is a cluster operator being a sum of operators responsible for the excitations. For the CCSD method: singles (T 1 -S) and doubles (T 2 -D). The o is a reference function, i.e. a Slater determinant constructed from the RHF orbitals. The first step in the calculations is solving the RHF equations and then the CCSD ones for the reference system. In our case -Li +2 2 . In the multireference approaches [28], the eigenvalue equation for the effective Hamiltonian operator, H eff , is solved: where and P is a projection operator and is the valence universal wave operator: where S is a FS wave operator defined in a universal manner for all sectors. Thus, the exact eigenvalues are obtained by diagonalisation of the H eff operator within the model space. The model space in the FS formalism [28][29][30][31][32][33][34][35][36] is the configurational subspace of the FS obtained as a direct sum of the particular sectors (see Ref. [1] for details). The active orbitals form the one-electron valence space composed of two parts: particle subspace and hole subspace. The model space is indicated by giving the top sector (k, l) and the size of the active-particle and active-hole spaces.
The double electron attached states correspond to the (2,0) sector. In this case, the active space includes only particle levels hence providing this number is enough to define the active space. In the (2,0) sector, the model space is spanned by the configurations αβ formed by the distributing the two added electrons in all possible ways among the active-virtual levels. The hierarchical nature of the FS solutions requires that in order to solve the FS equations for the (k, l) sector, all the solutions for the lower rank sectors (i, j) with i ࣘ k and j ࣘ l must be known. In the current case, this requirement enforces solving the FS equations for the (0,0) and (1,0) sectors. The (0,0) sector corresponds to the single-reference solution for the reference system (Li +2 2 ion). In the (1,0) sector, the model space consists of single-electron attached configurations, α , with the additional electron placed on the active-particle levels.
The energy values of the double electron attached states are obtained by the diagonalisation of the effective Hamiltonian within the αβ configurational space: The S (0, 0) ( ࣕ T), S (1,0) and S (2,0) are cluster operators for sector indicated by the superscripts. These operators have a cluster structure. In our model, they are limited to singles and doubles: The S (2,0) 1 operator does not contribute since it cannot be defined within the DEA formalism.
The solutions for the S (1, 0) sector can be obtained using EA-EOM-CC formalism [37,38], i.e. the EOM scheme applied to the single electron attached states [39,40]. This is known that the energy values obtained with the EA-EOM-CC scheme are identical to those of FS-CC for (1,0) sector while the corresponding eigenvectors differ only by straightforward transformation. Once the S amplitudes for the (0,0) and (1,0) sectors are found (using above scheme in order to omit convergence problems in one-valence sector), we can solve the FS equations for S (2,0) and after getting these amplitudes we construct the H eff operator the diagonalisation of which provides the sought energies of the double electron attached states.
It should be mentioned that in order to avoid convergence problems (connected to intruder states) in the (2,0) sector in the practical realisation, we applied intermediate Hamiltonian (IH) variant of the FS-CCSD method [39][40][41][42][43][44][45] described in detail in our paper devoted to the FS-CCSD (2,0) method [1]. This approach reduces the diagonalisation of theH operator within the subspace spanned by the ab configurations obtained by all possible distributions of two attached electrons among virtual levels. TheH is obtained by adding so-called dressing [1] to the standardH = e −T He T operator, employed in the EOM Table . Potential energies, E(R), at different internuclear separations, R, for the   g state of Li  at the FS-CCSD (,) method using ANO-RCC+ basis set.
theory. To obtain singlet solutions, required in the current study, the spin free equations need to be solved, see Table 1 of [14] (first line). This formalism allows to use large model spaces and basis set. In summary, we should mention one important feature of this approach, i.e. its rigorous size-extensivity which is particularly useful in the studies of the dissociation process. Moreover, a scaling is no worse than n 5 for the IH-FS-CCSD model for the DEA and EA parts (without counting the active orbitals whose number is limited). The underlying CCSD ground state scales as n 6 . Thus, the IH-FS-CCSD (2,0) approach offers a robust first principle size-extensive method which provides excellent results without any additional fixed parameters.

Results and discussion
All calculations are done using the ACES II [46] program system supplemented with the FS-CCSD (2,0) module [1]. For all 95 points included in calculations (see Table 1), the RHF function for the Li +2 2 cation has been adopted as a reference in accordance with the overall computational strategy presented in previous sections. The basis set used was obtained by supplementing the uncontracted ANO-RCC (relativistic Atomic Natural Orbital) [47] basis set with additional two diffuse functions for s, p, d, f shells with exponents given in [2], which resulted in the total size of the basis set equal to 114 functions per atom. The size of the active space for the FS-CCSD (2,0) part has been set to 86 which means that the same number of roots (and corresponding eigenvectors) has been found when solving the EA-EOM-CCSD equations. We observed that the size of the active space reached maximum and going beyond that limit did not affect significantly the results.
In Table 2, we compare the energy of the 1 1 g state at the infinite separation of atoms obtained with the FS-CCSD (2,0) with the sum of the energies of the separated Li atoms: one in the ground state and the other in the 3 P state, both obtained with the EA-EOM-CCSD method (equivalent to the FS-CCSD (1,0) scheme). For size-extensive methods, the two energy values should be identically equal and this situation occurs in the current case, comparing the energy values in the third and fourth column of Table 2. Note that in all cases the process of DEA engages both involved centres since the two added electrons are placed each in one atom. This observation is important since it means that the size intensive methods (like EOM DEA scheme) are not going to give the size extensive results contrary to the current FS-MRCC approach. In the range from 3.6 to 17.5 a.u., the energy calculations were done for points separated by 0.2 a.u., while in the range 18-40 a.u. with larger steps (0.5, 1.0 and 2.0). Due to the fact that the reference system (Li +2 2 ) dissociates into the closed shell fragments, both SCF (self-consistent field) and CC calculations easily converged even for large interatomic distances (for the largest one, the R was set to 200 a.u., see Table 1). In Figure 1, the computed PEC is shown (in the range 4.8-40 a.u.) and compared with the experimental curve. Available experimental points covered the range from 5.5 to 31.2 a.u. [48].
We also made an attempt to generate the PEC for the studied state with the MR(multireference)-DA(double attachment)-CCSD scheme introduced in [14], also known under DEA-EOM-CCSDT acronym [1], which relies on the CCSD ground state and EOM equations involving the R 2 (2p0h, p − particle, h − hole) and R 3 (3p1h) operators. Within the distances where we were able to converge the EOM equations, the curves closely overlap. However, in the range of 9.6-18.0 a.u., critical for the reproduction of the energies of the vibrational levels, the EOM solution failed to converge, see Figure 2.
The quality of the calculated PEC can be appreciated by comparison of the computed vibrational terms with the experimental data. In Table 3, we listed term values of the first 32 vibrational levels (v = 0 to v = 31) (obtained with help of the LEVEL-8 program of Le Roy [49]) and compared with the same number of the available experimental data. The agreement between two sets of data is surprisingly good. For majority of levels, the deviation from experiment is below 1 cm −1 and the value of mean absolute error is equal to 0.79 cm −1 . The only set of other theoretical term values for the 1 1 g state is reported in [50], see fourth and fifth column of Table 3, and the average deviation from experiment is larger nearly by an order of magnitude. Our estimate of the energy of the vibrational levels for the curves reported in other theoretical papers [51][52][53][54][55] indicates that the present results are very reliable. The excellent quality of the present theoretical results is demonstrated also in Figure 3 where the theoretical and experimental levels are represented by the solid and broken horizontal lines, respectively. As we see, they fully overlap which tells us that the applied computational scheme correctly reproduces the experimental data.
In Table 4, we compare inertial rotational constants B v obtained in this work with the experimental values. In most cases, the differences occur at the fourth decimal place which confirms the reliability of the computed PEC. To make the characteristics of the studied state more complete, we calculated also the set of the centrifugal distortion coefficients. They are included as the supplemental material. A summary of the calculation is given in Table 5, where the basic spectroscopic constants for 1 g state are presented. The deviations from the experimental values are really small, 0.003Å for the equilibrium bond length, three wavenumbers with respect to the dissociation energy and 0.004 eV for the adiabatic excitation energy. Compared to other theoretical works, the FS-CCSD (2,0) results look good; a similar quality results were obtained by Jasik et al. [51] using local potential approach. Note that DEA-EOM-CCSDT scheme turned out to be equally good as the FS approach providing the identical (R e and T e ) or nearly identical (D e = 1419 vs. 1426 cm −1 ) values of the studied constants.

Conclusions
The first principle quantum chemical calculations are applied in the theoretical study of PEC and vibrational energy levels for the 1 1 g state of the Li 2 molecule. Within the whole range of interatomic distances, the calculations were done for the RHF reference function with six electrons correlated.
The results remain in very good agreement with experiment both for the spectroscopic parameters as well as when the energy of vibrational levels is considered. A crucial thing in obtaining high accuracy results is an adoption of a proper reference system which is an Li +2 2 ion dissociating into a closed shell fragment in whole range of interatomic distances. This made it possible (1) to adopt a convenient RHF reference and (2)to use the FS scheme formulated in (2,0) sector to recover the original neutral Li 2 structure upon attaching two electrons to the Li +2 2 ion. The latter method being rigorously size extensive seems to be particularly useful in the study of the dissociation of the alkali metal diatomics.