Improvement of Ab Initio Ligand Field Theory by Means of Multistate Perturbation Theory.

Over the last few years, ab initio ligand field theory (AILFT) has evolved into an important tool for the extraction of ligand field models from ab initio calculations. The inclusion of dynamic correlation on top of complete active space self-consistent field (CASSCF) reference functions, which is important for accurate results, was so far realized at the level of second-order N-electron valence state perturbation theory (NEVPT2). In this work, we introduce two alternative methods for the inclusion of dynamic correlation into AILFT calculations, the second-order dynamic correlation dressed complete active space method (DCD-CAS(2)) and the Hermitian quasi-degenerate NEVPT2 (HQD-NEVPT2). These methods belong to the class of multistate perturbation theory approaches, which allow for the mixing of CASSCF states under the effect of dynamic correlation (state-mixing). The two new versions of AILFT were tested for a diverse set of transition-metal complexes. It was found that the multistate methods have, compared to NEVPT2, an AILFT fit with smaller root mean square deviations (rmsds) between ab initio and AILFT energies. A comparison of AILFT excitation energies with the experiment shows that for some systems, the agreement gets better at the multistate level because of the smaller rmsds. However, for some systems, the agreement gets worse, which could be attributed to a cancellation of errors at the NEVPT2 level that is partly removed at the multistate level. An investigation of trends in the extracted ligand field parameters shows that at the multistate level, the ligand field splitting Δ gets larger, whereas the Racah parameters B and C get smaller and larger, respectively. An investigation of the reasons for the observed improvement for octahedral CrIII halide complexes shows that the possibility of state-mixing relaxes constraints that are present at the NEVPT2 level and that keep Δ and B from following their individual preferences.

Introduction

Ligand field theory (LFT) is a powerful tool for the rationalization of the properties of transition-metal (TM) complexes. Traditionally, its parameters are obtained by a fit to experimental properties such as excitation energies, thermochemical data, electron paramagnetic resonance spectra, or magnetization curves, among many others. There have been some noticeable early successes of this procedure, for example, the explanation of the trends in the heats of hydration of the first-row TMs using ligand field splittings deduced from absorption spectroscopy.[1−3] In fact, ligand field parameters are invaluable and intuitively appealing guides to a host of chemical and physical trends of TM containing compounds. However, in general—and in particular for low-symmetry situations—fits of the LFT model to experimental data are often severely underdetermined. Consequently, the obtained model parameters may have only limited physical content.

Given the severity of the parameterization problem, it seems logical to turn to quantum chemistry in the hope that it can provide first-principles predictions of ligand field parameters. However, the ligand field parameters, being of a semiempirical nature, do not have a precise theoretical definition that would allow their computation in a straightforward manner from quantum chemistry. Hence, indirect procedures are necessary and have been proposed in the literature.[4−6] These approaches typically suffer either from a lack of generality or a lack of uniqueness in the reconstruction of the ligand field parameters from actual quantum chemical calculations. However, one of the early successful and general procedures is the ligand field density functional theory (LFDFT) developed by Atanasov, Daul, and Rauzy.[7,8]

In recent years, powerful wavefunction-based methods like the complete active space self-consistent field (CASSCF) method[9−11] combined with second-order N-electron valence perturbation theory (NEVPT2)[12−14] became applicable to large molecular systems and have surpassed Kohn–Sham density functional theory (DFT) as the method of choice for the treatment of the magnetic properties of TM complexes. Energies from these methods were used in a number of studies to fit ligand field parameters.[15−17] Around the same time, the idea of an ab initio LFT (AILFT) was first proposed.[15] Here, not the energies of ligand field states, but all matrix elements of the complete ligand field Hamiltonian (the ligand field “full-CI” matrix) are fitted simultaneously. As the latter complete ligand field matrix is linear in all ligand field parameters, the fit is unique and the parameters are well-defined. This is an important advantage over the fitting of only a few energy eigenvalues. Hence, AILFT gives access to the Racah B and C parameters, the complete 5 × 5 ligand field matrix (7 × 7 for f-elements), and the spin–orbit coupling constant ζ. A related approach was recently suggested for the analysis of lanthanide complexes, which also leads to a unique extraction of parameters from the matrix elements of an effective Hamiltonian.[18] The two approaches were compared recently.[19] Up to now, AILFT was parametrized using CASSCF or NEVPT2. The latter includes dynamic correlation, a physical effect that is important for quantitative agreement with experiments. Dynamic correlation arises from the contribution of Slater determinants (SDs) with nonvalence excitations to the wavefunction and energies.

AILFT was used with excellent success in a number of studies, for example, in the analysis of magnetostructural correlations in cobalt complexes,[20] the analysis of the ligand field of the azide ligand,[21] and the analysis of experimental results on a cobalt single ion magnet.[22] Furthermore, AILFT was used to investigate periodic trends in lanthanide[23,24] and actinide[24] ions and complexes. The approach was discussed in two recent review articles.[25,26]

We recently introduced a new multistate multireference perturbation theory (MS-MRPT) method called the second-order dynamic correlation dressed complete active space method (DCD-CAS(2)).[27,28] Already in our initial work on that method,[27] we observed for one specific complex, [CrF6]3–, that this new method might provide a better starting point for the AILFT parametrization than the NEVPT2 method that is predominantly used so far to include dynamic correlation. In the present work, we investigate if this finding is part of a more general trend by comparing the performance of DCD-CAS(2) and another recently introduced MS-MRPT (second-order Hermitian quasi-degenerate N-electron valence state perturbation theory, HQD-NEVPT2)[29] when they are combined with AILFT.

Methodology

Ab Initio LFT

In the following, we give an overview of the basic ideas of AILFT as described in the original article introducing the idea[15] and a recent review paper.[26]

At the heart of LFT is the finding that the ground state and many of the low-lying excited states of TM complexes or lanthanide and actinide complexes can be qualitatively described as linear combinations of SDs that only differ in the occupation of a set of 5 (for TM complexes) or 7 (for lanthanide and actinide complexes) molecular orbitals (MOs) that show similarities with the sets of d and f orbitals of free atoms and ions. For simplicity, one also refers to these MOs as d and f orbitals, although they also have covalent contributions from ligand orbitals. The Hamiltonian in the basis of these SDs is, apart from an irrelevant overall energy shift, a function of a 5 × 5 or 7 × 7 ligand field matrix hLFT and of the two-electron repulsion integrals involving the d or f orbitals. With one additional assumption, namely, that the active MOs transform like pure d or f orbitals among each other under rotations, the latter integrals can be written in terms of three (four) parameters A, B, C (F0, F2, F4, F6) for TM complexes (lanthanide or actinide complexes).[30] This is the most widespread parametrization of LFT. It follows from this discussion that the LFT model can exactly parametrize the complete active space configuration interaction (CASCI) matrix for free atoms or ions, where the assumption of spherical d or f orbitals is exactly fulfilled.

AILFT is based on the observation that the LFT model Hamiltonianis linear in the parameters p, collected here in a vector. The parameters p include the 15 (for a d configuration) or 28 (for an f configuration) elements of the one-electron ligand field Hamiltonian hLFT, and the electron repulsion parameters B and C (for a d configuration) or F2, F4, F6 (for an f configuration). Note that more than one electron or hole is necessary to define electron repulsion parameters, and that C is nonredundant only if more than one multiplicity block is considered. Note that in general there is more than one possible multiplicity for a given d configuration, meaning that the Hamiltonian matrix can have an additional index next to IJ that signifies the multiplicity. In the following, these indices are treated as a compound index. One can then write eq in the form of the matrix–vector equationHere, HLFT(p) is the vectorized form of the LFT Hamiltonian and the matrix A is defined by A = HLFT,. The LFT Hamiltonian is then identified with an ab initio-derived effective Hamiltonian Heff, for example, the CASCI Hamiltonian, describing the same part of the electronic spectrum. The optimal parameters describing the ab initio Hamiltonian are obtained by least-squares fitting to the model Hamiltonian, which for linear models has a unique solution given by[31]where A+ is the Moore–Penrose pseudoinverse of the matrix A. If the parameters are nonredundant and there are more matrix elements than parameters (meaning that the system is not underdetermined), one can express the pseudoinverse as

This allows to relate eq to the equations given in the original AILFT reference[15]

Choices for the Ab Initio Effective Hamiltonian

The most obvious choice for the ab initio effective Hamiltonian is a CASCI Hamiltonian with the corresponding d-like or f-like MOs chosen as active orbitals. One can write the CASCI matrix via a spectral resolution aswhere ECASCI is a diagonal matrix containing the CASCI energies and CCASCI is the CASCI coefficient matrix with respect to the configuration state function (CSF) basis. As we noted above, the only approximation when describing eq via the LFT model is the parametrization of the two-electron integrals. For a real complex, with a symmetry-lowering ligand environment, the active orbitals have lower than full spherical symmetry. For example, the d-orbitals in an octahedral complex split into distinct eg and t2g sets. If covalency is not too large, the CASCI Hamiltonian is usually nevertheless very well described by the LFT model. Apart from extreme cases, the root mean square deviations (rmsd) between ab initio and LFT state energies at the CASCI level are typically not much larger than 0.1 eV.[26]

One problem with this version of AILFT is that dynamic correlation is missing from the CASCI method, which also limits the accuracy of the derived AILFT model. A way to incorporate dynamic correlation in a computationally reasonably efficient way is second-order MRPT, a popular variant of which is the N-electron valence state perturbation theory (NEVPT2).[12−14] In this method, the CASCI-part of the wavefunction is unchanged, whereas the total energies are corrected for dynamic correlation. One can define an effective NEVPT2 effective Hamiltonian in analogy to eq via back-transformation of the NEVPT2 energies with the CASCI coefficientsHere, ENEVPT2 is a diagonal matrix containing the NEVPT2 energies. This choice of effective Hamiltonian often leads to better agreement with the experiment than the CASCI Hamiltonian. However, past experience has shown that the LFT model provides a significantly inferior parametrization for this Hamiltonian than for the CASCI Hamiltonian, as is evident from a larger rmsd between ab initio and LFT energies, which is often increased by an order of magnitude or more.[26] This large discrepancy can partially absorb the gain in accuracy of NEVPT2 over CASCI.

Recently, we proposed two MRPT methods of the perturb-then-diagonalize kind, DCD-CAS(2) and HQD-NEVPT2, which do not fix the CASCI parts of the wavefunctions, but allow them to mix and relax under the effect of dynamic correlation. This potentially provides a more balanced modification of the CASCI matrix under the effect of dynamic correlation, compared to the NEVPT2 approach where only the diagonal energies are modified but the CASCI coefficients stay identical.

In DCD-CAS(2), one defines the effective Hamiltonian asHere, the energies EDCD-CAS(2) contain the first-order bias correction. This is the default option for DCD-CAS(2) and was found to be essential for obtaining reliable excitation energies.[27]

In HQD-NEVPT2, the effective Hamiltonian in the CSF basis must also be reconstructed from energies and CI coefficientsas the effective Hamiltonian is usually formulated in the basis of eigenstates of a prior CASCI calculation. The effective Hamiltonians eqs and 11 have been newly implemented in a development version of the ORCA electronic structure program[32] for the present work.

Equations –11 represent the different choices for ab initio effective Hamiltonians including the effect of dynamic correlation that we investigate in this work. By fitting the LFT model to these effective Hamiltonians, one can then obtain “renormalized” parameters in the sense of Gerloch et al.[33]

Test Set and Computational Details

To test the different approaches, we compiled a diverse set of spectroscopically well-documented TM complexes. It includes homoleptic octahedral 3d3 complexes (CrIII with ligands F–, Cl–, Br–, I–, CN–, NH3), the 4d3 complexes [MoCl6]3– and [TcF6]2–, and the 5d3 complexes IrF6 and [ReX6]2– (with X = F, Cl, Br). In order to have different d occupations present, we also include the series of tetrahedral divalent chloride complexes of metals from the first transition row from Ti to Ni. The 3d3 Cr(acac)3 and 3d4 Mn(acac)3 complexes are included as examples for chelate complexes, and the 3d2 [FeO4]2– and [MnO4]3– complexes as examples for TM complexes with high +VI and +V oxidation states, respectively. Systems with only one electron or one hole in the d manifold were not included, as electron–electron repulsion does not play a role for them and they can always be exactly parametrized by the LFT model.

We performed all calculations with a development version of the ORCA electronic structure program.[32] The DKH2 scalar relativistic Hamiltonian[34,35] was used in all calculations. As there is no basis set optimized for a relativistic Hamiltonian for all elements present in the test set, we used DKH-def2-TZVP[36] for elements lighter than Kr (including 3d TMs), Sapporo-DK-TZP[37] for elements lighter than Xe (including 4d TMs), and SARC-DKH-TZVP[36,38−40] for everything heavier. The reason for this choice is that we want to treat as many elements as possible consistently, and the DKH-def2-TZVP was specifically recontracted for use together with the SARC basis sets.[36] Only for 4d elements, where neither of these two basis sets is defined, was the Sapporo basis set used. Auxiliary basis sets, if needed, were constructed using AutoAux.[41] Geometry optimizations were performed with DFT using the BP86[42,43] exchange–correlation functional with D3BJ[44,45] dispersion correction, grid5 integration grids, and C-PCM[46] (with infinite permittivity) in order to approximately account for environment effects. Unless otherwise mentioned, the multireference calculations were performed on top of reference states and orbitals obtained from a CASSCF(N,5) calculation averaged over all possible roots of all multiplicities. The strongly contracted[13,14] version of (HQD-)NEVPT2 was used throughout.

Results

Validation of the Computational Protocol

TM complexes always exist in some environment (solution, crystal) that influences their properties. An often applied computationally inexpensive approach to account for such environment effects is the use of an implicit solvation model like the conductor-like polarizable continuum model (C-PCM).[46] We start by investigating its effect for a set of octahedral CrIII complexes with different ligands. Table shows the bond lengths obtained in a geometry optimization with or without C-PCM. Table shows the quartet excitation energies obtained from CASSCF and NEVPT2 calculations on top of those geometries, with two different state-averaging protocols.

Table 1

Cr–X Bond Lengths (in Å) for the [CrX6] Series with and without C-PCM, Compared to the Experiment

X	expt[26]	BP86-D3	BP86-D3/C-PCM
F–	1.901 (CrF3)	1.998	1.941
	1.933 (K2NaCrF6)
	1.913 (Cs2NaCrF6)
Cl–	2.347 (CrCl3)	2.446	2.378
Br–	2.524 (CrBr3)	2.609	2.539
I–	...	2.817	2.750
C (CN–)	2.078 (K3Cr(CN)6)	2.063	2.017
N (NH3)	2.074 ([Cr(NH3)6](ClO4)3	2.140	2.081

Table 2

Quartet Excitation Energies (in eV) for the CrX6 Series Obtained with Different Geometriesa

	gas phase structure		C-PCM structure
	CASSCF	NEVPT2	CASSCF	NEVPT2	expt
[CrF6]3–	1.42	1.64	1.64	1.88	2.00b
	2.32	2.58	2.64	2.88	2.90b
	3.93	3.94	4.29	4.34	...
[CrCl6]3–	1.18	1.52	1.36	1.73	1.70c
	1.97	2.41	2.23	2.66	2.38c
	3.50	3.63	3.76	3.96	...
[CrBr6]3–	1.10	1.51	1.25	1.69	1.66d
	1.85	2.40	2.07	2.61	2.16d
	3.37	3.58	3.58	3.86	...
[CrI6]3–	1.09	1.64	1.21	1.80	...
	1.84	2.56	2.02	2.74	...
	3.33	3.72	3.51	3.99	...
[Cr(CN)6]3–	3.01	3.69	3.37	4.12	3.29e
	4.18	4.71	4.56	5.13	4.02e
	6.60	7.64	7.27	8.48	...
[Cr(NH3)6]3+	2.08	2.52	2.41	2.89	2.67f
	3.17	3.58	3.56	3.96	3.53f
	4.98	5.47	5.56	6.17	...

The orbitals were obtained from CASSCF with state averaging over all quartet and doublet roots (with equal weight of 0.5 for both multiplicities), yielding a weight of 0.5/10 = 0.05 for each quartet root and a weight of 0.5/40 = 0.0125 for each doublet root.

K2NaCrF6.[47] This system has an experimental metal-ligand distance that is closest to the C-PCM geometry among the systems in Table .

CrCl3.[47]

CrBr3.[47]

K3Cr(CN)6.[48]

[Cr(NH3)6](ClO4)3 in H2O.[49,50]

It can be observed that the structures optimized with C-PCM have smaller metal–ligand bond lengths, which agree better with the experimental bond lengths in all cases except for the CN– ligand, where even the gas phase optimized bond length is too short. Taken together, this suggests that the inclusion of C-PCM will in general be beneficial for accurate structures. For excitation energies the situation is less clear cut. Using C-PCM geometries, the first excitation energy is usually well represented, whereas the second one is typically overestimated. Using the gas-phase geometries, this situation is reversed: the first excitation energy is usually underestimated, whereas the second excitation energy is usually closer to experiment. There is a general trend that the NEVPT2 method with triple-zeta basis sets overestimates excitation energies for d–d transitions, as exemplified by our recent study of the excitation energies of free ions in the gas phase,[28] where discrepancies because of an incorrect description of the environment can be excluded. This overestimation of excitation energies was also discussed in the recent AILFT review article.[26] This means that the good results for second excitation energies using NEVPT2 and the gas-phase geometries rely on error cancellation between the too long bond distances (leading to smaller excitation energies) and errors in the NEVPT2 method (leading to larger excitation energies). We therefore prefer the results based on the C-PCM-optimized geometries, as they do not rely on such error cancellation. The data suggest that a smaller metal–ligand distance leads to larger excitation energies, which is also intuitively obvious as in this case the ligand field is stronger. This explains that the calculation of excitation energies based on the gas-phase geometry (which is closer to the experimental geometry) agrees better with the experiment for the CN– ligand.

Apart from state averaging over all roots (quartet and doublet), we also investigated results obtained with state averaging over only the quartet roots; see the Supporting Information. The results are very similar in the two state-averaging protocols. This implies that here the exact choice is not too important for the results. One should note however that the quality of the AILFT fit can of course suffer if several multiplicity blocks of the effective Hamiltonian are parametrized simultaneously.

rmsds and Comparison with Experimental Excitation Energies

Before comparing the rmsd between LFT and ab initio energies for the different methods investigated in this work, it is in order to discuss some fundamental limitations of the most simple LFT model used here. The parametrization of the electron–electron repulsion in terms of the three Racah parameters assumes that the partially occupied orbitals are pure spherically symmetric d-orbitals. This is only exactly fulfilled for free atoms and ions. In complexes, this assumption is not exactly fulfilled because of covalency, which renders the electron–electron repulsion anisotropic. In the most general case, the 120 permutationally distinct two-electron integrals between the d-orbitals are all unique and should be considered as independent parameters in the model. For systems with some residual symmetry, the number of independent parameters is of course smaller. For example, in an octahedral complex there are 10 independent parameters which completely capture the anisotropy of the electron–electron repulsion at the CASSCF level.[8,30,51] Neglecting this anisotropy will inevitably lead to disagreement between the model energies and the ab initio energies, which correctly include this anisotropy, and hence a larger rmsd. One should however note that for Werner-type complexes with limited covalency, the assumption of isotropic electron–electron repulsion is often a very good one and then rmsds at the CASSCF level are small. Whereas in this work we focus on the pure LFT model with isotropic electron–electron repulsion, we plan to investigate models that explicitly include anisotropy in future work.

Figure shows the rmsd between the LFT energies and the ab initio energies for all investigated methods and complexes. The exact numbers are given in the Supporting Information. Some general trends can be observed. Even at the CASSCF level, the rmsd varies between 0.01 and 0.32 eV, correlating well with the covalency of the d-like active MOs.

Figure 1

Total rmsds between AILFT and ab initio state energies for different methods and complexes.

This is expected, as overlap with ligand orbitals will reduce the symmetry of the orbitals from purely spherical and therefore make the approximation of the two-electron integrals in terms of the three Racah parameters worse. The NEVPT2 rmsd is always significantly larger than the CASSCF rmsd, and they correlate well: if the CASSCF rmsd is large for a particular compound, then usually also the NEVPT2 one will be large compared to other compounds. The rmsds for the multistate methods HQD-NEVPT2 and DCD-CAS(2) improve significantly on the NEVPT2 results. DCD-CAS(2) performs best for all molecules, but its performance is usually closely matched by HQD-NEVPT2. Only for a few cases (the octahedral 5d complexes, [MnCl4]2– and [FeO4]2–) is the HQD-NEVPT2 fit slightly worse than the NEVPT2 one.

The effect of these rmsds is seen in Table , which shows excitation energies calculated with the ab initio methods and the LFT models derived from them, compared to experimental band maxima. This comparison is obviously inflicted with a number of uncertainties as the band maxima do not correspond to vertical excitation energies and are also influenced by a variety of environment effects that we do not attempt to model here. Nevertheless, as shown below, there is still reasonably good agreement. It can be observed in Table that in cases where NEVPT2 gives results that are in closer agreement with the experiment, like for the CrIII halide complexes, the large rmsd of the corresponding AILFT fit (see the Supporting Information for exact numbers) leads to AILFT(NEVPT2) energies that are often not better than the CASSCF values. HQD-NEVPT2 and DCD-CAS(2) on the other hand give, thanks to the much better fit as reflected by the smaller rmsds, LFT models that have a closer agreement with the experiment than the one fitted to NEVPT2. Also, for many other complexes the smaller rmsd at the HQD-NEVPT2 and DCD-CAS(2) levels leads to LFT models that have a better agreement with the experiment than the models that are derived from state-specific NEVPT2. In some cases, one can observe that the large rmsd for NEVPT2 leads to fortuitous cancellation of errors when the NEVPT2 ab initio energies already deviate strongly from the experimental values. This is, for example, the case for IrF6. Here, NEVPT2 and DCD-CAS(2) both predict a first excitation energy of 1.25 eV, which is larger than the experimental 1.09 eV. Because of the larger rmsd, the AILFT(NEVPT2) value of 1.07 then agrees better than the corresponding AILFT(DCD-CAS(2)) value of 1.15 eV. Also, for the two complexes with metals in high oxidation states, [FeO4]2– and [MnO4]3–, the AILFT values, especially the DCD-CAS(2) ones, benefit from such error cancellation. One should be careful in using such a model that has a large rmsd, as the parameters might have limited physical significance.

Table 3

Excitation Energies (in eV) Compared with the Experiment

	(approx.) point group	term symbol	CASSCF	AILFT(CASSCF)	NEVPT2	AILFT(NEVPT2)	HQD-NEVPT2	AILFT(HQD-NEVPT2)	DCD-CAS(2)	AILFT(DCD-CAS(2)	expt
[CrF6]3–	Oh	4T2g	1.64	1.64	1.88	1.65	1.88	1.84	1.88	1.84	2.00[26]
		4T1g	2.64	2.63	2.88	2.61	2.84	2.80	2.83	2.80	2.90[26]
[CrCl6]3–	Oh	4T2g	1.36	1.30	1.73	1.39	1.73	1.61	1.68	1.56	1.70[26]
		4T1g	2.23	2.13	2.66	2.22	2.60	2.46	2.52	2.41	2.38[26]
[CrBr6]3–	Oh	4T2g	1.25	1.18	1.69	1.30	1.69	1.54	1.64	1.49	1.66[26]
		4T1g	2.07	1.95	2.61	2.08	2.54	2.36	2.45	2.30	2.16[26]
[Cr(CN)6]3–	Oh	4T2g	3.37	3.32	4.12	3.84	4.12	4.00	4.05	3.97	3.29[26]
		4T1g	4.56	4.49	5.13	4.96	5.09	4.98	5.00	4.96	4.02[26]
[Cr(NH3)6]3+	Oh	4T2g	2.41	2.36	2.89	2.54	2.89	2.78	2.85	2.75	2.67[26]
		4T1g	3.56	3.50	3.96	3.70	3.91	3.82	3.85	3.79	3.53[26]
[MoCl6]3–	Oh	4T2g	2.22	2.15	2.45	2.22	2.45	2.33	2.43	2.34	2.37[52]
		4T1g	3.05	2.97	3.10	2.94	3.07	2.95	3.03	2.96	2.96[52]
		2Eg	1.67	1.59	1.36	1.30	1.36	1.26	1.35	1.28	1.20[52]
		2T2g	2.56	2.42	2.18	2.03	2.16	1.98	2.15	2.01	1.83[52]
[TcF6]2–	Oh	4T2g	3.51	3.46	3.76	3.56	3.76	3.63	3.78	3.69	3.52[52]
		4T1g	4.53	4.47	4.52	4.41	4.50	4.37	4.50	4.44	4.27[52]
		2Eg	1.90	1.84	1.58	1.52	1.58	1.48	1.57	1.51	1.38[52]
		2T2g	3.00	2.89	2.57	2.42	2.55	2.37	2.55	2.41	2.21[52]
IrF6	Oh	2Eg	1.46	1.56	1.25	1.07	1.25	1.05	1.25	1.15	1.09[52]
		2T2g	2.59	2.54	2.00	1.75	2.00	1.72	1.99	1.87	1.75[52]
[ReF6]2–	Oh	4T2g	4.03	3.98	4.17	4.01	4.17	4.05	4.21	4.13	4.07[52]
		2Eg	1.79	1.74	1.58	1.48	1.58	1.45	1.58	1.47	1.38[52]
		2T2g	2.86	2.77	2.52	2.38	2.51	2.32	2.51	2.37	2.27[52]
[ReCl6]2–	Oh	2Eg	1.59	1.47	1.40	1.20	1.40	1.16	1.39	1.19	1.18[52]
		2T2g	2.55	2.32	2.24	1.93	2.22	1.86	2.20	1.91	1.84[52]
[ReBr6]2–	Oh	2Eg	1.56	1.41	1.38	1.13	1.38	1.08	1.36	1.11	1.14[52]
		2T2g	2.49	2.22	2.20	1.80	2.18	1.73	2.15	1.79	1.75[52]
[MnCl4]2–	Td	4T1	3.62	3.61	3.16	3.09	3.14	3.00	3.14	3.02	2.63[52]
		4T2	3.71	3.72	3.32	3.18	3.31	3.15	3.31	3.17	2.78[52]
		4A1	3.74	3.74	3.37	3.19	3.37	3.18	3.36	3.20	2.88[52]
		4E	3.75	3.74	3.38	3.19	3.38	3.18	3.37	3.20	2.88[52]
[CoCl4]2–	Td	4T1	0.55	0.58	0.75	0.58	0.75	0.77	0.74	0.77	0.64[52]
		4T1	2.57	2.60	2.33	2.23	2.34	2.38	2.33	2.38	1.82[52]
		2E	2.43	2.43	2.28	2.16	2.27	2.17	2.26	2.17	1.98[52]
[NiCl4]2–	Td	3T2	0.34	0.38	0.54	0.44	0.54	0.56	0.55	0.57	0.52[53]
		3A2	0.61	0.67	0.95	0.77	0.95	1.00	0.96	1.01	0.91[53]
		3T1	2.61	2.64	2.41	2.43	2.42	2.51	2.41	2.51	1.82[53]
		1T2	2.09	2.11	1.90	1.90	1.89	1.92	1.88	1.93	1.50[53]
		1E	2.23	2.26	2.11	2.07	2.11	2.15	2.10	2.16	1.50[53]
		1T2	3.17	3.19	3.01	2.88	3.01	2.94	3.00	2.95	2.29[53]
Cr(acac)3	D3 (Oh)	4A1	2.04	2.03	2.38	2.09	2.38	2.31	2.36	2.30	2.19[54]
		4E	2.07	2.07	2.43	2.12	2.43	2.37	2.41	2.36	2.29[54]
		4E	3.14	3.13	3.44	3.18	3.39	3.34	3.36	3.34	2.90[54]
Mn(acac)3	D4h (Oh)	5A1g	1.14	1.12	1.29	1.06	1.29	1.23	1.30	1.24	1.18[55]
		5B2g	2.13	2.07	2.59	2.10	2.59	2.44	2.58	2.42	2.22[55]
		5Eg	2.34	2.27	2.81	2.28	2.81	2.65	2.81	2.64	2.67[55]
[FeO4]2–	Td	3T2	1.91	2.11	3.81	2.13	3.81	2.31	3.05	1.70	1.60[56]
		3T1	2.59	2.93	4.48	3.09	4.23	2.83	3.46	2.30	2.29[56]
		1E	1.48	1.42	1.38	0.66	1.35	0.48	1.35	0.59	0.77[56]
		1A1	2.72	2.51	2.55	1.36	2.49	0.95	2.45	1.14	1.13[56]
[MnO4]3–	Td	3T2	1.71	1.69	2.47	1.83	2.47	2.01	2.25	1.81	1.40[57]
		3T1	2.55	2.56	3.30	2.76	3.15	2.68	2.91	2.52	1.80[57]
		1E	1.67	1.56	1.48	1.17	1.47	1.04	1.47	1.09	1.04[57]
		1A1	2.94	2.70	2.66	2.19	2.60	1.89	2.58	1.97	1.64[57]

We also display the correlation between the calculated and experimental excitation energies from Table in Figure for the exact ab initio numbers and in Figure for the numbers obtained from the different AILFT models. The complexes [FeO4]2– and [MnO4]3–, for which LFT was not expected to work well in the first place, were excluded from these figures as they would significantly distort the results.

Figure 2

Correlation between experimental excitation energies and excitation energies calculated with different ab initio methods for the data shown in Table (excluding the two complexes with high oxidation states). The gray dashed line denotes perfect agreement.

Figure 3

Correlation between experimental excitation energies and energies of the AILFT models derived from different ab initio methods for the data shown in Table (excluding the two complexes with high oxidation states). The gray dashed line denotes perfect agreement.

Table gives the slopes of the linear least-squares fits for the relationship between experimental and calculated excitation energies. For the methods including dynamic correlation (NEVPT2, HQD-NEVPT2, DCD-CAS(2)), it can be observed that the slopes of the ab initio numbers are usually larger than the slopes of the AILFT models. For these, the slopes are much closer to the ideal value of 1.0, which would mean that there is no systematic under- or overestimation. This behavior is another manifestation of the abovementioned error cancellation between the overestimation of excitation energies at the level of second-order perturbation theory and the systematic underestimation of the energies through the AILFT fits compared to experimental energies. As mentioned above, this error cancellation is partially removed by the better fit and hence lower rmsds of the AILFT models derived from the multistate theories HQD-NEVPT2 and DCD-CAS(2). Therefore, the AILFT slopes of these two methods are slightly larger than the corresponding AILFT(NEVPT2) slope. Comparing the Pearson correlation coefficients r also shown in Table , one can observe that the value is closer to the ideal value of 1.0 for the methods including dynamic correlation than for CASSCF, whereas the AILFT models have a slightly smaller correlation coefficient, which can be attributed to the fitting errors. One can also see that AILFT based on the multistate methods HQD-NEVPT2 and DCD-CAS(2) has a slightly improved correlation coefficient compared to the state-specific NEVPT2 method. The mean absolute deviations (MADs) shown in Table are around 0.3 eV for the correlated methods, smaller than the deviations at the CASSCF level and certainly in reasonable agreement considering the intrinsic accuracy of second-order perturbation theory[58] as well as the neglect of explicit environmental effects and experimental uncertainties. It can be observed that the multistate methods have a slightly smaller MAD than state-specific NEVPT2, whereas the AILFT models have the smallest MADs, again because of the abovementioned error cancellation.

Table 4

Slopes m of the Linear Regression Line, Pearson Correlation Coefficient r, and MAD between Theoretical and Experimental Excitation Energies (Excluding the Two Complexes with High Oxidation States)a

	CASSCF		NEVPT2		HQD-NEVPT2		DCD-CAS(2)
	AI	AILFT	AI	AILFT	AI	AILFT	AI	AILFT
m	1.11	1.09	1.14	1.05	1.13	1.08	1.13	1.09
r	0.90	0.89	0.98	0.96	0.98	0.97	0.98	0.97
MAD/eV	0.42	0.41	0.31	0.23	0.30	0.21	0.29	0.22

AI denotes the pure ab initio prediction, whereas AILFT denotes the prediction from the extracted LFT model.

Trends in the Extracted Ligand Field Parameters

We now investigate the behavior of the ligand field parameters extracted via AILFT from the different ab initio effective Hamiltonians. For octahedral and tetrahedral complexes, the one-electron ligand field matrix can be parametrized by a single number, the ligand field splitting Δ. This quantity is the difference between the ligand field orbital energies of the t2(g) and e(g) set. Figure shows the ligand field splitting for all complexes in the test set that are at least approximately tetrahedral or octahedral. For the complexes which do not have this symmetry exactly, orbital energies were averaged to derive the t2(g) and e(g) orbital energies.

Figure 4

Ligand field splittings Δ for all complexes in the test set that are approximately octahedral or tetrahedral. All ligand field orbital energies belonging to degenerate sets in the pseudo-symmetry group were averaged.

One can observe that there is a clear trend toward slightly larger ligand field splittings when using multistate perturbation theory methods compared to the state-specific NEVPT2. Only for the extremely covalent complexes [FeO4]2– and [MnO4]3– is there an exception to this finding at the DCD-CAS(2) level. Figures and 6 show the values of the extracted Racah parameters B and C, whereas Figure shows their ratio C/B.

Figure 5

AILFT Racah parameter B derived from different ab initio methods for all complexes in the test set.

Figure 6

AILFT Racah parameter C derived from different ab initio methods for all complexes in the test set.

Figure 7

AILFT Racah parameter ratio C/B derived from different ab initio methods for all complexes in the test set.

It can be seen that in almost all cases the extracted B is smaller, whereas C is larger for the multistate methods than for NEVPT2. For [FeO4]2–, the C parameter is even negative at the NEVPT2 level, which is clearly unphysical, and a result of the bad fit via the LFT model as a consequence of the large covalency and anisotropy in this system. At the level of the two multistate methods, the value of C for this complex is positive, but remains unreasonably small in magnitude. As a consequence of the larger values of C and smaller values of B at the multistate level, also the ratio C/B increases in most cases. This results in a value that is closer to the often cited estimate of C/B ≈ 4 for 3d TM complexes, which can be derived under the assumption that the radial parts of the d orbitals are Slater functions.[26]

Analysis of State-Mixing Effects in [CrX6]3–

In order to better understand the observations reported above, we investigate the effect of state-mixing in the MS-MRPT methods for the example of the chromium(III) halide complexes, in particular [CrF6]3–. All ab initio and AILFT energies of this system calculated at the state-specific NEVPT2 and at the HQD-NEVPT2 levels are depicted in Figure . At first sight, there is no significant difference between the energy levels for the different methods. However, upon closer inspection, one observes that the AILFT(NEVPT2) values differ much more strongly from the NEVPT2 values than the AILFT(HQD-NEVPT2) values from the HQD-NEVPT2 values. This finding is in agreement with the larger rmsd discussed above. When looking at the energy-level diagrams, this behavior is quite surprising as the results of NEVPT2 and HQD-NEVPT2 seem very similar. Hence, state-mixing seems to have a very minor effect on the energies. Also, the general trend that the imperfect AILFT fit leads to an under- rather than overestimation of the excitation energies at the AILFT level compared to the ab initio energies is nicely illustrated by Figure .

Figure 8

All ligand field energies of [CrF6]3– relative to the ground state calculated with NEVPT2 and HQD-NEVPT2 together with the corresponding AILFT fits. The quartet energy levels are shown in red and the doublet energy levels in blue. The vertical axis is cut at 5 eV, i.e., some higher-lying doublet states are calculated but not shown, to make the differences between the different methods better visible.

The ligand field energies calculated with the state-specific NEVPT2 and HQD-NEVPT2 are also shown in Table . One can clearly see that some state-mixing between states of the same symmetry takes place at the HQD-NEVPT2 level, which repels the energy levels compared to NEVPT2. For example 4T1g(1) and 4T1g(2), 2T2g(3) and 2T2g(4), and 2Eg(3) and 2Eg(4) (highlighted with bold font in Table ) are such pairs of levels whose energies are repelled by letting them mix under the influence of dynamic correlation.

Table 5

Energies (in eV, Relative to the Ground State) of All States of [CrF6]3– for Different Methodsa

level	NEVPT2	HQD-NEVPT2	difference
4A2g	0.000	0.000	0.000
4T2g	1.877	1.877	0.000
4T1g(1)	2.879	2.838	–0.041
4T1g(2)	4.345	4.385	0.040
2Eg(1)	2.181	2.173	–0.008
2T1g(1)	2.316	2.310	–0.006
2T2g(1)	3.135	3.114	–0.021
2A1g	3.699	3.699	0.000
2T2g(2)	4.060	4.058	–0.002
2T1g(2)	4.164	4.161	–0.003
2Eg(2)	4.354	4.350	–0.004
2T1g(3)	4.853	4.844	–0.009
2T2g(3)	5.651	5.567	–0.084
2T1g(4)	5.951	5.953	0.002
2A2g	5.988	5.989	0.001
2T2g(4)	6.163	6.250	0.087
2Eg(3)	6.854	6.762	–0.092
2T1g(5)	7.052	7.069	0.017
2Eg(4)	8.962	9.064	0.102
2T2g(5)	9.224	9.242	0.018

Pairs of states in which clearly a repulsion of energies levels because of state-mixing happens are printed in bold.

For simplicity, we now focus on only the quartet roots and try to understand why state-mixing leads to a better parametrization via the LFT model. As shown in the Supporting Information, one component of the 4T2g level can be chosen as the xy → x2 – y2 excited SD, whereas the xy → z2 singly excited SD and the xz, yz → z2, x2 – y2 doubly excited SD correspond to the same component of the two 4T1g CSFs, which will mix to form the final 4T1g energy levels. In order to make sure that these excitations stay pure in the calculations, that is, not mixed with the other excited CSFs of the same energy, we very slightly tetragonally distorted the complex along the z direction from perfect octahedral symmetry, such that the effective point group is D4. More explicitly, we elongated the bonds in the z direction by 10–5 Å. A graphical depiction of the CSFs is shown in Figure .

Figure 9

Quartet CSFs for an octahedral d3 complex. Only one representative for each multidimensional irreducible representation is shown.

The ligand field Hamiltonian for the quartet states of an octahedral d3 system (setting the energy of the lowest state equal to 0) can be written in the basis of the quartet CSFs Φ1 = |xy, xz, yz⟩, Φ2 = |x2 – y2, xz, yz⟩, Φ3 = |z2, xz, yz⟩, and Φ4 = |xy, z2, x2 – y2⟩ as

One can see that in this model, the first excitation energy is exactly identical to H22LFT = Δ. The difference between the diagonal values H22LFT and H33LFT (the energies of the two singly excited SDs) is 12B, which is exactly twice the off-diagonal element of the 4T1g block. The last diagonal element is equal to H44LFT = 2H22LFT + 3B. It is also important to note that in the LFT model eq , the ratio of ligand field splitting Δ and Racah parameter B can be written as a function of the coefficients of the lower-energy 4T1g state, that is,where C1 and C2 are the coefficients of the first and second 4T1g CSF, respectively. The derivation of this equation is shown in the Supporting Information.

The ab initio Hamiltonians for the fluoride, chloride, and bromide complexes are given in Table .

Table 6

Quartet Blocks of the Ab Initio Effective Hamiltonian Matrices for Complexes [CrX6]3– with Different Halide Ligands X–a

X–	HCASCIeff	HNEVPT2eff	HHQD-NEVPT2eff
F–
Cl–
Br–

The matrix elements are given in eV. The DCD-CAS(2) effective Hamiltonians are similar to the HQD-NEVPT2 ones and can be found in the Supporting Information.

When going from CASCI to NEVPT2, one can observe in all three cases that the first excitation energy (corresponding to Δ in the LFT model) increases, whereas the off-diagonal matrix element of the 4T1g block (proportional to B in the LFT model) decreases. At the same time, the eigenfunctions of the CASCI and NEVPT2 effective Hamiltonians are identical by definition. This means that, according to eq , the ratio Δ/B should stay constant if the LFT model provides a perfect fit. These two requirements are obviously in conflict, which explains why it is not possible to fit the NEVPT2 effective Hamiltonian as well as the CASCI Hamiltonian with the LFT model. For the DCD-CAS(2) and HQD-NEVPT2 effective Hamiltonians on the other hand, the diagonal elements of the 4T1g block have a much larger difference, corresponding to a larger Δ/B ratio. This leads to better fits of these effective Hamiltonians via the LFT model, with smaller rmsds.

That the HQD-NEVPT2 effective Hamiltonian can be better approximated by the LFT model than the NEVPT2 effective Hamiltonian can also be seen in a different way. From eq , it follows that, once the parameter Δ is set equal to the matrix element H22LFT, B can be determined from three different matrix elements, that is,

In the LFT model, these three equations must of course lead to exactly the same B value. One can however also apply these equations to the ab initio effective Hamiltonian instead of the LFT Hamiltonian, and then the values for B will in general be different. Their variance is then a measure for how good the effective Hamiltonian can be described with the LFT model. Table gives the values for different ab initio methods and ligands.

Table 7

Racah Parameter B Extracted from Different Matrix Elements of the Effective Hamiltonian after Fixing the Ligand Field Splitting

		CASCI	NEVPT2	HQD-NEVPT2
F–	B33	0.133	0.128	0.108
	B34	0.133	0.117	0.106
	B44	0.133	0.018	0.098
Cl–	B33	0.124	0.122	0.098
	B34	0.126	0.106	0.097
	B44	0.137	–0.005	0.091
Br–	B33	0.121	0.121	0.094
	B34	0.124	0.102	0.094
	B44	0.145	–0.015	0.091

It can be seen that for CASCI, the three values are almost identical for the quite ionic complex [CrF6]3–, whereas their difference grows when going to the more covalent chloride and bromide complexes. The reason is that the t2g and eg orbital sets differ more in the complexes with stronger covalency, which leads to more anisotropic electron–electron repulsion. For NEVPT2, the B33 and B34 values are similar for all three complexes, but B44 is very small or even negative, which is unphysical. This can again be seen as an effect of the too constrained form of the NEVPT2 effective Hamiltonian, which incorporates dynamic correlation only on the level of total energies, but not on the level of the wavefunctions. HQD-NEVPT2 corrects this shortcoming and the three resulting B values are quite similar, as is expected by the better AILFT fit for this level of theory.

When comparing the NEVPT2 and HQD-NEVPT2 effective Hamiltonians for [CrF6]3– in Table , one can also observe that the matrix elements change quite substantially, for example, by 0.24 eV for the diagonal matrix elements of the 4T1g block. At the same time, the effect of the state-mixing on the total energies is relatively small (about 0.04 eV, see Table ). This explains the question raised by Figure of why the fit is so much worse at the NEVPT2 level and shows that one should not judge the importance of state-mixing by only looking at total energies.

The fitted LFT parameters for [CrF6]3– are shown in Table together with the quartet energies. The nephelauxetic ratio β = Bcomplex/Bgaseous also given in Table is determined with respect to calculations on the free CrIII ion. It can be clearly seen that the better fit for DCD-CAS(2) and HQD-NEVPT2 improves the value of Δ, which corresponds to the first excitation energy and has an experimental estimate of 2.00 eV (see Table ). Allen et al., for example, obtained Δ = 1.88 eV and B = 0.092 eV by LFT-fitting to experimental energies, corresponding to a roughly 20% nephelauxetic reduction with β = 0.81.[59] It is interesting that the LFT parameters derived from the multistate methods are closer to those parameters than the prediction of the state-specific NEVPT2 method.

Table 8

Excitation Energies and LFT Parameters (All Quantities Except the Nephelauxetic Ratio in eV) for Octahedral [CrF6]3– with a Bond Length of 1.9408 Å (Averaged)

	CASCI	NEVPT2	HQD-NEVPT2	DCD-CAS(2)
4T2g	1.644	1.877	1.877	1.880
4T1g(1)	2.636	2.878	2.837	2.822
4T1g(2)	4.288	4.342	4.383	4.397
Δ	1.641	1.654	1.839	1.835
B	0.133	0.121	0.109	0.111
β	0.92	0.97	0.88	0.89

Conclusions

In this work, we introduced two new ab initio effective Hamiltonians for use in AILFT analyses. They are defined in terms of the recently introduced multistate perturbation theory methods DCD-CAS(2) and HQD-NEVPT2.

We tested the different AILFT versions on a diverse test set of TM complexes with metals from different periods and groups of the periodic table, different ligands, different coordination environments, and different oxidation states. Compared to NEVPT2, which is so far the standard choice for including the effects of dynamic correlation into AILFT, the multistate methods were shown to yield better fits of the ligand field model to the ab initio effective Hamiltonians. We expect that other multistate methods like XMS-CASPT2[60] or XMCQDPT2[61] would show a similar behavior. As they are not implemented in ORCA, we did not test their performance in the present work. For some systems (e.g., the octahedral CrIII halide series), the better rmsd between the ab initio energies and the energies predicted by the extracted LFT models led to better agreement between the AILFT energies and the experiment. However, in some cases the agreement got worse, for example, for IrF6. This could be explained by a cancellation of errors happening on the level of AILFT(NEVPT2) between an overestimation of excitation energies in the NEVPT2 method and an underestimation because of fitting errors. Because of that, we observed on average only minor improvement of the agreement between AILFT and experimental energies when going from NEVPT2 to the multistate methods. It should however be emphasized that the extracted model parameters do change and reflect more closely the physical picture described by the ab initio methods if the fit is better.

An investigation of the LFT parameters for the whole test set showed that for the multistate methods there is a clear trend for the ligand field splitting Δ of tetrahedral and octahedral complexes to increase compared to NEVPT2. The Racah parameter B usually decreases, whereas the value of C (and therefore also the ratio C/B) increases.

In order to understand the mechanism for the smaller rmsds at the multistate levels, we investigated in detail the case of the quartet states of the [CrF6]3– complex. We found that after inclusion of dynamic correlation, there is a tendency to increase the ligand field splitting Δ, which can be interpreted in terms of increased metal–ligand covalence, and to decrease at the same time the magnitude of the Racah parameter B. At the NEVPT2 level, these parameters cannot change independently from their CASSCF values as the wavefunctions are required to stay the same. At the multistate level on the other hand, the CASSCF wavefunctions can mix under the effect of dynamic correlation, which allows Δ and B to vary independently. We think that a similar mechanism for the reduction of rmsds when going from NEVPT2 to the multistate methods is occurring for the whole test set, that is, that only after the possibility of state-mixing are the parameters Δ and B independent of each other and can follow their individual preferences for change after dynamic correlation. This could explain the general trends observed for the extracted LFT parameters.

The present work focused on the most straightforward LFT parametrization, where the electron–electron repulsion is described in an isotropic way via only three parameters. A more exact model would incorporate the possibility of anisotropic electron–electron repulsion. For example, for an octahedral complex, where the orbitals with variable occupation are of t2g and eg symmetry, the 120 electron–electron repulsion integrals can be written in terms of only 10 parameters because of symmetry.[30,51] Future work will focus on extensions of AILFT in this direction.

Our work also showed that methods that allow for state-mixing can lead to better fits via the LFT model. The accuracy is however limited by some systematic errors of the perturbation theory methods. A natural extension would therefore be the use of more accurate electron correlation treatments that include state-mixing, for example, using the MR-EOM-CC method.[62−64] This can potentially have the same benefits of a better fit between the LFT model and ab initio effective Hamiltonian as the new methods introduced in the present work, while not suffering from the limitations of low-order perturbation theory.