Crystal Field Splittings in Lanthanide Complexes: Inclusion of

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Crystal Field Splittings in Lanthanide Complexes: Inclusion of Correlation Effects beyond Second Order Perturbation Theory P. P. Hallmen,*,†,‡ G. Rauhut,‡ H. Stoll,‡ A. O. Mitrushchenkov,§ and J. van Slageren*,† †

Institute of Physical Chemistry, University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany Institute of Theoretical Chemistry, University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany § Laboratoire de Modélisation et Simulation Multi Echelle, Université Paris-Est, MSME UMR 8208, 5 bd Descartes, 77454 Marne-la-Vallée, France Downloaded via LA TROBE UNIV on July 5, 2018 at 01:47:02 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: State-averaged complete active space self-consistent field (CASSCF) calculations and a subsequent spin−orbit calculation mixing the CASSCF wave functions (CASSCF/state-interaction with spin−orbit coupling) is the conventional approach used for ab initio calculations of crystal-field splittings and magnetic properties of lanthanide complexes. However, this approach neglects dynamical correlation. Complete active space second-order perturbation theory (CASPT2) can be used to account for dynamical correlation but suffers from the well-known problems of multireference perturbation theory, e.g., intruder state problems. Variational multireference configuration interaction (MRCI) calculations do not show these problems but are usually not feasible due to the large size of real lanthanide complexes. Here, we present a quasi-local projected internally contracted MRCI approach which makes MRCI calculations of lanthanide complexes feasible and allows assessing the influence of dynamical correlation beyond second-order perturbation theory. We apply the method to two well-studied molecules, namely, [Er{N(SiMe3)2}3] and {C(NH2)3}5[Er(CO3)4]·11H2O.

I. INTRODUCTION Single-molecule magnets (SMMs) have been investigated intensively during the last years due to their potential applications in high-density data storage, molecular spintronics, and quantum computers.1−7 These compounds are complexes of transition metal, lanthanide, or actinide ions, surrounded by ligands, and show characteristic features such as slow relaxation and quantum tunneling of the magnetization. Paramagnetic lanthanide ions have been shown to be superior to transition metal ions regarding the slow relaxation of the magnetization and the magnetization blocking temperature. This is due to the large total angular momentum of their ground state and their intrinsic anisotropy, which can generate a high effective energy barrier toward relaxation of the magnetic moment.8−10 Very recently, extraordinarily high blocking temperatures (∼60 K) were reported for dysprosiumbased SMMs.11,12 A detailed understanding of the correlation between the magnetic properties and the electronic and geometric structure is necessary to rationally improve the properties of these materials and thus enable their application in devices. In the case of lanthanide SMMs, the electronic structure is essentially determined by their low-lying crystalfield (CF) levels which can be regarded as spin−orbit multiplets of the total angular momentum J of their ground state split by the crystal field generated by the ligands. The total angular momentum has its origin in the unpaired 4felectrons of the compounds, and the strong spin−orbit © XXXX American Chemical Society

coupling in lanthanides makes their magnetism strongly anisotropic. The direction and the magnitude of the anisotropy depends on the nature and arrangement of the ligands surrounding the trivalent lanthanide ion. In order to achieve a better understanding of the relation between geometric and electronic structure and the magnetic properties, a combined use of experiment and ab initio calculations is necessary. The common approach for ab initio calculations of crystalfield splittings in lanthanide complexes is the state-averaged CASSCF/SI-SO3,13 method (complete active space selfconsistent field and state interaction with spin−orbit coupling), also known as CASSCF/RASSI-SO. It is implemented, for example, in the quantum chemistry software Molcas.14,15 This method is able to provide energy levels and together with the SINGLE_ANISO module14,16 also g-factors, crystal-field parameters, and other relevant quantities and was successfully applied to many SMM systems during the past decade.3,11,16−25 However, the CASSCF/SI-SO method covers only static correlation. In order to increase the accuracy of the CASSCF/SI-SO calculations by including dynamical correlation, the multistate complete active space perturbation theory of the second order (CASPT2)26−29 method can be employed to improve the diagonal elements of the spin−orbit matrix,30 but this method often suffers from problems if the minimal Received: February 19, 2018 Published: June 15, 2018 A

DOI: 10.1021/acs.jctc.8b00184 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation active space is chosen (e.g., intruder state problems, “overshooting effect”)27,30,31 or is not feasible for large active spaces and large molecules. Alternatives to state-averaged CASSCF/ SI-SO have been proposed in the literature in order to speed up the time-consuming optimization of the state-averaged orbitals of CASSCF, namely, configuration-averaged Hartree− Fock (CAHF)32 and local-density fitted CAHF (LDFCAHF).33 These methods are more time-efficient than CASSCF/SI-SO calculations but also do not include dynamical correlation. The multistate MRCI method34−37 is a method to variationally include dynamical correlation. It does not suffer from the standard problems of multistate CASPT2 and can be used to correct both diagonal and nondiagonal elements of the spin−orbit matrix by including dynamical correlation. The main problem of this method is that it is not feasible even for small lanthanide complexes, especially if the calculation of many states/roots is required. Here, we present a solution to this problem in the form of a quasi-local projected internally contracted MRCI approach. Here, we make use of localized orbitals to only include the most important effects of dynamical correlation on the crystalfield splitting. In addition, the required states are calculated sequentially instead of simultaneously, which is computationally much more efficient. We test this approach on two complexes for which extensive experimental data is available, namely, [Er{N(SiMe3)2}3] and {C(NH2)3}5[Er(CO3)4]· 11H2O.

Orbitals occurring in the reference configurations, which are not closed-shell orbitals, are called active orbitals, and orbitals which are not occupied in any reference configuration R are called external or virtual orbitals. These CASSCF wave functions are usually obtained by multiconfigurational self-consistent field39−41 (MCSCF) calculations where the spatial orbitals and the CI coefficients are optimized simultaneously by alternating optimization steps. In the case of lanthanide complexes, the active space is most often chosen to consist of the seven 4f-orbitals. Scalar relativistic effects are included by using a pseudopotential for the lanthanide ion or by employing the (second order) Douglas−Kroll−Hess scalar operator42 in the one-electron part of the Hamiltonian. The resulting CASSCF wave functions are used to build the spin−orbit matrix using a convenient SO-operator such as the no-pair Douglas−Kroll− Hess, the Breit−Pauli, or a pseudopotential SO-operator.13,15,43,44 Diagonalization of the spin−orbit matrix yields the crystal field (CF) energy levels, and the eigenvectors can be used to calculate different properties, e.g., crystal field parameters or g-tensors. As S is not a good quantum number in the case of strong spin−orbit coupling, often all possible states of all spin-manifolds are used to build the spin−orbit matrix. Optimizing the orbitals for different states and different spin-manifolds individually would make the calculations unfeasible not only because of the enormous required computational effort but also because this would lead to nonorthogonal (within one spin-manifold) CASSCF wave functions. Therefore, the CASSCF calculations are performed as state-averaged (SA) calculations (SA-CASSCF),41 where one single set of orbitals is calculated by minimizing the average energy of the states within the spin-manifolds of interest. It was demonstrated32 that these SA-CASSCF calculations for a full averaging over all roots of all spinmanifolds to obtain the averaged orbitals for lanthanide complexes can be replaced by configuration-averaged Hartree− Fock (CAHF) calculations. Very recently, we presented a modification of the CAHF method, namely, local-density fitted CAHF (LDF-CAHF),33 that, with a subsequent complete active space configuration interaction (CASCI) step (LDFCAHF+CASCI), can be used to calculate SA-CASSCFequivalent wave functions in a much more time-efficient way. These multiconfigurational wave functions, constructed by the LDF-CAHF+CASCI method, can then be used as a reference for methods including dynamical correlation, such as CASPT2 or multireference configuration interaction (MRCI).34−36,38,45−47 To our knowledge, so far only CASPT2 has been employed for accounting of dynamical correlation effects in lanthanide complexes and within these calculations only the diagonal elements of the spin−orbit matrix have been replaced by CASPT2 energies.30,48,49 These calculations often resulted in an overestimation (“overshooting effect”) of the crystal field splitting, if the minimal active space (seven 4f-orbitals) was chosen.30 Calculations using larger active spaces are usually not feasible for large lanthanide complexes. The reason that so far only CASPT2 has been employed are computational limitations. We want to go beyond a perturbational treatment of dynamical correlation and therefore aim at variational MRCI calculations. II.B. Internally Contracted MRCI (icMRCI). To achieve accurate results, the uncontracted MRCI method can be employed, where the correlated wave function |ΨMRCI⟩ is

II. THEORY We will briefly outline the conventional approach to ab initio calculations in the molecular nanomagnetism community, CASSCF/SI-SO and the recently published (LDF-)CAHF +CASCI/SI-SO method, in Section II.A, because these methods are used to generate the reference wave functions for subsequent calculations including dynamical correlation. In Section II.B, we will review the multistate internally contracted MRCI (icMRCI) method as a method to variationally include correlation effects beyond CASPT2 and explain the problems arising when it is used for lanthanide complexes. Then, we will summarize the projection operator icMRCI method for excited states as a methodology to circumvent the problems of multistate icMRCI if many states have to be calculated. Finally, we will present our new approach of introducing local restrictions in the wave function ansatz of projected icMRCI to render variational icMRCI calculations of crystal-field splittings in lanthanide complexes feasible. II.A. Fast Generation of the Reference Functions. Conventional calculations of crystal field splittings in lanthanide complexes often employ the CASSCF/SI-SO method3 and thus consider only static correlation. Dynamical correlation is treated at most at the CASPT226,27,29 level of theory. The CASSCF wave functions |Φ⟩ are expanded38 in terms of spin-adapted configuration state functions (CSFs) |ΦRμ⟩ |Φ⟩ =

∑ C Rμ|ΦRμ⟩ Rμ

(1)

where R labels the spatial molecular orbital configuration and μ the spin-eigenfunctions of the desired total spin and CRμ are the configuration interaction (CI) expansion coefficients. The spatial orbitals, which are doubly occupied in all reference configurations R, are called closed-shell or inactive orbitals. B


Article

Journal of Chemical Theory and Computation

have to be performed which further increases the size of the space of ICCs. This is likely the reason why MRCI calculations of crystal field splittings of lanthanide complexes including spin−orbit coupling have not been carried out yet. In this work, we present a strategy to include dynamical correlation at the MRCI level of theory and still calculate sufficient states in order to account for the strong spin−orbit coupling. This is achieved by combining projection operator techniques with the exploitation of locality as explained in the following sections. The reference wave functions for this methodology can be efficiently obtained by our previously described LDFCAHF+CASCI method.33 II.C. Projection Operator Technique for Excited States. A particular problem of conventional multistate icMRCI is that, as already mentioned, if several excited states are computed, different contracted configurations are required for each state of interest. Therefore, the number of doubly external ICCs generated from these reference wave functions increases linearly with the number of states, as the calculation of k states requires the calculation of the first k eigenstates of the Hamiltonian matrix in this compound contracted basis.36 This leads to a scaling behavior as k2 to k3 of the whole calculation.36 This becomes a serious problem for calculating the low-energy electronic structure of lanthanide complexes, as many states have to be calculated to account for the strong spin−orbit coupling. A solution for this unfavorable scaling is found in the combined use of a projection operator technique and an ICC basis for state k which is generated only from the kth reference function. In this approach developed by Werner and Knowles,36 the desired kth state is obtained by calculating the ground state of a modified Hamiltonian matrix H(k) iteratively by a direct CI procedure with the (k − 1) lower eigensolutions shifted away (their eigenvalues are shifted to zero). Therefore, the methodology consists of a linear sequence of calculations, starting with the determination of the ground state of H ≡ H(1), using the solution to calculate H(2) and subsequently finding the lowest root of H(2) and then use the obtained eigenvectors to calculate H(3) and so on, up to the desired state k. Thus, the computational effort scales linearly in the number of desired states k, because of the successive instead of simultaneous calculation of the desired states and the fact that the ICC basis for each state is generated only from its own reference function. The calculations of excited states using this projection operator technique are faster and the memory requirements are smaller than for conventional multistate icMRCI calculations. This is due to the linear scaling in the number of states k and the use of an individual ICC basis for each state. This also makes it easier to parallelize the calculations (since each core needs less RAM). Nevertheless, the large number of closed orbitals to be correlated for most lanthanide complexes significantly slows down the calculation and inhibits the use of icMRCI. We present a remedy for this problem below, by combining projected icMRCI calculations with the exploitation of locality to include only the most important effects on the crystal field splittings. II.D. Exploiting Locality in Projected icMRCI. The exploitation of locality was previously used to determine the reference wave function in the LDF-CAHF+CASCI method.33 There was also a great deal of research in the past decade in the field of local correlation methods including dynamical correlation, e.g., local coupled cluster or local CASPT228,55−65

constructed from a subset of CSFs generated by applying single and double excitations to each reference configuration: |ΨMRCI⟩ =

∑ C Iμ|ΦIμ⟩ + ∑ CaSμ|ΦaSμ⟩ + ∑ Iμ

Sμ , a

Pμ ab Cab |Φ Pμ⟩

Pμ , ab

(2)

The index I refers to the internal space, i.e., to the space of all CSFs which contain only closed-shell and active orbitals. The index S refers to the singles space, i.e., to the space of CSFs which contain one electron in some external orbital and the index P to the pair space consisting of all CSFs with two electrons in external orbitals. The indices a and b refer to external/virtual orbitals. We will omit the spin-label μ from now on for convenience. This ansatz is inserted into the molecular electronic Schrödinger equation which in turn is solved variationally. The problem of uncontracted MRCI is that the number of variational parameters (the CI coefficients) quickly becomes so large as to render the method computationally unfeasible. A solution for this problem is the use of internally contracted configurations (ICCs)34,45,47 |Φem⟩ = Ê em|Φ,⟩ and |Φefmn⟩ = Ê efmn|Φ⟩ which are excited configurations generated by applying spin-free (or spin-summed) excitation operators50−53 Ê em and Ê efmn to a fixed linear combination |Φ⟩ (the reference CASSCF function) of the reference CFSs |ΦR⟩, where the indices m and n refer to internal orbitals and the indices e and f to active or external orbitals. The internal contraction leads to a significant reduction in the number of variational parameters and therefore in the computational effort. Nevertheless, the use of ICCs leads to various difficulties due to nonorthogonality, linear dependencies, and the necessity of calculating high-order reduced density matrices.38 Therefore, the implementation of icMRCI by Werner and Knowles34 in the Molpro54 program suite employs a mixed excitation basis which consists of ICCs and CSFs where only the doubly external CSFs are replaced by doubly external ICCs. The wave function ansatz for this case38 is |ΨicMRCI⟩ =

∑ C R|ΦR ⟩ + ∑ C I|ΦI ⟩ + ∑ CaS|ΦaS⟩ R

1 + 2

I

∑ Cabmn|Φabmn⟩ mnab

Sa

(3)

where |ΦI⟩ are CSFs which are not contained in the reference space. Two problems arise when using icMRCI to calculate crystal field splittings in lanthanide complexes: the required time and random access memory (RAM) limitations. These problems arise for different reasons. For most of the lanthanide complexes, the number of closed-shell orbitals is very large and it is not feasible to correlate all of them, because the number of excited CSFs and ICCs grows rapidly with the number of closed-shell orbitals. Moreover, the calculation of crystal field splittings in these complexes requires the calculation of a large number of states belonging to different spin-manifolds in order to account for the spin−orbit coupling. In conventional multistate icMRCI calculations one typically uses a basis which is the union of all contracted configurations generated from separate reference functions for each state up to and including the highest state of interest. Since the number of states of interest is usually large for lanthanide complexes, excitations out of many reference functions for different states C


Article

Journal of Chemical Theory and Computation or embedding approaches.66 Nevertheless, these methods either do not go beyond the level of perturbation theory (local CASPT2) or are not yet applicable for multireference cases (local coupled cluster), or multistate options have not yet been implemented/developed for the calculation of the large number of required states (local coupled cluster and local CASPT2). A local multistate icMRCI method would be highly desirable but has not been reported yet, to the best of our knowledge. Our approach is to use the icMRCI method in a quasi-local manner for the closed-shell orbital space. Our aim was to reduce the computation time and memory requirements to make icMRCI calculations of many states including dynamical correlation feasible. This allows increasing the accuracy of crystal field splitting calculations of lanthanide complexes compared to LDF-CAHF+CASCI in a nonperturbative way (beyond the usually used CASPT2 approach). Therefore, we combined the projection operator technique for excited states with a physically reasonable subspace of closed orbitals to be correlated, which we identify by localization. It is expected that the most important effects on the crystal field splitting are due to the interaction of the forbitals of the central ion with the spatially closest ligand orbitals, e.g., the orbitals localized on the ligand atoms of the first coordination sphere of the lanthanide ion. Since the conventional canonical orbitals are normally delocalized over the whole molecule, we first perform a localization within the space of closed-shell orbitals, e.g., Foster−Boys67 localization, and group the localized orbitals according to their distance (using the center of charge as measure) to the central ion. In the next step, we identify the space of localized closed-shell orbitals at or spatially closest to (e.g., the lone pairs of the directly coordinating atoms) the central ion, which we will call L in the following. Within this subspace L, we make the chosen closed-shell orbitals pseudocanonical again (for faster convergence in the icMRCI iterations), while leaving the remaining localized closed-shell orbitals unchanged, which are then treated as core (uncorrelated) orbitals in the subsequent icMRCI calculation. The icMRCI wave function then has the following form in our approach |ΨicMRCI − loc⟩ =

∑ C R|ΦR ⟩ + ∑ R

+

∑

by a spatial distance criterion. Moreover, this approach can be systematically improved, as it is possible to increase the size of coordination spheres included in L. In combination with the projection operator technique and the efficiently obtained reference functions by LDF-CAHF+CASCI, this finally makes it feasible to include dynamical correlation in a more accurate and variational way in calculations of crystal field splittings in lanthanide complexes and we will call this procedure quasilocal projected internally contracted MRCI (qlp-icMRCI).

III. COMPUTATIONAL DETAILS In this section, the technical details of the calculations of the systems investigated in this work are presented, i.e., of the free Er(III) ion, the Er(III) surrounded by three singly negative point charges in C2v symmetry, and the lanthanide complexes {C(NH2)3}5[Er(CO3)4]·11H2O and [Er{N(SiMe3)2}3]. The structure of the Er(III) ion surrounded by point charges was obtained by placing the three point charges at a distance of 2.1914 Å (the N−Er distance in [Er{N(SiMe3)2}3]68) around the central ion with one bond angle of 90° and two of 135°. All calculations were carried out with the Molpro54,69 suite of ab initio programs. For the all-electron SA-CASSCF calculations (and all subsequent calculations including dynamical correlation) of the free Er(III) ion and the Er(III) surrounded by point charges, the ANO-RCC70 basis set reduced to triple-ζ ([8s7p4d3f2g1h]) was used for erbium. Scalar-relativistic effects were accounted for by the second order Douglas− Kroll−Hess transformation, and the Breit−Pauli operator was used for the SO-coupling step. The LDF-CAHF calculations (and all subsequent calculations including dynamical correlation) were in all cases performed using the 28 electron pseudopotential/effective core potential (ECP) ECP28MWB43 for erbium, and an appropriate ECP-SO operator43 was used for the SO-coupling step. In case of all LDF-CAHF calculations (and the subsequent correlated calculations) the def2-TZVPP/def2-ATZVPP basis sets and the corresponding def2-TZVPP/def2-ATZVPP-JKFIT auxiliary basis sets for density fitting were employed for erbium. For the molecular calculations, the cc-pVDZ/cc-pVTZ/aug-cc-pVTZ71−73 basis sets were used for the remaining elements. In case of [Er{N(SiMe3)2}3], a minimal basis (MINAO74) was used for hydrogen, as no significant effect of these atoms on the crystal field splitting is expected for this molecule. The active space was always chosen to be CAS(11,7) for Er(III) and the corresponding erbium containing systems. The weighting factors in the all-electron SA-CASSCF calculations were chosen to be the same for each spin-manifold. In case of the CASPT2 calculations, only the diagonal elements of the spin− orbit matrix were replaced by CASPT2 energies, while the rest of the matrix was calculated using SA-CASSCF/LDF-CAHF +CASCI wave functions. The full spin−orbit matrix was calculated using icMRCI wave functions in the case of the icMRCI calculations, if not mentioned otherwise. The icMRCI energies used for the diagonal elements of the SO-matrix in this work were obtained without the use of the Davidson correction, which can be employed to approximately account for the size-consistency error of MRCI. However, the deviation of the energy differences obtained by calculations with and without Davidson correction is about 10−5 to 10−6 hartree for the systems investigated in this work and hence the Davidson correction has a negligible effect on the crystal field splitting. For the calculation of {C(NH2)3}5[Er(CO3)4]·11H2O, only the [Er(CO3)4]5− anion was explicitly treated quantum

m

CmnREn̂ |ΦR ⟩

R n ∈{L}, m ∈{A}

∑ ∑

a

CanREn̂ |ΦR ⟩

R n ∈{L ∪ A}, a

1 + 2

∑ mn ∈{L ∪ A}, ab

mn ab |Φmn⟩ Cab

(4)

where A is the space of active orbitals and where we suppose that the reference space is a CAS. Note that the internal configurations not contained in the reference space (second term on the right-hand side (r.h.s.) of the equation) and the singly external configurations (third term in the r.h.s. of the equation) are uncontracted. This means that all ICCs and all CSFs not contained in the reference space are generated by excitations from a restricted set of orbitals, i.e., the union of L and A ({L∪A}), such that only the active orbitals and the closed-shell orbitals in the vicinity of the central ion are correlated. This approach leads to a drastic reduction in the number of correlated closed-shell orbitals but is still expected to cover the most important effects of dynamical correlation on the crystal field splitting, since the closed-shell orbitals are selected in a physically reasonable way D


Article

Journal of Chemical Theory and Computation Table 1. Computational Details for the Different Systems free Er3+ ion

Er3+ + point charges

Basis set

def2-TZVPP

def2-TZVPP

Auxiliary basis set

def2-TZVPP/JKFIT

def2-TZVPP/JKFIT

Pseudopotenial/all-electron calculation? Pseudopotential for Er Structure

Pseudopotential + allelectron ECP28MWB -

Energy convergence threshold (hartree) Number of electrons

10−9

Pseudopotential + allelectron ECP28MWB See section Computational Details 10−9

65

65

a

[Er(CO3)4]5−

[Er{N(SiMe3)2}3] Er: def2-TZVPP C, N, Si: cc-pVDZ H: MINAO Er: def2-TZVPP/JKFIT C, N, Si, H: cc-pVDZ/ JKFIT Pseudopotential

See Table 2

Pseudopotential

ECP28MWB Crystallographicb

ECP28MWB Crystallographica

10−7

10−7

335

193

JKFIT basis sets corresponding to basis in Table 2

b

From ref 75. From ref 68.

converged crystal field energies of the ground multiplet than to obtain a converged energy of the (split) first excited multiplet. This analysis is important, because it is crucial to restrict the space of calculated states as much as possible, especially in view of the demanding icMRCI calculations of crystal field splittings for larger molecules. Moreover, the influence of dynamical correlation on the crystal field splitting employing the point charge model is investigated in Section IV.B using different levels of theory to assess the qualitative and quantitative differences of multistate CASPT2 and multistate icMRCI as well as the performance of the projected icMRCI method compared to multistate icMRCI. Section IV.C covers the application of the presented methodology to molecular systems to demonstrate its applicability to lanthanide complexes. IV.A. Free Er3+ Ion. IV.A.1. Inclusion of Different Numbers of States/Spin-Manifolds. We performed stateaveraged (all-electron) CASSCF and also LDF-CAHF+CASCI calculations with subsequent SI-SO calculations with different numbers of states included for each spin-manifold and changed the number of states always by at least one full degenerate set of scalar-relativistic states. For Er3+, quartets and doublets are possible. In the SA-CASSCF calculations, the orbitals are optimized for the number of states to be mixed in the SI-SO step while the LDF-CAHF orbitals are equivalent to orbitals obtained from a SA-CASSCF calculation averaging over all possible states of all spin-manifolds within the chosen active space (35 quartets and 112 doublets for Er3+), giving the roots of each spin-manifold a weighting factor corresponding to their relative MS-degeneracy (4 for quartets, 2 for doublets). The results are displayed in Table 3. Note that the spin-state mixture 4 includes all possible states of all possible spinmanifolds within the chosen active space. The results of the CAHF and SA-CASSCF calculations are slightly different because of the different state-averaging, but also due to the fact that the LDF-CAHF calculations employ a pseudopotential for the scalar-relativistic part while the SA-CASSCF calculations are all-electron calculations and both types of calculations use a different SO-operator (pseudopotential vs Breit-Pauli). Table 3 shows a large gap of more than 500 cm−1 between spin-state mixtures 1 and 2, while the difference between mixtures 2 and 4 (the mixture including all possible states) is only about 30 cm−1 for SA-CASSCF and ∼20 cm−1 for LDFCAHF+CASCI. Therefore, it is necessary to include at least 35 quartets and 40 doublets (mixture 2) for the free Er3+. We will

chemically. The effects of the counterions and the surrounding crystal were described by putting a +1-point charge at the crystallographic position75 of the central C of the [C(NH2)3]1+ cation and a −5-point charge at the position of the central Er of all [Er(CO3)4]5− anions, besides the one treated quantummechanically, for all ions within 4 crystallographic shells of unit cells in all directions of the crystal. The specific details for the different systems are shown in Table 1. For the LDF-CAHF +CASCI calculations of {C(NH2)3}5[Er(CO3)4]·11H2O different basis sets were employed. In order to distinguish these basis sets, they are listed in Table 2 and numbered from 1 to 3. Table 2. Different Basis Sets Used for the Calculation of [Er(CO3)4]5−a Er O (close) O (distant) C

Basis 1

Basis 2

Basis 3

def2-ATZVPP aug-cc-pVTZ aug-cc-pVTZ aug-cc-pVTZ

def2-TZVPP cc-pVTZ cc-pVTZ cc-pVTZ

def2-TZVPP cc-pVTZ cc-pVDZ cc-pVDZ

a

Close/distant refers to the distance to the central Er3+ ion.

IV. RESULTS This paper focuses on Er containing systems as exemplary systems for our investigations. First, the effect of the inclusion of different numbers of states of different spin-manifolds on the first spin−orbit splitting (i.e., the energy of the first excited Jmultiplet) is analyzed for the free Er3+ ion on SA-CASSCF/ LDF-CAHF+CASCI level of theory to assess how many roots are needed to obtain a converged reference. Moreover, we will investigate the influence of correcting the nondiagonal elements of the SO-matrix for the free ion by dynamical correlation in Section IV.A.2, to analyze the difference of using correlated icMRCI wave functions to build up the SO-matrix or only correcting the diagonal elements by CASPT2 or icMRCI energies. In Section IV.B, the effect of the number of included states on the crystal field splitting of the ion surrounded by three singly negatively charged point charges is investigated at the SA-CASSCF/LDF-CAHF+CASCI level of theory. It is important to assess the minimum number of states needed to achieve converged energies for the reference wave functions for both the free ion and the ion surrounded by point charges and to identify the differences of both cases, because it is expected that fewer states are required to obtain E


Article

Journal of Chemical Theory and Computation Table 3. Spin-Orbit Energies of First Excited J-Multiplet for Er3+ for Different Spin-State Mixtures Included in the SI-SO Calculation Using SA-CASSCF Orbitals (optimized for the corresponding number of states) or LDF-CAHF Orbitalsa mixture pure spin 3/2 1 2 3 4

no. of S = 3/2 states

no. of S = 1/2 states

SACASSCF

LDF-CAHF +CASCI

35

0

6208.05

5893.45

35 35 35 35

20 40 73 112

6196.26 6700.38 6675.40 6670.10

5893.45 6370.68 6348.67 6348.67

the full SO-matrix instead of the simple replacement of the diagonal elements IV.B. Effect of Dynamical Correlation on the CrystalField Splitting Using a Point Charge Model. We performed SA-CASSCF, LDF-CAHF+CASCI, CASPT2, and icMRCI calculations of the crystal field splitting in an Er3+ surrounded by singly negatively charged point charges in C2v symmetry (we did not use D3h symmetry, because Molpro can exploit only Abelian point groups) as a model system. The effect of the spin-state mixture and the level of theory (CASPT2 vs icMRCI) as well as the performance of the projected icMRCI method compared to multistate icMRCI is discussed. First, the influence of the inclusion of different numbers of states of different spin-manifolds to build up the SO-matrix on the crystal field splitting is analyzed at the SACASSCF and LDF-CAHF+CASCI level of theory. This is done to assess how many states are required for each method to achieve a converged reference for the subsequent calculations including dynamical correlation where it is very important to reduce the number of states as much as possible to make the calculations feasible. In the SA-CASSCF calculations, the orbitals are again optimized for the number of states to be mixed in the SI-SO step while the LDF-CAHF orbitals are equivalent to orbitals obtained from a SA-CASSCF calculation averaging over all possible states of all spinmanifolds within the chosen active space (35 quartets and 112 doublets), giving the quartets a weighting factor of 4 and the doublets a weighting factor of 2 (corresponding to their relative MS-degeneracy). The results are displayed in Tables 5 and 6.

a

SA-CASSCF: All-electron calculations. LDF-CAHF: Pseudopotential calculations. Energies are given in cm−1.

investigate below if this large number of required states for the free ion to achieve converged SO-energies is also necessary to obtain converged results for the crystal-field splitting for the free ion surrounded by point charges or in molecules. IV.A.2. Influence of Dynamical Correlation on the SOSplitting. In the following we will consider the influence of dynamical correlation on the SO-splitting energy, i.e., the energy of the first excited J-multiplet, using all-electron SACASSCF wave functions as reference. This will be done using two methods: Multistate CASPT2 and multistate icMRCI. In case of the CASPT2 calculations, only the diagonal elements of the SO-matrix are replaced by the CASPT2 correlation energies and the reference CASSCF wave functions are used to calculate the nondiagonal elements. In the multistate icMRCI calculations, it is possible to calculate the whole SOmatrix using correlated MRCI wave functions and therefore to investigate the effect of correcting also the nondiagonal elements by dynamical correlation. The energies of the first excited multiplet for CASSCF, CASPT2 and icMRCI correcting only the diagonal elements of the SO-matrix and for a calculation building the full SO-matrix using icMRCI wave functions are shown in Table 4.

Table 5. Crystal-Field Splitting of the Ground J-Multiplet of Er3+ Surrounded by Three Singly Negative Point Charges in C2v Symmetry for Different Spin-State Mixtures Included in the SA-CASSCF/SI-SO Calculationa

Table 4. Effect of Dynamical Correlation on the Energy of the First Excited J-Multiplet for the Free Er3+ Ion Using SACASSCF Wave Functions as Referencea level of theory

Er3+ SO-splitting

SA-CASSCF ms-CASPT2 ms-icMRCI diag ms-icMRCI full

6700.38 6692.44 6725.41 6957.80

a

a

CASSCF state-averaging SI-SO step: 35 quartets/40 doublets. Diag: Only the diagonal elements of the SO-matrix are corrected by dynamical correlation. Full: The whole SO-matrix is calculated using correlated wave functions. ms: multi-state. Energies are given in cm−1.

doublet

13/0

35/20

35/40

35/73

35/112

1 2 3 4 5 6 7 8

0.00 127.11 226.35 306.82 391.56 503.15 640.72 811.17

0.00 121.52 216.52 293.73 376.33 485.01 618.76 784.61

0.00 135.12 241.36 327.00 417.13 535.96 682.27 862.29

0.00 130.28 233.02 316.26 403.07 517.51 658.81 833.46

0.00 122.21 218.61 297.24 379.80 488.67 623.08 789.27

Energies are given in cm−1.

Table 6. Crystal-Field Splitting of the Ground J-Multiplet of Er3+ Surrounded by Three Singly Negative Point Charges in C2v Symmetry for Different Spin-State Mixtures Included in the CASCI/SI-SO Calculation Using CAHF Orbitalsa

It is apparent that correcting the nondiagonal elements by dynamical correlation has an important effect (>200 cm−1) on the SO-energies. The CASPT2 and icMRCI calculations replacing only the diagonal elements shift the SO-splitting energy approximately by 10 and 30 cm−1 respectively, while the calculation of the full SO-matrix using icMRCI wave functions leads to an energy change of more than 250 cm−1 compared to CASSCF. Therefore, the effect of dynamical correlation on the SO-splitting in the free ion is quite large and should be treated in a nonperturbative way and by correcting

a

F

doublet

13/0

35/20

35/40

35/73

35/112

1 2 3 4 5 6 7 8

0.00 118.51 212.76 288.86 367.27 470.72 598.63 756.59

0.00 118.22 212.26 287.97 366.75 470.61 598.73 756.90

0.00 122.20 219.95 298.79 379.70 486.53 618.60 781.10

0.00 122.48 220.46 299.48 380.49 487.45 619.72 782.56

0.00 122.48 220.46 299.47 380.49 487.45 619.72 782.57

Energies are given in cm−1. DOI: 10.1021/acs.jctc.8b00184 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article


configuration space is larger for all states. Multistate icMRCI calculations for spin-state mixtures including more states (than for mixture 13/0) only converged employing symmetry. Those calculations can still be used to analyze the effect of dynamical correlation and to show the equivalence of the multistate and projection method, despite the small differences to the calculations not employing symmetry. The molecular calculations will be all performed without the use of symmetry. We present the results for the spin-state mixture 35/20 in Table 8, where we also compare the icMRCI results using the multistate method with the energies obtained by the projection method to justify the use of the latter.

Tables 5 and 6 clearly show that the crystal-field splitting is less sensitive to the spin-state mixture than the SO-splitting, i.e., the energy of the first excited J-multiplet. For the SACASSCF calculations, the mixture 35/20 and 35/112 (full averaging) yield very similar results (Table 5). In the case of the LDF-CAHF+CASCI calculations, there is a slight gap between mixture 35/20 and mixture 35/40, but the average relative difference between 35/20 and the mixture 35/112 including all states is still only ∼3.5% (and the root-meansquare deviation of both mixtures only ∼16 cm−1). Therefore, we will restrict the demanding icMRCI calculations for molecules to a spin-state mixture of 35 quartets and 20 doublets in order to render them feasible. To investigate the effect of dynamical correlation in the point charge model, we performed multistate icMRCI and CASPT2 calculations. We started with the spin-state mixture 13/0, due to the large required computational effort for the multistate icMRCI calculations, to assess which effect dynamical correlation has on the crystal-field splitting using different levels of theory. The results are displayed in Table 7.

Table 8. Influence of Dynamical Correlation on the Crystal Field Splitting for Er3+ Surrounded by Three Singly Negative Point Charges in C2v Symmetry for Different Levels of Theory Calculating 35 Quartet and 20 Doublet States Using SA-CASSCF Wave Functions as Reference and Mixing Them in the SI-SO Stepa

Table 7. Influence of Dynamical Correlation on the Crystal Field Splitting for Er3+ Surrounded by Three Singly Negative Point Charges in C2v Symmetry for Different Levels of Theory Calculating 13 Quartet States Using SACASSCF Wave Functions as Reference and Mixing Them in the SI-SO Stepa doublet

13/0 CASSCF

13/0 msCASPT2 no sym

13/0 msicMRCI no sym

13/0 msicMRCI sym

1 2 3 4 5 6 7 8

0.00 127.11 226.35 306.82 391.56 503.15 640.72 811.17

0.00 137.12 240.57 325.04 417.81 540.47 691.98 882.66

0.00 128.63 226.94 306.98 394.31 509.16 650.18 825.38

0.00 131.98 229.55 309.32 397.92 512.67 654.03 825.48

doublet

35/20 CASSCF

35/20 msCASPT2 no sym

35/20 msicMRCI sym

35/20 picMRCI sym

1 2 3 4 5 6 7 8

0.00 121.52 216.52 293.73 376.33 485.01 618.76 784.61

0.00 136.01 238.47 322.28 415.61 538.93 691.14 882.98

0.00 134.59 230.22 307.40 396.95 511.51 650.74 815.23

0.00 129.48 227.43 305.69 393.84 508.27 646.24 816.65

a

ms: multi-state. p: projection-method. sym/no sym: Symmetry was employed/not employed in the calculations. Energies are given in cm−1.

The CASPT2 method yields very similar results for both spin-state mixtures and overestimates the crystal field splitting (sometimes known as “overshooting effect”30) for both spinstate mixtures compared to MRCI. This is a well-known deficiency of the CASPT2 method and can, for example, be caused by intruder state problems.27,31 Increasing the size of the active space by another set of seven metal-based f-orbitals in the SA-CASSCF calculations can reduce this effect, as some of the dynamical correlation is already included at the CASSCF level.30 On the other hand, increasing the size of the active space often renders subsequent calculations including dynamical correlation unfeasible for larger systems. Table 8 shows that the stretching of the crystal field splitting pattern to larger energies by icMRCI is even increased for the 35/20 mixture (compared to the 13/0 mixture). This is consistent with the empirical use of scaling-factors23 to match the experimental splitting and the energies obtained by the CASSCF/SI-SO approach which systematically neglects dynamical correlation. The projected icMRCI method yields results which differ at most by ∼5 cm−1 from multistate icMRCI which justifies its use for the calculation of crystal-field splittings. The differences stem from the use of different ICC bases: While multistate icMRCI employs a compound ICC basis generated from internally contracted double excitations applied to all reference functions for all states, projected icMRCI calculations employ an ICC basis for state k which is generated only from the kth reference function. This is an essential reason for the reduction of memory and computation time in the projected icMRCI method.

a

ms: multi-state. sym/no sym: Symmetry was employed/not employed in the calculations. Energies are given in cm−1.

Table 7 shows that dynamical correlation included by multistate CASPT2 and multistate icMRCI for the 13/0 mixture leads to a stretching of the crystal field splitting pattern to larger energies. The stretching is significantly smaller for icMRCI than for CASPT2. The CASPT2 method overestimates the crystal field splitting compared to icMRCI, which is much closer to CASSCF than to CASPT2. Moreover, we performed the multistate icMRCI calculations employing symmetry. Table 7 shows that there are small differences between the calculations using and not using symmetry. This is due to the fact that the configuration spaces for calculations with and without symmetry are different. This effect occurs because the wave function ansatz is constructed by applying excitations (see eq 4) to the reference wave functions. If an excitation operator is applied to a reference wave function belonging to a certain irreducible representation, the resulting configuration or ICC is discarded if it does not belong to the same irreducible representation as the reference wave function, but the discarded configurations/ICCs can in principle contribute to the wave function of a state belonging to a different irreducible representation. In calculations without symmetry all configurations/ICCs are kept and therefore the G


Article


Table 9. Crystal-Field Splitting Pattern of [Er{N(SiMe3)2}3] for LDF-CAHF+CASCI/SI-SO, SA-CASSCF/SI-SO, and LDFCAHF+CASCI/Quasi-Local Projected icMRCI (qlp-icMRCI)/SI-SO Level of Theory and Their Comparison with Experimental Values Obtained by Luminescence and the CF Fit to Magnetic Data and the Computation Time for the Full Calculation Including the SI-SO Stepa doublet

LDF-CAHF

CASSCF

qlp-icMRCI

exptl. CF energyb

CF fit energyb

1 2 3 4 5 6 7 8 RMSD expt. full calculation time

0.00 107.53 191.11 262.11 320.74 424.35 480.97 517.94 15.99 1.98 hc

0.00 111.98 199.12 272.83 333.26 439.86 497.96 535.95 15.03 49.41 hc

0.00 129.89 212.72 278.21 338.42 464.70 530.26 572.95 21.16 396.97 hc

0 110 190 245 327 455 0.00 -

0 111 173 228 280 450 495 528 23.72 -

a

Seven closed (spatially closest to the central ion) + seven active orbitals were correlated in the qlp-icMRCI calculation. 35 quartet and 20 doublet states were calculated in the ab initio calculations and mixed in the SI-SO step. RMSD: Root mean square deviation. Energies are given in cm−1. b From ref 76. cCalculations were performed in a parallelized fashion on 8 CPU cores (AMD Opteron 6140, 2.6 GHz).

IV.C. Application to Molecules. IV.C.1. Application to [Er{N(SiMe3)2}3]. The first molecular system chosen was [Er{N(SiMe3)2}3]33,68,76 because its low lying energy level structure is well-known from experiment,76 and it was already used as a benchmark system for the LDF-CAHF+CASCI methodology33 which is used to generate the reference wave functions for the methodology described in this work. It was shown33 that on LDF-CAHF+CASCI/SI-SO level of theory energetically converged results are obtained if one employs a def2-TZVPP77 basis for erbium and a cc-pVDZ71−73 basis for the remaining elements and that is it sufficient to calculate 35 quartet and 20 doublet states and mix them in the subsequent SI-SO step. The RMSD of the LDF-CAHF+CASCI/SI-SO calculation including only static correlation to the experimental energies obtained by luminescence76 was shown33 to be smaller than the RMSD of the CF fit76 to these experimental energies and only 16 cm−1. We employed the quasi-local projected icMRCI procedure and localized the LDF-CAHF orbitals and identified 7 closed-shell orbitals to be correlated, i.e., the 5s and 5p orbitals at the erbium site and the 3 lone pairs of the directly coordinating nitrogen atoms pointing toward the erbium ion, calculated 35 quartets and 20 doublets, and mixed them in the SI-SO step. The corresponding crystal field energy levels of the LDF-CAHF+CASCI/SI-SO, SACASSCF/SI-SO, and LDF-CAHF+CASCI/qlp-icMRCI/SISO calculation and the full calculation time including the SISO step, as well as the experimental energies obtained by luminescence and the CF fit,75 are displayed in Table 9. We tried to perform multistate CASPT2 calculations for this system employing the same basis, the same orbitals to be correlated, and the same spin-state mixture, but these calculations failed to converge (we performed calculations without level shift and with a level shift of 10 hartree31). Inclusion of dynamical correlation by the qlp-icMRCI procedure leads to a stretching of the crystal field splitting pattern to larger energies as previously observed in the point charge model. For doublets 2, 3, and 4 this leads to a larger deviation from the experimental CF energy compared to LDFCAHF+CASCI. For doublet 5 the deviation stays more or less the same and doublet 6 is closer to the experimental energy. We attribute the increased deviation of the lower doublets from the experimental energies (and therefore the resulting increased root-mean-square deviation/RMSD) to the effect of

error compensation for LDF-CAHF+CASCI (and SACASSCF), that coincidentally leads to the exceptional agreement between experimental and ab initio values even though the method completely neglects dynamical correlation. The reason why the qlp-icMRCI method yields slightly worse results compared to the experiment is that dynamical correlation is not yet probably fully captured. This is due to the several restrictions in the wave function ansatz, e.g., the restriction to have at most double excitations and limitations in the basis set size. In the full CI limit (if also relativistic effects, crystal environment effects, etc. are treated properly and very large basis sets are used), this deviation should vanish. The increased deviation from the experiment is therefore a measure for the incompleteness of the correlation space and basis set. The variational inclusion of dynamical correlation by qlpicMRCI leads to a stretching of the energy level pattern with values that are in reasonable agreement with the experimental values obtained by luminescence and no large overshooting occurs as it was observed for CASPT2 calculations of lanthanide complexes using small active spaces (seven 4forbitals).30 We will investigate a second example where LDFCAHF+CASCI or SA-CASSCF calculations fail to predict the energy levels properly, namely, {C(NH2)3}5[Er(CO3)4]· 11H2O, to assess if the qlp-icMRCI methodology can improve the results in cases where the neglect of dynamical correlation leads to large deviations from experimental data. IV.C.2. Application to {C(NH2)3}5[Er(CO3)4]·11H2O. We applied the described methodology to the single ion magnet (SIM) {C(NH2)3}5[Er(CO3)4]·11H2O.75 We chose this system because of the highly negatively charged anion (−5) that increases the influence of dynamical correlation on the crystal field energy levels and because of the availability of extensive experimental data of the crystal field splitting.75 The high charge and large influence of dynamical correlation lead to a large deviation of ab initio calculations considering only static correlation from the experimental energies and makes the system interesting to assess the performance of the methodology described in this work to improve the results. We first performed LDF-CAHF+CASCI/SI-SO calculations to assess the influence of the basis set size and the inclusion of different numbers of states of different spin-manifolds to determine the smallest energetically converged LDF-CAHF+CASCI reference for our quasi-local projected icMRCI methodology H


Article


calculation was performed for 35 quartet and 20 doublet states. The results, the full calculation time including the SI-SO step, and the comparison with the LDF-CAHF+CASCI/SI-SO level of theory as well as the experimental energies obtained by far-infrared spectroscopy (FIR)75 and the crystal-field fit to magnetic data75 are shown in Table 11. We moreover performed SA-CASSCF calculations and also corresponding multistate CASPT2 calculations (using the LDF-CAHF orbitals) for this system employing the same basis, the same orbitals to be correlated and the same spin-state mixture to assess the differences in accuracy and calculation time. The inclusion of dynamical correlation by CASPT2 or qlpicMRCI leads again to a stretching of the crystal field splitting pattern compared to LDF-CAHF+CASCI and SA-CASSCF. CASPT2 overestimates the total splitting and qlp-icMRCI slightly underestimates it compared to the CF fit. The first excited energy level obtained by qlp-icMRCI perfectly coincides with the FIR measurement and the deviation of the second excited state from the FIR measurement is also reduced to less than 9 cm−1. The root-mean-square deviation (RMSD) and average relative error of the energy levels obtained by multistate CASPT2 calculation is worse than for qlp-icMRCI for these lowest doublets, for which FIR data are available, and lies in the range of the RMSD and average relative error of the SA-CASSCF and LDF-CAHF calculations. The energy value of doublet 4 obtained by FIR spectroscopy deviates more from the MRCI and CASPT2 calculations including dynamical correlation than from the LDF-CAHF +CASCI or SA-CASSCF calculation including only static correlation, but it is not certain if this is actually a real peak in the FIR spectrum or some background noise. It was included in the publication by Rechkemmer et al.75 because a quite close energy value resulted from the CF fit to the magnetic data. The RMSD to the CF fit is reduced by ∼40% if dynamical correlation is included by the qlp-icMRCI methodology described in this work and the average relative error of the

regarding the basis set size and the number of included states. The results are shown in Table 10. Table 10. Crystal-Field Splitting of {C(NH2)3}5[Er(CO3)4]· 11H2O for Different Basis Sets and Spin-State Mixtures Included in the CASCI/SI-SO Calculation for LDF-CAHF +CASCI/SI-SO Level of Theorya doublet 1 2 3 4 5 6 7 8

35/112 Basis 1 35/112 Basis 2 35/112 Basis 3 35/20 Basis 3 0.00 42.49 67.61 151.42 205.41 233.36 283.07 366.04

0.00 43.72 66.98 149.35 206.95 234.38 285.02 365.82

0.00 44.04 68.49 153.78 208.15 236.16 285.80 370.10

0.00 45.07 66.68 149.93 209.55 236.93 288.39 371.34

Basis set numbers: See Table 2. Energies are given in cm−1.

a

Since the energy differences between basis 1 and mixture 35/112 and basis 3 and mixture 35/20 are at most approximately 5 cm−1, it is sufficient to use the latter basis and number of states to obtain an energetically converged LDF-CAHF+CASCI reference for reasonable subsequent icMRCI calculations, for which it is crucial to restrict both the basis and the number of states to be calculated and mixed in the SI-SO step. Nevertheless, it should be noted that the convergence of icMRCI with respect to the basis set size is slower than the convergence of LDF-CAHF/CASCI. As a next step, we employed the previously described procedure and localized the LDF-CAHF orbitals obtained using basis 3 and identified 12 closed-shell orbitals to be correlated, i.e., the 5s and 5p orbitals at the erbium site and the 8 lone pairs of the directly coordinating oxygens pointing toward the erbium ion. After making the chosen closed-shell orbitals pseudocanonical within their subspace, a subsequent projected icMRCI

Table 11. Crystal-Field Splitting Pattern of {C(NH2)3}5[Er(CO3)4]·11H2O for LDF-CAHF+CASCI/SI-SO, SA-CASSCF/SISO, LDF-CAHF+CASCI/CASPT2/SI-SO, and LDF-CAHF+CASCI/Quasi-Local Projected icMRCI (qlp-icMRCI)/SI-SO Level of Theory and Their Comparison with Experimental Values from FIR Spectroscopy and the CF Fit to Magnetic Data and the Computation Time for the Full Calculation Including the SI-SO stepa doublet

LDF-CAHF

CASSCF

CASPT2

qlp-icMRCI

FIRb

CF fitb

1 2 3 4 5 6 7 8 RMSD FIR data to ab initio relative error FIR data to ab initio RMSD CF fit to ab initio relative error CF fit to ab initio full calculation time

0.00 45.07 66.68 149.93 209.55 236.93 288.39 371.34 13.19 20.7% 80.23 30.1% 0.6 hd

0.00 46.42 68.17 153.27 214.20 241.20 294.43 380.27 11.87 17.6% 76.15 28.7% 5.16 hd

0.00 66.77 93.29 213.23 313.04 354.05 436.01 549.20 12.34 16.0% 53.87 17.0% 727.6 he

0.00 51.32 75.39 169.72 249.02 280.97 341.61 435.15 6.11 6.4% 48.33 18.7% 520.8 hd

0 52 84 (105)c -

0 44 91 112 280 325 437 462 -

a

Twelve closed (spatially closest to the central ion) + seven active orbitals obtained by LDF-CAHF were correlated in the multi-state CASPT2 and qlp-icMRCI calculations; 35 quartet and 20 doublet states were calculated in the ab initio calculations and mixed in the SI-SO step, and basis 3 (see Table 2) was employed. RMSD: Root mean square deviation. Energies are given in cm−1. bFrom ref 75. cValue in brackets because it was not clear if it is a real peak; it was included in the publication because the CF fit resulted in an energy level which is close to the uncertain band (doublet 4 at 112 cm−1 in the CF fit). Therefore, this value is excluded in the determination of the RMSD and relative error. dCalculations were performed in a parallelized fashion on 8 CPU cores (AMD Opteron 6140, 2.6 GHz). eCalculations were performed in a parallelized fashion on 8 CPU cores (Intel Xeon CPU X5690, 3.47 GHz). I


Article

Journal of Chemical Theory and Computation CF fit to the ab initio energies is also reduced from over 30% to less than 19%. The CASPT2 calculation leads to similar results regarding the RMSD and relative error to the CF fit. It is expected that this RMSD/average relative error can be further decreased mainly by increasing the basis set size (due to the slow convergence of the dynamical correlation energy with respect to the size of the basis) but also by including more (groups of) closed-shell orbitals to be correlated (for comparison: correlating only 4 closed-shell orbitals at the erbium site resulted in a RMSD of ∼55 cm−1 and an average relative error of ∼21%). Moreover, the results can be improved by calculating more states to be mixed in the SI-SO step, using a more sophisticated point charge lattice to describe the Madelung potential of the crystal environment and by increasing the quantum mechanically treated part of the system, e.g., by also treating the cation quantum mechanically, and finally by also including higher excitations than double excitations in the wave function ansatz. Furthermore, no overshooting effect, as often observed for CASPT2, occurs, where the (total) CF splitting is often significantly overestimated. The time required for the full calculation including the SI-SO step is significantly less than for multistate CASPT2, even though the CASPT2 calculation was performed on a faster computer. The qlp-icMRCI method proves to be robust and is able to systematically improve the results in a case, where conventional SA-CASSCF or LDF-CAHF+CASCI calculations fail to reproduce the experiment and where multistate CASPT2 calculations do not converge/are not feasible.

known problems of multireference perturbation theory (as, e.g., intruder state problems or “overshooting effect”). Moreover, it turned out to be more time-efficient than multistate CASPT2 for the systems (and the chosen active space, basis set, etc.) investigated in this work. As any other (MR)CI method, the method is systematically improvable by increasing the basis set size and the number of (groups of) correlated closed-shell orbitals while keeping the minimal size of the active space (seven 4f-orbitals). This is especially interesting for systems with large correlation effects, where conventional SA-CASSCF and LDF-CAHF+CASCI calculations, that neglect dynamical correlation, strongly deviate from the experiment. The reference wave functions for the qlpicMRCI method can be time-efficiently constructed by our recently presented LDF-CAHF+CASCI method. In order to further increase the applicability to large systems using large basis sets and large numbers of correlated orbitals, the local restrictions can in principle be extended from the closed space to the virtual space, i.e., restricting the excitations from the closed/active orbitals to localized virtual orbitals which are spatially close to the orbitals from which the electrons are excited. However, a fully local icMRCI method is a nontrivial task to implement but would be highly desirable.

■

AUTHOR INFORMATION

Corresponding Authors

*(J.v.S.) E-mail: [email protected]. *(P.P.H.) E-mail: [email protected]. ORCID

V. CONCLUSIONS We have analyzed the influence of dynamical correlation beyond the second order perturbation theory on the energy levels of free lanthanide ions and lanthanide ions surrounded by point charges as model systems for single-molecule magnets. We showed that correcting the nondiagonal elements of the spin−orbit matrix by dynamical correlation for the free Er3+ has a large influence (over 250 cm−1) on the energy levels. Moreover, we performed SA-CASSCF, LDF+CAHF+CASCI, multistate CASPT2, and variational multistate icMRCI calculations with subsequent spin−orbit coupling for an Er3+ ion surrounded by three point charges. The results showed that the crystal-field splitting of the ground multiplet is significantly less sensitive to the spin-state mixture than the energy levels of the free ion on SA-CASSCF/LDF-CAHF +CASCI level of theory, i.e., fewer states are required to obtain converged results. Dynamical correlation included by CASPT2 and icMRCI stretched the crystal-field splitting pattern, and replacing the diagonal elements of the spin−orbit matrix by CASPT2 energies overestimates the splitting compared to the calculation of the full spin−orbit matrix using icMRCI wave functions. We showed that projected icMRCI calculations can reliably reproduce the crystal field energies obtained by the multistate icMRCI method. Furthermore, we have presented a methodology for the inclusion of dynamical correlation effects beyond the second order perturbation theory in the calculation of crystal-field splittings of molecular lanthanide complexes consisting of the combination of techniques exploiting locality for the closed-shell/active orbital space with projected internally contracted MRCI (qlp-icMRCI). This method makes the variational inclusion of dynamical correlation feasible for molecules, even for large numbers of states to be mixed in the SI-SO step, and does not suffer from the well-

P. P. Hallmen: 0000-0002-2874-308X G. Rauhut: 0000-0003-0623-3254 Funding

This work was funded by the Vector Foundation and the Landesgraduiertenförderung Baden-Württemberg (scholarship to P.P.H.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

We kindly thank Professor H.-J. Werner (University of Stuttgart) for fruitful discussions.

(1) Pedersen, K. S.; Ariciu, A.-M.; McAdams, S.; Weihe, H.; Bendix, J.; Tuna, F.; Piligkos, S. Toward Molecular 4f Single-Ion Magnet Qubits. J. Am. Chem. Soc. 2016, 138, 5801−5804. (2) McAdams, S. G.; Ariciu, A.-M.; Kostopoulos, A. K.; Walsh, J. P.S.; Tuna, F. Molecular single-ion magnets based on lanthanides and actinides: Design considerations and new advances in the context of quantum technologies. Coord. Chem. Rev. 2017, 346, 216−239. (3) Lanthanides and actinides in molecular magnetism; Layfield, R., Murugesu, M., Eds.; Wiley-VCH: Weinheim, 2015. (4) Urdampilleta, M.; Klyatskaya, S.; Cleuziou, J.-P.; Ruben, M.; Wernsdorfer, W. Supramolecular spin valves. Nat. Mater. 2011, 10, 502−506. (5) Katoh, K.; Isshiki, H.; Komeda, T.; Yamashita, M. Molecular spintronics based on single-molecule magnets composed of multipledecker phthalocyaninato terbium(III) complex. Chem. - Asian J. 2012, 7, 1154−1169. (6) Vincent, R.; Klyatskaya, S.; Ruben, M.; Wernsdorfer, W.; Balestro, F. Electronic read-out of a single nuclear spin using a molecular spin transistor. Nature 2012, 488, 357−360. J


Article

Journal of Chemical Theory and Computation (7) Thiele, S.; Balestro, F.; Ballou, R.; Klyatskaya, S.; Ruben, M.; Wernsdorfer, W. Electrically driven nuclear spin resonance in singlemolecule magnets. Science 2014, 344, 1135−1138. (8) Ishikawa, N.; Sugita, M.; Ishikawa, T.; Koshihara, S.-Y.; Kaizu, Y. Lanthanide double-decker complexes functioning as magnets at the single-molecular level. J. Am. Chem. Soc. 2003, 125, 8694−8695. (9) Woodruff, D. N.; Winpenny, R. E. P.; Layfield, R. A. Lanthanide single-molecule magnets. Chem. Rev. 2013, 113, 5110−5148. (10) Sessoli, R.; Powell, A. K. Strategies towards single molecule magnets based on lanthanide ions. Coord. Chem. Rev. 2009, 253, 2328−2341. (11) Guo, F.-S.; Day, B. M.; Chen, Y.-C.; Tong, M.-L.; Mansikkamäki, A.; Layfield, R. A. A Dysprosium Metallocene Single-Molecule Magnet Functioning at the Axial Limit. Angew. Chem., Int. Ed. 2017, 56, 11445−11449. (12) Goodwin, C. A. P.; Ortu, F.; Reta, D.; Chilton, N. F.; Mills, D. P. Molecular magnetic hysteresis at 60 K in dysprosocenium. Nature 2017, 548, 439−442. (13) Berning, A.; Schweizer, M.; Werner, H.-J.; Knowles, P. J.; Palmieri, P. Spin-orbit matrix elements for internally contracted multireference configuration interaction wavefunctions. Mol. Phys. 2000, 98, 1823−1833. (14) Aquilante, F.; Autschbach, J.; Carlson, R. K.; Chibotaru, L. F.; Delcey, M. G.; De Vico, L.; Fdez. Galván, I.; Ferré, N.; Frutos, L. M.; Gagliardi, L.; Garavelli, M.; Giussani, A.; Hoyer, C. E.; Li Manni, G.; Lischka, H.; Ma, D.; Malmqvist, P. Å.; Müller, T.; Nenov, A.; Olivucci, M.; Pedersen, T. B.; Peng, D.; Plasser, F.; Pritchard, B.; Reiher, M.; Rivalta, I.; Schapiro, I.; Segarra-Martí, J.; Stenrup, M.; Truhlar, D. G.; Ungur, L.; Valentini, A.; Vancoillie, S.; Veryazov, V.; Vysotskiy, V. P.; Weingart, O.; Zapata, F.; Lindh, R. Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table. J. Comput. Chem. 2016, 37, 506−541. (15) Malmqvist, P. -Å.; Roos, B. O.; Schimmelpfennig, B. The restricted active space (RAS) state interaction approach with spinorbit coupling. Chem. Phys. Lett. 2002, 357, 230−240. (16) Chibotaru, L. F.; Ungur, L. Ab initio calculation of anisotropic magnetic properties of complexes. I. Unique definition of pseudospin Hamiltonians and their derivation. J. Chem. Phys. 2012, 137, 064112. (17) Chibotaru, L. F.; Ungur, L.; Soncini, A. The origin of nonmagnetic Kramers doublets in the ground state of dysprosium triangles: evidence for a toroidal magnetic moment. Angew. Chem., Int. Ed. 2008, 47, 4126−4129. (18) Bernot, K.; Luzon, J.; Bogani, L.; Etienne, M.; Sangregorio, C.; Shanmugam, M.; Caneschi, A.; Sessoli, R.; Gatteschi, D. Magnetic anisotropy of dysprosium(III) in a low-symmetry environment: a theoretical and experimental investigation. J. Am. Chem. Soc. 2009, 131, 5573−5579. (19) Cucinotta, G.; Perfetti, M.; Luzon, J.; Etienne, M.; Car, P.-E.; Caneschi, A.; Calvez, G.; Bernot, K.; Sessoli, R. Magnetic anisotropy in a dysprosium/DOTA single-molecule magnet: beyond simple magneto-structural correlations. Angew. Chem., Int. Ed. 2012, 51, 1606−1610. (20) Boulon, M.-E.; Cucinotta, G.; Liu, S.-S.; Jiang, S.-D.; Ungur, L.; Chibotaru, L. F.; Gao, S.; Sessoli, R. Angular-resolved magnetometry beyond triclinic crystals: out-of-equilibrium studies of Cp*ErCOT single-molecule magnet. Chem. - Eur. J. 2013, 19, 13726−13731. (21) Boulon, M.-E.; Cucinotta, G.; Luzon, J.; Degl’Innocenti, C.; Perfetti, M.; Bernot, K.; Calvez, G.; Caneschi, A.; Sessoli, R. Magnetic anisotropy and spin-parity effect along the series of lanthanide complexes with DOTA. Angew. Chem., Int. Ed. 2013, 52, 350−354. (22) Jung, J.; Cador, O.; Bernot, K.; Pointillart, F.; Luzon, J.; Le Guennic, B. Influence of the supramolecular architecture on the magnetic properties of a Dy(III) single-molecule magnet: an ab initio investigation. Beilstein J. Nanotechnol. 2014, 5, 2267−2274. (23) Marx, R.; Moro, F.; Dörfel, M.; Ungur, L.; Waters, M.; Jiang, S. D.; Orlita, M.; Taylor, J.; Frey, W.; Chibotaru, L. F.; van Slageren, J. Chem. Sci. 2014, 5, 3287. (24) Ding, Y.-S.; Chilton, N. F.; Winpenny, R. E. P.; Zheng, Y.-Z. On Approaching the Limit of Molecular Magnetic Anisotropy: A

Near-Perfect Pentagonal Bipyramidal Dysprosium(III) Single-Molecule Magnet. Angew. Chem., Int. Ed. 2016, 55, 16071−16074. (25) Vonci, M.; Giansiracusa, M. J.; Gable, R. W.; van den Heuvel, W.; Latham, K.; Moubaraki, B.; Murray, K. S.; Yu, D.; Mole, R. A.; Soncini, A.; Boskovic, C. Ab initio calculations as a quantitative tool in the inelastic neutron scattering study of a single-molecule magnet analogue. Chem. Commun. (Cambridge, U. K.) 2016, 52, 2091−2094. (26) Finley, J.; Malmqvist, P.-Å.; Roos, B. O.; Serrano-Andrés, L. The multi-state CASPT2 method. Chem. Phys. Lett. 1998, 288, 299− 306. (27) Ghigo, G.; Roos, B. O.; Malmqvist, P.-Å. A modified definition of the zeroth-order Hamiltonian in multiconfigurational perturbation theory (CASPT2). Chem. Phys. Lett. 2004, 396, 142−149. (28) Menezes, F.; Kats, D.; Werner, H.-J. Local complete active space second-order perturbation theory using pair natural orbitals (PNO-CASPT2). J. Chem. Phys. 2016, 145, 124115. (29) Shiozaki, T.; Werner, H.-J. Communication: Second-order multireference perturbation theory with explicit correlation: CASPT2F12. J. Chem. Phys. 2010, 133, 141103. (30) Ungur, L.; Chibotaru, L. F. Ab Initio Crystal Field for Lanthanides. Chem. - Eur. J. 2017, 23, 3708−3718. (31) Roos, B. O.; Andersson, K. Multiconfigurational perturbation theory with level shift  the Cr2 potential revisited. Chem. Phys. Lett. 1995, 245, 215−223. (32) van den Heuvel, W.; Calvello, S.; Soncini, A. Configurationaveraged 4f orbitals in ab initio calculations of low-lying crystal field levels in lanthanide(iii) complexes. Phys. Chem. Chem. Phys. 2016, 18, 15807−15814. (33) Hallmen, P. P.; Köppl, C.; Rauhut, G.; Stoll, H.; van Slageren, J. Fast and reliable ab initio calculation of crystal field splittings in lanthanide complexes. J. Chem. Phys. 2017, 147, 164101. (34) Werner, H.-J.; Knowles, P. J. An Efficient Internally Contracted Multiconfiguration Reference CI Method. J. Chem. Phys. 1988, 89, 5803−5814. (35) Werner, H.-J.; Knowles, P. J. A comparison of variational and non-variational internally contracted multiconfiguration-reference configuration interaction calculations. Theoret. Chim. Acta 1991, 78, 175−187. (36) Knowles, P. J.; Werner, H.-J. Internally contracted multiconfiguration-reference configuration interaction calculations for excited states. Theoret. Chim. Acta 1992, 84, 95−103. (37) Mitrushchenkov, A.; Werner, H.-J. Calculation of transition moments between internally contracted MRCI wave functions with non-orthogonal orbitals. Mol. Phys. 2007, 105, 1239−1249. (38) Shamasundar, K. R.; Knizia, G.; Werner, H.-J. A new internally contracted multi-reference configuration interaction method. J. Chem. Phys. 2011, 135, 054101. (39) Knowles, P. J.; Werner, H.-J. An Efficient Second Order MCSCF Method for Long Configuration Expansions. Chem. Phys. Lett. 1985, 115, 259−267. (40) Werner, H.-J.; Knowles, P. J. A Second Order MCSCF Method with Optimum Convergence. J. Chem. Phys. 1985, 82, 5053. (41) Werner, H.-J.; Meyer, W. A quadratically convergent MCSCF method for the simultaneous optimization of several states. J. Chem. Phys. 1981, 74, 5794−5801. (42) Wolf, A.; Reiher, M.; Hess, B. A. The generalized Douglas− Kroll transformation. J. Chem. Phys. 2002, 117, 9215−9226. (43) Dolg, M.; Stoll, H.; Preuss, H. Energy-adjusted pseudopotentials for the rare earth elements. J. Chem. Phys. 1989, 90, 1730−1734. (44) Heß, B. A.; Marian, C. M.; Wahlgren, U.; Gropen, O. A meanfield spin-orbit method applicable to correlated wavefunctions. Chem. Phys. Lett. 1996, 251, 365−371. (45) Werner, H.-J.; Reinsch, E.-A. The self-consistent electron pairs method for multiconfiguration reference state functions. J. Chem. Phys. 1982, 76, 3144−3156. (46) Shiozaki, T.; Werner, H.-J. Explicitly correlated multireference configuration interaction with multiple reference functions: Avoided crossings and conical intersections. J. Chem. Phys. 2011, 134, 184104. K


Article

Journal of Chemical Theory and Computation (47) Knowles, P. J.; Werner, H.-J. An Efficient Method for the Evaluation of Coupling Coefficients in Configuration Interaction Calculations. Chem. Phys. Lett. 1988, 145, 514−522. (48) Baldoví, J. J.; Duan, Y.; Morales, R.; Gaita-Ariño, A.; Ruiz, E.; Coronado, E. Rational Design of Lanthanoid Single-Ion Magnets: Predictive Power of the Theoretical Models. Chem. - Eur. J. 2016, 22, 13532−13539. (49) Dolg, M. Computational Methods in Lanthanide and Actinide Chemistry; John Wiley & Sons Ltd: Chichester, U.K., 2015. (50) Kutzelnigg, W. Quantum chemistry in Fock space. I. The universal wave and energy operators. J. Chem. Phys. 1982, 77, 3081− 3097. (51) Kutzelnigg, W. Quantum chemistry in Fock space. III. Particlehole formalism. J. Chem. Phys. 1984, 80, 822−830. (52) Kutzelnigg, W.; Koch, S. Quantum chemistry in Fock space. II. Effective Hamiltonians in Fock space. J. Chem. Phys. 1983, 79, 4315− 4335. (53) Kutzelnigg, W.; Mukherjee, D. Normal order and extended Wick theorem for a multiconfiguration reference wave function. J. Chem. Phys. 1997, 107, 432−449. (54) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Györffy, W.; Kats, D.; Korona, T.; Lindh, R.; et al. MOLPRO, a package of ab initio programs, version 2015.1; 2012. See: http://www.molpro.net. (55) Schwilk, M.; Ma, Q.; Köppl, C.; Werner, H.-J. Scalable Electron Correlation Methods. 3. Efficient and Accurate Parallel Local Coupled Cluster with Pair Natural Orbitals (PNO-LCCSD). J. Chem. Theory Comput. 2017, 13, 3650−3675. (56) Schwilk, M.; Usvyat, D.; Werner, H.-J. Communication: Improved pair approximations in local coupled-cluster methods. J. Chem. Phys. 2015, 142, 121102. (57) Werner, H.-J. Communication: Multipole approximations of distant pair energies in local correlation methods with pair natural orbitals. J. Chem. Phys. 2016, 145, 201101. (58) Krause, C.; Werner, H.-J. Comparison of explicitly correlated local coupled-cluster methods with various choices of virtual orbitals. Phys. Chem. Chem. Phys. 2012, 14, 7591−7604. (59) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (60) Saitow, M.; Becker, U.; Riplinger, C.; Valeev, E. F.; Neese, F. A new near-linear scaling, efficient and accurate, open-shell domainbased local pair natural orbital coupled cluster singles and doubles theory. J. Chem. Phys. 2017, 146, 164105. (61) Adler, T. B.; Werner, H.-J. Local explicitly correlated coupledcluster methods: Efficient removal of the basis set incompleteness and domain errors. J. Chem. Phys. 2009, 130, 241101. (62) Adler, T. B.; Werner, H.-J.; Manby, F. R. Local explicitly correlated second-order perturbation theory for the accurate treatment of large molecules. J. Chem. Phys. 2009, 130, 054106. (63) Mata, R. A.; Werner, H.-J. Calculation of smooth potential energy surfaces using local electron correlation methods. J. Chem. Phys. 2006, 125, 184110. (64) Mata, R. A.; Werner, H.-J. Local correlation methods with a natural localized molecular orbital basis. Mol. Phys. 2007, 105, 2753− 2761. (65) Goll, E.; Leininger, T.; Manby, F. R.; Mitrushchenkov, A.; Werner, H.-J.; Stoll, H. Local and density fitting approximations within the short-range/long-range hybrid scheme: Application to large non-bonded complexes. Phys. Chem. Chem. Phys. 2008, 10, 3353− 3357. (66) Coughtrie, D. J.; Giereth, R.; Kats, D.; Werner, H.-J.; Köhn, A. Embedded Multireference Coupled Cluster Theory. J. Chem. Theory Comput. 2018, 14, 693−709. (67) Foster, J. M.; Boys, S. F. Canonical Configurational Interaction Procedure. Rev. Mod. Phys. 1960, 32, 300−302. (68) Zhang, P.; Zhang, L.; Wang, C.; Xue, S.; Lin, S.-Y.; Tang, J. Equatorially coordinated lanthanide single ion magnets. J. Am. Chem. Soc. 2014, 136, 4484−4487.

(69) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: A general-purpose quantum chemistry program package. WIREs Comput. Mol. Sci. 2012, 2, 242−253. (70) Roos, B. O.; Lindh, R.; Malmqvist, P.-A.; Veryazov, V.; Widmark, P.-O.; Borin, A. C. New relativistic atomic natural orbital basis sets for lanthanide atoms with applications to the Ce diatom and LuF3. J. Phys. Chem. A 2008, 112, 11431−11435. (71) Wilson, A. K.; Woon, D. E.; Peterson, K. A.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. IX. The atoms gallium through krypton. J. Chem. Phys. 1999, 110, 7667− 7676. (72) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358−1371. (73) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties. J. Chem. Phys. 1994, 100, 2975−2988. (74) Knizia, G. Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory and Chemical Concepts. J. Chem. Theory Comput. 2013, 9, 4834−4843. (75) Rechkemmer, Y.; Fischer, J. E.; Marx, R.; Dörfel, M.; Neugebauer, P.; Horvath, S.; Gysler, M.; Brock-Nannestad, T.; Frey, W.; Reid, M. F.; van Slageren, J. Comprehensive Spectroscopic Determination of the Crystal Field Splitting in an Erbium Single-Ion Magnet. J. Am. Chem. Soc. 2015, 137, 13114−13120. (76) Jank, S.; Amberger, H.-D.; Edelstein, N. M. Electronic structures of highly-symmetrical compounds of f elements. Spectrochim. Acta, Part A 1998, 54, 1645−1650. (77) Gulde, R.; Pollak, P.; Weigend, F. Error-Balanced Segmented Contracted Basis Sets of Double-ζ to Quadruple-ζ Valence Quality for the Lanthanides. J. Chem. Theory Comput. 2012, 8, 4062−4068.

L


Crystal Field Splittings in Lanthanide Complexes: Inclusion of

Recommend Documents