A comprehensive approach for elucidating the interplay between 4fⁿ⁺¹ and 4fⁿ5d¹ configurations in Ln²⁺ complexes

Abstract

Lanthanides (Ln) are typically found in the +3 oxidation state. However, in recent decades, their chemistry has been expanded to include the less stable +2 oxidation state across the entire series except promethium (Pm), facilitated by the coordination of ligands such as trimethylsilylcyclopentadienyl, C₅H₄SiMe₃ (Cp′). The complexes have been the workhorse for the synthesis and theoretical study of the fundamental aspects of divalent lanthanide chemistry, where experimental and computational evidence have suggested the existence of different ground state (GS) configurations, 4fⁿ⁺¹ or 4fⁿ5d¹, depending on the specific metal. Standard reduction potentials and 4fⁿ⁺¹ to 4fⁿ5d¹ promotion energies have been two factors usually considered to rationalize the occurrence of these variable GS configurations, however the driving force behind this phenomenon is still not clear. In this work we present a comprehensive theoretical approach to shed light on this matter using the [LnCp₃]⁻ model systems. We begin by calculating 4fⁿ⁺¹ to 4fⁿ5d¹ promotion energies and successfully correlate them with existing experimental data. Furthermore, we analyze how changes in the GS charge distribution between the Ln ions, LnCp₃ and the reduced [LnCp₃]⁻ complexes (Ln = La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu) correlate with experimental trends in redox potentials and the calculated promotion energies. For this purpose, a comprehensive theoretical work that includes relativistic ligand field density functional theory (LFDFT) and relativistic ab initio wavefunction methods was performed. This study will help the rational design of suitable environments to tune the different GS configurations as well as modulating the spectroscopic properties of new Ln²⁺ complexes.

---

1. Introduction

Lanthanide (Ln) chemistry is usually governed by the +3 oxidation state. Before the 1990s, just a small number of molecules containing metal ions in the +2 oxidation state were known for Sm, Eu, and Yb.¹ In these elements, reduction results in stable electron configurations where the 4f orbitals are nearly half-filled (Sm, 4f⁶), half-filled (Eu, 4f⁷) and fully filled (Yb, 4f¹⁴). Consequently, they exhibit the most accessible Ln³⁺/Ln²⁺ reduction potentials of the lanthanide series.² Over the last few decades, lanthanide (Ln) chemistry has been extended to the less stable +2 oxidation state in molecular compounds for the entire series by coordinating the metal center to ligands such as cryptates,³ cryptands,⁴ aryloxides,^5,6 crown ethers,^7–9 and cyclopentadienyl ([C₅H₅]⁻, Cp) derivatives.^8,10–17 The unique electronic configurations and chemical properties displayed by these reduced systems have been explored from a fundamental perspective leading to promising applications in the field of quantum information science.^18,19 However, further studies on divalent lanthanide complexes are essential to understand and optimize their electronic properties.

Among the diversity of Ln²⁺ complexes that have been reported to date, (Cp′ = trimethylsilylcyclopentadienyl, C₅H₄SiMe₃) systems are the most widely explored primarily because of the well-established protocols that allow for their synthesis and characterization. Experimental and theoretical evidence has shown that these compounds exhibit variable ground state (GS) configurations across the lanthanides series, which results in different spectroscopic features and variable bond lengths.^11–13 For example, a 4fⁿ⁺¹ GS configuration has been reported for Sm²⁺, Eu²⁺, Tm²⁺ and Yb²⁺ whereas the 4fⁿ5d¹ configuration has been seen for La²⁺, Ce²⁺, Pr²⁺, Gd²⁺, Tb²⁺, Ho²⁺, Er²⁺ and Lu²⁺.^11,12 Interestingly, Nd²⁺ and Dy²⁺ can switch configurations depending on the coordination environment for which they are considered “crossover ions”.^13,20 Density functional theory (DFT) has been used to rationalize the origin of these 4fⁿ5d¹ configurations, suggesting that the tris(cyclopentadienyl) ligand field, stabilizes the 5d_z² orbital relative to the 4f orbital shell making their population feasible.^11–13 Nonetheless, the driving force that determines the GS configuration is still not clear. In general, two main factors have been considered to get a better understanding of this phenomenon: the standard reduction potentials and the 4fⁿ⁺¹ to 4fⁿ5d¹ promotion energies for the Ln²⁺ free ions.

The standard reduction potentials for the Ln³⁺/Ln²⁺ couple have usually been the first approximation to rationalize the accessibility to the different Ln²⁺ ions.^1,2,21 Accordingly, several reduction potentials have been reported in literature (Fig. S1†). Over five decades ago, Nugent et al. extrapolated the values of the unknown Ln²⁺ ions using the experimental redox potentials of Eu²⁺, Yb²⁺, and Sm²⁺. Their redox potentials were correlated through their effective nuclear charge, Racah crystal field, spin–orbit coupling parameters, and 4f–5d absorption energies.^21–23 At the time, another set of reduction potentials derived from thermodynamic cycles were reported by Johnson in 1973.²⁴ Nearly a decade later, potentials based on correlations obtained from molten salt chemistry were also reported by Mikheev, Kamenskaya, and coworkers.^25–29 Despite the differences among these different redox potentials, they agree on the Eu²⁺, Yb²⁺, Sm²⁺, Tm²⁺, Dy²⁺, and Nd²⁺ values. This is not surprising because these divalent ions were characterized in the solid-state more than fifty years ago.

Some authors have looked for a correlation between the Ln³⁺/Ln²⁺ reduction potentials and 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies employing 4f–5d absorption energies for the Ln²⁺ ions.³⁰ A pioneering study published in 1963 by McClure and Kiss paved the way for this pursuit.²³ Firstly, they reported the absorption spectra for the Ln²⁺ series in cubic crystals of CaF₂, which was a valuable source of information since all the data was obtained through a systematic procedure. Then, they used this spectroscopic data, along with some theoretical approximations to calculate the 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies as the difference between the Russell Saunders ground states of the 4fⁿ⁺¹ and 4fⁿ⁻¹5d¹ configurations. Since significant deviations were detected in some lanthanides such as La²⁺, Ce²⁺, Gd²⁺ and Tb²⁺, new adjustments were incorporated to get more reasonable values. However, a deficient correspondence with the standard reduction potential curves was observed.

Several 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energy curves have been reported using a similar methodology as Kiss and McClure and have been employed for comparative purposes.^{10,29,31–33} If the energy values reported by Dorenbos et al. (Fig. S2†) are considered to analyze the energetics of the GS configurations in systems, only La²⁺ and Gd²⁺ exhibit a 4fⁿ5d¹ GS which suggests that there is an energy gap of around 20 000 cm⁻¹ that must be overcome in order to allow the 4fⁿ5d¹ GS configuration.³⁰ As mentioned by Fieser et al., this presumably relates to the crystal field stabilization provided by the tris(cyclopentadienyl) ligand set.¹³ This scheme also shows that Ho²⁺ and Er²⁺ are closer to be considered crossover ions than Dy²⁺ and Nd²⁺. However, experimental, and theoretical evidence does not support this observation which suggests that more factors should be considered to solve this puzzle.

A third parameter yet to be explored is the GS charge distribution of the lanthanide ions. The charge distribution is determined by the composition of the wave function in terms of projections (M_J) of the total angular momentum J that can be associated to a prolate, oblate, or spherical shape, based on the quadrupole moment Q. Therefore, Q > 0 indicates a prolate shape whereas Q < 0 an oblate charge distribution.³⁴ Under the assumption that the lanthanide ion is in its largest M_j state, equatorial ligand field environments such as those provided by Cp derivatives in LnCp₃/[LnCp₃]⁻ systems should destabilize the fⁿ⁺¹-states for oblate-shaped ions such as cerium or terbium. This agrees with the rules provided by Long and Rinehart, where this type of charge distributions minimizes the M_J value of their GS which means that a prolate shape is favored.³¹ The opposite is expected for prolate-shaped densities because an equatorial coordination environment is energetically more favorable. From this perspective, it is reasonable to think that variations in the preferred Ln charge distribution shape can be an important piece to explain the intricate behavior of the electronic structure of Ln²⁺ complexes.

In this work, we propose a general approach to elucidate the variable nature of the GS in divalent lanthanide complexes. We start revisiting well-known concepts in the field of Ln²⁺ chemistry such as 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies and its correlation with some standard reduction potentials. Then we analyze the Ln charge distribution shapes in terms of the M_J projections of J to determine if this is the missing piece that allows to explain the interplay between the two configurations observed in systems. For this purpose, an extensive theoretical work including the relativistic ligand field density functional theory (LFDFT) and relativistic ab initio methods was performed on Ln³⁺ free ions as well as in LnCp₃ and [LnCp₃]⁻ model systems (Fig. 1). The impact of this study will enable the rational design of suitable environments to tune the different GS configurations as well as modulating the spectroscopic properties of new Ln²⁺ complexes.

Fig. 1 (a) Structures of the complex and (b) the truncated model where the SiMe₃ groups were replaced with hydrogen atoms. See Table S3† for a comparison in the 4fⁿ⁺¹ to 4fⁿ5d¹ promotion energies when using both structures for Tm and Eu systems.

2. Results and discussion

2.1 Theoretical approach

There are several reports accounting for the 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies in Ln²⁺ systems and most of them were derived from lanthanide-doped inorganic matrices through experimental fittings and theoretical extrapolations.^23,32,35,36 To the best of our knowledge, there are no reports of 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies for the Ln²⁺ series obtained through systematic procedures. From the theoretical and computational perspective, this is understandable since it is well known that the electronic structure of Ln²⁺ complexes is difficult to predict. Although single-reference DFT calculations have contributed to the general understanding of the spectroscopic properties of some divalent lanthanide complexes, they usually fail in the determination of their electronic structure due to their lack of static correlation. Methods that account for this type of correlation, such as LFDFT and CASSCF, become more appropriate to analyze these elements.^37,38

The LFDFT approximation is a promising tool in the field of Ln²⁺ chemistry. It has proven accurate and efficient from a computational perspective, and it has been employed successfully in previous studies.^3,7,39,40 Nonetheless, due to the over delocalization of the molecular electron density intrinsic to DFT that significantly affects the description of the more polarizable and diffuse 5d shell in lanthanides, it requires significant fine tuning to obtain reliable results. On the other hand, the CASSCF method is a robust tool to analyze GS of multiconfigurational systems and it has proven valuable to analyze several Ln²⁺ complexes.^3,7,41 However, a couple of issues have been identified regarding the calculation of excited states. Firstly, a double f shell is needed to account for the so-called “radial correlation” that is particularly relevant in divalent lanthanides within this CASSCF formalism.³⁵ Additionally, the lack of dynamic correlation can introduce non-negligible errors in the prediction of the 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies. For this reason, second-order perturbation (PT2) methods are employed to solve this problem, but they are usually impractical to calculate molecular systems.^42,43 An alternative method that overcome this issue at a lower computational cost is the multiconfigurational pair-density DFT method (MC-pDFT).⁴⁴ This approximation combines the multiconfigurational approach with on-top density functionals to obtain electronic energies and has been successfully applied by Fleischauer et al. to predict the magnetic circular dichroism (MCD) spectra of some Ln²⁺ complexes.⁴⁵

In the first part of this work, we use the LFDFT and CASSCF in combination with CASPT2/MC-pDFT approximations to calculate, in a systematic way, the 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies of both Ln²⁺ free ions and [LnCp₃]⁻ complexes (Ln = La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb) to evaluate the correlation between the resulting curves and the aforementioned standard reduction potentials alongside previously reported 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies curves.

2.2 LFDFT and CASSCF/CASPT2/MC-pDFT 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies in Ln²⁺ free ions

LFDFT and CASSCF/CASPT2/MC-pDFT approximations show the same trend in the promotion energies through the lanthanide series which results in similar curves, as observed in Fig. 2a. As expected, Gd²⁺ is the only ion that possesses a 4f⁷5d¹ configuration, which is predicted to be around 2000–3000 cm⁻¹ more stabilized than the 4f⁸ configuration. Lanthanum and lutetium are also expected to exhibit 4f⁰5d¹ and 4f¹⁴5d¹ GS configurations, respectively. In La³⁺ the 5d orbital occupation is favored since the 5d shell is more stabilized than the 4f shell at the beginning of the lanthanide series while in Lu³⁺, a full-filled 4f shell is the responsible. However, neither method was able to recover a reasonable difference of energy between both configurations. With respect to the promotion energies reported by Dorenbos and coworkers (Fig. S2†).³⁰ resemblances are observed in the positions of Sm²⁺, Eu²⁺, Tm²⁺ and Yb²⁺. Nevertheless, significant differences are evidenced in Nd²⁺, Dy²⁺, Ho²⁺ and Er²⁺ which correspond to those ions close to the “crossover” region in the systems. According to these authors, the promotion energy decreases following the sequence Ho²⁺ > Er²⁺ > Dy²⁺ > Nd²⁺, contrary to our calculations where the promotion energy decreases following the sequence Dy²⁺ > Nd²⁺ > Ho²⁺ > Er²⁺. Interestingly, this order agrees with the curves of reduction potentials shown in Fig. 2a and S1,† where the best correspondence is with those values derived from thermodynamic data.

Fig. 2 (a) Correlation between the theoretical 4fⁿ⁺¹ to 4fⁿ5d¹ promotion energy curves and standard reduction potentials calculated from thermodynamic data. (b) Theoretical CASSCF/MC-pDFT 4fⁿ⁺¹ to 4fⁿ5d¹ promotion energies for the Ln²⁺ free ions and the [LnCp₃]⁻ complexes. The black dashed line indicates the crossover between both configurations.

2.3 CASSCF/CASPT2/MC-pDFT 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies in [LnCp₃]⁻

Once the ions are in the tris(cyclopentadienyl) coordination environment, the promotion energies are significantly lower. While the overall shape of the curve is retained, not all the Ln²⁺ ions are stabilized equally (Fig. 2b). In the case of the ions with a 4fⁿ⁺¹ GS, the stabilization is ca. 14 500 cm⁻¹. Conversely, ions possessing a 4fⁿ5d¹ GS show are stabilized by ca. 30 000 cm⁻¹ (vide infra). Below the black dashed-line that indicates the crossover between GS configurations, the promotion energy decreases following the same sequence as in the Ln²⁺ free ions: Dy²⁺ > Nd²⁺ > Ho²⁺ > Er²⁺. This result agrees with experimental evidence that classifies the Dy²⁺ and Nd²⁺ as crossover ions depending on the ligand environment surrounding the metal.

Since the experimental 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies have not been measured within this molecular framework, the Eu²⁺ ion was used as a reference. Given its 4f⁷ ground state, Eu²⁺ is the ideal candidate to experimentally measure the promotion energy between configurations since its value is equivalent to the lowest electronic transition, easily observable in an emission spectrum. This suggests that the is the most suitable system to measure the position of this band and to validate our theoretical approach. In Fig. 3, a comparison between the experimental and simulated emission spectrum of [EuCp₃]⁻ is shown, where a prominent band centered at 449 nm (22 271 cm⁻¹) is observed and ascribed to a 4f⁶5d¹ → 4f⁷ transition. Overall, the agreement between this and the simulated spectra is good, with errors in the range of 26 nm (1257.7 cm⁻¹) and 33 nm (1766.8 cm⁻¹) for CASSCF/MC-pDFT and CASSCF/CASPT2, respectively. Once we validated our calculations against this system, the rest of the series was analyzed.

Fig. 3 (a) Comparison between experimental and simulated emission spectra for the [EuCp₃]⁻ system calculated with CASSCF/MC-pDFT/SO (solid green line) and CASSCF/CASPT2/SO (dashed blue line). (b) To better compare the shape of both simulated and experimental spectra, the CASSCF/MC-pDFT/SO spectrum was blue-shifted ∼26 nm whereas the CASSCF/CASPT2/SO was red-shifted ∼33 nm to match the experimental peak.

It is noteworthy that historically Dy(II) and Nd(II) were thought to be 4fⁿ⁺¹ ions because they were the most stable in the +2 oxidation state along with Sm²⁺, Eu²⁺, Tm²⁺ and Yb²⁺. However, experimental evidence supporting a 4fⁿ5d¹ configuration for Nd²⁺ and Dy²⁺ has been available since 1976.^46–49 These precedents were based on structural comparisons, where small metal–ligand bond length differences going from trivalent to divalent species were observed, which is usually attributed to the localization of the additional electron in a d orbital. In the case of the metals coordinated to Cp derivatives, both and complexes display 4f³5d¹ and a 4f⁹5d¹ GS configurations, respectively. However, in 2018 Jenkins et al. reported the isolation of [LnCp^tet₃]⁻ (Ln = La, Ce, Pr, Nd, Sm, Gd, Tb, Dy; Cp^tet = C₅Me₄H) complexes across the lanthanide series, where all the systems excepting [SmCp^tet₃]⁻ and [DyCp^tet₃]⁻ exhibited a 4fⁿ5d¹ GS. This confers additional experimental validation to the curves shown in Fig. 2b as it shows that after the four well-known 4fⁿ⁺¹ ions, dysprosium is the next to display this GS configuration in a molecular environment. Furthermore, this suggests that the Cp′ can stabilize the 5d_z² orbital and/or destabilize the f-orbitals better than the Cp^tet ligand, thus making it more energetically available to accept the additional electron in Ln²⁺ ions closer to the crossover line. The reason behind this observation may be attributed to an increased energy of the 5d_z² orbital due to its potential interaction with occluding methyl groups of the Cp^tet with respect to the Cp′ ligands. This increases the electronic repulsion produced by the methyl groups of the Cp^tet rings, which results in a probable destabilization of the 5d_z² orbital thus making it less accessible and in consequence more difficult to populate.

The theoretical calculations of 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies for both the Ln²⁺ free ions and the [LnCp₃]⁻ complexes provide a consistent source of information that supports several experimental observations. Interestingly, the differential stabilization of the promotion energies observed upon coordination of the Cp₃ ligands correlates with the observed differences in the GS configurations in Ln²⁺ systems. However, it does not explain the fundamental reason on why lanthanides such as Dy and Nd are able to display both 4fⁿ⁺¹and 4fⁿ5d¹ configurations, depending on the coordinating ligands. To enhance our understanding of this phenomenon, we investigated how variations in charge distributions among [LnCp₃]⁻, LnCp₃, and free Ln ions correlate with the trends observed in their 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies.

2.4 Charge distributions of Ln single ions, [LnCp₃]⁻, and LnCp₃ model systems

The charge distribution of lanthanide GS is typically employed in magnetism to quantify the anisotropy of these ions and it is defined by prolate (axially elongated), oblate (equatorially expanded) or spherical shapes, based on its quadrupole moment. In 1981, Sievers reported the angular dependence of the charge distribution based on pure ± M_j states to visualize the asphericity of the 4f shells in their Russell Saunders ground states.⁵⁰ Under these conditions, the Ln(III) ions of Ce, Pr, Nd, Tb, Dy and Ho exhibited oblate shapes while the Ln(III) ions of Sm, Er, Tm and Yb exhibited prolate shapes for their maximal M_j projections. Gd(III) and Lu(III) are the only ions displaying a spherical shape given their 4f⁷ and 4f¹⁴ GS configurations. Following this classification, Long and Rinehart³¹ proposed that a judicious choice of the crystal field that surrounds the metal center can be done to maximize the GS anisotropy and minimize the electron repulsion between the intrinsic metal electron density and the ligand field. In consequence, two general ligand architectures can be used. In the case of ions with oblate-shaped densities for their maximal M_j projections, the anisotropy can be maximized by using axial ligands (above and below the xy plane, assuming that z is axial in the frame of reference). Conversely, for ions with prolate-shaped densities, equatorial ligands should be preferred.

This concept can be extended to further explore the complex interplay between GS configurations in divalent lanthanides systems. If we consider the simplest yet chemically intuitive model, it is reasonable to assume an ideal ground state where the anisotropy of the metal center is maximized (largest M_j value) in a favorable crystal field environment.^51,52 Though this is not necessarily the case for every system, it helps to simplify the discussion. By using this ideal model, it is possible to analyze sequential differences in charge density changes between Ln³⁺, Ln²⁺, LnCp₃ (model trivalent complex), and [LnCp₃]⁻ (model divalent complex). Interestingly, this sequence suggests that [LnCp₃]⁻ complexes displaying 4fⁿ⁺¹ GS configurations (Ln = Sm, Eu, Tm, Yb) retain a spherical (or nearly spherical) density. It is possible to infer that these Ln²⁺ ions show small variations with respect to the Ln³⁺ free-ion in their largest M_j state and the analogous trivalent systems. Conversely, ions displaying more pronounced variations in the shape of the metal density are consistent with a 4fⁿ5d¹ GS (Ln = La, Ce, Pr, Nd, Gd, Tb, Dy, Ho, Er, Lu). Although, Nd²⁺ and Dy²⁺, which can display both GS configurations depending on the ligand environment, should be considered intermediate cases. The results are rationalized in terms of the composition of the wavefunction, occupation of the natural orbitals, and the Long–Rinehart model described above (the details of these calculations can be found in the Computational methods section and ESI†).

2.4.1 Charge distribution of Ln ions. As previously mentioned, the Ln³⁺ free ions serve as a natural starting point to validate our theoretical approach and to analyze changes in the charge distribution (from the Hund's rule ground state in their maximal M_j state) upon coordination of the Cp₃ ligands (LnCp₃) and further chemical reduction ([LnCp₃]⁻). Initially, the ions were subjected to state-average CASSCF calculations, enforcing a spherical symmetry to obtain pure f-orbitals, where the occupation of the natural orbitals (NO) equivalent to the spin-free state with the largest M_j was examined (Table S4†).

Although the spin-free analysis is consistent with the model proposed by Sievers,⁵⁰ the effect of the spin–orbit (SO) coupling was also incorporated by state interactions between CASSCF wavefunctions to obtain the natural spin-orbitals (NSO). At this level, the mixing of spin-free natural orbitals driven by SO coupling can further deviate occupations from 0 or 1, as observed in Table 1 and Table S4.† This is shown by the Ce³⁺ ion where the electron in the largest M_j state now resides in an NSO that mixes both f_{y(3x²−y²)} and f_{x(x²−3y²)} spin-free orbitals featuring fractional occupations (0.43) with a maximum value of M_J = 5/2. Similarly, the rest of the series shows the same behavior of maximizing M_J in the absence of a ligand field (Table 1). Since the distribution of the electron density is defined by the nature of the occupied orbitals, its density can be determined and ascribed to oblate, prolate or spherical shapes. As observed in Fig. 4, Ln ions having an electron in the f_{y(3x²−y²)} or f_{x(x²−3y²)} orbitals (L = 3), such as Ce or Tb, display the most pronounced oblate shapes. Conversely, ions such as Sm or Tm have a more prolate shape due to the occupation of the f_(yz²), f_(xz²), and f_(z³) orbitals, which are axially elongated. Finally, in case of Eu, Gd and Yb, the electron density is more spherical, consistent with their almost half-filled, and full filled 4f shells. Overall, these results support reports by Sievers on the Ln³⁺ series.

Table 1 Ground state Natural Spin Orbital (NSO) occupations for Ln³⁺ free ions, LnCp₃ and [LnCp₃]⁻ complexes from the MC-pDFT calculations. The total angular momentum quantum number (J) along with its predominant M_j are specified. In Fig. S5 a depiction of the f-orbitals and their labels are shown

Ln	J	Mj	f3−	f2−	f1−	f0	f1+	f2+	f3+	d0
Ce^3+	5/2	5/2	0.43	0.08	0	0	0	0.08	0.43	0
CeCp3	5/2	1/2	0	0	0.33	0.34	0.33	0	0	0
[CeCp3]^−			0	0	0.06	0	0.06	0	0.92	0.98
Pr^3+	4	4	0.49	0.41	0.09	0.01	0.09	0.41	0.49	0
PrCp3	4	0	0.03	0.36	0.53	0.15	0.53	0.36	0.03	0
[PrCp3]^−	7/2	5/2	0.25	0.07	0.35	0.3	0.35	0.07	0.62	1
Nd^3+	9/2	9/2	0.50	0.49	0.41	0.21	0.41	0.49	0.50	0
NdCp3	9/2	5/2	0.15	0.31	0.56	0.36	0.56	0.31	0.77	0
[NdCp3]^−	4	3	0.11	0.31	0.49	0.44	0.5	0.33	0.82	1
Sm^3+	5/2	5/2	0.54	0.69	0.83	0.88	0.83	0.69	0.54	0
SmCp3	5/2	5/2	0.4	0.64	0.87	0.93	0.87	0.64	0.66	0
[SmCp3]^−	0	0	0.62	0.88	0.94	0.94	0.94	0.88	0.8	0
Eu^3+	0	0	0.86	0.86	0.86	0.86	0.86	0.86	0.86	0
EuCp3	0	0	0.61	0.88	0.95	0.94	0.95	0.88	0.79	0
[EuCp3]^−	—		1	1	1	1	1	1	1	0
Gd^3+	—	—	1	1	1	1	1	1	1	0
GdCp3	—	—	1	1	1	1	1	1	1	0
[GdCp3]^−	—	—	1	1	1	1	1	1	1	1
Tb^3+	6	6	1.49	1	1	1	1	1	1.49	0
TbCp3	6	0	0.99	1.04	1.2	1.43	1.3	1.04	0.99	0
[TbCp3]^−	13/2	13/2	1.45	1	1	1	1	1	1.55	1
Dy^3+	15/2	15/2	1.5	1.5	1.01	1	1.01	1.5	1.5	0
DyCp3	15/2	1/2	1	1.36	1.48	1.2	1.48	1.3	1.13	0
[DyCp3]^−	8	8	1.46	1.49	1.01	1	1.01	1.49	1.52	1
Ho^3+	8	8	1.5	1.5	1.5	1	1.5	1.5	1.5	0
HoCp3	8	2	1.11	1.36	1.6	1.29	1.6	1.36	1.69	0
[HoCp3]^−	17/2	17/2	1.39	1.5	1.5	1	1.5	1.5	1.61	1
Er^3+	15/2	15/2	1.5	1.50	1.5	2.00	1.5	1.50	1.50	0
ErCp3	15/2	11/2	1.36	1.5	1.81	1.27	1.86	1.5	1.7	0
[ErCp3]^−	8	3	1.27	1.61	1.48	1.52	1.81	1.4	1.91	1
Tm^3+	6	6	1.5	1.50	2	2	2	1.50	1.50	0
TmCp3	6	6	1.46	1.5	2	2	2	1.5	1.54	0
[TmCp3]^−	7/2	7/2	1.38	2	2	2	2	2	1.62	0
Yb^3+	7/2	7/2	1.54	1.96	2	2	2	1.96	1.54	0
YbCp3	7/2	7/2	1.7	1.8	2	2	2	1.8	1.7	0
[YbCp3]^−	—	—	2	2	2	2	2	2	2	0

---

Fig. 4 Classification of Ln ions based on their charge density distribution when going from the largest M_j state of the Ln³⁺ free ions to the trivalent and divalent complexes. The underlined metals Ce, Eu and Nd were taken as reference to plot the densities. Isovalue = 0.03.

From a convenient standpoint, we compared differences in 4f asphericities between Ln³⁺ and Ln²⁺ in their maximal M_j states (Table S5†). Quadrupole moments (Q₂) have shown useful to discuss differences in the shapes of the 4f electron densities. Q₂ values have been previously reported for trivalent ions but not for divalent lanthanides possibly due to the interplay of the 4fⁿ⁺¹ and 4fⁿ5d¹ configurations.^51,52 For comparisons, we have calculated Q₂ values of the Ln²⁺ for their maximal M_j states assuming a 4fⁿ⁺¹ configuration. Fig. S5† shows how lanthanide single ion shapes vary as they move from a trivalent to divalent ion. Interestingly, ions displaying highly anisotropic 4f densities correlate with 4fⁿ5d¹ GS configurations such as Ce, Gd, and Tb (except Tm). As expected, ions with spherical densities match the 4fⁿ⁺¹ GS configuration leaving the ions with decreased or low anisotropies as intermediate or crossover cases. The exception of Tm is discussed below.

2.4.2 Charge distribution of complexes. The perturbation of the electron distribution on the Ln³⁺ ions induced by the Cp₃ coordination can be rationalized through the rules proposed by Long and Rinehart.³¹ Because the Cp ligands coordinate equatorially imposing a pseudo-trigonal symmetry around the metal center (Fig. 5), it forces Ln³⁺ ions become prolate to reduce the electronic repulsion in the equatorial region (Fig. S3,† and Table 1). However, since not all the Ln display the same intrinsic density distribution, the energy “cost” involved to become prolate varies according to the ion. For instance, ions such as Ce or Tb experience the most pronounced variations in terms of shape, going from oblate to prolate. Unlike the free ion analogs, where the f_{y(3x²−y²)} or f_{x(x²−3y²)} orbitals (L = 3) were preferentially occupied and the M_j was maximized (5/2 and 6 for Ce and Tb, respectively), in the CeCp₃ and TbCp₃ complexes, the electrons are mostly distributed across the f_xz² and f_yz² (L = 1), and f_z³ (L = 0) orbitals thus minimizing the M_j value (1/2 and 0 for Ce and Tb, respectively) favoring a prolate shape (see Table S4†). This is not the case for ions such as Sm, Eu, Gd, Er, Tm and Yb which possess intrinsic prolate or spherical shapes. In these cases, the “cost” of coordinating the Cp₃ ligand set is lower. Finally, for ions such as Pr, Nd, Dy, Ho, and Er the GS shapes fall outside the general classification as oblate or prolate owing to a more complex angular dependence. Thus, these Ln³⁺ ions can be interpreted as intermediate cases.

Fig. 5 Qualitative comparison between spherical and pseudo-trigonal symmetries.

2.4.3 Charge distribution of [LnCp₃]⁻ complexes. When the LnCp₃ complexes are reduced, the additional electron may either occupy the 4f shell, resulting in a 4fⁿ⁺¹ GS configuration, or populate the 5d shell, leading to a 4fⁿ5d¹ GS configuration. Given the pseudo-trigonal symmetry imposed by the Cp₃ ligand set, the 5d_z² orbital is the most suitable to allocate the reducing electron. If the energy penalty of populating the 5d_z² orbital is larger than the 4f-intraelectron repulsion (destabilization of the 4f shell), the resulting GS should correspond to a 4fⁿ⁺¹. If the opposite is true, then the 4fⁿ5d¹ configuration will be the most stable. Therefore, to rationalize the occurrence of either configuration, the analysis is organized in groups of ions, based on their density distribution and promotion energies.
2.4.3.1 Sm, Eu, Tm, Yb. In their trivalent oxidation state, these ions possess 4f⁵, 4f⁶, 4f¹² and 4f¹³ GS configurations, respectively. Upon reduction, it would be reasonable to think that the additional electron can either occupy the low-lying 5d_z² orbital in cases like Sm and Yb or populate the 4f shell in the case of Eu or Yb. For example, and according to Fig. 2b, the promotion energies favoring the 4fⁿ⁺¹ configuration over the 4fⁿ5d¹ are ∼23 325 cm⁻¹ and 20 941 cm⁻¹ in Eu and Yb, respectively. This case is probably the simplest to rationalize as achieving a half- and full-filled shell is more energetically favorable than populating the d orbital. Conversely, while Sm and Tm do not achieve a half-filled shell, their 4f density seems to be nearly spherical in the [LnCp₃]⁻ system (Table 1 and Fig. S4†). This suggests that ions able to acquire a spherical or nearly spherical density shapes upon reduction will display a 4fⁿ⁺¹ configuration.
2.4.3.2 Ce, Tb, Gd. According to Fig. 2b, these ions possess the most stable 4fⁿ5d¹ GS configurations with promotion energies of ca. −15 524 cm⁻¹, −10 415 cm⁻¹, and −22 250 cm⁻¹ for Ce, Tb and Gd, respectively. In case of Ce and Tb, we observe that this correlates with the pronounced changes in the density distribution when going from the maximal M_j state (oblate) to (prolate) complexes (Table 1, Fig. 4 and S4†). For Gd, the reason is rather simple as the energy penalty to break the half-filled 4f⁷ configuration is greater than the energy required to populate the 5d_z² orbital. This is also why Gd is the Ln with the most stable 4fⁿ5d¹ configuration in the series (Fig. 2).
2.4.3.3 Pr, Nd, Dy, Ho, Er. In terms of 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies, these ions are close to the crossover line (Fig. 2b). Their 4fⁿ5d¹ GS configurations are found to be ∼2300–6300 cm⁻¹ more stable than the 4fⁿ⁺¹, suggesting that their GS configurations can be altered by modulating the ligand field. This correlates with their density distribution because these ions do not fall under the oblate/prolate classification. If contrasted with the single ion anisotropies, we could infer that because the less anisotropic densities in their divalent oxidation state the electron repulsion is reduced in a the crystal field compared to other more anisotropic Ln ions. They can be considered intermediate cases lying in between the rest of the Ln series. Therefore, upon coordination and reduction they do not undergo pronounced variations in terms of density shape.

3. Conclusions

In conclusion, we have presented an extensive and systematic theoretical study to rationalize experimental observations regarding the complex interplay of GS configurations of Ln²⁺ ions in a tris(cyclopentadienyl) molecular framework. Using ab initio wavefunction methods, we have calculated the 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies across the entire series of Ln²⁺ free ions along with the [LnCp₃]⁻ model complexes. Our results show a good correlation with Ln³⁺/Ln²⁺ reduction potentials along with available structural and spectroscopic data that assigns 4fⁿ⁺¹ GS to Eu, Sm, Tm and Yb, and 4fⁿ5d¹ GS to La, Ce, Pr, Nd, Gd, Tb, Dy, Ho and Er. Interestingly, unlike previous studies of 4fⁿ⁺¹ → 4fⁿ5d¹ promotion energies based on semi-empirical approximations, our results find Nd and Dy ions next in line to display the 4fⁿ⁺¹ configuration, consistent with previous reports that describe these ions as “crossover” due to its ability to display both GS configurations depending on the coordinating ligands.

We finally correlated the electronic structure of the Ln series with the spatial shape of the metal electron density in both Ln single ions as well as the LnCp₃/[LnCp₃]⁻ complexes to provide further explanation on the intricate interplay between GS configurations, particularly for ions close to the crossover line. Building upon the work of Sievers and Long,^31,50 we investigated changes in the shape of the density resulting from the coordination of the Cp ligands to the Ln³⁺ center (LnCp₃) as well as changes induced by the reduction of these systems ([LnCp₃]⁻). We demonstrate that Ln centers that retain near-spherical densities are more prompt to access 4fⁿ⁺¹ GS configurations, which is supported by a lower energy penalty in the 4f-intraelectron repulsion. Conversely, ions that display pronounced variation in their shapes, particularly going from oblate to prolate, prefer 4fⁿ5d¹ GS configurations. The case of Nd(II) and Dy(II) and other ions near the crossover line can be considered as intermediate, where the near-spherical density shape is not retained, but it remains accessible since they do not exhibit marked changes in their density. This leads us to think that Pr(II), Ho(II), and Er(II) could display 4fⁿ⁺¹ configurations if using a suitable crystal field that aligns with their shapes.

These calculations mark an initial step toward a deeper understanding of Ln²⁺ chemistry. To strengthen these findings, it is essential to extend the analysis to linear geometries and additional coordination frameworks, such as aryloxides or crown ethers, that can stabilize this oxidation state. Ultimately, this study will provide valuable insights into the rational design of environments to optimize various ground state configurations and modulate the spectroscopic properties of new Ln²⁺ complexes, thereby guiding and enhancing future experimental efforts.

4. Computational methods

Below a summary of the computational methods employed for these calculations is found. More details are given in the ESI.†

4.1 Geometry optimizations

The experimental coordinates of and (Ln = La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm and Yb) complexes were obtained from the experimental structures¹³ and optimized in ADF version 2019,⁵³ with no symmetry or geometry constraints at M06-2X/TZP level of theory.⁵⁴ The relativistic effects were incorporated through the zeroth-order regular approximation (ZORA) of the Dirac equation.⁵⁵ Finally, analytical frequency calculations were performed on the optimized structures to confirm them as true minima. Since these optimized structures are being subjected to expensive ab initio wavefunction calculations, the SiMe₃ groups were replaced with hydrogen atoms as shown in Fig. 1.

4.2 Wavefunction calculations

Both LnCp₃ and [LnCp₃]⁻ complexes along with the Ln²⁺ and Ln³⁺ free ions were subjected to ab initio wavefunction calculations using the OpenMolcas software.⁵⁶ The choice of theoretical approximation and active space was made depending on the lanthanide oxidation state. For Ln³⁺ and LnCp₃ systems, a minimum active space, nf-electrons in 7 orbitals, was calculated through the Complete Active Space Self-Consistent Field (CASSCF) method³⁸ (see details of the multiplicities and roots in Table S1†). Conversely, for the Ln²⁺ an extended active space including n electrons in 19 orbitals (4f shell + 5f shell + 5d shell orbitals) was calculated through the Restricted Active Space Self-Consistent Field (RASSCF)⁵⁷ to properly account for the radial correlation as defined by Seijo and Barandiarán (Table S02†).⁴¹ Finally in [LnCp₃]⁻ systems the active space consisted of n + 2 electrons in 17 orbitals (4f shell + 5f shell + 6s orbital + 5d_z² orbital). Dynamic correlation effects were calculated through MC-pDFT and CASPT2 approximations.^43,44 In a final step, the spin–orbit coupling (SOC) was included via state interactions between the CASSCF/RASSCF wavefunctions through the restricted active space state interaction (RASSI) method.⁵⁸

4.3 Ligand field DFT calculations

This approach was considered to obtain the 4fⁿ⁺¹ to 4fⁿ5d¹ promotion energies corresponding to the Ln³⁺ free ions. The calculations were performed using the ADF package version 2019.⁵³ The level of theory consisted of the PBE functional with Grimme corrections in its D3 formulation to account for the dispersion forces along with TZP basis sets for all the elements.⁵⁹ The relativistic effects were incorporated through the ZORA Hamiltonian. The multiplet manifolds of each Ln were calculated by diagonalizing the matrix elements of the LFDFT Hamiltonian.

Data availability

The data supporting this article have been included as part of the ESI.†

Author contributions

The study was conceived by M. J. B. L., C. C. B., T. E. A., and W. J. E. The theoretical calculations and analysis of the results was performed by M. J. B. L. and C. C. B., W. N. G. M., T. F. J., and W. J. E. contributed with the experimental portion. All the authors contributed to the writing and proofreading the manuscript.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, Biosciences and Heavy Elements Chemistry Program, under award number DE-SC0023693 (T. E. A.) and the U.S. National Science Foundation under award CHE-2154255 (W. J. E). M. J. B.-L. acknowledge the postdoctoral support provided by the Glenn T. Seaborg Institute. The authors thank Professor Zoila Barandiaran for her valuable input on the theoretical treatment of these systems and Valentina Villegas for her help designing the TOC. This work was partially executed at Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (contract no. 89233218CNA000001).

References

1. D. H. Woen and W. J. Evans, in Handbook on the Physics and Chemistry of Rare Earths, ed. J.-C. G. Bünzli and V. K. Pecharsky, Elsevier, 2016, vol. 50, pp. 337–394.
2. L. R. Morss, Chem. Rev., 1976, 76, 827–841.
3. T. N. Poe, M. J. Beltrán-Leiva, C. Celis-Barros, W. L. Nelson, J. M. Sperling, R. E. Baumbach, H. Ramanantoanina, M. Speldrich and T. E. Albrecht-Schönzart, Inorg. Chem., 2021, 60, 7815–7826.
4. S. S. Bokouende, D. N. Kulasekara, S. A. Worku, C. L. Ward, A. B. Kajjam, J. C. Lutter and M. J. Allen, Inorg. Chem., 2024, 63, 9434–9450.
5. L. M. Anderson-Sanchez, J. M. Yu, J. W. Ziller, F. Furche and W. J. Evans, Inorg. Chem., 2023, 62, 706–714.
6. S. A. Moehring, M. J. Beltrán-Leiva, D. Páez-Hernández, R. Arratia-Pérez, J. W. Ziller and W. J. Evans, Chem. – Eur. J., 2018, 24, 18059–18067.
7. T. N. Poe, S. Molinari, M. J. Beltran-Leiva, C. Celis-Barros, H. Ramanantoanina and T. E. Albrecht-Schönzart, Inorg. Chem., 2021, 60, 15196–15207.
8. D. N. Huh, J. W. Ziller and W. J. Evans, Dalton Trans., 2018, 47, 17285–17290.
9. T. N. Poe, S. Molinari, S. Justiniano, G. M. McLeod and T. E. Albrecht-Schönzart, Cryst. Growth Des., 2022, 22, 842–852.
10. J. K. Peterson, M. R. MacDonald, J. W. Ziller and W. J. Evans, Organometallics, 2013, 32, 2625–2631.
11. M. R. MacDonald, J. E. Bates, J. W. Ziller, F. Furche and W. J. Evans, J. Am. Chem. Soc., 2013, 135, 9857–9868.
12. M. R. MacDonald, J. E. Bates, M. E. Fieser, J. W. Ziller, F. Furche and W. J. Evans, J. Am. Chem. Soc., 2012, 134, 8420–8423.
13. M. E. Fieser, M. R. MacDonald, B. T. Krull, J. E. Bates, J. W. Ziller, F. Furche and W. J. Evans, J. Am. Chem. Soc., 2015, 137, 369–382.
14. P. B. Hitchcock, M. F. Lappert, L. Maron and A. V. Protchenko, Angew. Chem., Int. Ed., 2008, 47, 1488–1491.
15. R. E. MacKenzie, T. Hajdu, J. A. Seed, G. F. S. Whitehead, R. W. Adams, N. F. Chilton, D. Collison, E. J. L. McInnes and C. A. P. Goodwin, Chem. Sci., 2024, 15, 15160–15169.
16. K. R. McClain, C. A. Gould, D. A. Marchiori, H. Kwon, T. T. Nguyen, K. E. Rosenkoetter, D. Kuzmina, F. Tuna, R. D. Britt, J. R. Long and B. G. Harvey, J. Am. Chem. Soc., 2022, 144, 22193–22201.
17. P.-B. Jin, Q.-C. Luo, G. K. Gransbury, I. J. Vitorica-Yrezabal, T. Hajdu, I. Strashnov, E. J. L. McInnes, R. E. P. Winpenny, N. F. Chilton, D. P. Mills and Y.-Z. Zheng, J. Am. Chem. Soc., 2023, 145, 27993–28009.
18. K. Kundu, J. R. K. White, S. A. Moehring, J. M. Yu, J. W. Ziller, F. Furche, W. J. Evans and S. Hill, Nat. Chem., 2022, 14, 392–397.
19. P. W. Smith, J. Hrubý, W. J. Evans, S. Hill and S. G. Minasian, J. Am. Chem. Soc., 2024, 146, 5781–5785.
20. T. F. Jenkins, D. H. Woen, L. N. Mohanam, J. W. Ziller, F. Furche and W. J. Evans, Organometallics, 2018, 37, 3863–3873.
21. L. J. Nugent, R. D. Baybarz, J. L. Burnett and J. L. Ryan, J. Phys. Chem., 1973, 77, 1528–1539.
22. L. J. Nugent, in Lanthanides and Actinides, University Park Press, Baltimore, 2nd edn, 1975, vol. 7, pp. 195–219.
23. D. S. McClure and Z. Kiss, J. Chem. Phys., 1963, 39, 3251–3257.
24. D. A. Johnson, J. Chem. Soc., Dalton Trans., 1974, 1671–1675.
25. A. N. Kamenskaya, Russ. J. Inorg. Chem., 1984, 29, 251–258.
26. N. B. Mikheev and A. N. Kamenskaya, Coord. Chem. Rev., 1991, 109, 1–59.
27. N. B. Mikheev, L. N. Auerman, I. A. Rumer, A. N. Kamenskaya and M. Z. Kazakevich, Russ. Chem. Rev., 1992, 76, 827–841.
28. N. B. Mikheev, Inorg. Chim. Acta, 1984, 94, 241–248.
29. N. B. Mikheev, Russ. J. Inorg. Chem., 1984, 29, 258–263.
30. P. Dorenbos, J. Phys.: Condens. Matter, 2003, 15, 575.
31. J. D. Rinehart and J. R. Long, Chem. Sci., 2011, 2, 2078–2085.
32. L. Brewer, J. Opt. Soc. Am., 1971, 61, 1666–1682.
33. P. Dorenbos, J. Lumin., 2000, 91, 91–106.
34. J. Tang and P. Zhang, in Lanthanide Single Molecule Magnets, Springer HeidelbergNew York Dordrecht London, 2015, pp. 2–11.
35. E. Loh, Phys. Rev., 1968, 175, 533–536.
36. J. Rubio O., J. Phys. Chem. Solids, 1991, 52, 101–174.
37. H. Ramanantoanina, M. Sahnoun, A. Barbiero, M. Ferbinteanu and F. Cimpoesu, Phys. Chem. Chem. Phys., 2015, 17, 18547–18557.
38. B. O. Roos, P. R. Taylor and P. E. M. Sigbahn, Chem. Phys., 1980, 48, 157–173.
39. A. García-Fuente, F. Cimpoesu, H. Ramanantoanina, B. Herden, C. Daul, M. Suta, C. Wickleder and W. Urland, Chem. Phys. Lett., 2015, 622, 120–123.
40. H. Ramanantoanina, F. Cimpoesu, C. Göttel, M. Sahnoun, B. Herden, M. Suta, C. Wickleder, W. Urland and C. Daul, Inorg. Chem., 2015, 54, 8319–8326.
41. Z. Barandiarán and L. Seijo, J. Chem. Phys., 2013, 138, 074102.
42. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J.-P. Malrieu, J. Chem. Phys., 2001, 114, 10252–10264.
43. J. Finley, P.-Å. Malmqvist, B. O. Roos and L. Serrano-Andrés, Chem. Phys. Lett., 1998, 288, 299–306.
44. L. Gagliardi, D. G. Truhlar, G. Li Manni, R. K. Carlson, C. E. Hoyer and J. L. Bao, Acc. Chem. Res., 2017, 50, 66–73.
45. V. E. Fleischauer, G. Ganguly, D. H. Woen, N. J. Wolford, W. J. Evans, J. Autschbach and M. L. Neidig, Organometallics, 2019, 38, 3124–3131.
46. H. P. Beck, Z. Naturforsch. B, 1976, 31, 1548–1549.
47. F. Jaroschik, F. Nief, X.-F. Le Goff and L. Ricard, Organometallics, 2007, 26, 1123–1125.
48. F. Jaroschik, A. Momin, F. Nief, X.-F. Le Goff, G. B. Deacon and P. C. Junk, Angew. Chem., Int. Ed., 2009, 48, 1117–1121.
49. M. D. Walter, D. Bentz, F. Weber, O. Schmitt, G. Wolmershäuser and H. Sitzmann, New J. Chem., 2007, 31, 305–318.
50. J. Sievers, Z. Phys. B: Condens. Matter, 1982, 45, 289–296.
51. J. Tang and P. Zhang, Lanthanide Single Molecule Magnets, Springer Berlin, Heidelberg, 1st edn, 2015.
52. R. Skomski and D. J. Sellmyer, J. Rare Earths, 2009, 27, 675–679.
53. G. te Velde, F. M. Bickelhaupt, E. J. Baerend, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders and T. Ziegler, J. Comput. Chem., 2001, 22, 931–967.
54. E. Van Lenthe and E. J. Baerends, J. Comput. Chem., 2003, 24, 1142–1156.
55. E. van Lenthe, E. J. Baerends and J. G. Snijders, J. Chem. Phys., 1993, 99, 4597–4610.
56. G. Li Manni, I. F. Galván, A. Alavi, F. Aleotti, F. Aquilante, J. Autschbach, D. Avagliano, A. Baiardi, J. J. Bao, S. Battaglia, L. Birnoschi, A. Blanco-González, S. I. Bokarev, R. Broer, R. Cacciari, P. B. Calio, R. K. Carlson, R. Carvalho Couto, L. Cerdán, L. F. Chibotaru, N. F. Chilton, J. R. Church, I. Conti, S. Coriani, J. Cuéllar-Zuquin, R. E. Daoud, N. Dattani, P. Decleva, C. de Graaf, M. G. Delcey, L. De Vico, W. Dobrautz, S. S. Dong, R. Feng, N. Ferré, M. Filatov(Gulak), L. Gagliardi, M. Garavelli, L. González, Y. Guan, M. Guo, M. R. Hennefarth, M. R. Hermes, C. E. Hoyer, M. Huix-Rotllant, V. K. Jaiswal, A. Kaiser, D. S. Kaliakin, M. Khamesian, D. S. King, V. Kochetov, M. Krośnicki, A. A. Kumaar, E. D. Larsson, S. Lehtola, M.-B. Lepetit, H. Lischka, P. López Ríos, M. Lundberg, D. Ma, S. Mai, P. Marquetand, I. C. D. Merritt, F. Montorsi, M. Mörchen, A. Nenov, V. H. A. Nguyen, Y. Nishimoto, M. S. Oakley, M. Olivucci, M. Oppel, D. Padula, R. Pandharkar, Q. M. Phung, F. Plasser, G. Raggi, E. Rebolini, M. Reiher, I. Rivalta, D. Roca-Sanjuán, T. Romig, A. A. Safari, A. Sánchez-Mansilla, A. M. Sand, I. Schapiro, T. R. Scott, J. Segarra-Martí, F. Segatta, D.-C. Sergentu, P. Sharma, R. Shepard, Y. Shu, J. K. Staab, T. P. Straatsma, L. K. Sørensen, B. N. C. Tenorio, D. G. Truhlar, L. Ungur, M. Vacher, V. Veryazov, T. A. Voß, O. Weser, D. Wu, X. Yang, D. Yarkony, C. Zhou, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2023, 19, 6933–6991.
57. P. A. Malmqvist, A. Rendell and B. O. Roos, J. Phys. Chem., 1990, 94, 5477–5482.
58. P. Å. Malmqvist, B. O. Roos and B. Schimmelpfennig, Chem. Phys. Lett., 2002, 357, 230–240.
59. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sc05438e

---