Ab initio non-covalent crystal field theory for lanthanide complexes: a multiconfigurational non-orthogonal group function approach

Abstract

We present a multiconfigurational ab initio methodology based on non-orthogonal fragments for the calculation of crystal field energy levels and magnetic properties of lanthanide complexes, implementing a systematic description of non-covalent contributions to metal–ligand bonding. The approach consists of two steps. In the first step, appropriate ab initio wave functions for the various ionic fragments (lanthanide ions and coordinating ligands) are optimized separately, accounting for the influence of the surrounding environment within various approximations. In the second and final step, the scalar relativistic (DKH2) electrostatic Hamiltonian of the whole molecule is represented on the basis of the optimized metal–ligand multiconfigurational non-orthogonal group functions (MC-NOGFs), and reduced to an effective (2J + 1)-dimensional non-orthogonal configuration interaction (CI) problem via Löwdin-partitioning. Within the proposed formalism, the projected non-orthogonal CI Hamiltonian can be expanded to any desired order of perturbation theory in the fragment-localised excitations out of the degenerate space, and its eigenvalues and eigenfunctions provide systematic approximations to the crystal field energies and wave functions. We present here a preliminary implementation of the proposed MC-NOGF method developed for first-order degenerate perturbation theory within our own ab initio code CERES, and compare its performance both with the simpler non-covalent orthogonal ab initio approach, Fragment Ab Initio Model Potential (FAIMP) approximation, and the full CAHF/CASCI-SO method, accounting for metal–ligand covalency in a mean-field manner. We found that the energies and magnetic properties of 44 complexes obtained via an iteratively optimized version of our MC-NOGF first-order non-covalent method compare remarkably well with those obtained using the full CAHF/CASCI-SO method including metal–ligand covalency, thus exposing the predominantly electrostatic character of the metal–ligand interactions, and are superior to those obtained using the FAIMP approach, which in its iteratively optimised variant was believed to date to be the best ab initio description of non-covalent metal–ligand interactions.

---

Introduction

Trivalent lanthanide ions are extensively employed in many technological applications, ranging from energy production to life sciences.¹ While they are certainly fundamental components in many optical applications,^2–4 these ions play a very special role in magnetism, thanks to their large magnetic moments and magnetic anisotropy.^5–7 At the root of their enhanced magnetic properties are the intrinsically strong electronic correlation and spin–orbit coupling interactions and comparatively weak interactions between the 4f electron shell and its coordination environment, mainly as a result of the 4f shell's radially compact nature, and its consequent efficient shielding from the surrounding electrostatic and chemical environment by the outer 5s²5p⁶ orbital shells.⁵ The simplest and earliest approach to describe the splitting of a ground (2J + 1)-fold degenerate ground multiplet produced by the surrounding chemical environment is called crystal field theory (CFT).⁸ CFT describes the splitting of the 4f many-electron wave functions associated with the (2J + 1) ground spin–orbit multiplet of the pristine atomic 4f-shell as operated solely by the static electric fields induced by the charge distribution on the ligands. CFT entails a chemical bonding description which is purely ionic, and it is in fact referred to also as the ionic model.

While clearly CFT is an oversimplified description for quantitative purposes, which even at the phenomenological level is riddled by an over-parameterization problem (in low symmetry the CFT potential consists of 27 free parameters), its heuristic value should not be underestimated, even today within the field of molecular magnetism.

For instance, Rinehart and Long⁹, discussing the asphericities of the electron density of the magnetic states in trivalent lanthanide ions calculated by Sievers some twenty years earlier¹⁰, pointed out how a proper axial (equatorial) distribution of the charge density in the coordinating ligands could be harnessed to synthesise Ln complexes with large magnetic anisotropy, by stabilizing the appropriate oblate (prolate) charge density distribution in that specific ion. This simple but insightful observation, based on the assumption of purely ionic bonding interactions between the metal and the ligand, has proven rather useful in guiding synthetic chemists towards ever more efficient Ln-based single-molecule magnets (SMMs).

These simple electrostatics arguments were also successfully translated into a more quantitative classical electrostatics repulsion energy minimization procedure, which was shown to be capable of semi-quantitatively predicting the specific direction of the magnetic anisotropy axis in many Dy(III) single-molecule magnets,¹¹ a finding that was recently confirmed via X-ray diffraction experiments.¹²

Despite these useful results, it is known since the 1950s that quantitative descriptions of lanthanide–ligand bonding cannot be achieved without accounting for non-covalent interactions other than pure electrostatics (e.g. induction, contact or dispersion interactions), and for various flavors of so-called covalent contributions (charge transfer excitations of various kinds).^13–22

A way to systematically account for all possible metal–ligand electronic interactions consists of ab initio calculations, typically based on molecular orbital (MO) theory, which in recent years has proved to be a reliable tool for the description of low-lying multiplets of lanthanide complexes with good accuracy.²³

However, in these methods the interactions between the electrons in the environment and the ion are not treated in a perturbative manner, leading not only to an increased computational cost, but also to a difficult interpretation of the different physical mechanisms involved in the description of metal–ligand bonding.

A recent work presented a study on the role of the electrostatic interactions in the prediction of crystal field levels in lanthanide complexes.²¹ In that paper, the authors analyzed the ligand field generated by both different atomic charge distributions and a more rigorous approach based on the Fragment Ab Initio Model Potential (FAIMP) method.^24–27 Besides the electrostatic interaction, the FAIMP approach allows inclusion of (i) the exchange interaction between the electrons of the ligand and the metal, and (ii) Pauli's repulsion contribution originating from the requirement of orthogonality between the metal and ligand molecular orbitals, which prevents over-delocalization of the electrons of one fragment in the core regions of the other. The authors' conclusions are that the electrostatic approach, even at the ab initio level via the FAIMP approach, does not suffice for an accurate description of the crystal field levels in Ln(III) complexes, and hence that covalency effects must play an important role.

In this work we investigate fragment ab initio modelling of crystal field levels in Ln(III) systems from an innovative point of view. Because it is considered one of the most rigorous non-covalent embedding methods of a metal ion in an ab initio potential created by the ligand, we start by briefly reviewing the FAIMP model as an ab initio electrostatic description. Then we introduce a novel ab initio approach based on the theory of multi-configurational non-orthogonal group functions (MC-NOGF), with the aim of modelling the strongly ionic metal–ligand interactions in terms of a systematic ab initio description based on purely intermolecular non-covalent interactions. This approach allows us to treat the metal–ligand interactions as a perturbation acting on the wave functions of the isolated fragments, which accounts in principle for both purely non-covalent (electrostatics, induction and dispersion) interactions and overlap effects which arise from the non-orthogonal theoretical framework. While in other methods focused on intermolecular interactions, such as Symmetry Adapted Perturbation Theory (SAPT),^28–30 the expansions of the intermolecular interaction operator and of the overlap integrals are truncated to some order, the electronic Hamiltonian in MC-NOGF is the full molecular Hamiltonian, and the overlap effects are fully accounted for. Moreover, in the developments the preservation of the orbital localization on the fragments allows for a clear identification and separation of the different molecular interactions, which is not achievable in the conventional non-orthonormal configuration interaction (CI) methods, where usually the non-orthogonality is partially or totally removed before the evaluation of the Hamiltonian operator matrix elements by appropriate orbital mixing (e.g. biorthogonal basis sets).^31–34 After presenting the equations for an ab initio multiconfiguration development involving matrix elements between the non-orthogonal sets of orbitals, as we implemented in our own quantum chemistry code CERES³⁵, we study how different approximations to the ligands' wave function perform within the MC-NOGF method in the simulation of the crystal field level and magnetic properties of a set of Ln(III) complexes.

Hamiltonian and product functions

Within the Born–Oppenheimer approximation, the non-relativistic electronic Hamiltonian of N_TOT interacting electrons can be defined in atomic units aswithwhere ∇²(i) = ∂²/∂x_i² + ∂²/∂y_i² + ∂²/∂z_i², V(i) is the electron-nucleus potential, and r_ij = r_i − r_j. Within the group function (GF) approximation,^36–43 the total Hamiltonian [script letter H] can be represented on the basis of the antisymmetrized products of the wave functions for the single parts composing the system aswith ψ^R_r describing the group R, which collects N_R electrons in the state r, and [scr A, script letter A] is the partial antisymmetrizer operator, which exchanges the electron coordinates x between two, or more, different groups. The group functions ψ^R_r are individually antisymmetric with respect to the electron exchange

ψ^R_r(…,x_i,…,x_j,…) = −ψ^R_r(…,x_j,…,x_i,…) (4)

and ortho-normalizedwhere δ_ij is the Kronecker delta. The product function Ψ_κ is fully characterized once the set of states κ ≡ (r, s,…) is made explicit. No inter-molecular correlation is explicitly included in a wave function of the form in eqn (3), since electrons belonging to different groups are assumed to move independently, but a large part of the intra-molecular correlation may be accounted for by admitting some degree of configuration interaction within each electron group.

Strongly orthogonal group functions and the FAIMP approach

A convenient way to handle fragment approaches based on group functions requires the group functions ψ^R_r to fulfil a so-called strong-orthogonality condition³⁸:so that integrating over any one variable x₁ common to two different group functions makes the integral vanish identically for all values of the other variables x_i. This is a much stronger condition than the orthonormality requirement in eqn (5) which is assumed for the different group functions within each group.

It is well known in the literature that when a wave function has the form given in eqn (3) and satisfies the conditions in eqn (4), (5), and (6), the system is group-separable³⁸. This means that the total energy of the system E can be written as (after integration over the spin variables)

where

H^R = 〈 ψ^R_r|[script letter H]^R|ψ^R_r〉 (8)

where P^R_r is the spinless one-particle density matrix of the function ψ^R_r. [script letter H]^R in eqn (8) is the electronic Hamiltonian for the group R:where V(i) in h(i) is the potential energy of electron i in the field of the nuclei of the total system.

Let us assume now that ψ^R_r are functions of the sets of mono-electronic spin-orbital functions {ϕ^R_i}, where ϕ^R_i(x) = χ^R_i(r)η^R_i with spatial part χ^R_i(r) and spin η^R_i. Then eqn (6) is automatically satisfied if we impose the orthonormality conditions:

Starting from the total wave function in eqn (3), and defining the groups R, S,… as doubly occupied closed-shell singlet systems described by the functions

where [scr A, script letter A]^R is the antisymmetrizer operator which exchanges all the electron spin–spatial coordinates x within the group R, and det|…| indicates the Slater determinant. Huzinaga and coworkers found that an optimum {ϕ^R_i} set, which minimizes the energy E while holding fixed the {ϕ^S_i} sets with S ≠ R, and imposing the conditions in eqn (12) and (13), is given by^24–26[scr F, script letter F] is the Fock operator defined bywith h_eff(i) given bywhere [scr J, script letter J] and [scr K, script letter K] are the Coulomb and exchange operators, respectively, describing the effective field at point i due to the electrons in the groups S:Note that [scr F, script letter F] is formally identical to the Fock operator for the total wave function Ψ_κ. The projector [scr P, script letter P]^S originates from the orthonormality conditions, and reads^24–26Eqn (15) can be used as the self-consistent field (SCF) equation to determine the best orbitals for the group R interacting with the closed-shell singlets included as one or more groups S, enforcing the orthogonality between the orbitals of all groups. Eqn (15) is the starting point of all the developments that go under the name FAIMP.^27,44

Let us now introduce a group M which is better described by a multi-configurational wave function ψ^M_m defined as

where C_mm′ is the weight of the m′-th Slater determinant as given in eqn (14) built from the set of othonormal spin–orbitals {ϕ^M_i}. As usual when dealing with a multi-configurational function, we can distribute the N_M electrons into two sets of orbitals: (i) one set collecting N^inact_M orbitals that are doubly occupied in all the configurations , called inactive orbitals, and (ii) a second set collecting N^act_M orbitals which allow for an occupation number from 0 to 2, called active orbitals. Introducing ψ^M_m in the product function Ψ_κ defined by eqn (3) we obtainwhere now λ ≡ (m′, r, s,…) and κ ≡ (m, r, s,…). The problem of finding an optimal set {ϕ^R_i} for the closed-shell singlet group R, still described by a single Slater determinant wave function ψ^R_r, leads to the Fock operator as defined in eqn (16), but now h_eff(i) has the formwhere [scr J, script letter J]^S_j(i)/[scr J, script letter J]^M_j(i) and [scr K, script letter K]^S_j(i)/[scr K, script letter K]^M_j(i) are defined in eqn (18) and (19), andwhere ρ^m_jk is the one-electron transition density matrix for ψ^M_m, given bywhich is averaged on N_m states in eqn (23) to account for the energy degeneracy of ψ^M_m functions. In this work the projector operator [scr P, script letter P]^M is defined as

Non-orthogonal group functions

FAIMP is an ab initio rigorous way to introduce the electrostatic interactions and Pauli repulsion between the metal and the environment; however, perhaps not surprisingly, it has shown poor performance in the accurate modelling of the electronic structure of Ln(III) complexes.²¹ To improve on this model, while keeping on pursuing a fragment ab initio approach to describe the metal–ligand interactions essentially in terms of intermolecular forces, let us start again from eqn (3), fix the sets of spin–orbitals {ϕ^R_i} for all groups to some optimal orbitals obtained for the isolated fragments, and, crucially, relax the strong-orthogonality constraints. Thus we define the wave function for the whole system aswhere Ψ_κ is a product function as defined in eqn (3), including a few closed-shell singlet groups described by ψ^R_r wave functions given in eqn (14), i.e. the ligands, and one center characterized by the multi-configurational wave function ψ^M_m shown in eqn (21), i.e. the metal. We fix the spin-orbitals entering the latter functions, and the C_mm′ coefficients entering eqn (21), to the sets obtained by energy minimization on the isolated fragments. Under this requirement, the orthogonality relationships in eqn (12) and (13) readWe are then interested in defining the optimal coefficients D_κ entering eqn (28), which can be estimated by variational minimization of the total energyor, equivalently, in a matrix form solving the eigenvalue problem

HD = EMD (32)

where

M_κκ′ = 〈Ψ_κ|Ψ_κ′〉 (33)

H_κκ′ = 〈Ψ_κ|[script letter H]|Ψ_κ′〉 (34)

The dimension of the variational problem in eqn (32) is typically very large. Hence, to reduce the dimensionality, and also to set up an effective ab initio crystal field splitting theory i.e. a degenerate perturbation theory within the lanthanide ion's degenerate ground level, we apply the Löwdin partitioning procedure. This is achieved by partitioning the full space spanned by H, D and M in the two subspaces A (here the (2J + 1)-degenerate space of the lanthanide ground spin–orbit multiplet) and B (the complement excitation space)³⁸Thus an effective Hamiltonian H̃^AA can be written for the subspace A as

H̃^AAD^A = EM^AAD^A (36)

where

H̃^AA = H^AA − (H^AB − EM^AB)(H^BB − EM^BB)⁻¹(H^BA − EM^BA) (37)

Note that H̃^AA is itself the function of the unknown energy E. Using the matrix identity

(X − Y)⁻¹ = X⁻¹ + X⁻¹Y(X − Y)⁻¹ (38)

and using the fact that, upon iteration of eqn (38), the inverse of the (X − Y) term can be written as

(X − Y)⁻¹ = X⁻¹ + X⁻¹YX⁻¹+… (39)

the matrix elements of H̃^AA have the formwhere eqn (39) has been used together with

X = H^BB_diag − EM^BB_diag (41)

Y = H^BB_off-diag − EM^BB_off-diag (42)

to make the different expansion terms explicit. This partitioning scheme projects the problem to the NOGFs formed from the metal-based multielectron wave functions belonging to the (2J + 1)-degenerate ground spin–orbit multiplet, defining subspace A, whereas the interaction with the metal-based and ligand-based excitations belonging to subspace B can be accounted for in a perturbative manner, as per eqn (40).

Here we focus our attention on terms up to first-order in the perturbative expansion eqn (40) (i.e. the very first term on the right-hand side of eqn (40)). The energies are then determined by eqn (36) with H̃_κκ′ = H_κκ′, and the variational optimization involves the diagonalization of a matrix of dimension (2J + 1) × (2J + 1) once the integrals in eqn (43)–(45) are made explicit.

Implementation in the code CERES: cofactor versus inverse matrix formalism

We describe here the implementation strategy pursued within our own ab initio code CERES^35,45 of the non-orthogonal CI problem arising from eqn (36) and (37) represented on the basis of the multiconfigurational non-orthogonal group functions, which is at the heart of the proposed MC-NOGF fragment approach.

From now on we drop the indices (r, s,…) ≡ (0, 0,…) in λ and κ indices for simplification. The explicit form for the normalization factor M_κκ′ on the Ψ_κ basis set is given by

In H_κκ′, the one-electron terms have the formand two-electron terms

Defining the overlap matrix S as

it is well known in the literature that

〈Ψ_n|Ψ_n′〉 = det(S) (47)

andwhere det(S) indicates the determinant of S, cof[S](ij) indicates the first-order cofactor arising from S when the row i and the column j are removed, and cof[S](ij|kl) indicates the second-order cofactor produced by further deleting the row k and the column l. If S⁻¹ exists, the cofactors can be expressed using the Jacobi ratio theorem

cof[S](ij) = det(S)[S⁻¹]_ji (50)

cof[S](ij|kl) = det(S)([S⁻¹]_ji[S⁻¹]_lk − [S⁻¹]_li[S⁻¹]_jk) (51)

Although this theorem is a powerful tool that makes the calculation of cofactors simpler, the use of eqn (50) and (51) requires a matrix inversion for every couple of determinants in eqn (43)–(45), and can be a non-valid development when the matrix to be inverted is singular and the inverse does not exist.

Löwdin first presented a formula for the matrix elements of the Hamiltonian between Slater determinants expressed in terms of non-orthogonal spin-orbitals,^46–49 and after this different methods have been suggested to avoid the need for calculating the inverse, allowing for a rapid evaluation of the cofactors,^50–52 but in general they do not explicitly consider that the determinants from the set employed in a calculation may differ only in few spin-orbitals from a reference determinant. Hayes and Stone are the first to emphasize the advantages coming from properly taking this peculiarity into account,⁵³ followed by Figari and Magnasco, who generalized Hayes' results.⁵⁴ Starting from the latter developments, we present here a slightly modified approach, which allows for both computational efficiency and generality.

We start reordering the spin-orbital in S to obtain the block partitioning

where A is a symmetric matrix of dimension N_inact × N_inact, which includes the overlap elements between the spin-orbitals that are common to all the product functions, which are the inactive set of spin-orbitals for the metal and the spin-orbitals included in the ψ^R₀ group functions. The C block collects the overlap terms between the N_act orbitals that can vary in the different product functions, which in this paper coincide with the active set of spin–orbitals for the metal. Note that N_act ≪ N_inact. The U and V blocks include the overlap elements between the active and inactive sets of spin-orbitals. While U, V, and C depend on a specific product function pair, A never changes. Then, S⁻¹ can be written by the use of the Woodbury matrix identity as⁵⁵where Z = (C − VA⁻¹U), X = A⁻¹U, and Y = VA⁻¹. To avoid problems when Z is singular, we can apply eqn (50) back to the Z matrix's elements:

[Z⁻¹]_ij = cof[Z](ji)det(Z)⁻¹ (54)

or, in the matrix form

Z⁻¹ = cof[Z]^T det(Z)⁻¹ (55)

where cof[Z]^T indicates the matrix collecting all the first-order cofactors of Z transpose. Defining the two matricesandand remembering that

det(S) = det(A)det(Z) (58)

eqn (50) and (51) can be written by the use of eqn (53) and (55) as

cof[S](ij)det(A)⁻¹ = det(Z)T_ji + W_ji (59)

where K_ijkl is defined as (all the possible cases under the conditions i < k and j < l in eqn (49))with the p, q, r, and s indices running on the Z block. The relationship

cof[Z](ij|kl)det(Z) = cof[Z](ij)cof[Z](kl) − cof[Z](il)cof[Z](kj) (62)

has been used to cancel out the 1/det(Z) factor from K_ijkl terms, and the permutational symmetry relationships

cof[Z](ij|kl) = −cof[Z](ji|kl) = −cof[Z](ij|lk) = cof[Z](ji|lk) (63)

have been used for simplification. We noted that, defining the two matriceswhere I is the identity matrix with rank(I) = rank(C), Z_ijkl can be written by a unique expression asEqn (59) and (60) are not affected anymore by singularities, and they do not require any conditions but the existence of A⁻¹. Moreover, the need for the calculation of first- and second-order cofactors of Z is not a big computational issue, because of the generally small dimension of this matrix with respect to the dimension of S. Defining det(C) = 1 when rank(C) = 0, note that if rank(C) = 2 only one second-order cofactor for Z exists, and it is 1, while the K_ijkl term disappears when rank(C) < 2. Finally, eqn (59) and (60) reduce to just the T terms when rank(C) = 0, as expected. Some more simplifications can be done when making the spin partitioning of S explicit, and by the use of the Cholesky representation of the electron repulsion integrals,⁵⁶ see the ESI.† The equations for the evaluation of the magnetic moment and g-tensor within the MC-NOGF approach are also reported in the ESI.†

Computational details

Calculations were carried out on a set of Ln(III) complexes to analyze how the FAIMP and NOGF approximations to the molecular wave function affect the energy gaps within the lowest energy spin–orbit multiplet and the magnetic properties. The results were compared with the energies obtained at the complete active space configuration interaction with spin–orbit (CASCI-SO) level coupled with the configuration averaged Hartree–Fock (CAHF) method.^45,57

In the CAHF/CASCI-SO strategy, the set of molecular orbitals are generated at the CAHF level minimizing the average-energy functional represented on the basis of all possible Slater determinants, of any M_S quantum number, built up allowing N_act active electrons to be distributed in all the possible ways in M_act active orbitals. Then, CAHF orbitals are used in the CASCI-SO step to construct the representation of the total Hamiltonian of the system, which includes both the Born–Oppenheimer electrostatic and spin–orbit Hamiltonian, still on the basis of all possible Slater determinants in the CAS space of N_act electrons in M_act orbitals. Finally, the total Hamiltonian is diagonalized to obtain the energy levels.

The number M_act = 7 active orbitals has been used for all the Ln(III) ions studied (i.e. 4f atomic shell), and the number N_act = 8, 9, 10, and 11 active electrons has been used for Tb(III), Dy(III), Ho(III), and Er(III) ions, respectively.

Within the FAIMP approach, the crystal field levels and magnetic proprieties are estimated from the wave function of the metal ψ^M_m, optimized at the CAHF/CASCI-SO level in the presence of the ligands described by the FAIMP model as

where ψ^R₀, ψ^S₀, and ψ^T₀ indicate the ground state closed shell single Slater determinants of the separated molecules of the ligands optimized at the HF level and ψ^W₀ indicates the ground state closed shell single Slater determinant of the whole group of the ligands optimized at the HF level. indicates that the functions on the left-hand side affect the function on the right-hand side by the FAIMP model, and indicates that the functions on both the left- and right-hand sides are adjusted to each other by a self-consistent procedure.

Within the NOGF theory presented in this paper, the crystal field levels and magnetic proprieties are estimated by the variational optimization of the pure electrostatic Hamiltonian on the basis of the product functions in eqn (28), which includes (i) the wave functions of the isolated metal ψ^M_m optimized at the CAHF/CASCI-SO level and (ii) the ground state closed shell single Slater determinants optimized at the HF level describing the ligands as

Note that in NOGF3 the functions ψ^M_m included in the product functions are the wave functions of the isolated metal without any influence from the ligands, whereas ψ^M_m used in the self-consistent procedure to adjust ψ^W₀ is discharged after the optimization.

The geometries of the Ln complexes analyzed have been fixed to their experimental crystallographic structures reported in the literature.^58–66 We indicate the ligands with acac = acetylacetonate, dppz = dipyridophenazine, dpq = dipyridoquinoxaline, phen = 1,10-phenanthroline, hfac = hexafluoroacetylacetonate, glyme = dimethoxyethane, paah = N-(2-pyridyl)-ketoacetamide, tta = 2-thenoyltrifluoroacetonate, bipy = 2,2′-bipyridine, meOH = methanol, H₃ trensal = 2,2′,2′′-tris(salicylideneimido)trimethylamine, and tfpb = 4,4,4-trifluoro-1-phenylbutane-1,3-dionato, see Fig. 1.

Fig. 1 Sketch, mark and color code for the molecules studied in this work.

All atoms were described by the ANO-RCC basis set,⁶⁷ with the contraction [8s7p5d3f2g1h] for Ln atoms, [4s3p2d1f] for C, N, and O atoms coordinating the metal, [3s2p] for C, N, and O atoms not coordinating the metal, and [2s] for H.

Scalar relativistic terms were included in the one-electron part of the electrostatic Hamiltonian in all the HF, CAHF, and CASCI-SO computations, within the second order Douglas–Kroll–Hess (DKH2) approximation.⁶⁸ In CASCI-SO evaluations, the spin–orbit interaction was included by the atomic approximation to the Breit–Pauli Hamiltonian, in which the two-electron spin–orbit integrals involving the atomic basis functions on multiple centers are discharged. The Cholesky representation of the electron repulsion integrals was used to speed up the calculations, with δ = 10⁻⁸.⁵⁶

All the calculations were performed using the software package CERES, an ab initio quantum chemistry package specifically designed for the calculation of the electronic structure and magnetic properties of lanthanide complexes.

The errors affecting the energy gaps ΔE^APP_i = E^APP_i − E^APP₀, calculated within a given spin–orbit multiplet M = (2S + 1)L_J by the use of the approximation APP, were analyzed by the statistical indicators

where ΔE^BENC_i = E^BENC_i − E^BENC₀ are the energy gaps calculated by the use of the true H^BP_SO, and N is the total number of states in the multiplet (minus one for the ground state in the ground spin–orbit multiplet). The results were graphically represented by the use of the probability density of the normal distribution:

Results

Fig. 2–5 present the crystal field levels within the ground spin–orbit multiplets of 44 different Tb(III)/Dy(III)/Ho(III)/Er(III) complexes estimated by the use of the different methods presented in the previous sections. In the figures, the energy gaps calculated at FAIMP or NOGF levels are reported on the x-axis, while the results obtained from CAHF/CASCI-SO simulations on the integral systems are reported on the y-axis for comparison (see also Tables S1–S4 in the ESI†).

Fig. 2 Graphical representation of the crystal field levels within the ground spin–orbit multiplet ⁷F₆ of different Tb(III) complexes, estimated at the CAHF/CASCI-SO level on the whole system (y axis) versus the different approximations presented in this work (x axis). The mean μ and standard deviation σ for the %errors affecting the energy gaps are represented by the normal distribution function.

Fig. 3 Graphical representation of the crystal field levels within the ground spin–orbit multiplet ⁶H_15/2 of different Dy(III) complexes, estimated at the CAHF/CASCI-SO level on the whole system (y axis) versus the different approximations presented in this work (x axis). The mean μ and standard deviation σ for the %errors affecting the energy gaps are represented by the normal distribution function.

Fig. 4 Graphical representation of the crystal field levels within the ground spin–orbit multiplet ⁵I₈ of different Ho(III) complexes, estimated at the CAHF/CASCI-SO level on the whole system (y axis) versus the different approximations presented in this work (x axis). The mean μ and standard deviation σ for the %errors affecting the energy gaps are represented by the normal distribution function.

Fig. 5 Graphical representation of the crystal field levels within the ground spin–orbit multiplet ⁴I_15/2 of different Er(III) complexes, estimated at the CAHF/CASCI-SO level on the whole system (y axis) versus the different approximations presented in this work (x axis). The mean μ and standard deviation σ for the %errors affecting the energy gaps are represented by the normal distribution function.

In FAIMP1 and NOGF1 strategies the metal interacts with the electron densities of the separated molecules, which are included in the environment after separate optimisation. For FAIMP1 simulations, the graphs show a general large underestimation of the energy gaps with respect to the CAHF/CASCI-SO results for all the studied ions, with mean percentage errors of about −80%/−90%. These results can be attributed to a poor description of the ligands' density, which, as expected, cannot be described as the simple sum of the densities of the isolated molecules, but the latter must adjust to each other within the complex structure. However, it is noteworthy that generally slightly better results are obtained when the densities of the isolated molecules are treated within the NOGF method (see NOGF1 results). This improvement is especially observed for the Tb-pc₂⁻ complex, which is formed by two highly π conjugated pc⁻² ligands, where the NOGF1 methodology reproduces the CAHF/CASCI-SO results with a mean %error smaller than −20% and a small standard deviation of about 4%.

A general improvement in the results is obtained when the ligands are described by a single wave function optimized on the whole set of atoms present in the environment, and the presence of a FAIMP model potential describing the metal ion. This approach gives rise to the FAIMP2 and NOGF2 methods. In this case, while the percentage errors affecting the crystal field levels calculated by the FAIMP2 method are still rather large for all the ions, the NOGF2 strategy performs significantly better. In NOGF2 in fact the wavefunction of the ligands and that of the isolated metal are combined together in a single antisymmetrised non-orthogonal product function, which leads to results in rather good agreement with those obtained at the full CAHF/CASCI-SO level on the whole set of systems explored here. This not only proves the importance of including an accurate description for the electron density of the ligands, but also seems to indicate that the rigid interaction (i.e. no orbital relaxation) between the metal density optimized on the isolated ion and a suitable density which accurately describes the atoms in the environment already suffices for the modelling of the crystal field splittings in these systems.

Finally, the third family of approximations explored here is named FAIMP3 and NOGF3, where the densities of the metal and the ligands are iteratively optimized to each other, to fully capture the polarization effects of both ligands and metal electrons, and, consequently, a refined description of electrostatic interactions between metals and ligands. The only family of non-covalent metal–ligand interactions that cannot be captured here is dispersion interactions, since the dynamical correlation is currently not included in our wave function models. We note however that in our proposed NOGF strategy such dispersion interactions can indeed be captured by going in higher order expansion of eqn (40), while this is not possible in the FAIMP approximation.

Despite the significant improvement observed in both FAIMP3 and NOGF3 approaches, note that the FAIMP3 strategy still underestimates the energy gaps. On the other hand, the NOGF3 method shows crystal field gaps almost coincident with the CAHF/CASCI-SO results for the majority of the Tb(III), Dy(III), and Er(III) complexes in all the range of frequencies. Slightly larger discrepancies are shown by paah₂-2NO₃-meOH ligands, where an average overestimation of the energy gaps of about +40% and +17% is found for Tb(III) and Dy(III) ions, respectively, and an underestimation of about −26% for the Er(III) metal. Differently, the NOGF3 strategy does not lead to a significant improvement of the crystal field energies of all Ho(III) complexes.

We also note that Dy-trensal and Er-trensal show very different behaviours. On the one hand, the Dy(III) system is accurately reproduced by the NOGF3 method both in the low and high energy ranges. On the other hand, NOGF3 for Er-trensal shows quite large discrepancies with the CAHF/CASCI-SO energies, especially for the levels above 400 cm⁻¹, leading to a large underestimation of the NOGF3 crystal field gaps when compared with the fully interacting CAHF/CASCI-SO results.

Besides the crystal field energy spectrum, we also tested the proposed methods for the calculation of the magnetic properties. In particular, in Tables 1 and 2 we report the principal values of the g tensor for the ground and first excited pseudo-Kramers/Kramers doubles (KD) for Tb(III)/Dy(III) complexes, estimated via different approximations. From the tables, the discrepancies between the NOGF3 and CAHF/CASCI-SO results for the ground KDs never exceed 1.5% for all the Tb(III) and Dy(III) systems. A very good agreement is also found in the orientation of the main magnetic axis, where the angle between the axes estimated by the NOGF3 and CAHF/CASCI-SO methods is in all cases smaller than 8 degrees, see Tables 3 and 4. Larger discrepancies are observed for the first excited KDs, where now the g-tensor calculated at the NOGF3 level is affected by errors of about 5%/10%, and the angles of the main magnetic axes are between 5 and 15 degrees, with the exception of acac₃-dpq (22°), acac₃-dppz (24°), tta₃-phen (39°), and hfac₃-glyme (21°) ligands for Tb(III), and the tta3bipy ligand for Dy(III).

Table 1 g tensor values for the two lowest lying pseudo-Kramers doublets (KD) within the ground spin–orbit multiplet ⁷F₆ of different Tb(III) complexes

		CAHF/CASCI-SO	FAIMP1	FAIMP2	FAIMP3	NOGF1	NOGF2	NOGF3
acac3-2H2 O	KD1	17.578	+0.143	−0.575	+0.054	+0.006	−0.567	+0.100
	KD2	17.061	−2.928	−5.403	−3.722	−3.454	−3.404	−2.489
acac3-dpq	KD1	17.335	−0.356	+0.321	+0.455	−0.358	−0.201	+0.072
	KD2	13.189	−0.573	+0.687	+1.089	−0.862	−0.102	+0.664
acac3-dppz	KD1	17.620	−0.890	+0.077	+0.251	−2.756	−0.061	+0.143
	KD2	17.109	−4.348	−3.384	−2.556	−6.072	−2.372	−1.655
acac3-phen	KD1	17.463	−0.109	+0.357	+0.392	−0.129	+0.055	+0.187
	KD2	14.755	−1.748	−0.304	−0.228	−1.401	−0.574	−0.286
pc2^−	KD1	17.918	−0.472	+0.020	−0.038	+0.019	+0.020	+0.002
	KD2	14.546	+1.543	+0.310	+0.180	+0.358	+0.407	+0.210
paah2-2NO3-meOH	KD1	17.772	−0.316	+0.146	+0.157	+0.039	+0.137	+0.149
	KD2	15.199	+0.746	−0.392	−0.382	−0.167	−0.144	−0.197
tta3-bipy	KD1	17.759	−1.175	+0.128	+0.173	−0.822	+0.027	+0.095
	KD2	15.569	−3.282	−0.595	−0.661	−0.583	−1.355	−0.902
tta3-phen	KD1	17.608	−0.217	+0.138	+0.255	−0.645	−0.039	+0.097
	KD2	13.736	+0.636	+0.424	+0.858	−0.973	−0.421	+0.448
tfpb3-dppz	KD1	17.484	−0.285	+0.057	+0.341	−1.459	−0.147	+0.100
	KD2	12.229	+1.537	+1.191	+2.290	+0.038	−0.394	+0.996
hfac3-glyme	KD1	17.630	−1.596	−0.061	+0.168	−1.774	−0.144	+0.035
	KD2	13.774	−2.114	−0.427	+0.574	−2.462	−1.070	−0.029

---

Table 2 Main g tensor values for the two lowest lying Kramers doublets (KD) within the ground spin–orbit multiplet ⁶H_15/2 of different Dy(III) complexes

		CAHF/CASCI-SO	FAIMP1	FAIMP2	FAIMP3	NOGF1	NOGF2	NOGF3
acac3-2H2 O	KD1	19.471	+0.380	−0.276	+0.054	−0.627	−0.127	+0.009
	KD2	15.549	+1.647	−0.467	+0.465	−0.141	−0.564	+0.123
acac3-dpq	KD1	19.219	−0.038	+0.459	+0.517	+0.229	+0.152	+0.275
	KD2	16.230	+0.288	+0.096	+0.354	+0.000	−0.011	+0.267
acac3-dppz	KD1	19.366	+0.311	+0.131	+0.381	−0.817	−0.045	+0.159
	KD2	15.369	+1.292	+0.784	+1.460	−2.917	−0.063	+0.690
acac3-phen	KD1	19.344	+0.242	+0.310	+0.390	−0.203	+0.053	+0.192
	KD2	15.701	+0.308	+0.809	+1.034	+0.247	+0.414	+0.669
pc2^−	KD1	17.490	+2.413	+2.411	+2.410	+0.993	+0.133	+0.283
	KD2	14.690	+2.465	+2.504	+2.496	+5.126	+1.678	+1.657
paah2-2NO3-meOH	KD1	19.592	−0.391	+0.273	+0.292	+0.162	+0.138	+0.193
	KD2	16.490	+0.636	+0.648	+0.680	+0.584	+0.338	+0.501
tta3-bipy	KD1	19.576	−0.224	+0.131	+0.279	−0.616	+0.073	+0.200
	KD2	15.102	+3.354	+1.860	+2.120	−0.869	+2.047	+2.495
tta3-phen	KD1	19.485	+0.296	+0.263	+0.317	−0.079	+0.130	+0.224
	KD2	15.674	+2.944	+1.249	+1.305	+1.226	+1.038	+1.158
tfpb3-dppz	KD1	19.296	+0.190	+0.363	+0.504	−0.270	+0.031	+0.207
	KD2	14.549	+4.845	+2.015	+2.388	−1.588	+1.190	+1.883
hfac3-glyme	KD1	19.450	+0.336	+0.160	+0.293	−0.278	+0.086	+0.230
	KD2	15.930	−3.308	+0.548	+0.889	+0.532	+0.522	+0.789

---

Table 3 Main magnetic axes for the two lowest lying pseudo-Kramers doublets (KD) within the ground spin–orbit multiplet ⁷F₆ of different Tb(III) complexes: angle (in degrees) between the axes estimated by the different methods proposed in this paper and the ab initio computation on the whole systems at the CAHF/CASCI-SO level

		FAIMP1	FAIMP2	FAIMP3	NOGF1	NOGF2	NOGF3
acac3-2H2O	KD1	75	36	23	81	6	8
	KD2	59	27	32	43	37	22
acac3-dpq	KD1	45	14	19	26	4	7
	KD2	40	5	6	16	10	6
acac3-dppz	KD1	89	18	20	66	3	6
	KD2	84	39	40	57	23	24
acac3-phen	KD1	57	12	13	24	2	4
	KD2	31	20	20	30	15	15
pc2^−	KD1	17	5	9	3	0	2
	KD2	17	5	4	7	0	9
paah2-2NO3-meOH	KD1	10	9	9	7	5	6
	KD2	86	11	12	1	7	9
tta3-bipy	KD1	81	8	8	72	1	1
	KD2	38	35	41	61	32	39
tta3-phen	KD1	42	7	9	19	2	3
	KD2	58	16	11	50	15	9
tfpb3-dppz	KD1	62	8	10	12	3	3
	KD2	78	3	9	53	5	5
hfac3-glyme	KD1	89	6	8	46	2	3
	KD2	90	21	20	41	24	21

---

Table 4 Main magnetic axes for the two lowest lying Kramers doublets (KD) within the ground spin–orbit multiplet ⁶H_15/2 of different Dy(III) complexes: angle (in degrees) between the axes estimated by the different methods proposed in this paper and the ab initio computation on the whole systems at the CAHF/CASCI-SO level

		FAIMP1	FAIMP2	FAIMP3	NOGF1	NOGF2	NOGF3
acac3-2H2O	KD1	75	11	5	28	2	2
	KD2	73	8	2	49	2	5
acac3-dpq	KD1	51	15	15	18	6	6
	KD2	43	12	11	13	10	11
acac3-dppz	KD1	86	10	11	11	2	3
	KD2	83	7	9	16	4	8
acac3-phen	KD1	70	14	13	19	6	5
	KD2	23	15	18	9	8	13
pc2^−	KD1	12	9	9	21	2	4
	KD2	4	9	10	16	5	6
paah2-2NO3-meOH	KD1	37	4	3	2	2	2
	KD2	28	1	1	14	0	0
tta3-bipy	KD1	85	15	19	9	4	7
	KD2	79	36	32	62	24	25
tta3-phen	KD1	78	19	16	20	6	6
	KD2	45	23	23	24	13	15
tfpb3-dppz	KD1	89	16	18	10	4	5
	KD2	50	22	21	30	9	12
hfac3-glyme	KD1	73	14	16	19	6	8
	KD2	39	16	22	14	10	16

---

All results presented here clearly show that our proposed NOGF methodologies for the description of non-covalent crystal field interactions in lanthanide complexes perform significantly better than the FAIMP methods of equivalent quality, and in fact, quite remarkably, they often manage to reproduce the fully interacting results obtained via CAHF/CASCI-SO to an unexpected degree of accuracy. To appreciate the meaning of these results, it is important to remember that the CAHF/CASCI-SO and CASSCF/SI-SO methods, while neglecting the dynamical correlation between metal and ligand electrons, do provide a mean-field description of metal–ligand charge-transfer excitations, and hence of covalency effects, via the refinement of non-redundant active-to-virtual and inactive-to-active orbital rotations carried out during the orbital optimisation process (in other words, via the optimisation of orbitals containing a small degree of metal–ligand hybridisation). On the other hand, FAIMP and NOGF by definition do not formally account for charge transfer processes between different groups; hence they are formally incapable of describing any metal–ligand interactions other than non-covalent interactions. Thus two points of discussion immediately arise from these results.

First, the observation that NOGF performs consistently better than the FAIMP version of equivalent quality would seem to suggest that in NOGF, the process of building antisymmetrised product functions from non-orthogonal wave functions associated with the single groups, and of producing linear combinations of such product functions to form the eigenstates of the first-order degenerate perturbation theory interaction Hamiltonian, can lead to a recovery of some degree of effective metal–ligand hybridisation, as much as it is necessary to ensure orthogonality between different orbitals in an effective orbital-delocalised description leading to equivalent results. The minimal hybridisation required between localised non-orthogonal orbitals within an orthogonalisation process can be considered as a component of covalency as described by an orthogonal and hybridised basis, accounting for some charge transfer processes between metal-centered and ligand-centered orbitals. One should also consider that a rather significant difference between NOGF and FAIMP methods is not fully unexpected, given that non-orthogonality makes the matrix elements of the electrostatic Hamiltonian between group product functions significantly more complex and rich in information about group–group interactions with respect to the effective matrix elements of the FAIMP method, where the strong orthogonality constraint makes the Hamiltonian describing the interaction between metal and ligands eqn (23) significantly simpler than the full non-orthogonal results of the NOGF method, eqn (48) and (49). After all, FAIMP can be considered as an approximation to the full NOGF method, hence bound to lead to a poorer performance.

Second, the observation that, generally speaking, the NOGF3 approach appears to reproduce to a high degree of accuracy the fully interacting CAHF/CASCI-SO results suggests that fully-interacting multiconfigurational ab initio methods ignoring the dynamical correlation, such as CAHF/CASCI-SO or CASSCF/SI-SO, lead to a description of the Ln(III)–ligand interactions that is strongly dominated by intermolecular interactions, and thus probably tend to underestimate covalency effects, so much so that for most systems explored here the mean field covalency effects captured by the molecular orbital optimisation process are captured also by our non-covalent NOGF methods.

However, we have also exposed here a few cases where even our best NOGF approach displays some difficulties in recovering the fully interacting CAHF/CASCI-SO results. We can identify here a major systematic trend in these more difficult cases. We note in fact that, generally speaking, the NOGF methodology tends to perform much better for complexes of Dy(III) and Tb(III) ions, rather than for equivalent complexes of Er(III) and Ho(III) ions. We provide here a rationalisation of this behaviour in terms of the multipolar expansion of the 4f-charge distributions associated with the M_J states of a spin–orbit multiplet with a total angular momentum quantum number J. As shown in fact first by Sievers,¹⁰ the asphericity of the 4f charge density distribution, ρ_J,MJ(θ, ϕ), can be expanded in terms of just three spherical harmonics C_kq(θ, ϕ), corresponding to the 4f multipolar expansion of the charge density associated with a given state M_J as in eqn (70)

ρ_J,MJ(θ,ϕ) = A_2,MJC₂₀(θ, ϕ) + A_4,MJC₄₀(θ,ϕ) + A_6,MJC₆₀(θ, ϕ), (70)

where A_2,MJ is the quadrupole coefficient for the state M_J, while A_4,MJ and A_6,MJ are the hexadecapole and hexacontadecapole coefficients, respectively, for the same state. In the approximation that the J-multiplet is dominated by a single Russell–Saunders term with specific values of the L and S total angular momentum and spin quantum numbers (known to be a very good approximation for complexes of trivalent lanthanide ions), the three multipolar coefficients can be exactly calculated asEqn (71) shows that the magnitude of a given multipolar term describing the charge density of different crystal field states for a given lanthanide, aside from a simple Wigner 3j symbol differentiating different M_J states within the same multiplet, is essentially determined by the reduced matrix element 〈LSJ||C^k||LSJ〉, for which an expression can be found e.g. in ref. 10. The direct calculation of this reduced matrix element determining the magnitude of the quadrupolar component of the charge density (k = 2), for ions Tb (L = 3, S = 3, J = 6), Dy (L = 5, S = 5/2, J = 15/2), Ho (L = 6, S = 2, J = 8), and Er (L = 6, S = 3/2), gives

〈LSJ‖C²‖LSJ〉 = −0.42(Tb); −0.40(Dy); −0.16(Ho); 0.16(Er). (72)

Furthermore, within a simple point charge model, we can easily show that the strength of the electrostatic interactions between metal and ligands in lanthanides is strongly dominated by the quadrupolar term. This can be shown in the following manner. If a ligand charge is at a distance R from the lanthanide ion, the strength of the electrostatic potential generated at distance R by the 4f-multipole of order k associated with the 4f charge density distribution can be estimated as . For a typical metal–ligand distance of R = 5a₀(2.65 Å), and using the radial multipoles tabulated for all trivalent lanthanide ions by Edvardsson and Klintenberg,⁶⁹ this leads to rather small ratios between the hexadecapolar (k = 4) and quadrupolar (k = 2) crystal field strengths of , and to even smaller ratios between hexacontadecapolar (k = 6) and quadrupolar terms for all four ions: , which is evidence that metal–ligand electrostatics are dominated by the k = 2 quadrupolar interactions (quadrupole–charge interactions in this simple example).

This simple argument thus shows that intermolecular electrostatic interactions between the ligands and the lanthanide ions in a complex are dominated by quadrupolar terms, and thus are much stronger in Dy and Tb than in Ho and Er, where the strength of the quadrupolar terms of the 4f charge density expansion calculated in eqn (72) is found to be less than half the magnitude of those found for Dy and Tb. The fact that in Er and Ho the interactions between metal and ligands will be less strongly dominated by pure electrostatics provides a rationalisation of our finding that ab initio non-covalent approaches like NOGF have more difficulties in providing a satisfactory description of metal–ligand interactions with respect to methods like CAHF/CASCI-SO, in which covalent ligand–metal interactions can be partially captured via optimisation of ligand–metal orbital hybridisation.

At this time, we were not able to rationalise the relatively poor performance of the NOGF methods for Dy(III) and Tb(III) for the single case of the paah₂-2NO₃-meOH ligands.

Conclusions

We presented a novel multiconfigurational non-orthogonal group function (MC-NOGF) approach for the calculation of the crystal field energies and magnetic properties of trivalent lanthanide complexes. The novel MC-NOGF approach has been implemented in the ab initio program Computational Emulator of Rare Earth Systems (CERES).^35,45 By direct comparison of two different fragment approaches of the crystal field energies of 44 complexes, namely the FAIMP and the here proposed MC-NOGF, we showed that the MC-NOGF model offers a better performing approach to an ab initio non-covalent description of metal–ligand interactions in trivalent lanthanide complexes. This suggests that many bonding situations in lanthanide complexes, especially involving lanthanide ions whose 4f charge density has a strong quadrupolar component such as Tb(III) and Dy(III), can be satisfactorily described in terms of “intermolecular interactions” between ligands and 4f electrons, a fact that could be used not only to better understand the chemistry of 4f coordination related to different ligands, but also to design more efficient and localised ab initio methods for the inclusion of dynamical correlation effects, which could focus on the description of those excitations that are by definition absent from NOGF strategies.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

A. S. thankfully acknowledges support from the Australian Research Council, Future Fellowship FT180100519, and Discovery Project Grant DP210103208. This research was supported by the Nectar Research Cloud (https://nectar.org.au), a collaborative Australian research platform supported by the National Collaborative Research Infrastructure Strategy (NCRIS).

References

1. R. M. Pallares and R. J. Abergel, Nanoscale, 2020, 12, 1339–1348.
2. I. Georgieva, T. Zahariev, A. J. A. Aquino, N. Trendafilova and H. Lischa, Spectrochim. Acta, Part A, 2020, 240, 118591.
3. A. G. Cosby, J. J. Woods, P. Nawrocki, T. J. Sørensen, J. J. Wilson and E. Boros, Chem. Sci., 2021, 12, 9442–9451.
4. S. E. Bodman, C. Breen, S. Kirkland, S. Wheeler, E. Robertson, F. Plasser and S. J. Butler, Chem. Sci., 2022, 13, 3386–3394.
5. J. Tang and P. Zhang, Lanthanide Single Molecule Magnets, Springer-Verlag Berlin Heidelberg, 1st edn, 2015.
6. C. A. P. Goodwin, F. Ortu, D. Reta, N. F. Chilton and D. P. Mills, Nature, 2017, 548, 439–442.
7. F. S. Guo, B. M. Day, Y. C. Chen, M. L. Tong, A. Mansikkamaki and R. A. Layfield, Science, 2018, 362, 1400.
8. H. A. Bethe, Splitting of terms in crystals: Complete English translation, Consultants Bureau, 1929.
9. J. D. Rinehart and J. R. Long, Chem. Sci., 2011, 2, 2078–2085.
10. J. Sievers, Z. Phys. B, 1982, 45, 289.
11. N. F. Chilton, D. Collison, E. J. L. McInnes, R. E. P. Winpenny and A. Soncini, Nat. Commun., 2013, 4, 2551.
12. C. Gao, A. Genoni, S. Gao, S. Jiang, A. Soncini and J. Overgaard, Nat. Chem., 2020, 12, 213–219.
13. M. T. Hutchings and D. K. Ray, Proc. Phys. Soc., 1963, 81, 663–676.
14. D. K. Ray, Proc. Phys. Soc., 1963, 82, 47–57.
15. C. J. Lenander and E. Y. Wong, J. Chem. Phys., 1963, 38, 2750–2752.
16. R. E. Watson and A. J. Freeman, Phys. Rev., 1964, 134, A1526–A1546.
17. K. Rajnak and B. G. Wybourne, J. Chem. Phys., 1964, 41, 565–569.
18. D. Newman, Adv. Phys., 1971, 20, 197–256.
19. C. A. Morrison and R. P. Leavitt, Handbook on the Physics and Chemistry of Rare Earths. vol. 5, North Holland, 1982, vol. 5, pp. 461–692.
20. O. Malta, S. Ribeiro, M. Faucher and P. Porcher, J. Phys. Chem. Solids, 1991, 52, 587–593.
21. L. Ungur and L. F. Chibotaru, Chem. – Eur. J., 2017, 23, 3708–3718.
22. R. Alessandri, H. Zulfikri, J. Autschbach and H. Bolvin, Chem. – Eur. J., 2018, 24, 5538–5550.
23. L. Ungur, Lanthanide-Based Multifunctional Materials, Elsevier, 2018, pp. 1–58.
24. J. D. Weeks and S. A. Rice, J. Chem. Phys., 1968, 49, 2741–2755.
25. S. Huzinaga and A. A. Cantu, J. Chem. Phys., 1971, 55, 5543–5549.
26. S. Huzinaga, D. Mcwilliams and A. A. Cantu, Adv. Quantum Chem., 1973, vol. 7, 187–220.
27. B. Swerts, L. F. Chibotaru, R. Lindh, L. Seijo, Z. Barandiaran, S. Clima, K. Pierloot and M. F. A. Hendrickx, J. Chem. Theory Comput., 2008, 4, 586–594.
28. B. Jeziorski, R. Moszynski and K. Szalewicz, Chem. Rev., 1994, 94, 1887–1930.
29. K. Szalewicz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 254–272.
30. K. Patkowski, P. S. Zuchowski and D. G. A. Smith, J. Chem. Phys., 2018, 148, 164110.
31. P. W. Payne, J. Chem. Phys., 1982, 77, 5630–5638.
32. P. A. Malmqvist, Int. J. Quantum Chem., 1986, 30, 479–494.
33. P.-Å. Malmqvist and B. O. Roos, Chem. Phys. Lett., 1989, 155, 189–194.
34. A. Mitrushchenkov and H.-J. Werner, Mol. Phys., 2007, 105, 1239–1249.
35. A. Soncini, S. Calvello, M. Piccardo and S. V. Rao, Ceres, an ab initio quantum chemistry package for the electronic structure and magnetic properties of lanthanide complexes, 2018.
36. A. L. Tchougréeff, Hybrid Methods of Molecular Modeling, Springer, Netherlands, Dordrecht, 2008, vol. 17.
37. J. G. Ángyán and G. Náray-Szabó, in Theoretical Treatment of Large Molecules and Their Interactions, ed. Z. B. Maksic, Springer Berlin Heidelberg, Berlin, Heidelberg, 1991, pp. 1–49.
38. R. McWeeny and B. T. Sutcliffe, Methods of molecular quantum mechanics, Academic Press, San Diego, 1969.
39. E. L. Mehler, J. Math. Chem., 1992, 10, 57–91.
40. E. L. Mehler, J. Chem. Phys., 1981, 74, 6298.
41. E. L. Mehler, J. Chem. Phys., 1977, 67, 2728.
42. C. Amovilli and R. McWeeny, Chem. Phys., 1990, 140, 343–361.
43. R. McWeeny, Proc. R. Soc. A, 1959, 253, 242–259.
44. T. Shimazaki, K. Kitaura, D. G. Fedorov and T. Nakajima, J. Chem. Phys., 2017, 146, 084109.
45. S. Calvello, M. Piccardo, S. V. Rao and A. Soncini, J. Comput. Chem., 2018, 39, 328–337.
46. P.-O. Löwdin, J. Chem. Phys., 1950, 18, 365–375.
47. W. J. Carr, Phys. Rev., 1953, 92, 28–35.
48. P.-O. Löwdin, Phys. Rev., 1955, 97, 1490–1508.
49. P.-O. Löwdin, Phys. Rev., 1955, 97, 1474–1489.
50. F. Prosser and S. Hagstrom, Int. J. Quantum Chem., 1968, 2, 89–99.
51. S. C. Leasure and G. Balint-Kurti, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 2107–2113.
52. J. Verbeek and J. H. Van Lenthe, Theochem, 1991, 229, 115–137.
53. I. Hayes and A. Stone, Mol. Phys., 1984, 53, 69–82.
54. G. Figari and V. Magnasco, Mol. Phys., 1985, 55, 319–330.
55. W. W. Hager, SIAM Rev., 1989, 31, 221–239.
56. M. Piccardo and A. Soncini, J. Comput. Chem., 2017, 38, 2775–2783.
57. W. Van den Heuvel, S. Calvello and A. Soncini, Phys. Chem. Chem. Phys., 2016, 18, 15807–15814.
58. S.-D. Jiang, B.-W. Wang, G. Su, Z.-M. Wang and S. Gao, Angew. Chem., 2010, 122, 7610–7613.
59. G.-J. Chen, C.-Y. Gao, J.-L. Tian, J. Tang, W. Gu, X. Liu, S.-P. Yan, D.-Z. Liao and P. Cheng, Dalton Trans., 2011, 40, 5579.
60. G.-J. Chen, Y.-N. Guo, J.-L. Tian, J. Tang, W. Gu, X. Liu, S.-P. Yan, P. Cheng and D.-Z. Liao, Chem. – Eur. J., 2012, 18, 2484–2487.
61. E. M. Fatila, E. E. Hetherington, M. Jennings, A. J. Lough and K. E. Preuss, Dalton Trans., 2012, 41, 1352–1362.
62. N. F. Chilton, S. K. Langley, B. Moubaraki, A. Soncini, S. R. Batten and K. S. Murray, Chem. Sci., 2013, 4, 1719.
63. Z.-G. Wang, J. Lu, C.-Y. Gao, C. Wang, J.-L. Tian, W. Gu, X. Liu and S.-P. Yan, Inorg. Chem. Commun., 2013, 27, 127–130.
64. Y. Bi, Y.-N. Guo, L. Zhao, Y. Guo, S.-Y. Lin, S.-D. Jiang, J. Tang, B.-W. Wang and S. Gao, Chem. – Eur. J., 2011, 17, 12476–12481.
65. M. Habib, S. Sain, B. Das and S. Chandra, J. Indian Chem. Soc., 2011, 88, 1501–1508.
66. Q. Sun, P. Yan, W. Niu, W. Chu, X. Yao, G. An and G. Li, RSC Adv., 2015, 5, 65856–65861.
67. B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov, P.-O. Widmark and A. C. Borin, J. Phys. Chem. A, 2008, 112, 11431–11435.
68. M. Reiher and A. Wolf, Relativistic quantum chemistry: the fundamental theory of molecular science, Wiley-VCH, 2015.
69. S. Edvardsson and M. Klintenberg, J. Alloys Compd., 2013, 275–277, 230–233.

---

Footnotes

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

‡ These authors contributed equally to this work.

---