Assessing Quantum Calculation Methods for the Account of Ligand Field in Lanthanide Compounds

Abstract

We obtained thorough insight into the capabilities of various computational methods to account for the ligand field (LF) regime in lanthanide compounds, namely, a weakly perturbed ionic body and quasidegenerate orbital multiplets. The LF version of the angular overlap model (AOM) was considered. We intentionally took very simple idealized systems, the hypothetical [TbF]²⁺, [TbF₂]⁺ and [Tb(O₂NO)]²⁺, in order to explore the details overlooked in applications on complex realistic systems. We examined the 4f and 5d orbital functions in connection to f–f and f–d transitions in the frame of the two large classes of quantum chemical methods: wave function theory (WFT) and density functional theory (DFT). WFT methods are better suited to the LF paradigm. In lanthanide compounds, DFT faces intrinsic limitations because of the frequent occurrence of quasidegenerate ground states. Such difficulties can be partly encompassed by the nonstandard control of orbital occupation schemes. Surprisingly, we found that the simplest crystal field electrostatic approximation, reconsidered with modern basis sets, works well for LF parameters in ionic lanthanide systems. We debated the largely overlooked holohedrization effect that inserts artificial inversion symmetry into standard LF Hamiltonians.

1. Introduction

In this study, we considered computational experiments in different technical settings, investigating the title problem and retaining the conceptual beacon ideas of ligand field (LF) theories. The LF represents a class of phenomenological models, proposing parameters for the split of atomic orbitals (d for transition metal compounds [1] or f-type for lanthanide ions [2]) in specific environments from molecules and crystals. Initially, in historical versions, due to Bethe and van Vleck [3,4], it was thought that perturbation can be estimated from first principles if it is considered to be mostly electrostatic. Practice, acquired mostly from the spectroscopy of transition metal ions [5], showed that the electrostatic simple hypothesis did not work well, but the LF models were fruitful after making all parameters adjustable, fitted from experimental spectra. The term ‘ligand field’ is also encountered under another name, ‘crystal field’ (CF) [2]. Here, we would like to use the term ‘CF’ for situations when the parameters are taken by electrostatic models, similar to early ideas.

Long before the computer era of molecular or solid-state electronic structure codes, the LF paradigm was the sole route to the examination of the quantum effects of d or f electrons. Currently, the LF is used as a tool for interpreting and validating the results of brute-force calculations [6,7].

The LF Hamiltonian in the Wybourne formulation [8] is as follows:

$${\widehat{V}}_{LF}={\displaystyle \sum }_{k=2,4\dots }^{2l}\sqrt{\frac{4\pi }{2k+1}}{B}_0^k{Y}_{k,0}+{\displaystyle \sum }_{k=2,4\dots }^{2l}\sqrt{\frac{4\pi }{2k+1}}{\displaystyle \sum }_{q=1}^k\left[B_q^k\left(Y_{k,q}+{\left(-1\right)}^q{Y}_{k,-q}\right)+i{B}_{-q}^k\left(Y_{k,-q}-{\left(-1\right)}^q{Y}_{k,q}\right)\right]$$

(1)

Here, Y_k,q = Y_k,q(θ, φ) are spherical harmonics, the explicit dependence on polar angles being dropped in the above description. The LF potential operator takes a matrix representation on the basis of the considered atomic orbitals (i.e., l = 2 for transition metals, or l = 3 for lanthanides):

$$H_{m,m\prime }^{\left[l\right]}=\langle Y_{l,m}|{\widehat{V}}_{LF}|Y_{l,m^{\prime }}\rangle ={{\displaystyle \int }}_{\theta =0}^{\pi }{{\displaystyle \int }}_{\phi =0}^{2\pi }{Y}_{l,m}^{*}\left(\theta ,\phi \right){\widehat{V}}_{LF}\left(\theta ,\phi \right)Y_{l,m^{\prime }}\left(\theta ,\phi \right)\sin \left(\theta \right)d\theta d\phi$$

(2)

Instead of complex spherical harmonics, we can work with the real angular parts of the orbitals, converting the $Y_{l,m}^{}$ and $Y_{l,m}^{*}$ conjugate couples into the real sine and cosine components, $Y_{l,m}^s$ and $Y_{l,m}^c$, respectively.

Designs with spherical harmonics are convenient because these Hamiltonian matrices must be integrated with atomic orbitals (AOs) that pertain to a shell with a given secondary quantum number, l. Then, the integrals from Equation (2) contain Y^*_l,m·Y_k,q·Y_l,m’ products, with the left- and right-side factors coming from AO definitions (the “bra” and “ket” functions); the middle factor comes from the Hamiltonian matrix. The use of even k indices, k = 0, 2, … 2l is determined via algebraic rules in order to have nonvanishing integrals [9]. The k = 0 component, which causes an overall shift in the diagonal of the LF Hamiltonian matrix and the resulted eigenvalues, is omitted from Equation (1), where the summation starts at k = 2. This means that the Hamiltonian representation is traceless, with the sum of diagonal elements and eigenvalues being zero, i.e., having the energy spectrum shifted in its barycenter, as a matter of convention.

The somewhat complicated expansion of Equation (1) pays tribute to the historical development of this field. Otherwise, the above formulation is quite trivial because it generally proposes just as many parameters as in a Hamiltonian matrix written for a given l shell. For instance, with f orbitals, the LF matrix has 7 × 7 dimensions and 7·8/2 = 28 general elements (i.e., taking the diagonal and half of the matrix from below or above the diagonal because the transposed elements are the complex conjugate of the selected set). Making the matrix traceless reduces the general degrees of freedom to 27, in the absence of any symmetry (C₁ point group). In Equation (1), the k = 2, 4 and 6 sums comprise 2k + 1 elements, having the 5, 9, and 13 lengths, respectively, leading to a total of 27 general $B_q^k$ parameters. When the problem has symmetry, the number of effective parameters is reduced, imposing certain ratios between the $B_q^k$ values and determining the vanishing of the other ones.

We do not discuss the subtleties of Wybourne-type LF Hamiltonians here but simply state that the number of expected independent parameters is equal to the number of totally symmetric representations selected from the reduction of k = 2, … 2l spherical harmonic sets in the given point group. This is a basic requirement because a well-designed effective Hamiltonian should remain invariant to group operations. As an example, let us consider the O_h point group. Here, k = 2 (quadrupole) behaves as E_g + T_2g, yielding no contribution to the LF operator because no A₁ is included. k = 4 (hexadecapole) yields A_1g + E_g + T_1g + T_2g; k = 6 (64-pole) corresponds to A_1g + A_2g + E_g + T_1g + 2T_2g. These two sets contribute equally, with an independent parameter [10]. Here, there is an established proportionality between $B_q^4$ parameters such that only $B_0^4$ can be taken as an independent term. Similarly, the k = 6 set ends with only two LF parameters for the f shell in octahedral symmetry, $B_0^4$ and $B_0^6$. For the d shell, which is confined to the k = 4 set, there is only one LF parameter, $B_0^4$, proportional to the well-known 10Dq gap. In the absence of symmetry, the active k-type spherical harmonic operators are trivially decomposed in (2k + 1)·A₁ totally symmetric elements, making each $B_q^k$ an independent parameter.

There are several other versions of the LF beyond the above-described Wybourne parameterization, a very intuitive one being the angular overlap model (AOM) [11,12,13,14,15,16,17,18,19]. The AOM was enthusiastically embraced by chemists dealing with d-type coordination chemistry and spectroscopy given its transparency virtues. For the f electrons, it was pioneered by W. Urland [20], which remained lesser known because of the more complicated algebraic apparatus and confinement to the somewhat exotic lanthanide chemistry. The AOM does not make the Hamiltonian operator explicit, directly parameterizing the LF matrix with the help of ${\mathrm{e}}_{\lambda }^L$ parameters assigned to the ligand L and bonding type (λ = σ, π, …):

$$H_{m,m\prime }^{\left[l\right]}={\displaystyle \sum }_L^{ligands}{\displaystyle \sum }_{\lambda =\sigma ,{\pi }_{\perp },{\pi }_{\left|\right|}}^k{D}_{m\lambda }^{\left[l\right]}\left({\theta }_L,{\varphi }_L\right)D_{m\prime \lambda }^{f\left[l\right]}{\left({\theta }_L,{\varphi }_L\right)}_{}{\mathrm{e}}_{\lambda }^L$$

(3)

where the D values are functions of φ_L and θ_L polar coordinates of the coordinated ligand (in a frame centered on the metal ion), having well-defined expressions for the l-type shell to which the LF problem refers. An alternative model was advanced by Malta [21], mixing AOM-like terms assigned to chemical bonding with electrostatic terms. Recent works and reviews [22,23] have proven the current interest in AOM on lanthanide compounds in connection with modern applications of these materials in new light-emission technologies [24,25,26].

Because we revisit the LF ideas with modern tools, the following part of the Introductory is devoted to state-of-the-art computational chemistry. There are two distinct classes of theoretical approaches: wave function theory (WFT) [27] and density functional theory (DFT) [28]. For the following discussion, the former is closer to the physics of LF split, because the spectral terms encountered in spectroscopy and magnetism have a multiconfigurational and quasidegenerate nature, tractable with complete active space self-consistent field (CASSCF) [29] (sometimes followed by second-order perturbations’ correction (PT2) [30,31,32]) and spin orbit (SO) [33,34] coupling. On the contrary, DFT is limited, by the grounding principles, to nondegenerate states and, via Kohn–Sham [35,36] techniques, to single-determinant functions. Then, we must adopt certain techniques to circumvent the limitations of the DFT.

Although conceptually simple, because mononuclear lanthanide compounds, having one ion with fⁿ configuration, can be treated with a CASSCF(n,7) keyword, meaning an active space of seven orbitals lodging n electrons, the technical realization may be difficult because the calculations must be initiated with well-controlled guess orbitals, having a pure or an almost-pure atomic orbital (AO) f characteristic [37]. Converged molecular orbitals (MOs) are expected to retain quasiatomic f preponderance. In this study, e applied the guidelines outlined in Chapters 5 and 7 of the structural chemistry textbook by Putz, Ferbinteanu and Cimpoesu [38] and in debated studies [39,40], where these authors pioneered the use of CASSCF for realistic lanthanide compounds. In detail, the strategy consists of performing separate calculations of lanthanide fragment (CASSCF on naked ion) and single-determinant procedures on ligands, merging the resulting orbitals in a common block, initially placing zero nondiagonal values for the inter-fragment mixing coefficients. It is advisable to perform the LF calculations in state-averaged mode, including a number of states equal to the degeneracy of the ground state of the free lanthanide ion. In the case of polynuclear systems, the WFT procedures are more complicated but still similarly possible. In a previous instance, one of us applied this computational protocol to magnetochemical problems in a series of lanthanide complexes, now revisiting the methodology in its basic aspects in this study [41].

Remarkable results in the WFT electronic structure of lanthanide compounds were reported by Neese and Atanasov [42] in their code, Orca [43,44], which has a set of keywords dedicated to the interpretation of ligand field theory [45]. Because the final WFT solutions retain a small mixing of lanthanide f-type orbitals, the guess from merged fragments is well suited as a start. Otherwise, the common practice would be to trigger a CASSCF calculation by the orbitals obtained, in the preamble, through lower-level single-determinant calculations. However, this routine is faulty in lanthanide complexes, with such attempts facing severe convergence problems.

Although DFT calculations on lanthanide compounds can converge in routine mode using an unrestricted frame (i.e., with independent sets of molecular orbitals for α and β electrons), the results may break the physical meaning, being particularly incompatible with the LF phenomenology [46]. In nonroutine modes, we may obtain meaningful results, with a rational method being the use of the fractional population n/7 on the seven molecular orbitals akin to the f-atomic functions (in the case of mononuclears with fⁿ lanthanide ions) [47,48,49,50,51,52,53,54,55,56]. This somewhat resembles the state-averaged procedures from CASSCF, emulating an artificial object having a well-defined single-determinant state in DFT, avoiding the spurious results obtained from the orbital optimization of a single-determinant configuration accidentally selected among the quasidegenerate configurations possible, in general, for a lanthanide compound. The population-averaged DFT also enables the use of the restricted mode (i.e., a unique orbital set for both α and β electrons), a situation more compatible with the LF-oriented interpretation.

The above-cited works produced good results in terms of electronic structure in both WFT [37,39,41] and DFT [47] branches, a directly approaching compounds with a realistic level of complexity, as demanded by the debated practical problems. Here, we revisit these methods on extremely simple systems, aiming to critically examine their capabilities to obtain systematic insight that, to the best of our knowledge, was not realized until now given the current hurry to obtain direct relevance for practical respects.

2. Methods

Multiconfigurational calculations were performed with General Atomic and Molecular Electronic Structure System (GAMESS) code [57], using the SARC-ZORA basis for lanthanides and cc-pVTZ for ligand atoms [58]. The DFT calculations were realized with Amsterdam Density Functional (ADF) code [59,60] using the built-in ZORA/TZ2P basis sets. The developed analyses and graphical representations were conducted with the original codes written in the Matlab-Octave environment [61,62]. Analytic formula derivation was realized with the help of Mathematica^TM symbolic algebra code [63,64].

3. Results and Discussion

3.1. The Multiconfigurational Account of Ligand Field Problems

We intentionally considered very simple model calculations in order to identify and clearly express the factors and problems in this domain. LF parameterization refers to one-electron quantities, partly due to two-electron interactions being borrowed from atomic theory (Slater-Condon [65,66] or Racah-type [67] parameters, adjusted by certain factors, with respect to the values resulting from the free-atom spectroscopy data). An LF model describes how the f-orbitals split in a given environment (as a consequence of effective one-electron factors), with the complete application to spectral states being the technical problem of applying the given LF scheme to the poly-electronic context. To effectively deal with the LF part only, ignoring the electron–electron interactions for acceptable reasons, we may conveniently choose systems showing, as the free-ion ground state, an F-type term, as is the case of Tb(III). In such a situation, the split of the poly-electronic ⁷F term is parallel with the conceivable pattern of f-type orbitals in the corresponding LF scheme.

Let us consider some Tb(III)-based model systems in CASSCF calculations, including eight electrons in seven orbitals, CASSCF(8,7). [TbF]⁺² does not correspond to any experimentally known species but offers a first glance into LF parameterization. The Tb–F distance is conventionally set to 2.35 Å, corresponding to the range encountered in various solid-state fluorides [68,69,70]. The CASSCF-computed states assignable to the ⁷Φ, ⁷Δ, ⁷Π, and ⁷Σ spectral terms, ordered in this sequence, have the relative energies 0.0, 497.5, 1008.7, and 1339.8 cm⁻¹, respectively. Figure 1 depicts the computed levels and suggests the relationship with the orbital LF scheme, depicting the canonical MOs, which is practically almost pure AOs, which can be described as the location of the β electron for the given spectral term. In this simple case, the terms are parallel to the itinerant f^β component from the f⁸ ≡ f^7αf^β configuration of Tb(III). The scheme also suggests the meaning of the AOM-type parameterization, namely, the relative magnitudes of the e^f_σ, e^f_π, and e^f_δ perturbations for f electrons, after conventionally setting the e^f_φ = 0 energy origin. The above energy levels were obtained in the primary active space of CASSCF(8,7), related to the f⁸ configuration. Extending the space to also comprise the five d-type virtual orbitals to CASSCF(8,12), the split of ⁷F term is slightly modified, with the results being interpreted as the following AOM parameters: e^f_σ(F) = 1175.0 cm⁻¹, e^f_π(F) = 910.2 cm⁻¹, and e^f_δ(F) = 466.5 cm⁻¹, in addition to the imposed e^f_φ = 0.

Note that the e_δ parameter is not negligible, as customary in the established AOM practice, presuming e^f_δ and e^f_φ are null, because the ligands are not expected to be able to exert δ and φ true orbitals overlap. Let us ignore this apparent discrepancy for the moment, debating the available details of the proposed computational experiments.

The Tb(III) ion is also well suited to discuss the d-type ligand field in lanthanides, analyzing the f⁸ → f⁷d orbital promotions (implying 5d empty AOs) in the CASSCF(8,12) setting. The f⁷ subconfiguration is spherically symmetric, the set of f⁷d excited states then has D representation, and the spectrum is parallel with the LF scheme of d electrons. In more detail, there is a ⁷D term arising from the f^7αd^β configuration and a ⁹D term arising from the spin swap to f^7αd^α. The high-spin states are lower in energy, benefiting from stabilization by the f–d inter-hell exchange coupling. In the following, as a measure of the d-type LF split, we consider the spin-septet terms, keeping the same multiplicity as in the ⁷F ground state. To be distinguished from the d-type LF encountered in transition metal complexes, here, we deal with an excited-state d orbital population. In linear symmetry, the ⁷D term splits into ⁷Δ, ⁷Π, and ⁷Σ levels. Shifting the ⁷Δ orbital doublet to zero, namely, imposing e^d_δ = 0, the relative gaps can be assigned to the d-type AOM parameters, resulting in e^d_σ(F) = 14,959.7 cm⁻¹ and e^d_π(F) = 4518.9 cm⁻¹. Figure 2 shows the discussed spectral terms and the assigned d-type orbitals from the f⁷d orbital promotions. In debating the effective one-electron parameters, the overall shift in the ⁷Δ is ignored because it is due to non-LF effects. As expected, the d-type LF split is sensibly larger than the f-type LF split because, having a larger radial extension, the 5d orbitals are closer to the ligand sites, undergoing a larger perturbation.

Considering that an LF scheme has the meaning of orbital energies in a given environment, there is a naïve expectation that computed canonical orbital levels may directly provide this scheme. Strictly speaking, this is not generally true, and is particularly not expected in the case of CASSCF techniques, where the orbitals are merely intermediate objects used to set the many-electron problem, their energies being subject to certain conventions (MO canonicalization [27]). However, we can verify that, in our case, the relative energies of f-type optimized orbitals, corresponding to φ, δ, π, and σ labels, {0.0, 482.8, 987.6, 1448.5} cm⁻¹, parallel the computed gaps between the ⁷Φ, ⁷Δ, ⁷Π, and ⁷Σ spectral terms. This regularity also occurs with respect to the d-type sequence, with the corresponding molecular orbitals having relative energy gaps similar to those of the D spectral terms: {0.0, 5464.9, 14,704.8} cm⁻¹. Recall that in the axial symmetry, the σ or Σ labels correspond to nondegenerate states, while the other ones represent doubly degenerate couples.

After conducting a thorough investigation on the grounds of simplified systems, we now consider the artificial [TbNe]⁺³ molecule to produce a measure of LF perturbation in the absence of charge. The neon atom is roughly similar as an electronic structure to the structure of a fluoride anion. The AOM parameters extracted from computed ⁷F spectral terms are e^f_σ(Ne) = 264.5 cm⁻¹, e^f_π(Ne) = 148.1 cm⁻¹, and e^f_δ(Ne) = 66.5 cm⁻¹, with imposed e^f_φ = 0.0; those of the ⁷D origin are e^d_σ(Ne) = 2342.9 cm⁻¹, e^d_π(Ne) = 661.3 cm⁻¹ while setting e^d_δ(Ne) = 0.0. The values are sensibly smaller than those of the fluoride ligand, indicating, in a semiquantitative sense, the important contribution of the electrostatic part.

We next considered linear [TbF₂]⁺ in the next computer experiment, observing similar spectral terms, because the molecule now spans the D_∞h symmetry, having the previous C_∞v of the asymmetric linear molecule as a subgroup. Then, for the same series of spectral terms, ⁷Φ, ⁷Δ, ⁷Π, and ⁷Σ, from the atomic ⁷F state, the calculation produces the 0.0, 1044.5 cm⁻¹, 2109.1 cm⁻¹, and 2926.9 cm⁻¹ gaps, namely, about double values, in comparison with those of [TbF]⁺². In this way, for f electrons, we reasonably retrieved the additivity hypothesis intrinsic to LF models. For the d shell, the split in [TbF₂]⁺ is 0.0, 7764.5 cm⁻¹ and 15,033.9 cm⁻¹ for ⁷Δ, ⁷Π, and ⁷Σ. Then, quite intriguingly, when the π-type perturbation in [TbF₂]⁺ is roughly twice that in [TbF]⁺², the d-type σ perturbations are similar in [TbF]⁺² and [TbF₂]⁺. We suppose that because of the stronger interaction between the 5d shell and the ligand, the simple premises of LF are not strictly followed, particularly not in the σ overlap, achieving an advanced intercept of the ligand near the maximum 5d radial zone. As in the previous case, the relative f-type orbital energies are quite similar to the spectral terms: {0.0, 1031.5, 2041.1, and 3006.8} cm⁻¹ for the {φ, δ, π, and σ} sequence of f-type orbitals, respectively. In turn, the orbital energies are not a good measure for the d-based LF sequence, with the relative values for {δ, π, and σ} being {0.0, 9986.1, and 20,455.0} cm⁻¹, respectively.

Given the good parallelism between the spectral terms corresponding to the f-based LF scheme and f-type orbital energies, we may suggest that these can be taken as approximations in situations where the relationship with the spectral term energies is not simply defined. Thus, for a Gd(III) free ion, having a nondegenerate ⁸S ground state, with the higher spectral terms being only spin sextets, the {φ, δ, π, and σ} relative orbital energies in [GdF]⁺² are {0.0, 482.8, 943.7, and 1316.8} cm⁻¹, respectively, quite similar to those of the [TbF]⁺² case. The same regularities hold for [GdF₂]⁺, which has {0.0, 1009.6, 1997.2, and 2875.1} cm⁻¹ f-type orbital gaps, approximately double the values of the single-ligand model molecule.

The following numeric experiment demonstrates the variation in the computed spectral terms in the [TbF₂]⁺ unit as function of the F–Tb–F angle, which we performed to check whether the geometry dependence follows the pattern of the angular overlap model version of ligand field theory. Figure 3a shows the computation data with open circles and the fitted AOM variation with a continuous line. The adjusted parameters, accounting for the whole variation, are e^f_σ(F) = 1341.7 cm⁻¹, e^f_π(F) = 1018.6 cm⁻¹, and e^f_δ(F) = 514.3 cm⁻¹. The CASSCF and AOM data were fitted in their barycenters for better visualization. The fit is quite good in the 90°–180° domain of the F–Tb–F angles, with the deviation at lower angles being pardonable because it enforces the strong ligand–ligand repulsion, and the small inter-ligand angles are avoided in real lanthanide chemistry. The conclusion that the AOM pattern is satisfactorily obtained with the calculation also holds for the d-type LF, represented in Figure 3b. The parameters giving an overall satisfactory fit to the d-LF angular variation are e^d_σ(F) = 8356.6 cm⁻¹ and e^d_π(F) = 3690.5 cm⁻¹, with a constrained e^d_δ(F) = 0.0.

Imposing F–Tb–F = 90° allows us to debate a subtle issue of ligand field models. One may observe that, at this point, two pairs of emulated AOM states are crossing, as if degeneracies are occurring. The TbF₂⁺ unit has C_2v symmetry in all geometries when it is nonlinear, not admitting degenerate states. Indeed, the computed points are not superposed in a state-crossing situation, although the corresponding pairs of values are mutually close. In this situation, we observe a hidden feature of LF models, called holohedrization. The term was coined by Schäffer [71,72], being a definite feature of LF models, but it is largely ignored in common practice. It arises because the expansion of LF potential in Equation (1) is based only on even spherical harmonics, k = 2, 4, and 6, in the case of an f shell. Thus, it cannot account for asymmetric LF maps. In other words, this creates, for every perturbation coming from a given ligand, a ghost image at its antipode. Then, the perturbation of a ligand L is actually divided into two halves, one generated from its real position, say at {x_L, y_L, z_L} Cartesian coordinates; and another from the {−x_L, −y_L, −z_L} inverted point. This artificial effect of LF models is tacitly inherited in AOM too. Then, in AOM as well, the perturbation due to the Tb–F contact is smeared as if produced by two-halves of fluorine ions in diametrically opposed placement, (F/2)–Tb–(F/2). By the same reason, a [TbF₂]⁺ unit with F–Tb–F = 90° behaves as a square of four ligand halves. Therefore, the LF potential and AOM matrix artificially acquire the D_4h symmetry of the square, where the representation of the f orbitals is a_2u + b_1u + b_2u + 2e_u. The two doubly degenerate e_u elements are the crossing points described in the above discussion. The d-type LF scheme, shown in Figure 3b, shows a single crossing due to holohedrization, because the representation of d orbitals in D_4h, a_1g + b_1g + b_2g + e_g, contains one doubly degenerate component.

An illustration of the holohedrization is represented using the color map of Equation (1). If the θ and φ variables are drawn as the parallels and meridians of a sphere, the variation in the V_LF(θ, φ) potential is qualitatively suggested by the color shading, taking deep blue for the lowest values, purple-red for maxima, and the orange-yellow-green palette for the in-between range. The sphere is taken on an arbitrary scale with respect to the molecular skeleton.

One can observe in Figure 4 that the action of the ligand is mirrored at its inversion center. For instance, in panel 4a, the ligand is placed at the north pole, but the caps of red coloring, marking the positive perturbation from the negative ligand against the negative electrons, show an equal size at the south pole. With the same factor and the holohedrization effect, the LF potential of the angular [TbF₂]⁺ unit at F–Tb–F = 90° has four-fold symmetry. The actual C_2v symmetry of the system has an artificial inversion center, generally becoming D_2h and particularly D_4h at the rectangular geometry. The artificial-square-like symmetry is clearly seen in the inset in the right-bottom corner.

In this intriguing situation, we proceeded to the examination of the ab initio computed electrostatic potential. Figure 5 shows the same systems with electrostatic potential drawn on a sphere with radius 0.5 Bohr, corresponding to the maximum of the radial profile of the f shell in the given atomic basis. To improve visibility, the spheres represented here are scaled in arbitrary mode. In more detail, we rendered the averaged electrostatic potential collected from the seven states of the state-averaged CASSCF calculations of the discussed systems.

In the upper half of Figure 5, somewhat surprisingly, a map suggesting the partial subsistence of the holohedrization effect appears in the case of the [TbF]⁺² diatomic molecule. A red area marking the axial penetration of the LF perturbation at the antipode of the Tb-F bond is clearly visible but is smaller than those appearing under the ligand. The map of [TbF₂]⁺ in the lower half of Figure 5 more substantially differs from the emulated LF map in Figure 4, showing only small spots at the antipodes of the Tb-F bonds and an absolute maximum between ligands along the line bisecting the F–Tb–F angle.

We may conclude that the holohedrization, although due to the limitations of the LF models is in part a true phenomenon for the f-shell compounds. This is an important result of this work, with the issue deserving further debate. As previously noted, holohedrization is largely overlooked in the current LF practices. Graphical LF mapping, as we show here, has not been exploited despite the relevance of this issue.

We could not check the above-discussed LF quantities against experimental data, because the selected case studies are artificial, given our aim of extreme simplicity to obtain academic insight. However, one may attempt a step toward more realistic objects considering the octahedral [TbF₆]⁻³ unit. As a stand-alone molecule, this molecule is not fully realistic, but similar units can be found as doped sites in elpasolites [73,74,75].

The AOM equates the octahedral LF eigenvalues as follows:

$$\epsilon \left(a_{2u}\right)=6{\mathrm{e}}_{\mathsf{\delta }}$$

(4)

$$\epsilon \left(t_{1u}\right)=\frac{5}{2}{\mathrm{e}}_{\mathsf{\phi }}+\frac{3}{2}{\mathrm{e}}_{\mathsf{\pi }}+2{\mathrm{e}}_{\mathsf{\sigma }}$$

(5)

$$\epsilon \left(t_{2u}\right)=\frac{3}{2}{\mathrm{e}}_{\mathsf{\phi }}+2{\mathrm{e}}_{\mathsf{\delta }}+\frac{5}{2}{\mathrm{e}}_{\mathsf{\pi }}$$

(6)

Inserting previously ascribed the e^f_λ(F) values from the CASSCF(8,12) calculation on [TbF]⁺², we predicted values of 916.3 cm⁻¹ and 409.5 cm⁻¹ for the ${\mathsf{\Delta }}_1=\epsilon \left(t_{1u}\right)-\epsilon \left(a_{2u}\right)$ and ${\mathsf{\Delta }}_2=\epsilon \left(t_{2u}\right)-\epsilon \left(a_{2u}\right)$ gaps, respectively. The calculation of the [TbF₆]⁻³ unit, in settings similar to those of the diatomic metal–ligand system, gives Δ₁ = 762.9 cm⁻¹ and Δ₂ = 291.9 cm⁻¹, which are roughly comparable to the result assuming the transferability of the AOM parameters. The discrepancy is explained by the different balances in the non-LF effects subsisting in the positively charged artificial diatomic system versus the negatively charged larger system. The AOM practice does not consider the octahedral gaps with the full set of {e_σ, e_π, e_δ, e_φ} parameters, being confined only to the {e_σ, e_π} couple, first by the empirical assumption that the effects of δ and φ are not chemically intuitive and then by the pragmatic reason that the available data do not allow obtaining a full parametric scheme. Thus, if we consider that the two gaps for the [TbF₆]⁻³ unit are fitted only in the σ-π scheme, we obtain e_σ = 293.9 cm⁻¹ and e_π = 116.7 cm⁻¹, which is in a range compatible with the experimental spectroscopy [42,76]. Although the chemical intuition may refute the φ and δ bonding effects, the electrostatics support their underground action. However, we leave this matter as a challenge open to further debate. We honestly admit that the calculations have certain limitations; from a different perspective, they can offer more details than are available from experiments.

Let us now take a polyatomic ligand, maintaining the simplicity; a relevant case study is the nitrate coordinated in chelate symmetric mode. The Tb-O distances are conventionally imposed at 2.35 Å. The CASSCF calculation results in the full split of the ⁷F term into seven distinct levels (though with very small spacing between first and second states) with the following relative energies: 0.0, 1.7, 451.1, 488.8, 499.6, 765.0, and 1184.5 cm⁻¹. Via the handling of the Hamiltonian matrix of the CASSCF calculations, these values are fitted by the $B_q^k$ parameters from Equation (1), as follows: $B_0^2$ = 1410.8 cm⁻¹, $B_0^4$ = −203.2 cm⁻¹, $B_0^6$ = −465.9 cm⁻¹, $B_2^2$ = 179.4 cm⁻¹, $B_2^4$ = 376.6 cm⁻¹, $B_2^6$ = 885.5 cm⁻¹, $B_4^4$ = 44.3 cm⁻¹, $B_4^6$ = −19.5 cm⁻¹, and $B_6^6$ = −1.9 cm⁻¹. We observe that in the given symmetry, only parameters with non-negative and even q indices do not vanish.

The determination of the parameters enabled the drawing of the V_LF(θ, φ) ligand field potential map, shown in the upper part in Figure 6. The LF maps undergo holohedrization, with the two maxima due to coordinated oxygen atoms being reflected exactly at the opposite pole. Interestingly, holohedrization is also visible in the ab initio electrostatic potential, averaged over the seven CASSCF states, as illustrated in the bottom half in Figure 6. The potential contours under the ligand and those recorded at antipodes are not perfectly similar but are comparable, again emphasizing that the identification of holohedrization in this numeric experiment is a nontrivial finding.

The previously outlined $B_q^k$ parameters are not transparent in their physical meaning. To attempt a more intuitive interpretation, we considered the AOM scheme. Here, the classical AOM knowledge indicates a limited set of parameters. One e^f_σ(O) is assignable to the lone-pair oriented toward the metal, with the π effects dichotomized in one due to the nominal π orbitals, with lobes perpendicular to the ligand plane e^f_π⊥(O), and another acting in plane e^f_π||(O), which may be negligible. The two oxygen atoms are equivalent, receiving equivalent parametric definitions.

Judging by the degrees of freedom available in the given C_2v symmetry, the parametric leverage proposed by intuition-driven AOM is too constrained, having three variables, while the approach in Equation (1) produces nine adjustable terms. To verify the available degrees of freedom, we counted the repetitions of the A₁-irreducible component in the representation of k = 2, 4, and 6 sets of spherical harmonics in C_2v. Thus, we picked two from k = 2 (2A₁ + A₂ + B₁ + B₂), three from k = 4 (3A₁ + 2A₂ + 2B₁ + 2B₂), and four from k = 6 (4A₁ + 3A₂ + 3B₁ + 3B₂), acquiring the nine-dimensional parametric variance.

The attempt to fit the simplest AOM scheme gives parameters that seem acceptable: e^f_σ(O) = 835.2 cm⁻¹, e^f_π⊥(O) = 438.4 cm⁻¹, and e^f_π||(O) = 264.0 cm⁻¹ (the in-plane π was not vanishing, as hoped) and a sequence of LF energies {0, 14.5, 352.4, 369.3, 507.4, 650.38, 1181.2} cm⁻¹, which very roughly reproduce the computed values. Allowing more parameters, such as e^f_δ⊥(O), e^f_δ||(O), e^f_φ⊥(O), and e^f_φ||(O), improves the match of the LF spectrum (although it does not make it perfect) at the cost of nonvanishing values for these terms, which do not belong to the classical canon of the AOM.

3.2. Density Functional Account of Ligand Field Problems

In the following, as case studies, we use the same objects as in the previous section, devoted to multiconfigurational calculations. Because the fⁿ configurations, except the f⁷, have degenerate ground states, and their compounds (within the expected weak LF) are quasidegenerate systems, DFT should not be applied without special caution and wise techniques. Considered in an unrestricted self-consistent version, i.e., allowing different molecular orbitals for α and β spins, one may obtain convergent solutions, which should be considered with precaution [46]. Unrestricted calculations are not well-suited because the energies of the f-type orbitals may drastically differ in α versus β series, which is a situation incompatible with the LF paradigm. On the other hand, restricted-type DFT calculations are prone to severe convergence problems, a situation that can be described as due to the non-Aufbau configuration of lanthanide ions in compounds [37]. As illustrated, we briefly mention the unrestricted DFT attempt (with BLYP functional) of the [TbF]⁺² system in the same GAMESS code used for the previous CASSCF calculations. The simplest unrestricted calculation with no educated guess does not reach convergence. Inserting a guess provided by the CASSCF calculation (with the same MOs initially doubled for the α and β sets), we obtained a convergent solution containing seven f-type α functions in the range −0.902 to −0.857 Hartrees, i.e., with a gap of 9690 cm⁻¹, which is too large for the LF scale. There are two β f-type orbitals instead of the expected one f^β component, whose energies, at −0.770 and −0.777 Hartree, are very distant from the range of f^α orbitals by about 29,000 cm⁻¹. The restricted calculation attempt was a total failure. It is clear then that the brute-force DFT approach is not satisfactory.

In order to assess the density functional theory (DFT) approach to LF, we employed the Amsterdam density functional (ADF) code [59,60], which allows a rather wide range of computational experiments based on the control of orbital populations, a feature usually not available in many codes.

The advanced ADF controls allow restricted-type calculations, enforcing fractional occupation with 8/7~1.1428 on each of the seven f-type orbitals expected to emulate the LF sequence. It is important to note that fractional occupations are conceptually allowed in DFT [77,78] and are technically enabled in the ADF environment. Averaging the population of the fⁿ configuration to n/7 for each f-type orbital, we can overcome the ban on using DFT for quasidegenerate ground states. Taking the concrete case of Tb(III) with an f⁸ configuration, where we must have a doubly occupied f-type orbital, it is expected that, in principle, small energy differences occur if different configurations are prepared with the doubly occupied MO moving along the seven f-type orbitals; this is outside the validity of DFT according to its basic theorems [79]. If, in turn, the content of the double-occupied orbital is smeared on the components of the LF-type set, this dilemma disappears because we consider an object having a singly defined, though artificially devised, state. From another perspective, the orbital averaged scheme is somewhat similar to the state-averaged calculations performed in multiconfigurational mode.

The population control in ADF is also helped by symmetry, making it possible to organize the orbitals by their point-group representation. Thus, in the C_∞v symmetry of the [TbF]⁺² molecule, there is a core due to the [Xe] and [Ne] configurations of the inner Tb(III) shells and fluoride anion, with occupation [Xe-Ne] = σ²⁸π²⁸δ⁸. The LF part has the fractional populations φ^16/7δ^16/7π^16/7σ^8/7.

We took the orbital energies of the restricted calculation with this definition, and we interpreted the gaps as akin to f-type AOM parameters, with e^f_φ = 0 fixed at zero. The results are shown in the first row in

Table 1. Numerical experiments on the [TbF]⁺² molecular model, assessing the extraction of AOM-type ligand field parameters from various settings in the frame of DFT, performed with specific input controls from the ADF code. All values are in reported in cm⁻¹.

[TbF]+2		efδ(F)	efπ(F)	efσ(F)
Restricted	MO energies	2984.3	5750.9	6944.6
	Frozen LF configurations	−1202.6	6991.4	4948.3
	Iterative LF configurations	−851.7	2710.1	5372.6
	Frozen density MOs	500.1	879.2	1032.4
Unrestricted	α MO energies	2201.9	4549.1	4315.2
	β MO energies	2387.5	3218.2	5258.9
	Frozen LF configurations	−1857.5	−738.0	5171.7
	Iterative LF configurations	−1424.4	−292.0	1554.3
	Frozen density α MOs	459.7	806.6	959.8
	Frozen density β MOs	540.4	951.8	1121.1

. In comparison with the results from CASSCF, it seems that this computational setting overestimates the LF parameters. We can explain this by the slightly faulted description of the f characteristic in the MOs determining the LF sequence. In the multiconfigurational calculations from the previous section, the canonical orbitals were obtained with almost 100% content in f AOs; in the above-described DFT calculation, this happens only for the δ and φ orbitals. The π functions were obtained with only 74% f and the σ with about 89% f, with the remainder being due to the ligand. This creates a nonphysical asymmetry in the shape of the f-type components, leading to the overestimation of the LF parameters. On the other hand, although improperly balanced, we comment that this situation is in the spirit of AOM, i.e., driven by weak covalency occurring on σ and π channels but not on the δ and φ components.

Here, we comment on the expected dependence on the choice of the density functional. For instance, conducting the calculation with the B3LYP hybrid form, we obtain e^f_σ(F) = 6121.9 cm⁻¹, e^f_π(F) = 1645.4 cm⁻¹, and e^f_δ(F) = 3008.5 cm⁻¹, noticing a nonphysical swap of π vs. δ orbital energies.

The leverages of the ADF code allow the emulation of excited orbital configurations. In the following numeric experiments, the doubly occupied MO from the f⁸ configuration was successively placed in different LF-type orbitals, i.e., among functions with a preponderant f characteristic. The orbitals were kept frozen as a result of the above-discussed averaged-population result. The second line in Table 1 lists the total energy of the φ³δ²π²σ¹ configuration as zero of the energy scale, with the relative energies of the φ²δ³π²σ¹, φ²δ²π³σ¹ and φ²δ²π²σ² states assimilated to the e^f_δ, e^f_π and e^f_σ perturbations, respectively. As such, we break the rule that DFT is not enabled for excited states. However, in this particular case, we may invoke a symmetry-sustained dispense. The above-outlined configurations span different representations of their Slater determinants, ⁷Φ, ⁷Δ, ⁷Π, and ⁷Σ. Then, we state that each configuration acts as a ground state on the different symmetry channels, providing validity to the population-controlled calculation in the DFT frame. The numeric values do not well reflect the expected LF pattern. The negative e^f_δ(F) means that φ²δ²π³σ¹ is a ground configuration, instead of the initial expectation of a φ²δ³π²σ¹ configuration. Although this numeric experiment is, in principle, more advanced than the assimilation of orbital energies with the LF diagram, it did not lead to an actual improvement. The third line in Table 1 illustrates a tentative refinement by taking iterative calculations on the tuned configurations, with the results remaining unsatisfactory.

The fourth line in Table 1 addresses a specialized handling implemented in ADF, called frozen density embedding (FDE) [80,81]. Although not specifically aimed at ligand field problems, the FDE approach serves this idea well because, in our case, the calculation is performed only on the electrons of the metal ion, with the ligand counting only as surrounding electron density. The approach is not confined only to the classical electrostatic field because it also accounts for the exchange and kinetic energy effects from the frozen electronic cloud of the ligands, taken as fragments described in the preamble from free anions. The FDE-DFT calculation is then performed only on the Tb(III) ion, selecting the energies of the perturbed f orbitals. The frozen ligand density result seems to be in the correct range, being close to the CASSCF data reported in the previous section.

In the following, we discuss the results obtained in unrestricted mode. Now, population averaging occurs only for one electron in β-type orbitals, while the α f-type orbitals are filled with seven electrons, i.e., [φ^2αδ^2απ^2ασ^α] φ^2β/7δ^2β/7π^2β/7σ^β/7, corresponding to the f^7αf^β/7 formulation of the f⁸ configuration of Tb(III). In this case, we selected two series of orbital energies for the two spin polarizations. We remark on the persistence of parametric overestimations and on a strange swap of e^f_σ vs. e^f_π ordering for the α-type orbitals. The numeric experiments with the imposed orbital occupation schemes considered the following configurations: [φ^2αδ^2απ^2ασ^α]φ^β, [φ^2αδ^2απ^2ασ^α]δ^β, [φ^2αδ^2απ^2ασ^α]π^β, and [φ^2αδ^2απ^2ασ^α]σ^β, with the energy of the first taken as a relative reference. The results of the configurations for the frozen MOs were catastrophic; it is hard to explain why. However, the results improved if orbital iterations were allowed for each configuration, with a quite good e_σ but poor e_π and e_δ values. The unrestricted FDE approach must consider two sets of orbital energies. As in the restricted version, this technique, which is an approximation, yields a result compatible with that of the LF paradigm. Particularly, the FDE β-type e^f_λ parameters show a resemblance to the CASSCF results.

Table 2. Numerical experiments on [TbF₂]⁺ molecular model, assessing the extraction of AOM-type ligand field parameters from various settings in the frame of DFT and performed with specific input controls from the ADF code. All values are reported in cm⁻¹.

[TbF2]+		efδ(F)	efπ(F)	efσ(F)
Restricted	MO energies	1770.4	3722.3	4920.1
	Frozen LF configurations	−382.3	2189.05	4441.8
	Iterative LF configurations	−128.25	2211.6	4218.4
	Frozen density MOs	556.5	988.1	1193.7
Unrestricted	α MO energies	1290.5	2355.2	2661.7
	β MO energies	1318.7	2621.4	3472.3
	Frozen LF configurations	−846.9	740.4	2123.3
	Iterative LF configurations	−572.3	831.6	1971.7
	Frozen density α MOs	512.2	907.4	1096.9
	Frozen density β MOs	604.9	1096.9	1330.8

shows the results of the same series of calculations applied to the linear [TbF₂]⁺ molecule. In this case, all computed energy gaps were divided by two with the idea of obtaining LF perturbations assignable to one fluoride ligand.

In the [TbF₂]⁺ system, although a general overestimation of LF parameters occurs, it has a smaller amplitude than that of the [TbF]⁺² system (except for the well-behaved FDE estimation). We tentatively explain this by the symmetric characteristic due to the F–Tb–F axis. Although this does not reduce the overestimated metal–ligand mixing, the orbital shapes with inversion symmetry better account for the effective LF phenomenology than the polarized bonds along the diatomic Tb–F line.

Let us briefly examine the DFT results on the [TbF₆]⁻³ unit that were obtained using the symmetry and occupation keywords. The energies of the restricted f-type MOs with eight electrons smeared on the seven orbitals spanning the a_2u + t_1u + t_2u composition gave gaps of Δ₁ = 1855.1 cm⁻¹ and Δ₂ = 709.8 cm⁻¹ (see notations for Equations (4)–(6)), which were about twice larger than those produced with the CASSCF estimation. The unrestricted calculation is puzzling, yielding spurious results for α-type MOs of Δ₁ = −2242.3 cm⁻¹ and Δ₂ = −685.6 cm⁻¹, while producing semiquantitatively acceptable results for β-type MOs of Δ₁ = 1371.2 cm⁻¹ and Δ₂ = 524.3 cm⁻¹. Connecting the points, it seems that DFT generally overestimates the LF perturbations, which is curiously attenuated in larger systems: severe in [TbF]⁺², large in [TbF₂]⁺, and partly acceptable in [TbF₆]⁻³. This can be understood by the fact that in larger systems, the central ion feels an averaged field closer to the quasispherical situation in which the LF paradigm is constructed. The quantitative lapses in the smaller systems deserve, however, to be kept under observation, signaling an unbalanced interplay of factors inside the black box of DFT computation. As a conclusion, we state the caveat about the general problems with DFT in describing the LF parameters, emphasizing however, the reliability of the FDE approach, at least for the situation of the f shell in ionic lanthanide complexes.

3.3. Revisiting the Crystal Field Electrostatic Approximation

Let us now consider the calculation results with standard LF Hamiltonian from Equation (1), taking the simplest case of the [TbF]⁺² model molecule, as described in Section 3.1. As outlined in the Introduction, Equation (1) should obey the symmetry of the problem, taking the combinations of the spherical harmonics that render total symmetric representation in the given point group. In the case of linear systems, such as [TbF]⁺², in the C_∞v point group, the total symmetric representations are simply the Y_k,0 components in each spherical harmonic set. In other words, only the $B_0^k$ parameters are retained, and all the other $B_q^k$ ones, with q ≠ 0, vanish. The ab initio results on [TbF]⁺² from Section 3.1 were fitted with $B_0^2$ = 2037.7 cm⁻¹, $B_0^4$ = 662.3 cm⁻¹, and $B_0^6$ = 232.4 cm⁻¹. According to convention, the eigenvalues are placed with their barycenter in the zero of the energy scale. If we add an overall shift, with $B_0^0$ = 386.9 cm⁻¹, the eigenvalues are the same as interpreted in the AOM parameterization, with a fixed e_φ = 0.0 cm⁻¹. We here relabel the parameters $B_0^k=b_k$. For the one-ligand case, in the C_∞v point group, the eigenvalues of the Hamiltonian from Equations (1) and (2) in terms of $b_k$ parameters, classified by (σ, π, δ, ϕ) symmetry labels, are:

$${\epsilon }_{\sigma }=b_0+\frac{4}{15}{b}_2+\frac{2}{11}{b}_4+\frac{100}{429}{b}_6$$

(7)

$${\epsilon }_{\pi }=b_0+\frac{1}{5}{b}_2+\frac{1}{33}{b}_4-\frac{25}{143}{b}_6$$

(8)

$${\epsilon }_{\delta }=b_0-\frac{7}{33}{b}_4+\frac{10}{143}{b}_6$$

(9)

$${\epsilon }_{\phi }=b_0-\frac{1}{3}{b}_2+\frac{1}{11}{b}_4-\frac{5}{429}{b}_6$$

(10)

Note that the sum of the eigenvalues is ${\epsilon }_{\sigma }+2{\epsilon }_{\pi }+2{\epsilon }_{\delta }+2{\epsilon }_{\phi }=7b_0$, with the standard LF scheme imposing b₀ = 0. In the AOM, we postulate ${\epsilon }_{\phi }=0$ or, equivalently, $b_0=({\epsilon }_{\sigma }+2{\epsilon }_{\pi }+2{\epsilon }_{\delta })/7$.

If the LF Hamiltonian is taken in the crudest electrostatic definition, determined by a q_L point charge at distance R_L from the lanthanide center, the b_k parameters are

$$b_k\left(L\right)=q_L{{\displaystyle \int }}_{r=0}^{\infty }{\left(R\left(r\right)\right)}^2\frac{\mathit{min}{(r,R_L)}^k}{\mathit{max}{(r,R_L)}^{k+1}}{r}^{2}dr$$

(11)

where r is the radial coordinate of the electron, and R(r) is the radial function of the atomic orbitals undergoing the LF perturbation. This expression comes from the multipolar expansion of the electrostatic potential [1,2]. Although it quickly became clear that this estimation was not working, tacitly referring to the d-type LF frame, we revisited this issue because, to the best of our knowledge, it has not been re-examined in the modern age of computation with a focus on the lanthanides case.

The evaluations with the above equation do not depend on a dedicated calculation method: they are determined merely by the basis set of the LF shell. Because we worked in the Gaussian-type orbital (GTO) based frame in Section 3.1 and we used a Slater-type orbital (STO) set, specific to the ADF infrastructure, in Section 3.2, we exploited the opportunity to check the idea of first-principles electrostatic approximation using these different definitions.

The radial part of the LF shell is taken as a linear combination:

$$R\left(r\right)={\displaystyle \sum }_{i=1}^{n_f}{c}_i\cdot N_i\cdot \chi \left(n_i, {\zeta }_i,r\right)$$

(12)

where n_i and ζ_i are the power and exponent parameters of the primitive χ functions defining the basis, respectively; c_i are the coefficients resulting from the self-consistent calculation; and N_i are the normalization factors.

In the GTO case, the primitives are

$${\chi }_{GTO}\left(n_i, {\zeta }_i,r\right)=r^{n_i-1}\cdot \exp \left(-{\zeta }_i\cdot r^2\right)$$

(13)

with the normalization factors

$$N_i^{\mathrm{GTO}}=\sqrt{2{\left(2{\zeta }_i\right)}^{\frac{2n_i+1}{2}}/\mathsf{\Gamma }\left(\frac{2n_i+1}{2},0\right)}$$

(14)

expressed in terms of the gamma incomplete function.

In the STO frame, the primitives and normalization factors are:

$${\chi }_{STO}\left(n_i, {\zeta }_i,r\right)=r^{n_i-1}\cdot \exp \left(-{\zeta }_i\cdot r\right)$$

(15)

$$N_i^{\mathrm{STO}}=\sqrt{{\left(2{\zeta }_i\right)}^{2n_i+1}/\left(2n_i\right)! }$$

(16)

In both cases, the general formulation of the pure electrostatic LF parameters is:

$$b_k={\displaystyle \sum }_{i=1}^{n_f}{\displaystyle \sum }_{j=1}^{n_f}{c}_i\cdot c_j\cdot N_i\cdot N_j\cdot \left(\frac{1}{R_L^{k+1}}\cdot P_{ij}+R_L^k\cdot Q_{ij}\right)$$

(17)

The P_ij and Q_ij terms have specific expansions in GTO vs. STO instances. With the help of Mathematica^TM [63,64] computer-algebra code, we were able to derive closed analytical formulas, being close in the GTO case:

$$P_{ij}^{\mathrm{GTO}}=\frac{1}{2}{\left({\zeta }_1+{\zeta }_2\right)}^{-\frac{1}{2}\left(n_1+n_2+k+1\right)}\left(\mathsf{\Gamma }\left[\frac{1}{2}\left(n_1+n_2+k+1\right),0\right]-\mathsf{\Gamma }\left[\frac{1}{2}\left(n_1+n_2+k+1\right),\left({\zeta }_1+{\zeta }_2\right)R_L^2\right]\right)$$

(18)

$$Q_{ij}^{\mathrm{GTO}}=\frac{1}{2}{\left({\zeta }_1+{\zeta }_2\right)}^{\frac{1}{2}\left(k-n_1-n_2\right)}\mathsf{\Gamma }\left[\frac{1}{2}\left(n_1+n_2-k\right),\left({\zeta }_1+{\zeta }_2\right)R_L^2\right]$$

(19)

For STO, we obtained:

$$P_{ij}^{\mathrm{STO}}={\left({\zeta }_1+{\zeta }_2\right)}^{-n_1-n_2-k-1}\left(\mathsf{\Gamma }\left[\left(n_1+n_2+k+1\right),0\right]-\mathsf{\Gamma }\left[\left(n_1+n_2+k+1\right),\left({\zeta }_1+{\zeta }_2\right)R_L\right]\right)$$

(20)

$$Q_{ij}^{\mathrm{STO}}={\left({\zeta }_1+{\zeta }_2\right)}^{k-n_1-n_2}\mathsf{\Gamma }\left[\left(n_1+n_2-k\right),\left({\zeta }_1+{\zeta }_2\right)R_L\right]$$

(21)

Let us apply the GTO formulas for the free Tb(III) ion, taken on the SARC-ZORA basis [82] with a ζ = {32.1811, 10.727, 3.57568, 1.19189, 0.397298, 0.132433} sequence of exponent parameters in Bohr⁻² units, and all n_i = 4. The state-averaged CASSCF calculation on the ⁷F ground state defines the set of c= {0.405816, 0.614349, 0.227727, −0.0382306} coefficients. Inserting these values into the above formulas, we obtain electrostatic parameters b₀ = 49,464.5 cm⁻¹, b₂ = 1964.1 cm⁻¹, b₄ = 189.2 cm⁻¹, and b₆= 49.7 cm⁻¹, which determine the relative eigenvalues e_σ = 1207.8 cm⁻¹, e_π = 1028 cm⁻¹, e_δ = 601.4 cm⁻¹, and e_φ = 0 (equivalent to the AOM definitions), respectively. The b₀ parameter, usually ignored in the LF expansion, is very close to the classic point-charge potential, 1/R_L = 0.225182 a.u. ≡ 49,421.7 cm⁻¹. The b_k parameters are in the same range as those fitted from the ab initio spectral terms, but they do not match perfectly. On the other hand, the relative eigenvalues are very similar in the two schemes.

In the following, we check the pure electrostatic approximation of the LF parameters for the STO-type ZORA/TZ2P basis used in the calculations in Section 3.2. This basis is defined by ζ = {10.2, 4.9, 2.15, 1.65} exponent parameters in Bohr⁻¹ and the n = {4,4,4,5} sequence of power factors. The BLYP calculation of the free Tb(III) ion provides c = {0.405816, 0.614349, 0.227727, −0.0382306} coefficients. The STO-based estimation yields b₀ = 49,421.5 cm⁻¹, b₂ = 2257.0 cm⁻¹, b₄ = 271.2 cm⁻¹, and b₆ = 67.3 cm⁻¹, which are just a bit larger than the results from GTO. With this, we derive the following parameters: e_σ = 1395.3 cm⁻¹, e_π = 1176.3 cm⁻¹, and e_δ = 675.64 cm⁻¹ if conventionally considering e_φ = 0. The results are fully comparable to those for GTO.

We obtained, via the crudest crystal field approximation, results that are compatible with those of high-level multiconfigurational calculations. This contradicts the long-standing belief about the failure of pure electrostatic estimation. Actually, the description of pure electrostatic paradigm was pronounced about d-type LF systems, where partial covalent bonds truly occur in coordination compounds. Sed contra, lanthanide compounds are more ionic and are therefore better suited for the original ideas of electrostatic ligand fields. This situation deserves further debate; here, the limitation is that the results may be due to the very simplified nature of the discussed cases, with the imposed idealization having yet insightful value and even practical connotations, suggesting that if proven that we can rely on a simpler level of approximation in applications.

4. Conclusions

We performed an unprecedented examination of the computational methods used for the LF-dedicated description. The deliberate use of extremely idealized case studies produced details usually hidden in calculations on large-scale systems, which deserve to be considered in further explorations and exploitations of current practical problems.

Comparing the methods pertaining to wave function theory (WFT) versus those of density functional theory (DFT), we clearly concluded that the former is better suited for LF problems, while adapting the rather unsystematic DFT to this aim may be a useful quest.

Strikingly, revisiting the early hypothesis of the electrostatic nature of LF perturbation, we found it suited for the ionic complexes of lanthanides. We produced very similar LF parameters using crude estimations with the point-charge perturbation of the atomic f shell via multiconfigurational WFT methods and by the DFT, confined to the frozen density embedding (FDE) approximation. These corroborated results hint at the possible design of simplified approaches that are usable in complicated situations, specific to property-design challenges in the search for new materials. Here, LF can be a valuable ancillary tool in luminescence problems.

Here, we confined our study to hard-base ionic ligands. The conclusions are valid as extrapolation for large classes of solid-state systems with such characteristics (e.g., fluorides, oxides, phosphates, and silicates). We acknowledge that the factors may change if we consider systems in the opposed extreme, i.e., neutral and soft-base ligands, creating a challenge open to debate for the outlined issues.