Ab Initio Non-Covalent Crystal Field Theory for Lanthanide Complexes: A Multiconfigurational Non- Orthogonal Group Functions Approach†

Lanthanide Complexes: A Multiconfigurational NonOrthogonal Group Functions Approach† Alessandro Soncini⇤,‡ and Matteo Piccardo‡ We present a non-orthogonal fragment ab initio methodology for the calculation of crystal field energy levels and magnetic properties in lanthanide complexes, implementing a systematic description of non-covalent contributions to metal-ligand bonding. The approach has two steps. In the first step, appropriate ab initio wavefunctions for the various ionic fragments (lanthanide ion and coordinating ligands) are separately optimized, accounting for the electrostatic 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-NOGF), and reduced to an effective (2J + 1)–dimensional non-orthogonal Configuration Interaction (CI) problem via Löwdin-partitioning. Within the proposed formalism, the projected Hamiltonian can be implemented to any desired order of perturbation theory in the fragment-localised excitations out of the degenerate space, and its eigenvalues and eigenfunctions are systematic approximations to the crystal field energies and wavefunctions. We present a preliminary implementation of the proposed MC-NOGF method to 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 with the full CAHF/CASCI-SO method, accounting for metal-ligand covalency in a mean-field manner. We find that energies and magnetic properties for 44 complexes obtained via an iteratively optimized version of our MC-NOGF first-order non-covalent method, compare remarkably well to the full CAHF/CASCI-SO method including metal-ligand covalency, and are superior to the best purely electrostatic results achieved via an iteratively optimized version of the FAIMP approach. sec:intro Introduction


Introduction
Trivalent lanthanide ions play a special role in magnetism, thanks to their large magnetic moments and magnetic anisotropy.At the root of their enhanced magnetic properties are the rather weak interactions between the 4f electron shell and its coordination environment, mainly as a result of the 4f-shell radially compact nature, and its consequent efficient shielding from the surrounding electrostatic and chemical environment by the outer 5s 2 5p 6 orbital shells 1 .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 School of Chemistry, The University of Melbourne, Australia.⇤ E-mail: asoncini@unimelb.edu.au‡ These authors contributed equally to this work.† Electronic Supplementary Information (ESI) available: [details of any supplementary information available should be included here].See DOI: 00.0000/00000000.Crystal Field Theory (CFT) 2 .CFT describes the splitting of the 4f many-electron wavefunctions associated to 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 euristic value should not be underestimated, even today within the field of molecualr magnetism.
For instance, Rinehart and Long 3 , discussing the ashpericities of the electron density of the magnetic states in trivalent lanthanide ions calculated by Sievers some 20 years earlier, pointed out how a proper axial (equatorial) distribution of the charge den-sity 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 metal and ligand, has proven rather useful in guiding synthetic chemists towards ever more efficient Ln-based SMMs.
These simple electrostatics arguments were also successfully translated into a more quantitative classical electrostatics repulsion energy minimization procedure, which showed to be capable of semi-quantitatively predicting the specific direction of the magnetic anisotropy axis in many Dy(III) single-molecule magnets.4 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 flavours of so-called covalent contributions (charge transfer excitations of various kinds).[5][6][7][8][9][10][11][12][13][14] .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 have proved to be a reliable tool for the description of low-lying multiplets of lanthanide complexes with good accuracy 15 .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 13 .In that paper, the authors analyzed the ligand field generated by both different atomic charges distributions, and a more rigorous approach based on the Fragment Ab Initio Model Potential (FAIMP) method [16][17][18][19] .Besides the electrostatic interaction, the FAIMP approach allows to include (i) the exchange interaction between the electrons of the ligand and the metal, and (ii) the 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 a FAIMP approach, does not suffice for an accurate description of the crystal field levels in Ln(III) complexes, 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 of a metal ion in an ab initio potential created by the ligand, we start 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-ligands 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 for overlap effects which arise from the non-orthogonal theoretical framework.After presenting the equations for an ab initio multiconfiguration development involving matrix elements between non-orthogonal sets of orbitals, as we implemented in our own quantum chemistry code CERES 20 , 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 proprieties of a set of Ln(III) complexes.

Hamiltonian and product functions
Within the Born-Oppenheimer approximation, the non-relativistic electronic Hamiltonian of a N TOT interacting electrons can be defined in atomic units as: eq:ham eq:ham with: eq:h eq:h where - is the electronnuclei potential, and r i j = r i r j .Within the group function (GF) approximation 21-28 , the total Hamiltonian H can be represent on the basis of the antisymmetrized products of the wave functions for the single parts composing the system as: eq:Psik eq:Psik with y R r describing the group R, which collects N R electrons in the state r, and A is the partial antisymmetrizer operator, which exchanges the electron coordinates x between two, or more, different groups.The group functions y R r are individually antisymmetric with respect to the electron exchange: y R r (...,x i ,...,x j ,...) = y R r (...,x j ,...,x i ,...) eq:anti eq:anti and ortho-normalized: eq:gfortho eq:gfortho with d i j Kronecker delta.The product function Y k is fully characterized once the set of states k ⌘ (r, s,...) is made explicit.No inter-molecular correlation is explicitly included in a wave function of the form in eq. ( 3), since electrons belonging to different groups are assumed to move independently, but a large part of intra-molecular correlation may be accounted for by admitting a some degree of configuration interaction within each electron group. sec:ortho

Strongly-orthogonal group functions and the FAIMP approach
A convenient way to handle fragment approaches based on group functions is to require that the group functions y R r fulfil a socalled strong-orthogonality condition 23 : eq:strorth eq:strorth so that integrating over any one variable x 1 common to two different group functions make the integral vanish identically for all values of the other variables x i .This is a much stronger condition than the orthonormality requirement in eq. ( 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 eq. ( 3), and satisfies the conditions in eqs.( 4), (5), and (6), the system is group-separable 23 .This means that the total energy of the system E can be written as (after integration over the spin variables): eq:entot eq:entot where: eq:HR eq:HR with P R r spinless one-particle density matrix of the function y R r .H R in eq. ( 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 the nuclei of the total system.
Let us assume now that y R r are functions of the sets of mono- and spin h R i .Then eq. ( 6) is automatically satisfied if we impose the orthonormality conditions: Starting from the total wave function in eq. ( 3), and defining the groups R, S,... as doubly occupied closed-shell singlet systems described by the functions: eq:psiR2 eq:psiR2 where 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 out that an optimum {f R i } set which minimize the energy E while holding fixed the {f S i } sets with S 6 = R, and imposing the conditions in eqs.( 12) and ( 13), is given by 16-18 : " eq:euler4 eq:euler4 F is the Fock operator defined by: eq:fock2 eq:fock2 with h eff (i) given by: J and K are the Coulomb and exchange operators, respectively, describing the effective field at point i due to the electrons in the groups S: eq:Coul2 eq:Coul2 eq:Exc2 eq:Exc2 Note that F is formally identical with the Fock operator for the total wave function Y k .The projector P S originates from the orthonormality conditions, and reads [16][17][18] : eq:prj eq:prj Eq. ( 15) can be used as self-consistent field (SCF) equation for equation to determine the best orbitals for the group R interacting with the closed-shell singlets included as one or more groups S, enforcing orthogonality between the orbitals of all groups.Eq. 15 is the starting point of all the developments that go under the name Fragment Ab Initio Model Potential (FAIMP) 19,29 .
Let us now introduce a group M which is better described by a multi-configurational wave function y M m defined as: eq:wfm2 eq:wfm2 where C mm 0 is the weight of the m 0 -th Slater determinant y M m 0 as given in eq. ( 14) built from the set of othonormal spin-orbitals {f M i }.As usual when dealing with a multi-configurational functions, 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 y M m 0 , called inactive orbitals, (ii) a second set collecting N act M orbitals which allow for an occupation number from 0 to 2, called active orbitals.Introducing y M m in the product function Y k defined by eq. ( 3) we obtain: where now l ⌘ (m 0 , r, s ...) and k ⌘ (m, r, s ...).The problem of finding on optimal set {f R i } for the closed-shell singlet group R, still described by a single Slater determinant wave function y R r , leads to the Fock operator as defined in eq. ( 16), but now h eff (i) has form: eq:heff2 eq:heff2 where J S j (i)/J M j (i) and K S j (i)/K M j (i) are defined in eqs.(18) and ( 19), and r m jk is the one-electron transition density matrix for y M m , given by: which is averaged on N m states in eq. ( 23) to account for the energy degeneracy of y M m functions.In this work the projector operator P M is defined as: sec:nogf

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 13 .To improve on this model, while keep pursuing a fragment ab initio approach to describe the metalligand interactions essentially in terms of intermolecular forces, let us start again from eq.( 3), fix the sets of spin-orbitals {f 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 as: eq:Psik3 eq:Psik3 where Y k is a product function as defined in eq. ( 3), including few closed-shell singlet groups described by y R r wave functions given in eq. ( 14), i.e. the ligands, and one center characterized by the multi-configurational wave function y M m shown in eq. ( 21), i.e. the metal.We fix the spin-orbitals entering the latter functions, and the C m m 0 coefficients entering eq. ( 21), to the sets obtained by energy mimization on the isolated fragments.Under this requirement, the orthogonality relations in eqs.( 12) and ( 13) read: : We are then interested in defining the optimal coefficients D k entering eq. ( 28), which can be estimated by variational minimization of the total energy: or, equivalently, in matrix form solving the eigenvalues problem: eq:schr eq:schr where: The dimension of the variational problem in eq. ( 32) is typically very large.Hence, to reduce dimensionality, but 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) 23 : Thus an effective Hamiltonian HAA can be written for the subspace A: eq:HeffVar eq:HeffVar where: eq:PartE eq:PartE Note that HAA is itself function of the unknown energy E. Using the matrix identity: eq:Mexp eq:Mexp and using the fact that, upon iteration of eq. ( 38), the inverse of (X Y) term can be written as: eq:InvDev eq:InvDev the matrix elements of HAA have form: eq:Htilde eq:Htilde where eq. ( 39) has been used together with: to make explicit the different expansion terms.This partitioning scheme projects the problem to the NOGFs formed from the metal-based multielectron wavefunctions belonging to the (2J + 1)-degenerate ground spin-orbit multiplet, defining subspace A, whereas the interaction with the metal-based and ligandbased excitations belonging to subspace B can be accounted for in a perturbative manner, as per eq.( 40).
Here we focus our attention on terms up to first-order in the perturbative expansion eq. ( 40) (i.e. the very first term on the right-hand side of eq. ( 40)).The energies are then determined by eq. ( 36) with Hkk 0 = H kk 0 , and the variational optimization involves the diagonalization of a matrix of dimension (2J + 1) ⇥ (2J + 1) once the integrals in eqs.( 43)-(45) are made explicit.

Implementation in the code CERES: Cofactors versus inverse matrices formalism
We describe here the implementation strategy pursued within our own ab initio code CERES 20,30 of the non-orthogonal CI problem arising from eqs. ( 36)-(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 l and k indices for simplification.The explicit form for the normalization factor M kk 0 on the Y k basis set is given by: eq:HMO0 eq:HMO0 In H kk 0 , the one-electron terms have form: eq:HMO1 eq:HMO1 and two-electron terms: eq:HMO2 eq:HMO2 Defining the overlap matrix S as: : it is well known in the literature that: eq:HMO0det eq:HMO0det eq:HMO1cof eq:HMO1cof and: eq:HMO2cof eq:HMO2cof where det(S) indicates the determinant of S, cof[S](i j) the firstorder cofactor arising from S when the row i and the column j are removed, and cof[S](i j|kl) the second-order cofactor produced by further deleting the row k and the column l.If S 1 exist, the cofactors can be expressed using the Jacobi ratio theorem: eq:jac1 eq:jac1 eq:jac2 eq:jac2 Although this theorem is a powerful tool that makes the calculation of cofactors simpler, the use of eqs.( 51) and (52) requires a matrix inversion for every couple of determinants in eqs.( 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 [31][32][33][34] , and after him different methods have been suggested to avoid the need of calculating the inverse, allowing for a rapid evaluation of the cofactors [35][36][37] , 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 who emphasized the advantages coming from prop-erly taking this peculiarity into account 38 , followed by Figari and Magnasco, who generalized the Hayes' results 39 .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: 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, that are the inactive set of spin-orbitals for the metal plus the spin-orbitals included in the y R 0 group functions.The C block collects the overlap terms between the N act orbitals that can vary in the different product functions, that in this paper coincide con with active set of spinorbitals for the metal.Note that N act ⌧ N inact .The U and V block includes the overlap elements between the active and inactive sets of spin-orbitals.While U, V, and C depends on a specific product functions pair, A never changes.Then, S 1 can be written by the use of the Woodbury matrix identity as 40 : eq:Sblocks eq:Sblocks where Z = (C VA 1 U), X = A 1 U, and Y = VA 1 .To avoid problems when Z is singular, we can apply eq.(51) back to the Z matrix's elements: or, in matrix form: eq:Z-1 eq:Z-1 where cof[Z] T indicates the matrix collecting all the first-order cofactors of Z transpose.Defining the two matrices: and: and remembering that: ) eqs. ( 51) and ( 52) can be written by use of eqs.( 54) and (56) as: eq:cof1 eq:cof1 eq:cof2 eq:cof2 where K i jkl is defined as (all the possible cases under the condition i < k and j < l in eq. ( 50)): eq:Kijkl eq:Kijkl with the p, q, r, s indices running on the Z block.The relation: has been used to cancel out the 1/ det(Z) factor from K i jkl terms, and the permutational symmetry relations: have been used for simplification.We noted that, defining the two matrices: where I is the identity matrix with rank(I) = rank(C), Z i jkl can be written by an unique expression as: Eqs. ( 60) and ( 61) are not affected anymore by singularities, and they do not require any conditions but the existence of A 1 .Moreover, the need of 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 K i jkl term disappears when rank(C) < 2. Finally, eqs.( 60) and ( 61) reduce to just the T terms when rank(C) = 0, as expected.Some more simplifications can be done when making explicit the spin partitioning of S, and by the use of the Cholesky representation of the electron repulsion integrals, see Supporting Information (SI).

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 CAHF/CASCI-SO level 30,41 .In the CAHF/CASCI-SO strategy, the set of molecular orbitals are generated at 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.
Within the FAIMP approach, the crystal field levels and magnetic proprieties are estimated from the wave function of the metal y M m , optimized at CAHF/CASCI-SO level in presence of the ligands described by FAIMP model as: Within the NOGF theory presented in this paper, the crystal field levels and magnetic proprieties are estimated by variational optimization of the pure electrostatic Hamiltonian on the basis of the product functions in eq. ( 28), which includes (i) the wave functions of the isolated metal y M m optimized at CAHF/CASCI-SO level, (ii) the ground state closed shell single Slater determinants optimized at HF level describing the ligands as: Note that in NOGF3 the functions y M m included in the product functions is the wave function of the isolated metal without any influence from the ligands, whereas y M m used in the self-consistent procedure to adjust y W 0 is discharged after the optimization.
The geometries of the Ln complexes analyzed have been fixed to their experimental X-ray structures taken from the literature [REFS].We indicate the ligands with acac = acetylacetonate, dppz = dipyridophenazine, dpq = dipyridoquinoxaline, phen = 1,10-phenanthroline, hfac = hexafluoroacetylacetonate, glyme = dimethoxyethane, see Figure (1) for the mark and colour code used in the paper.
All atoms were described by the ANO-RCC basis set 42 , with the contraction [8s7p5d3f2g1h] for Ln atoms, [4s3p2d1f] for C, N, O atoms coordinating the metal, [3s2p] for C, N, O atoms not coordinating the metal, [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 43 .In CASCI-SO evaluations, the spinorbit interaction were 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 d = 10 8 44  .
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.(%error eq:err1 eq:err1 eq:err2 eq:err2 eq:err3 eq:err3 where 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: sec:res

Results
In Figures ( 2)-( 5) are represented 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 SI).In FAIMP1 and NOGF1 strategies the metal interacts with the electron densities of the separated molecules included in the environment.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%.This 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 others within the complex structure.However, it is noteworthy that general 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 Tb-pc 2 complex, which is formed by two highly p conjugated pc 2 ligands, where 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.Here note that the electron density of the ligands has been optimized under the influence of the density of the free metal ion, described by FAIMP model to avoid an empty region of space at the metal position, which would lead to un-physical repulsive interactions between the electrons of the environment (see FAIMP2 and NOGF2 graphs).In this case, while the percentage errors affecting the crystal field levels calculated by FAIMP method are still large for all the ions, it is impressive how the NOGF2 strategy, where the wave function of the ligands is mixed together with the wave functions of the isolated metal within the product functions, show results in pretty good agreement with those obtained at CAHF/CASCI-SO level on the whole system.This can prove the importance of including an accurate description for the electron density of the ligands, but also it 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.Moreover, it indicates the overlap between the functions describing the different fragments to play an import role.
FAIMP3 strategy still underestimates the energy gaps, showing that also when the densities of the metal and the ligands are iteratively optimized to each others the pure electrostatic approach by FAIMP model is not good enough for accurate predictions.On the other hand, NOGF3 method, where the density of the ligands is adjusted to the metal density by the iterative procedure, 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 bigger discrepancies are shown by no 3 -paah 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 Er(III) metal.Differently, NOGF3 strategy do not lead to a significant improvement in the result for Ho(III) complexes.
It is noteworthy that Dy-trensal and Er-trensal show very different behaviours.The first system is accurately reproduced by NOGF3 method both in the low and high energy ranges.Differently, NOGF3 for Er-trensal shows quite large discrepancies with the CAHF/CASCI-SO energies, especially for the levels above 400 cm 1  , with a large underestimation of the gaps by NOGF3 approach versus CAHF/CASCI-SO results.
To show how the different methods proposed in this paper perform in simulating the magnetic properties, we reported in Table 1/2 the main values of the g tensor for the ground and first excited pseudo-Kramers/Kramers doubles (KD) for Tb(III)/Dy(III) complexes, estimated by the 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 axis estimeted by NOGF3 and CAHF/CASCI-SO methods are in all cases smaller than 8 degrees, see Tables 3 and 4. Larger discrepances are observed for the first excited KDs, where now the g-tensor calculated at NOGF3 level are 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 acac3-dpq ( 22), acac3-dppz (24 ), tta3phen (39 ), hfac3glyme (21 ) ligands for Tb(III), and tta3bipy ligand for Dy(III).

Conclusions
We presented a novel multiconfigurational non-orthogonal group functions (MC-NOGF) approach to 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) 20,30 .By direct comparison of two different fragments approaches of the crystal field energies of 44 complexes, namely the FAIMP and the here propsed 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.
0 , y S 0 , y T 0 indicate the ground state closed shell single Slater determinants of the separated molecules of the ligands optimized at HF level, y W 0 the ground state closed shell single Slater determinant of the whole group of the ligands optimized at HF level.FAIMP !indicates that the functions on the left-hand side affect the function on the right-hand side by FAIMP model, and FAIMP ! that the functions on both the left-and right-hand sides are adjusted to each other through a self-consistent procedure.
The errors affecting the energy gaps DE APP i = E APP i E APP 0 , 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: Fig. 1 fig:legend Mark and color code for the molecules studied in this work.

Fig. 2
Fig. 2 fig:tb Graphical representation of the crystal field levels within the ground spin-orbit multiplet 7 F 6 of different Tb(III) complexes, estimated at 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 s for the %errors affecting the energy gaps are represented by the normal distribution function.

Fig. 3
Fig. 3 fig:dy Graphical representation of the crystal field levels within the ground spin-orbit multiplet 6 H 15/2 of different Dy(III) complexes, estimated at 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 s for the %errors affecting the energy gaps are represented by the normal distribution function.

Fig. 4
Fig. 4 fig:ho Graphical representation of the crystal field levels within the ground spin-orbit multiplet 5 I 8 of different Ho(III) complexes, estimated at 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 s for the %errors affecting the energy gaps are represented by the normal distribution function.

Fig. 5
Fig. 5 fig:er Graphical representation of the crystal field levels within the ground spin-orbit multiplet 4 I 15/2 of different Er(III) complexes, estimated at 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 s for the %errors affecting the energy gaps are represented by the normal distribution function.

Table 1
tab:gttb g tensor values for the two lowest-laying pseudo-Kramers doublets (KD) within the ground spin-orbit multiplet 7 F 6 of different Tb(III) complexes.

Table 2
tab:gtdy Main g tensor values for the two lowest-laying Kramers doublets (KD) within the ground spin-orbit multiplet 6 H 15/2 of different Dy(III) complexes.

Table 4
tab:angdy Main magnetic axis for the two lowest-laying Kramers doublets (KD) within ground spin-orbit multiplet 6 H 15/2 of different Dy(III) complexes: angle (in degrees) between the axis estimated by the different methods proposed in this paper and the ab initio computation on the whole systems at CAHF/CASCI-SO level.