Multiconfiguration Pair-Density Functional Theory for Transition Metal Silicide Bond Dissociation Energies, Bond Lengths, and State Orderings

Transition metal silicides are promising materials for improved electronic devices, and this motivates achieving a better understanding of transition metal bonds to silicon. Here we model the ground and excited state bond dissociations of VSi, NbSi, and TaSi using a complete active space (CAS) wave function and a separated-pair (SP) wave function combined with two post-self-consistent field techniques: complete active space with perturbation theory at second order and multiconfiguration pair-density functional theory. The SP approximation is a multiconfiguration self-consistent field method with a selection of configurations based on generalized valence bond theory without the perfect pairing approximation. For both CAS and SP, the active-space composition corresponds to the nominal correlated-participating-orbital scheme. The ground state and low-lying excited states are explored to predict the state ordering for each molecule, and potential energy curves are calculated for the ground state to compare to experiment. The experimental bond dissociation energies of the three diatomic molecules are predicted with eight on-top pair-density functionals with a typical error of 0.2 eV for a CAS wave function and a typical error of 0.3 eV for the SP approximation. We also provide a survey of the accuracy achieved by the SP and extended separated-pair approximations for a broader set of 25 transition metal–ligand bond dissociation energies.


Introduction
Transition metal silicides are useful in silicon-based devices due to their high thermal stability, low electric resistivity, and low density [1,2]. As electronics become smaller, the metal-silicon bond becomes increasingly important [3]. Studies involving transition-metalencapsulating Si clusters show that transition metals (specifically Zr, Nb, Mo, W) terminate the dangling bonds through selective interactions between the d and p orbitals [4], which suggests that focusing on the M-Si interaction could provide detailed understanding to tune properties of devices with such bonds. The nature of metal-silicide bonds has been explored experimentally by measuring bond dissociation energies (BDEs). Early measurements of BDEs for diatomic transition metal, lanthanide, and actinide silicides were done by Knudsen effusion mass spectrometry, which produces BDEs with errors in the range 0.1-0.9 eV [5]. To obtain more accurate BDEs, the Morse group has utilized a precise predissociation-based resonant two-photon ionization method [6][7][8][9] that has an accuracy of approximately 0.004 eV. This substantial increase in experimental accuracy produces benchmark values against which quantum chemistry methods can now be tested, and the accuracy attainable can be a significant indicator of whether the bonding is modelled adequately. This kind of test of theory is particularly interesting because prediction of the

Computational Methods
The computations were carried out using CASSCF and SP wave functions. In the SP method [29], the singly-occupied orbitals and their correlating orbitals are placed in a GAS subspace, and each pair of orbitals consisting of a doubly-occupied active orbital and its correlating orbitals are placed in separate two-electron, two-orbital subspaces. In this scheme, all intra-subspace electron excitations are allowed, but inter-subspace electron excitations are not allowed.
Previous work showed that the nom-CPO active space [16,20,21] provides an accurate treatment of TiSi [29], and so that scheme is adopted here. The nom-CPO active space includes all bonding, antibonding, and unpaired electrons and orbitals. For VSi, NbSi, and TaSi this results in a (7,10) active space. The active space for each molecule studied here includes the 2σ and 2σ * orbitals, the singly occupied 3σ orbital and its corresponding antibonding orbital, two sets of π bonding and antibonding orbitals, and a set of δ bonding and antibonding orbitals; pictures of the natural orbitals, including their occupation numbers are shown in Figures S1-S3 in the supporting information. Since VSi, NbSi and TaSi are all isoelectronic, the electron configurations corresponding to specific states are the same. Table 1 gives the dominant electron configuration and number of configuration state functions of both the CASSCF and SP active spaces with respect to the states considered in this work. This (7,10) nom-CPO active space is built from four orbitals in the a 1 symmetry: 2σ, 2σ * , 3σ, 3σ * , two orbitals in the b 2 symmetry: 1π x and 1π * x , two orbitals in the a 2 symmetry: 1δ x and 1δ * x , and two orbitals in the b 1 symmetry: 1π y and 1π * y . Seven electrons are distributed in these orbitals based on the dominant configuration of the state in question, as shown in Table 1. The CASSCF calculations use this active space without partitioning it, and the SP approximation partitions the CAS into four or five subspaces, where only low-lying states were considered with this approximation ( 4 Π, 6 Σ + , 2 ∆, 2 Σ + ). Each of these states has a unique set of subspaces, and they are listed in Table 2, along with the number of CSFs for each CASSCF and SP wave function used in this work. The large reduction in the amount of CSFs is appealing for computing properties of molecules that require prohibitively large CASSCF active spaces.
Potential energy curves were calculated for all three molecules with the OpenMolcas package [33] (version 18.09 tag 573-g0dff1fb) with CASSCF and SP reference wave functions and using CAS-PDFT, SP-PDFT, and CASPT2 with the ANO-RCC-VTZP [34,35] basis set. When calculating potential energy curves the natural orbitals from the previous geometry were used as an initial orbital guess to create a smooth curve. Scalar relativistic effects were included in all calculations with the second-order Douglas-Kroll-Hess Hamiltonian. For the CASPT2 computations, the ionization-potential-electron-affinity (IPEA) parameter was set to 0.25 a.u., and an imaginary shift of 0.25 a.u. was used to mitigate the presence of intruder states. Additionally, the following orbitals were frozen in the PT2 step: 1s and 2s for Si, up to 3s for V, up to 3d for Nb, and up to 5s for Ta. The CAS-PDFT and SP-PDFT potential energy curves were computed with a number of on-top pair-density functionals: tPBE, ftPBE, trevPBE, ftrevPBE, tBLYP, ftBLYP, tOreLYP, and ftOreLYP. [25,26,29] An in-house version of the OpenMolcas code was used for the tOreLYP and ftOreLYP on-top functionals. All MC-PDFT computations were done with an ultrafine grid size.
Heteronuclear diatomics have C ∞v symmetry, but OpenMolcas is only capable of handling Abelian symmetries, so the potential energy curves of the low-lying states investigated here were computed in C 2v symmetry (shown in parentheses after the correct C ∞v state symmetry): Σ + (A 1 ), Φ/Π (B 1 or B 2 ), and ∆ (A 2 ). The C 2v symmetry designations do not distinguish between Φ and Π states, so they are both labeled here as Π.
The VIBROT program in OpenMolcas was used to numerically solve the rovibrational spectra using the Numerov-Cooley method [36] with the potential energy curves interpolated using cubic splines. In this way we obtain the equilibrium bond lengths, zero point energies, and BDEs for all states considered. The bond dissociation energies, D 0 , were computed as: where D e is calculated as the difference in spin-orbit-free potential energy between the molecule at the equilibrium bond distance (r e ) and at the dissociation asymptote of the potential energy curve, SOC is the spin-orbit lowering (a negative number), and ZPE is the zero point energy. Taking ZPE and SOC into account allows us to directly compare to the experimental bond dissociation energies measured by Morse et al. [6]. The spin-orbit splittings for individual atoms, SOC(A) and SOC(B), are taken from experimental data [37] whereas the molecular spin-orbit splittings, SOC(AB), were computed as a difference between CASSCF energies and spin-orbit coupled restricted active space state-interaction [38] (RASSI-SO) energies at r e where the RASSI-SO energies were computed using 2 roots each of 6 Σ + , 4 Σ + , 2 Σ + , 6 ∆, 4 ∆, 2 ∆, 6 Π, 4 Π, and 2 Π states. This approach yields SOC(AB) results for VSi that are similar to those for other VX diatomics [39]. The SOC(AB) value is zero for NbSi due to the Σ ground state symmetry. The experimental D 0 , SOC(A), and SOC(B) values and the calculated SOC(AB) values are reported in Table 3.

Results and Discussion
The CASSCF dominant configurations for the ground state of each molecule are reported in Tables S1-S3. In the case of VSi and TaSi, the dominant configuration has a weight of 0.38 and three other configurations have weights between 0.20 and 0. 10. For NbSi, the dominant configuration has a weight of only 0.14, with two other configurations with weights above 0. 10. This shows that the ground state wave function of NbSi is more multiconfigurational in nature than VSi and TaSi.

State Ordering
Because the transition metal silicides have many low-lying electronic states, uncovering the true state ordering is a challenge for both experiment and computation. Here, we report the states of each molecule that lie within 1 eV of the ground state, as computed with CASSCF, CAS-tPBE, CAS-ftPBE, and CASPT2 for VSi, NbSi and TaSi. The excited states are computed as vertical excitation energies from the equilibrium geometry of the lowest-energy state. Figure 1 shows the state ordering for VSi for all states within 1 eV of the ground state. CASPT2 and CAS-tPBE predict 2 ∆ as the ground state, and CASSCF and CAS-ftPBE predict 4 Π as the ground state. The CASSCF results place 4 Π and 2 ∆ within 0.1 eV of each other, whereas these states are much closer together (0.01-0.03 eV) in the CASPT2, CAS-tPBE, and CAS-ftPBE results. These small energy differences are within the typical error ranges of the methods used. CASSCF predicts the 2 ∆ state to be higher in energy than the 6 Σ + state, whereas the post-SCF methods all place the 6 Σ + state 0.2 eV higher in energy than the ground state. All post-SCF methods predict quasidegeneracy between 2 ∆ and 4 Π states, 6 Σ + and 2 Σ + states, and 6 ∆ and 2 Π states, and the relative energy with respect to the ground state varies by up to 0.1 eV. The state ordering for NbSi is shown in Figure 2 for all states within 1 eV of the ground state. The CASSCF results give 4 Π as the ground state, with 2 ∆ and 4 ∆ states within 0.3 eV, and the 2 Σ + is about 0.32 eV higher than the ground state. The 6 Σ + state lies at 0.55 eV in the CASSCF results, but it is stabilized considerably in all post-SCF results such that these methods all report this state as the ground state. The CASPT2 results have the 6 Σ + state as the ground state, and the 4 Π state lies 0.05 eV above it. The 2 Σ + , 4 ∆, and 2 ∆ states all lie between 0.2 and 0.3 eV for CASPT2. The CAS-tPBE and CAS-ftPBE results give the same state ordering, where the 6 Σ + is the ground state, but the 4 Π state is close in energy at less than 0.1 eV above the ground state, and the 2 ∆, 4 ∆, and 2 Σ + states are clustered around 0.15-0.25 eV.
The state ordering for TaSi is shown in Figure 3. The CASSCF results give 4 Π as the ground state for TaSi, with 2 ∆ and 6 Σ + states about 0.2 eV and 0.25 eV higher in energy, respectively. The ground state is the same ( 4 Π) for all methods in Figure 3; however the excited state ordering varies slightly for the post-SCF methods. In the case of CASPT2, the 2 Σ + state is 0.11 eV above the ground state, and the 2 ∆ and 6 Σ + states are higher at 0.15 eV and 0.19 eV, respectively. For CAS-tPBE, the 2 Σ + state is the first excited state at 0.1 eV, the 2 Π state lies at 0.18 eV, and the 2 ∆ and 6 Σ + are degenerate at 0.23 eV above the ground state. The results from CAS-ftPBE show the 2 Σ + state is degenerate with the 2 Π state just below 0.1 eV, and the 2 ∆ state and 6 Σ + states lie close together at 0.22 and 0.25 eV, respectively. The other states are all above 0.5 eV and the ordering is the same for all methods.  For TaSi, the 4 Π state is clearly predicted to be the ground state, but for VSi and NbSi, no definitive conclusion can be drawn at this stage. A post-SCF treatment gives quasidegenerate 4 Π and 2 ∆ states as potential ground states for VSi, and a 6 Σ + ground state for NbSi (with a low-lying 4 Π state below 0.1 eV). These states will be used in the comparison to experimental bond dissociation energies in following sections.

Vanadium Silicide
Tables 4-7 give the BDEs (D 0 ), equilibrium bond lengths (r e ), and zero point energies (in cm −1 ) for the 2 ∆ and 4 Π states of VSi. Both states are included in the discussion since they are almost degenerate. The computed D 0 are compared to the experimental VSi BDE of 2.234(2) eV; the difference between computed BDE and experiment is reported as ∆D 0 . The last rows of the tables give the mean unsigned error (MUE)in the MC-PDFT results, averaged over eight on-top functionals. Table 4 shows the CASSCF results for the 4 Π state. The CASSCF BDE underestimates experiment by almost 1.5 eV, but both CASPT2 and CAS-PDFT provide a BDE that is closer to that of experiment. All pair-density functionals have an error less than 0.2 eV; CASPT2 underestimates D 0 by 0.04 eV, whereas CAS-tPBE overestimates it by 0.04 eV. The best pair density functionals in this case, CAS-ftPBE and CAS-tOreLYP, are both within 0.01 eV of experiment. Previous results computed with B3LYP [6] provide a BDE that overestimates the experimental value by 0.31 eV, which is almost double the error we get with some post-SCF methods. Additionally, the equilibrium bond distance is much longer than what is predicted by CASPT2 and CAS-PDFT methods in Table 4. Similarly to Table 4, the SP-PDFT results are shown in Table 5 with the ∆D 0 value computed with respect to the experimental value. Like CASSCF, SP also underestimates D 0 , but the SP difference from experiment is smaller. The SP-tPBE are SP-ftPBE differ from experiment by less than 0.04 eV. This implies that the SP scheme for selecting SP gives similar results to the much larger CASSCF active space of the same electron and orbital composition at both the equilibrium region and at dissociation. The results for the 2 ∆ state are similar to those for the 4 Π state and are reported in Table 6 for the CASSCF reference and Table 7 for the SP reference. The CASPT2, CAS-tPBE, CAS-ftPBE, CAS-trevPBE and CAS-ftrevPBE bond dissociation energies are all within 0.1 eV of the experiment. The SP-based results are similar with CAS-tPBE and CAS-ftPBE being only 0.06 eV and 0.01 eV below experiment.
The potential energy curves for the 4 Π state of VSi are shown in Figure 4. The CASSCF curve does not give a good description of the potential energy curve, as determined by the large ∆D 0 value in Table 4. The description is improved when post-SCF treatments are used; both of the shown pair-density functionals give similarly shaped curves to the orange CASPT2 one. There is a small discontinuity in the CASSCF-based results near 2.7 Å but the SP-based results are smoother. Like CASSCF, the SP wave function does not describe this process well, but SP-PDFT corrects for the poor description and yields a result similar to that of CASPT2.   The potential energy curves for VSi in the 2 ∆ state are shown in Figure 5. The CAS-PDFT results give a deeper potential than CASPT2, generating a slightly larger D e value. The SP-tPBE and SP-ftPBE calculations are much closer to the CASPT2 curve than are the CAS-PDFT results.  Even though the 2 ∆ state is reported as lowest energy state at equilibrium for most methods, the computed 4 Π bond dissociation energies match more closely with experiment on average.

Niobium Silicide
The results for NbSi are given in Tables 8 and 9. Table 8 gives CASSCF results for the 4 Π state and post-CASSCF results for the 6 Σ + state (these are the predictions of the ground state, as discussed above). The BDE computed with CASSCF underestimates experiment by 1 eV, but all post-SCF methods are much more accurate. The CASPT2 results are 0. 23 Table 9 contains BDEs, equilibrium bond distances, and zero point energy computed with an SP reference wave function for the ground state of NbSi. The best functionals for computing the BDE for NbSi are tPBE and ftPBE, which underestimate the experimental BDE by 0.26 and 0.22 eV, respectively. The equilibrium bond distances are very similar to what is reported in Table 8 with CAS-PDFT; the results with different reference wave functions but the same on-top functional are within 0.02 Å for equilibrium bond lengths.
The potential energy curves of NbSi in the 4 Π state are shown in Figure 6. The CASSCF-based results (left) show the CAS-PDFT potential energy curves track closely to CASPT2 at longer bond distances, however the CASPT2 curve is slightly shallower near equilibrium. These potential energy curves show a small discontinuity near 3.5 Å, and similarly in the SP-based potential energy curves (right), there is a small discontinuity near the equilibrium bond distance that is not present in the CASSCF or CASPT2 curve. This discontinuity is likely due to low-lying excited states present in the dense electronic excited state spectrum. The SP-PDFT curve is shifted to a lower equilibrium bond distance than CASPT2, which is also the case for CAS-PDFT.  The NbSi potential energy curves in the 6 Σ + state are shown in Figure 7. These curves are not as smooth as the ones for the 4 Π state. The CASPT2 and CAS-PDFT curves all have a discontinuity near 4 Å. The SP-based results suffer from more severe discontinuities and would probably benefit from a state-averaged treatment in future work. Although the potential curves are not smooth, the BDEs computed with six of the eight CAS-PDFT methods and with SP-tPBE and SP-ftPBE are all within 0.26 eV of the experimental value.

Tantalum Silicide
The BDE, equilibrium bond length, and zero point energy for the 4 Π state of TaSi computed with a CASSCF reference wave function are reported in Table 10. CASPT2 yields a slightly shorter equilibrium bond length than any of the CAS-PDFT calculations, and the BDE is 0.15 eV below the experimental value. The BDE results closest to experiment are those computed with CAS-tPBE and CAS-ftPBE, with deviations from experiment of 0.07 eV and 0.10 eV; these methods are closer to experiment than CASPT2. Other CAS-PDFT functionals have higher deviations, with tBLYP and ftBLYP differing from the measured BDE by about a half eV.  Table 11 reports the SP-based results for the dissociation of TaSi. Like that of CASSCF, the SP bond dissociation energy is too small, but most of this error is removed by CAS-PDFT. This is especially the case with SP-tPBE and SP-ftPBE, where the error with respect to experiment is on the same order as CAS-tPBE and CAS-ftPBE (approximately 0.1 eV). In this case, the B3LYP equilibrium bond length is similar to those computed with CAS-PDFT and CASPT2, but the BDE error in B3LYP is larger. Potential energy curves of TaSi computed with CASPT2, CASSCF, CAS-tPBE, and CAS-ftPBE are shown on the left side of Figure 8, and those computed with SP reference wave functions are on the right (with CASPT2 as a reference). Neither reference wave function gives an accurate potential energy curve, which is apparent in the large ∆D 0 value. However, both CAS-PDFT and SP-PDFT with the tPBE and ftPBE on-top functionals give BDEs very close to that of CASPT2, although the equilibrium distance is slightly longer. This trend was also seen for VSi where CASPT2 predicts a slightly shorter r e than CAS-PDFT or SP-tPBE by 0.02 Å.  Overall, the typical error involved in calculating the ground state BDE of VSi, NbSi, or TaSi with a nom-CPO active space with CAS-PDFT is 0.2 eV and with SP-PDFT is 0.3 eV.

Bond Dissociation Energies
Studying diatomic transition metal compounds provides an opportunity to affordably benchmark the SP approximation against CASSCF for active spaces of the same size, as well as to compare both methods to experiment. Here, we compare the accuracy of using CAS-PDFT and SP-PDFT on a set of 25 diatomic molecules containing one or two transition metal atoms. This includes the three molecules studied in the present paper (VSi, NbSi, TaSi) and 22 studied previously [16,17,29]. All 25 data points in this comparison are transition metal diatomic BDEs calculated with a nom-CPO active space [16,17,29]. For two of the molecules the SP method was not adequate, and in those cases we used the extended separated-pair (ESP) method [29]. In summarizing the accuracy for 25 molecules, we use the results for the ESP method when ESP was applied and the results for SP when ESP was not performed. Previous work with the SP approximation includes CrH, MnH, FeH, CoH, FeC, ScN, VN, CrN, TiO, FeO, NiO, V 2 , Cr 2 , CrF, NiC, FeS, NiS, FeSe, NiSe,and TiSi [16,17,29], and previous work involving the ESP approximation includes TiC and WCl [29].
The results of the survey of 25 diatomic transition metal compounds are shown in Table 12 The good results obtained with the ftOreLYP functional are very satisfying since this functional is based on full translation [26] of the OptX exchange functional developed by Handy and Cohen [40] with the goal of including both nonlocal exchange and leftright correlation energy and the reLYP correlation functional developed by Thakkar and McCarthy [41] with the goal of making a functional more accurate for the correlation energies of atoms with atomic numbers 20-86 (the transition metals considered here have atomic numbers 23, 41, and 73).

Conclusions
In this study, we report calculations of the ground state state symmetry, bond dissociation energy, equilibrium internuclear distance, and zero point energy for three transition metal silicides. The SP-PDFT bond dissociation energies are very similar to CAS-PDFT results, showing that the separated-pair approximation is a good way to approximate a CASSCF wave function while reducing the number of configuration state functions. The MC-PDFT functionals tPBE and ftPBE give bond dissociation energies that are just as accurate, or better than, the CASPT2 ones while also keeping the computational cost low. Additionally, there is, on average, a twofold improvement over previous Kohn-Sham density functional results. As the transition metal becomes heavier, the MUE increases slightly, but overall, the CAS-PDFT results have a MUE of 0.2 eV and SP-PDFT results have a MUE of 0.3 eV compared to experimental bond dissociation energies.
Comparing calculated bond dissociation energies to experiment is insufficient to reliably determine which state is the ground state for the transition metal silicides studied here because there are low-lying excited states. In the case of VSi, where the ground state 2 ∆ state is nearly degenerate with the 4 Π state, the ordering is similar for all post-SCF methods (CASPT2, CAS-PDFT, and SP-PDFT). For NbSi, some methods yield 6 Σ + as the ground state, and others yield 4 Π. The state ordering for NbSi trend with post-SCF methods is: 6 Σ + < 4 Π < 4 ∆ ≈ 2 ∆ < 2 Π. All results for TaSi are in agreement that the ground state is 4 Π, and the closest excited state is near 0.1 eV, but the ordering of the excited states varies with the method used.
We also combined the present results for three transition metal diatomics with previous results for 22 other transition metal diatomics to get a broader assessment of the accuracy attainable with MC-PDFT. The mean unsigned errors results are very similar for CAS-PDFT and SP-PDFT, and they are very similar for simply translated functionals and for fully translated ones; if we therefore combine all these into a single data set, it has 400 data, and the mean unsigned error in transition metal bond dissociation energies for this data set is only 0.27 eV. The best method, SP-PDFT with the ftOreLYP functional, has a mean unsigned error of only 0.21 eV.