Low-Lying Electronic States of the Nickel Dimer

The generalized Van Vleck second order multireference perturbation theory (GVVPT2) method was used to investigate the low-lying electronic states of Ni2. Because the nickel atom has an excitation energy of only 0.025 eV to its first excited state (the least in the first row of transition elements), Ni2 has a particularly large number of low-lying states. Full potential energy curves (PECs) of more than a dozen low-lying electronic states of Ni2, resulting from the atomic combinations 3F4 + 3F4 and 3D3 + 3D3, were computed. In agreement with previous theoretical studies, we found the lowest lying states of Ni2 to correlate with the 3D3 + 3D3 dissociation limit, and the holes in the d-subshells were in the subspace of delta orbitals (i.e., the so-dubbed δδ-states). In particular, the ground state was determined as X 1Γg and had spectroscopic constants: bond length (R e) = 2.26 Å, harmonic frequency (ωe) = 276.0 cm−1, and binding energy (D e) = 1.75 eV; whereas the 1 1Σg + excited state (with spectroscopic constants: R e = 2.26 Å, ωe = 276.8 cm−1, and D e = 1.75) of the 3D3 + 3D3 dissociation channel lay at only 16.4 cm−1 (0.002 eV) above the ground state at the equilibrium geometry. Inclusion of scalar relativistic effects through the spin-free exact two component (sf-X2C) method reduced the bond lengths of both of these two states to 2.20 Å, and increased their binding energies to 1.95 eV and harmonic frequencies to 296.0 cm−1 for X 1Γg and 297.0 cm−1 for 1 1Σg +. These values are in good agreement with experimental values of R e = 2.1545 ± 0.0004 Å, ωe = 280 ± 20 cm−1, and D 0 = 2.042 ± 0.002 eV for the ground state. All states considered within the 3F4 + 3F4 dissociation channel proved to be energetically high-lying and van der Waals-like in nature. In contrast to most previous theoretical studies of Ni2, full PECs of all considered electronic states of the molecule were produced.


INTRODUCTION
Since Ni 2 has few holes in otherwise complete subshells, one might expect theoretical studies of Ni 2 to be less complicated than for other first row transition metal dimers, like Cr 2 where one has many more possibilities of distributing 12 electrons in 12 valence orbitals. Reported information on Ni 2 , however, proves the contrary. For example, the exact symmetry of the ground electronic state of Ni 2 is still debated: Different studies have reported different space and spin symmetries for the molecule's ground term.
Experimental data on Ni 2 is sparse and the true ground state of the molecule is not unequivocally accepted. From the analysis of electronic absorption bands of Ni 2 in the visible spectral region in argon matrices, De Vore et al. (De Vore et al., 1975) determined ω e 192 cm −1 , whereas a frequency of 380.9 cm −1 was found in solid argon matrix (Ahmed and Nixon, 1979). The latter result was later criticized by Rasanen et al. (Rasanen et al., 1987) In photoelectron spectroscopic studies of Ni 2 -, ω e 280 ± 20 cm −1 was determined for the lowest electronic state of Ni 2 (Ho et al., 1993). Second and third law analyses of information derived from a combination of Knudsen effusion and mass-spectrometric techniques led to a binding energy of D 0 2.03 ± 0.30 eV (second law result) and D 0 2.36 ± 0.22 eV (third law result) for ground state Ni 2 (Kant, 1964). By using time-delayed resonant two-photon ionization, Morse et al. (Morse et al., 1984) determined D 0 2.068 ± 0.010 eV and R e 2.200 ± 0.007 Å for the lowest state of Ni 2 , but assigned as either 3 Γ u or 1 Γ g . Also from two-photon ionization studies on supersonic jet-cooled Ni 2 in argon carrier gas, Pinegar et al. (Pinegar et al., 1995) determined D 0 2.042 ± 0.002 eV and R e 2.1545 ± 0.0004 Å for the lowest state of Ni 2 but were unable to ascertain the symmetry of this state.
Theoretical studies of the electronic states of Ni 2 are complicated by the fact that the first excited state of the Ni atom, 3 D 3 (3d 9 4s 1 ), is only 0.025 eV above the 3 F 4 (3d 8 4s 2 ) ground atomic state, which is the least excitation energy of any of the first row of transition elements (Harrison, 2000;Kramida et al., 2013). This low promotion energy supports the existence of several low-lying molecular states of the Ni 2 molecule resulting from the 3 F 4 + 3 F 4 and 3 D 3 + 3 D 3 atomic combinations. For example, limited configuration interaction (CI) calculations (Shim et al., 1979) found 84 states of Ni 2 , corresponding to the 3 F 4 + 3 F 4 dissociation limit, to lie within an energy range of only 300 K (0.026 eV), and 45 states, also within a narrow energy range, to correlate with the 3 D 3 + 3 D 3 dissociation asymptote. Melius et al. (Melius et al., 1976) also noted that the manifold of electronic states within 0.50 eV of the ground state of Ni 2 was dense and complex.
Since the fully filled 4s-subshell of the 3 F 4 (3d 8 4s 2 ) ground state of Ni discourages significant bonding interaction, bonding in low-lying states of Ni 2 results largely from the coupling of excited state Ni atoms. In particular, the lowest states of the Ni 2 molecule can be expected to correlate with the 3 D 3 (3d 9 4s 1 ) + 3 D 3 (3d 9 4s 1 ) dissociation channel. At the generalized valence bond (GVB) and polarization CI (POL-CI) level of theory, (Upton and Goddard III, 1978) 30 of 45 low-lying states of Ni 2 , within the 3 D 3 + 3 D 3 dissociation channel, were found (Shim et al., 1979) to be singlets and triplets ordered energetically as δδ (6 states) < πδ (8 states) < δσ (4 states) < ππ (6 states) < πσ (4 states) < σσ (2 states), with the six lowest states (of symmetries 1 Γ g , 1 Σ g + , 3 Σ g − , 1 Σ u − , 3 Γ u , and 3 Σ u + ) being virtually degenerate and having an average equilibrium bond length, R e , of 2.04 Å, and binding energy, D e , of 2.29 eV. The designations δδ, πδ, et cetera, specify the positions of the holes in the dominant configurations of the 3d-orbitals at each atomic center which has a 3d 9 4s 1 configuration (Shim et al., 1979). In brief, the correct energy spacing of the low-lying states of the Ni dimer is problematic for both theory and experiment.
The determination of energy ordering and spacing is no less complicated for the Ni atom, and provides insight into necessary levels of experimental and theoretical approaches. The 3 F 4 (3d 8 4s 2 ) → 3 D 3 (3d 9 4s 1 ) excitation energy has been determined by at least one experimental study (Moore, 1952) to be negative (-0.029 eV). This negative value was supported by ab initio wave function and density functional theory (DFT) calculations (Bauschlicher et al., 1982;Raghavachari and Trucks, 1989;Russo et al., 1994). A more recent study, (Schultz et al., 2005) that employed several functionals at all five rungs of Jacob's ladder of DFT functionals, predicted the ground state configuration of the Ni atom as 3d 9 4s 1 ( 3 D 3 ) with most of the functionals when using a triple-ζ quality basis set. On the other hand, multireference studies (Andersson and Roos, 1992;Murphy and Messmer, 1992) predicted a 3d 8 4s 2 ( 3 F 4 ) ground state configuration for the Ni atom. Illustrating additional complexity, Upton and Goddard (Upton and Goddard III, 1978) found that averaging over J components (where J is the sum of spin and orbital angular momenta of the atom) of each state places the 3 D 3 (3d 9 4s 1 ) state lower energetically than the 3 F 4 (3d 8 4s 2 ) state.
Using an effective core potential basis set specifically optimized for the Ni atom in the 3 D 3 state within the generalized valence bond CI (GVBCI) method, Noell et al. (Noell et al., 1980) found the splitting of the six δδ-hole states of Ni 2 to be quite small (≤0.1 eV), with the lowest states being 1 Γ g and 1 Σ g + . Inclusion of polarization configurations involving single and double excitations to the virtual space (POLSDCI) placed the triplet states ( 3 Σ g -, 3 Γ u , 3 Σ u + ) approximately 0.07 eV below the singlets. At the singles and doubles CI (SDCI) level of theory, these authors found the energy splitting of the six lowest δδ-hole states of Ni 2 to be less than 0.009 eV, with an average bond length of 2.26 Å and binding energy of 1.88 eV.
With a basis set similar to that used by Noell et al., (Noell et al., 1980) a contemporaneous restricted Hartree-Fock (RHF) and CI with single and double excitations (CISD) study predicted a 3 Σ u + ground state for Ni 2 , with spectroscopic data: R e 2.33 Å, ω e 211 cm −1 , and D e 1.43 eV (Basch et al., 1980). Calculations by these authors at the same levels of theory using an all electron basis set corroborated the prediction of the ground state symmetry as 3 Σ u + . On the other hand, a local spin density method (Harris and Jones, 1979) predicted a 3 Σ g ground state with R e 2.18 Å, ω e 320 cm −1 , and D e 2.70 eV. A slightly earlier CASSCF/CASPT2 study (Pou-Amérigo et al., 1994) that used an atomic natural orbital (ANO) type contraction of the (21s15p10d6f4g) primitive basis to [6s5p4d3f2g] for calculations without correlation of the semi-core 3s3p electrons, and a [10s9p8d3f2g] contracted basis for calculations involving the correlation of 3s3p electrons, likewise found the six lowest δδ-hole states of Ni 2 to lie within a particularly narrow energy gap (0.04 eV) with the triplet states higher in energy than the singlets. However, after inclusion of scalar relativistic effects, the ground term was predicted as 1 Γ g , with the 1 Σ g + term lying only 0.01 eV higher at the equilibrium geometry. Correlating the 3s3p electrons in these calculations predicted the 1 Γ g and 1 Σ g + states to be degenerate to the reported accuracy, with slightly improved spectroscopic constants relative to reference experimental values [i.e., R e 2.23 Å, D e 2.06 eV; and ω e 293 cm −1 for 1 Γ g versus ω e 294 cm −1 for 1 Σ g + compared with experimental values of R e 2.1545 ± 0.0004 Å, (Pinegar et al., 1995) D 0 2.042 ± 0.002 eV, (Pinegar et al., 1995) and ω e 280 ± 20 cm −1 (Ho et al., 1993)].
DFT studies of Ni 2 have also been inconclusive. Yanagisawa et al. (Yanagisawa et al., 2000) used various DFT functionals to study the 3 Σ g and 3 Σ u + states of Ni 2 and found B3LYP (Becke, 1993) to predict the 3 Σ u + state to lie lower than 3 Σ g whereas the rest of the functionals predicted the latter state to lie lower at the equilibrium geometry. However, the 3 Σ g state that they found had a configuration that corresponded to the ππ-hole manifold rather than the δδ. Gutsev et al. (Gutsev and Bauschlicher, 2003) also found a 3 Σ g ground state, with the same configuration as did Yanagisawa et al., (Yanagisawa et al., 2000) when using a variety of hybrid functionals. Contrarily, Diaconu et al. (Diaconu et al., 2004) found a singlet δδ-hole ground state (with a mixture of 1 Γ g and 1 Σ g + ) for Ni 2 when using B3LYP with the (14s11p6d3f)/[8s6p4d1f] basis set, whereas use of the Stuttgart RSC ECP basis set (Dolg et al., 1987) with the same functional gave a triplet δδ-hole (with a mixture of 3 Σ g and 3 Γ u symmetries) that lay 0.001 eV lower than the singlet δδ-hole state at the equilibrium geometry. Using functionals at all levels of Jacob's ladder of DFT functionals, Schultz et al. (Schultz et al., 2005) also found different functionals to predict different ground state symmetries for Ni 2 , with all local spin density approximation (LSDA) functionals predicting a 3 Σ g ground state and all generalized gradient approximation (GGA) and meta GGA functionals predicting a 3 Π u ground state; whereas hybrid GGA and hybrid meta GGA functionals found either 3 Σ u + or 3 Σ g to lie lowest energetically. Du et al. (Du et al., 2008) used various functionals to study the low-lying states of Ni 2 . Their results that agreed best with experiment were obtained when using BLYP, which predicted a triplet σδ-hole (3d z2 σ u * 1 3d x2-y2 δ u * 1 ) ground state. The space symmetry of this state was not reported. With the B3P86 functional, (Perdew, 1986;Becke, 1993) a quintet ground state was predicted for Ni 2 , (Shi-Ying and Zheng-He, 2008) although the space symmetry was not reported. The Perdew-Burke-Ernzerhof (PBE) exchange correlation functional (Perdew et al., 1996) predicted a 3 Σ g ground state for Ni 2 , (Kamal et al., 2012) with spectroscopic constants R e 2.93 Å, D e 3.09 eV, and ω e 334.08 cm −1 , which showed significant deviations from experimental values.
Some of the most recent wave function based calculations on Ni 2 include those due to Dong et al. (Dong et al., 2013) using the symmetry-adapted-cluster configuration interaction (SAC-CI) method (Nakatsuji, 1979) and Cheskidov et al. (Cheskidov et al., 2012) using the average coupled pair functional (ACPF), (Gdanitz and Ahlrichs, 1988) average quadratic coupled cluster (AQCC), (Szalay and Bartlett, 1993) internally contracted single and double multireference configuration interaction (MRCI or MRCI with Davidson corrections, i.e., MRCI + Q), (Werner and Knowles, 1990) and N-electron valence state second-order perturbation theory (NEVPT2) (Angeli et al., 2001) methods. The study by Dong et al. (Dong et al., 2013) predicted a 3 B 1u ground state (with R e 2.56 Å) for Ni 2 in D 2h symmetry which corresponds to 3 Σ u + , 3 Δ u or 3 Γ u in D ∞h . The study by Cheskidov et al. (Cheskidov et al., 2012) used the Dunning-type quadruple-ζ quality basis set, cc-pVQZ-DK (22s18p11d3f2g1h/ [8s7p6d3f2g1h]), (Balabanov and Peterson, 2005) and found the 1 Γ g and 1 Σ g + δδ-hole states to be quasidegenerate for all five methods with the 1 Σ g + state lying lower when using AQCC, MRCI + Q, and MRCI methods and the two states fully degenerate (to the reported accuracy level) at the ACPF and NEVPT2 levels. At the ACPF level, the predicted ground state was instead 1 Σ u -; inclusion of spin-orbit relativistic corrections within ACPF calculations led to an 0 + g ground state ( 1 Σ g + + 3 Σ g δδ-hole states), whereas the 0 − u term ( 1 Σ u -+ 3 Σ u + δδ-hole states) lay at only 0.009 ± 0.004 eV above the predicted ground state.
The above synopsis of previous work on Ni 2 shows the difficulties involved in studying not only the spectroscopic constants, but even the ordering of the low-lying electronic states of Ni 2 . Although wave function methods generally support δδ-hole states ( 1 Γ g , 1 Σ g ) as lying lowest energetically, the methods predict different ground state symmetries with some finding all six states to be degenerate. Similarly, experimental spectroscopic data have been obtained but most of the studies could not ascertain the ground state symmetry of the molecule. Our current study exploits the ability of the GVVPT2 method (Khait et al., 2002;Jiang et al., 2009) to describe well full PECs of ground-and excited-electronic states of complicated transition element dimers, such as has already been demonstrated on other problematic transition metal molecules [e.g., Cr 2 and Y 2 (Tamukong et al., 2012;Tamukong et al., 2014)]. It should be noted that of all previous theoretical work described above on electronic states of Ni 2 , only six of the articles reported full PECs of the states they investigated. Consequently, where other data is available, this study also provides further assessment of the capabilities of the GVVPT2 method for difficult transition metal dimers, including the bond breaking regions. We have constructed full PECs of 21 states of Ni 2 . All (nonrelativistic) calculations used the Dunning-type cc-pVTZ basis set (Balabanov and Peterson, 2005), and calculations were performed using D 2h symmetry. Furthermore, low-lying 1 Γ g and 1 Σ g + states were further studied with scalar relativistic effects included in GVVPT2 (Tamukong et al., 2014) through the spin-free exact two component (sf-X2C) method (Liu, 2010;Cheng and Gauss, 2011;Li et al., 2012;Li et al., 2013). Relativistic calculations used the cc-pVTZ-DK basis set (Balabanov and Peterson, 2005). The rest of the paper is organized as follows: Methods section briefly reviews key features of GVVPT2 and the spin-free exact two component (sf-X2C) methods, and describes computational details; the results are presented and discussed in Results and Discussion section; while conclusions are drawn in Conclusion section.

GVVPT2
The GVVPT2 method for electron correlation in molecules has been thoroughly described elsewhere, (Khait et al., 2002;Jiang et al., 2009) as has its application to some challenging systems (Mbote et al., 2010;Khait et al., 2012;Tamukong et al., 2012;Mokambe et al., 2013;Tamukong et al., 2014). Here, salient features of the GVVPT2 method relevant to the present study are reviewed. In GVVPT, the total Hilbert space of many-electron functions [e.g., configuration state functions (CSFs)] with appropriate molecular space and spin symmetry (L) is partitioned into a model space (L M ), and an external space (L Q ) whose electronic configurations are derived from configurations generated from the model space by single and double electron excitations into virtual orbitals, L L M ⊕ L Q . The model space is further partitioned into a primary subspace, spanned by a set of reference functions that are linear combinations of CSFs (typically the lowest MCSCF or CASSCF states), and a secondary subspace. States within the primary subspace are then perturbatively corrected through primary-external (P-Q) interactions whereas the secondary subspace serves as a buffer energetically separating the primary and external subspaces. This energy buffer circumvents most intruder state problems that plague many multireference perturbation theory techniques. An effective Hamiltonian matrix (H eff MM ) is constructed with the same dimension as the model space, and its diagonalization includes both perturbatively corrected primary and unperturbed secondary subspace states. To guarantee continuity and smoothness of CSF responses (and ultimately PECs) even in situations of quasidegeneracies between primary and external states, GVVPT2 uses a continuously varying nonlinear denominator shift that arises from a resolvent that is both degeneracy-corrected and contains a hyperbolic tangent function as a switching function from nondegenerate to degenerate regimes (Khait et al., 2002). This resolvent can be written where the C mi denote components of eigenvectors of the unperturbed model Hamiltonian; ε (0) i is the reference Møller-Plesset-type energy while ε (0) q is the state-specific zeroth-order energy of external CSF q. As with all studies of specific molecules with GVVPT2 since ca. 2005, semicanonical orbitals from the MCSCF are used. The mathematical robustness of GVVPT2 allows it to support both complete and incomplete model spaces. By construction, GVVPT2 is subspace-specific (N.B. excited states of the same symmetry as lower-lying states can be calculated) and is spin-adapted.

Spin-Free Exact Two Component Method
Matrix formulations of two-component relativistic methods have had significant successes in molecular electronic structure calculations, following Dyall's seminal Normalized Elimination of Small Component (NESC) paper (Dyall, 1997) in 1997. For a summary of development of matrix two-component methods, we refer the reader to Ref. (Liu, 2014). In this work, we adopt the spin-free version of exact two component (sf-X2C) (Liu, 2010;Cheng and Gauss, 2011;Li et al., 2012;Li et al., 2013) that was previously used with GVVPT2 (Tamukong et al., 2014). The following Hamiltonian, written in second quantization form, incorporates scalar relativistic effects through the sf-X2C approach where the first term is the one-electron spin-free part of the exact two-component (X2C) Hamiltonian, while the second is the unmodified Coulombic two-electron term. The sf-X2C Hamiltonian for positive energy states derives from a modified Dirac Hamiltonian, h D which is decomposed into spin-free (sf) and spin-dependent (sd) parts.
where V is the matrix representation of the external nuclear attraction potential operator; T is the matrix representation of the kinetic energy operator; α is the fine-structure constant; while W is the matrix representation of the operator The spin-free (sf) part of W, that is W sf , describes scalar relativistic effects whereas the spin-dependent (sd) part, W sd , incorporates spin-orbit coupling effects.

Computational Details
Unperturbed wavefunctions for GVVPT2 calculations were obtained from multiconfigurational self-consistent field (MCSCF) reference functions, using the local code "undmol." The macroconfiguration approach (Khait et al., 2004) was used in all MCSCF and GVVPT2 calculations. In the macroconfiguration scheme, active orbitals are partitioned into orbital groups based on physical and mathematical intuition, such as orbital type or energy, and electrons are assigned to active orbital groups. Each Frontiers in Chemistry | www.frontiersin.org May 2021 | Volume 9 | Article 678930 unique assignment of electrons to active orbital groups defines a macroconfiguration, κ(n). Orbital rotations within each active orbital group are redundant (in MCSCF calculations), as all possible distributions of electrons within each group leading to the desired total spin and space symmetry are allowed. Each κ(n) generates a unique set of configurations, and hence configuration state functions (CSFs) that are orthogonal to those of all other κ(n). Additionally, the macroconfiguration approach permits a large number of noninteracting electronic configuration pairs to be efficiently screened by recognizing that all matrix elements resulting from configurations that differ by more than two electrons must be identically zero. The macroconfiguration formalism also provides an efficient way of generating excited configurations.
In all calculations, the active space consisted of 3d and 4sderived molecular orbitals (MOs) of Ni 2 (Figure 1 shows representative 3d and 4s-derived MOs that were used to constitute the active spaces in this study). Depending on the specific state being investigated, some of the 3d-derived MOs and/or 4s-derived MOs were restricted to be doubly occupied in the MCSCF calculations, but were included with the 3s-and 3pderived MOs in the active core and correlated at the GVVPT2 level of theory. For example, in all MCSCF calculations of δδ-hole states the 3dσ and 3dπ MOs were placed in the active core and only correlated at the GVVPT2 level. Similarly, the 3dσ and 3dπ electrons were kept doubly occupied in MCSCF calculations of ππ-hole states while only the 3dσ electrons were kept doubly occupied in MCSCF calculations of δπ-hole states, whereas the 4sσ, or 4sσ + 3dπ, or 4sσ + 3dσ orbitals were kept doubly occupied in MCSCF calculations of states within the 3 F 4 (3d 8 4s 2 ) + 3 F 4 (3d 8 4s 2 ) dissociation channel. The remaining orbitals in the active space were partitioned into sets (or so-called orbital groups), leading to configurations that describe δδ-, δπ-, δσ-, ππ-, πσ-, or σσ-hole states from the 3 D 3 (3d 9 4s 1 ) + 3 D 3 (3d 9 4s 1 ) atomic combination or configurations that describe 3 D 3 (3d 9 4s 1 ) + 3 F 4 (3d 8 4s 2 ) and 3 F 4 (3d 8 4s 2 ) + 3 F 4 (3d 8 4s 2 ) atomic couplings. In the present study, we indicate the positions of holes within the subspace of active 3d-derived MOs by specifying the active orbital types that qualitatively describe where holes exist (e.g., δδ-, δπ-, ππ-) as did Shim et al. (Shim et al., 1979) and Noel et al. (Noell et al., 1980) All δδ-states were computed using a single reference macroconfiguration, where the superscripts denote the number of electrons in each orbital group. The semi-core 3s3p electrons were correlated together with those derived from 3d z2, 3d xz and 3d yz at the GVVPT2 level. For four of the computed δδ-states, additional calculations were also performed in which the 3s3p were kept doubly occupied throughout (i.e., at both the MCSCF and GVVPT2 levels). When using the cc-pVTZ basis set, (Balabanov and Peterson, 2005) reference κ 1 (n) generated 8 model and 27,891,120 total CSFs (in D 2h symmetry) for the 1 3 Σ u + and 1 3 Γ u states; 8 model and 15,290,666 total space CSFs for the 1 1 Σ u − , 1 1 Σ g − , and 1 1 Γ u states; 10 model and 27,982,592 all space CSFs for the 1 3 Σ g and 1 3 Σ u states; and 12 model versus 15,270,687 all space CSFs for the X 1 Γ g and 1 1 Σ g + states. Without correlating the 3s3p electrons at the GVVPT2 level, the dimensions were: 12 model space CSFs versus 3,593,707 total CSFs for the X 1 Γ g and 1 1 Σ g + states; and 8 model versus 6,434,550 all space CSFs for the 1 3 Σ u + and 1 3 Γ u states. Relativistic calculations on the X 1 Γ g and 1 1 Σ g + states utilized the same reference macroconfiguration, κ 1 (n). A δδ 3 Γ u state was computed with 4 active electrons in 4 orbitals using reference κ 2 (n) 3d x 2 −y 2 δ g 3d x 2 −y 2 δ p u 6 4sσ g 4sσ p u 2 (6) which gave rise to 4 model space and 7,518,688 all space CSFs, using the cc-pVTZ basis set. The ππ-states were computed using reference which is similar to κ 1 (n) but with delta replaced with pi orbitals. This macroconfiguration gave rise to 12 model space and 15,267,629 all space CSFs for the 1 1 Δ g and 2 1 Σ g + states when using the cc-pVTZ basis set and D 2h molecular space symmetry. The δπ-states were computed from κ 4 (n) 3dδ g 3dδ p u 3dδ g 3dδ p u 7 3dπ u 3dπ p g 3dπ u 3dπ p g 7 4sσ g 4sσ p u 2 (8) κ 4 (n) generated 16 model space and 27,178,852 total space CSFs for the 1 1 Φ g and 1 1 Π g states versus 20 model and 50,736,846 all space CSFs for the 1 3 Φ g and 1 3 Π g states, when using the cc-pVTZ basis set. Within the 3 F 4 + 3 F 4 manifold, the 3 1 Σ g + , 1 3 Σ g + , 2 3 Σ g + , 1 3 Δ u , and 2 3 Σ u + δπ-states were computed using κ 5 (n) 3dδ g 3dδ p u 3dδ g 3dδ p u 6 3dπ u 3dπ p g 3dπ u 3dπ p g 6 (9) In these calculations, the 3s, 3p, 3d z2 , and 4s electrons were kept doubly occupied at the MCSCF level but correlated at the GVVPT2 level of theory. Reference κ 5 (n) resulted in 40 model versus 55,053,638 total CSFs for the 3 1 Σ g + state; 36 model and 103,306,512 all space CSFs for the 1 3 Σ g + and 2 3 Σ g + states; and 40 model versus 103,312,902 all space CSFs for the 1 3 Δ u and 2 3 Σ u + states when using the cc-pVTZ basis set.
Two δπ-states, 2 1 Γ g and 2 1 Δ g , were computed within the 3 F 4 + 3 F 4 manifold using κ 6 (n) 3dδ g 3dδ p u 3dδ g 3dδ p u 6 3dσ g 3dσ p u 2 (10) Reference κ 6 (n) generated 12 model space and 15,270,687 all space CSFs for the computed states using the cc-pVTZ basis set. Two quintet δπσ-states (i.e., 1 5 Φ u and 1 5 Π u ) were also computed within the 3 F 4 + 3 F 4 manifold, using κ 7 (n) 3dδ g 3dδ p u 3dδ g 3dδ p u 6 3dπ u 3dπ p g 3dπ u 3dπ p g 7 3dσ g 3dσ p u 3 (11) This reference, κ 7 (n), led to 12 model versus 69,738,914 total CSFs for the computed quintet states. Lastly, a limited study was done on two quintet states within the 3 D 3 (3d 9 4s 1 ) + 3 F 4 (3d 8 4s 2 ) manifold at short bond lengths. While reference κ 8 (n) 3dδ g 3dδ p u 3dδ g 3dδ p u 7 3dπ u 3dπ p g 3dπ u 3dπ p can describe the 1 5 Δ g and 2 5 Δ g δπσ-states resulting from the 3 D 3 (3d 9 4s 1 ) + 3 F 4 (3d 8 4s 2 ) manifold of molecular states adequately at short internuclear distances, it cannot be expected to do so at long internuclear distances and our studies were restricted to the shorter bond lengths for these two states. Reference κ 8 (n) generated 12 model CSFs versus 69,740,135 total space CSFs. The reference macroconfigurations described above were used to define the active space, while all lower energy MOs were doubly occupied in MCSCF calculations. Initial MOs for MCSCF calculations were obtained from approximate natural orbitals of second-order restricted Møller−Plesset perturbation (RMP2) calculations from a closed-shell Hartree−Fock (HF) reference. At the GVVPT2 level, 3s, 3p, and all 3d and/or 4s electrons were correlated whether they were or were not at the MCSCF level. For comparison purposes, a few of the GVVPT2 calculations were performed without correlating the 3s and 3p electrons. Calculations that accounted for scalar relativistic effects employed the sf-X2C method described above. Finally, to aid in interpretation, the effective bond order (EBO) was computed, and used the following expression where η is the EBO, χ i is the EBO for the ith CSF, while c i 2 is its corresponding weight. For each CSF used to estimate EBO, χ i was determined as where n b and n ab are the numbers of bonding and antibonding electrons, respectively. Vibrational frequencies were obtained by 3-point finite differencing near the minima of the curves.

RESULTS AND DISCUSSION
Where indicated, the letter "R" in parentheses following a molecular term denotes that scalar relativistic effects were included in the calculations, while the expression "no 3s3p" in parentheses after a molecular term symbol denotes that 3s and 3p electrons were not correlated in GVVPT2 calculations.

The δδ-Hole States
PECs of the δδ-hole states are shown in Figure 2 and the data describing them are in Table 1. In agreement with results from other high level ab initio methods, the lowest states of Ni 2 were found to be δδ-hole states of the 3 D 3 (3d 9 4s 1 ) + 3 D 3 (3d 9 4s 1 ) manifold. In particular, the ground state was found to be X 1 Γ g , with the 1 1 Σ g + state lying only 16.40 cm −1 (0.002 eV) higher at the equilibrium geometry. After including scalar relativistic effects, the energy gap between these states slightly increased to 23.39 cm −1 at equilibrium, with the X 1 Γ g term having spectroscopic constants: R e 2.20 Å, D e 1.95 eV, and ω e 296 cm −1 . These results are in good agreement with experimental data (R e 2.1545 ± 0.0004 Å, (Cooper et al., 1972) D o 0 2.042 ± 0.002 eV, (Cooper et al., 1972) and ω e 280 ± 20 cm −1 ) (Rösch and Rhodin 1974) and with the relativistic CASSCF/CASPT2 results of Pou-Amérigo et al. (Pou-Amérigo et al., 1994) who also found the lowest 1 Γ g and 1 Σ g + terms to be quasidegenerate (with R e 2.23 Å, D e 2.06 eV, and ω e 293 cm −1 for the 1 Γ g term). These results are also in good agreement with the time-delayed resonant two-photon ionization study of Morse and co-workers (Pinegar et al., 1995) that predicted either a 3 Γ u or 1 Γ g state as the ground state of Ni 2 and the scalar relativistic calculations of Cheskidov et al. (Cheskidov et al., 2012) who found the lowest 1 Γ g and 1 Σ g + terms to be degenerate at the ACPF and NEVPT2 levels of theory and the FIGURE 2 | PECs of low-lying electronic states of Ni 2 computed at the GVVPT2 level of theory using the cc-pVTZ basis set. All energies are plotted relative to the lowest energy value of the ground X 1 Γ g term. For all states, the holes are in the 3d delta orbitals (δδ-holes) except for 1 1 Δ g and 2 1 Σ g + states which are ππ-hole states.  (Cooper et al., 1972) (the reported binding energy is for D 0 ). b Ref. (Rösch and Rhodin, 1974).
Frontiers in Chemistry | www.frontiersin.org May 2021 | Volume 9 | Article 678930 Σ g + term to lie very slightly lower than the 1 Γ g at the AQCC, MRCI, and MRCI + Q levels of theory. At 2.25 Å, we found the EBO to be 0.963 and 0.960 for the 1 Γ g and 1 Σ g + states, respectively (using the largest dozen CSFs). In relativistic GVVPT2 calculations, these EBOs increased slightly to 0.975 for 1 Γ g and 0.972 for 1 Σ g + at this geometry. The major configurations describing the X 1 Γ g and 1 1 Σ g + states involved two holes in the same δ-type; i.e., 0.445 × (3d xy δ g 2 3d xy δ u p2 3d x2-y2 δ g 2 4sσ g 2 + 3d xy δ g 2 3d x2-y2 δ g 2 3d x2-y2 δ u p2 FIGURE 3 | PECs of low-lying δδ-hole electronic states of Ni 2 computed at the GVVPT2 level of theory, with and without the correlation of 3s3p semi-core electrons, using the cc-pVTZ basis set. All energies are plotted relative to the lowest energy value of the ground X 1 Γ g term. FIGURE 4 | PECs of the lowest-lying δδ-hole X 1 Γ g and 1 1 Σ g + states of Ni 2 computed at the GVVPT2 level of theory, with and without scalar relativity included, using the cc-pVTZ (or cc-pVTZ-DK) basis set. Non-relativistic energies are plotted relative to the lowest energy value of the ground X 1 Γ g term while relativistic energies are plotted relative to the lowest energy of the X 1 Γ g (R) term.
The semi-core 3s3p electrons were found to be important in the description of low-lying states of Ni 2 . The inclusion of 3s3p electron correlation at the GVVPT2 level increased the binding energies by 0.09 eV for X 1 Γ g , 0.10 eV for 1 1 Σ g + , and 0.09 eV for 1 3 Σ u + and 1 3 Γ u states in non-relativistic calculations. As can be seen in Figure 3 and Table 1, the effects of the 3s3p electrons on the equilibrium bond lengths and harmonic frequencies for these states are minimal whereas inclusion of such core-valence correlation raises, for example, the binding energy of X 1 Γ g from 1.66 to 1.75 eV compared to a reference D 0 value of 2.042 ± 0.002 eV (Cooper et al., 1972). Scalar relativistic effects shortened the bond length of X 1 Γ g by 0.06 Å and further increased the bond energy by 0.20 to 1.95 eV which agreed even better with the reference experimental values (see Figure 4 and Table 1). The 3s3p electrons did not have any effect on the EBOs of X 1 Γ g and 1 1 Σ g + ; the EBOs were determined as 0.962 and 0.959 at 2.27 Å for X 1 Γ g and 1 1 Σ g + , respectively, when the 3s3p electrons were not correlated compared to 0.963 vs. 0.960 when the semi-core electrons (3s3p) were correlated. Note the quasidegeneracy in the X 1 Γ g and 1 1 Σ g + states. For example, in Figure 4, the blue and green curves for the X 1 Γ g and 1 1 Σ g + states, respectively, lie on top of each other (only the green is visible). Also, the black and red curves for the X 1 Γ g (R) and 1 1 Σ g + (R) states lie on each other (only the red curve is visible).
The 1 1 Σ u state, which was predicted as the ground state of Ni 2 at the ACPF level of theory (Cheskidov et al., 2012) and found to lie quite close to a 1 Σ g + ground state in a limited CI study (Ahmed and Nixon, 1979), was found at the GVVPT2 level to lie 91.09 cm −1 above the X 1 Γ g state at equilibrium. The 2 1 Σ u state, however, lay much higher energetically (18,575.76 cm −1 above the ground state at equilibrium).
As can be seen in Table 1, GVVPT2 predicted the triplet δδ-hole states, 1 3 Σ g -, 1 3 Σ u + , and 1 3 Γ u , to lie energetically in the order 1 3 Σ g -< 1 3 Σ u + < 1 3 Γ u . Cheskidov et al. (Cheskidov et al., 2012) found this same ordering at the ACPF, AQCC, MRCI and MRCI + Q levels of theory, whereas their NEVPT2 calculations predicted the ordering 1 3 Σ u + < 1 3 Γ u < 1 3 Σ g -, with the 1 3 Σ g state lying at least 139 cm −1 higher than the other two states. It should be noted that the vertical excitation energies in Ref. (Cheskidov et al., 2012) were not determined at the equilibrium geometries of the computed states. The 1 3 Σ u + state, which was predicted as the ground state of Ni 2 in some previous wavefunction (Wood et al., 1980;Schultz et al., 2005) and DFT (Diaconu et al., 2004;Kramida et al., 2013) studies, was found in our current study to lie 349.60 cm −1 above the X 1 Γ g state at the equilibrium geometry. Likewise, the 1 3 Σ g state reported in other studies (Murphy and Messmer, 1992;Diaconu et al., 2004;Kramida et al., 2013) as the ground state of Ni 2 lay 221.98 cm −1 higher at equilibrium. The 1 3 Σ u state was found to have a bond length and bond energy comparable to those of 1 3 Σ g -, 1 3 Σ u + , and 1 3 Γ u but lying at least 531.48 cm −1 higher in energy. The EBOs for these triplet states were 0.971 for 1 3 Σ g -, 0.933 for 1 3 Σ u + and 1 3 Γ u , and 0.923 for 1 3 Σ u in the vicinity of their equilibrium geometries (again computed with 10 leading CSFs). The major configurations for the δδ-hole triplet states involved a doubly occupied 4sσ g bonding orbital.
The 1 3 Γ u state, in which both δ-holes were in the 3d x2-y2 δ g and δ u p orbitals, was computed using reference κ 2 (n). As can be seen in Table 1, this state was found to have spectroscopic constants comparable to other δδ-hole triplet states but lay much higher energetically (2442.21 cm −1 above the ground state at equilibrium). The present results suggest that the 3dδ orbitals are indeed split in the bonding interaction. Since they are nondegenerate, the Aufbau principle suggests that energetically low-lying orbitals (the bonding 3dδ orbitals) be occupied before higher ones. Moreover, Hund's rule suggests that orbitals with similar energies (in this case, 3dδ u p orbitals) be singly occupied before electron pairing occurs. This seems to bring qualitative understanding into why the 1 3 Σ g -, 1 3 Σ u + , and 1 3 Γ u states in which the 3d x2-y2 δ u p and 3d xy δ u p are singly occupied lie energetically lower than the 2 3 Γ u state for which the 3d x2-y2 δ g and 3d x2-y2 δ u p are singly occupied.

The δπ-Hole and ππ-Hole States
The PECs of the computed ππ-hole states (1 1 Δ g and 2 1 Σ g + ) are shown in Figure 2, while those for the δπ-hole states (1 1 Φ g , 1 1 Π g , 1 3 Π g , and 1 3 Φ g ) are shown in Figure 5 and compared with the ground state PEC. The data describing these curves are in Table 2. GVVPT2 predicted the ππ-hole states to lie higher in energy than the δπ-hole states, in agreement with previous studies (Ahmed and Nixon, 1979;Rasanen et al., 1987;Russo et al., 1994). For all four δπ-hole states studied, the major CSFs involved a doubly occupied 4sσ g bonding orbital. Thus, the main configurations of the 1 1 Φ g and 1 1 Π g states involved an unpaired spin-increasing electron in the 3dδ subspace and an unpaired spin-decreasing electron in the 3dπ subspace, e.g., 3dδ g 2 3dδ u p2 3dδ g 2 3dδ u p(↑) 3dπ u 2 3dπ g p2 3dπ u 2 3dπ g p(↓) 4sσ g 2 , whereas the major configurations of the 1 3 Φ g and 1 3 Π g states were similar to those of the corresponding singlet states but with two unpaired spin-increasing electrons: one in each of the 3dδ and 3dπ subspaces, e.g., 3dδ g 2 3dδ u p2 3dδ g 2 3dδ u p(↑) 3dπ u 2 3dπ g p2 3dπ u 2 3dπ g p(↑) 4sσ g 2 . At the equilibrium geometries, the EBOs were 0.930 for the singlet 1 1 Φ g and 1 1 Π g states and 0.933 for the 1 3 Π g and 1 3 Φ g states. GVVPT2 predicted the four δπ-hole states to lie energetically in the order 1 1 Φ g < 1 1 Π g < 1 3 Π g < 1 3 Φ g , in agreement with the Ref. (Cheskidov et al., 2012) study at the scalar relativistic ACPF, AQCC, MRCI, and MRCI + Q levels of theory. However, our calculations found all three states considered in Ref. (Cheskidov et al., 2012) (i.e., 1 1 Π g , 1 3 Π g , and 1 3 Φ g ) to lie some 500 cm −1 closer to the ground state; e.g., at the scalar relativistic MRCI + Q level, the 1 3 Φ g state was reported (Cheskidov et al., 2012) as lying 1238 cm −1 above the ground state at 2.5 Å while non-relativistic GVVPT2 calculations predicted this state to lie 546.76 cm −1 above the ground state at equilibrium. Based on our observation that including scalar relativistic effects increased the energy gap between the 1 1 Σ g + and X 1 Γ g states, it is likely that including such effects in GVVPT2 calculations on the δπ-hole states might lead to increases in corresponding adiabatic transition energies. It is not anticipated, however, that such effects would lead to any change in the energy ordering of the states.
Although lying higher in energy than the δπ-hole states, the ππ-hole states were found to have slightly shorter bond lengths and larger bond strengths than the δπ-hole states. The 1 1 Δ g state was 0.06 Å shorter while the 2 1 Σ g + state was 0.01 Å shorter in bond length than the δπ-hole states. At 2.24 Å, the EBOs of 1 1 Δ g and 2 1 Σ g + were 1.108 and 1.084 respectively (computed with a dozen CSFs); these were slightly higher than the EBOs of all computed δδ-hole and δπ-hole Ni 2 states. Near equilibrium, the major configurations of these ππ-hole states involved a doubly occupied 4sσ g bonding orbital and a configuration of the 3dπ subspace that had the two π-holes in the same π-orbital; e.g., 3dπ u 2 3dπ g p2 3dπ u 2 3dπ g p0 4sσ g 2 .
States of the 3 F 4 + 3 F 4 and 3 F 4 + 3 D 3 Manifolds Figure 6 contains PECs of states belonging to the 3 F 4 + 3 F 4 manifold. The data describing these curves are in Table 2.
Irrespective of how the model space was partitioned into macroconfigurations [i.e., κ i (n)], all such states were found to be van der Waals-like with interaction energies ≤0.04 eV. For example, near the equilibrium geometry (i.e., 3.77 Å), the 2 1 Γ g and 3 1 Σ g + states had EBO of only 0.005, while the 1 5 Φ u and 1 5 Π u states had EBOs of 0.003 and 0.00, respectively, at 3.84 Å. These latter quintet states were computed using reference κ 7 (n) and found to have the lowest total energies among the computed states of the 3 F 4 + 3 F 4 manifold; the 1 5 Π u state is 3.312 cm −1 less stable than the 1 5 Φ u state at equilibrium. The rest of the PECs in Figure 6 were plotted relative to the 1 5 Φ u PEC.
Lastly, the 1 5 Δ g and 2 5 Δ g states of the 3 F 4 + 3 D 3 manifold were investigated, at short bond lengths only, using reference FIGURE 5 | PECs of low-lying δπ-hole electronic states of Ni 2 computed at the GVVPT2 level of theory using the cc-pVTZ basis set compared with the PEC of the X 1 Γ g ground term. All energies are plotted relative to the lowest energy value of the X 1 Γ g ground term.

CONCLUSION
The GVVPT2 method was used to study low-lying electronic states of Ni 2 . The results indicate, in general, that bonding in these states involves predominantly the doubly occupied 4sσ g bonding orbital with the 3d-3d electrons pairwise spin coupled (e.g., consistent with antiferromagnetism). This statement is corroborated by EBOs that were found to be approximately 1.0 for most states studied and, moreover, states that did not allow this type of interaction [e.g., belonging to the 3 F 4 (3d 8 4s 2 ) + 3 F 4 (3d 8 4s 2 ) manifold] were found to be bound only by weak polarization forces with bond orders close to zero. For computed states of the 3 D 3 (3d 9 4s 1 ) + 3 D 3 (3d 9 4s 1 ) dissociation limit, all major configurations involved a doubly occupied 4sσ g bonding orbital and a vacant 4sσ u p antibonding orbital. The energy ordering of the computed states of the 3 D 3 (3d 9 4s 1 ) + 3 D 3 (3d 9 4s 1 ) manifold is in agreement with previous studies that found the δδ-hole states to lie lowest in energy followed by the δπ-hole and then the ππ-hole states (Rasanen et al., 1987). For the investigated δδ-hole states, the singlets were more stable than the triplet states at the GVVPT2 level of theory. As expected, based on previous GVVPT2 studies of transition metal dimers, all computed PECs are smooth and without unphysical artifacts (e.g., wiggles).
The ground state of Ni 2 was predicted as X 1 Γ g in agreement with previous results from other high level ab initio methods (Cheskidov et al., 2012;Pou-Amérigo et al., 1994). The determined equilibrium spectroscopic constants of the X 1 Γ g state, using a cc-pVTZ (or cc-pVTZ-DK) basis, were within the uncertainties of experimental results. The lowest 1 Σ g + state was found to be only 0.002 eV (16.40 cm −1 ) higher than the ground state at the equilibrium geometry. Core-valence correlation was found to be important in the description of low-lying states of Ni 2 where the inclusion of 3s3p electron correlation at the GVVPT2 level was shown to improve harmonic frequencies and bond energies (e.g., by an increase of 7.5 cm −1 in frequency and 0.09 eV in bond energy when 3s3p electron correlation was included in GVVPT2 calculations of the X 1 Γ g state). Scalar relativistic effects were also shown to be important, especially for dissociation energy, where spectroscopic constants from relativistic calculations were predicted to agree better with reference data (e.g., a decrease of 0.06 Å in bond length, increase of 0.20 eV in bond energy, and an increase of 20.0 cm −1 in harmonic frequency when including scalar relativistic effects in calculations of the X 1 Γ g state). In our previous study (Tamukong et al., 2014) on electronic states of Y 2 , Mn 2 , and Tc 2 , we did not find scalar relativistic effects to be as important for the Mn 2 molecule as they have proven to be in the present study. The inclusion of spin-orbit coupling effects was previously found (Cheskidov et al., 2012;Rasanen et al., 1987;Pou-Amérigo et al., 1994) to mix the low-lying states of Ni 2 , FIGURE 6 | PECs of electronic states of Ni 2 from the 3 F 4 (3d 8 4s 2 ) + 3 F 4 (3d 8 4s 2 ) manifold, computed at the GVVPT2 level of theory using the cc-pVTZ basis set and reference macroconfigurations κ 5 (n) to κ 7 (n). All energies are plotted relative to the lowest energy value of the 1 5 Φ u term.
Frontiers in Chemistry | www.frontiersin.org May 2021 | Volume 9 | Article 678930 leading to a 0 g + [ 1 Σ g + (δδ) + 3 Σ g -(δδ)] ground state. Since the two curves are close and nearly parallel for all internuclear separations in our study, we expect that including such effects after our scalar relativistic GVVPT2 treatment should lead to a similar mixing of states (i.e., and without change in spectroscopic constants).
The states investigated within the 3 F 4 (3d 8 4s 2 ) + 3 D 3 (3d 9 4s 1 ) manifold suggested significant bonding interaction, giving large harmonic frequencies and short bond lengths in comparison with states correlating with the 3 F 4 (3d 8 4s 2 ) + 3 F 4 (3d 8 4s 2 ) dissociation limit. Further work on Ni 2 should explore the 3 F 4 (3d 8 4s 2 ) + 3 D 3 (3d 9 4s 1 ) manifold more thoroughly, including possible expansion of the active space. It should be noted, however, that in the present study, we did not observe any significant electron excitations from the valence orbitals to 4p-dominated virtual orbitals. Also, while spin-free relativistic effects are noticeable (especially for bond dissociation energies), the parallelity and small energy separations for states in these manifolds are very similar to those in the 3 D 3 (3d 9 4s 1 ) + 3 D 3 (3d 9 4s 1 ) dissociation channel and, again, we do not expect that spin-dependent relativistic effects would change PECs appreciably, although current capabilities of our GVVPT2 code did not let us explore this.
In summary, the present study showed that Ni 2 does not form strong bonds with atomic configurations in which the 4s subshell is fully filled. This observation is consistent with other studies, by us and by others, of transition elements of the first row. Bonding in these systems is generally favored by atomic configurations that involve at least one of the participating atoms in an excited state (3d n+1 4s 1 ). For example, in our previous study (Tamukong et al., 2012) of the low-lying electronic states of Sc 2 , Cr 2 , and Mn 2 , the lowest states of Sc 2 were shown to correlate with the 2 D (3d 1 4s 2 ) + 4 F (3d 2 4s 1 ) dissociation asymptote, while those of Cr 2 correlated with the 7 S (3d 5 4s 1 ) + 7 S (3d 5 4s 1 ) dissociation limit. However, bonding in transition metal dimers is subtle and the most stable states of Mn 2 were obtained from weakly coupled 5 D (3d 5 4s 2 ) + 5 D (3d 5 4s 2 ) ground state Mn atoms, similar to the 3 F 4 (3d 8 4s 2 ) + 3 F 4 (3d 8 4s 2 ) coupling of ground state Ni atoms. Our results provide further evidence that GVVPT2, with sf-X2C treatment of relativistic effects, predict electronic excitation and bond energy trends in the first row transition metals consistent with experiment and the highest level ab initio calculations. In the present case, the quasidegeneracy of the 3 F and 3 D states of the Ni atom demonstrates that GVVPT2 can successfully be used for a system with an extraordinarily dense manifold of states, which generally requires computationally intensive variational treatments.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
MH and PT contributed to conception and design of the study. PT performed numerical calculations. MH and PT interpreted data and organized presentation. PT wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.