Impact of the KohnâSham Delocalization Error on ... - ACS Publications

Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Article

Impact of the Kohn-Sham Delocalization Error on the 4f Shell Localization and Population in Lanthanide Complexes Thomas Duignan, and Jochen Autschbach J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00238 • Publication Date (Web): 25 May 2016 Downloaded from http://pubs.acs.org on May 31, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Impact of the Kohn-Sham Delocalization Error on the 4f Shell Localization and Population in Lanthanide Complexes Thomas Duignana and Jochen Autschbacha,∗ a Department

of Chemistry University at Buffalo State University of New York Buffalo, NY 14260-3000, USA email: [email protected] Abstract: The extent of ligand to metal donation bonding and mixing of 4f (and 5d) orbitals with ligand orbitals is studied by Kohn-Sham (KS) calculations for LaX3 (X = F, Cl, Br, I), GdX3 , and LuX3 model complexes, CeCl62– , YbCp3 , and selected lanthanide complexes with larger ligands. The KS delocalization error (DE) is quantified via the curvature of the energy for non-integer electron numbers. The extent of donation bonding and 4f-ligand mixing correlates well with the DE. For Lu complexes, the DE also correlates with the extent of mixing of ligand and 4f orbitals in the canonical molecular orbitals (MOs). However, the localized set of MOs and population analyses indicate that the closed 4f shell is localized. Attempts to create situations where mixing of 4f and ligand orbitals occurs due to a degeneracy of fragment orbitals were unsuccessful. For La(III) and in particular for Ce(IV), Hartree-Fock, KS, and coupled-cluster singles and doubles calculations are in agreement in that excess 4f populations arise from ligand donation, along with donation into the 5d shell. Likewise, KS calculations for all systems with incompletely filled 4f shells, even those with ‘optimally tuned’ functionals affording a small DE, produce varying degrees of excess 4f populations which may be only partially attributed to 5d polarization.

1 Introduction The widespread application of complexes of lanthanide (Ln) elements begs a better understanding of their electronic structure.1, 2 In this context, the question of whether the 4f orbitals of Ln participate in any bonding with its ligands, or to which degree this may occur, has been the subject of theoretical as well as experimental research. The spatial extent of 4f orbitals is small, limiting the possibility of overlap with ligand orbitals.3 For paramagnetic species, the spin density is likewise considered to be localized at the 4f center.4–6 However, recent computational and experimental research has indicated that for early lanthanides the 4f shell may act as an acceptor for ligand → metal (L → M) dative (donation) bonding.7–9 The 4f shell radially contracts with increasing atomic number within the lanthanide series. But even for late lanthanides, the assumption of the 4f shell not participating in bonding has been challenged. For instance, a significant covalency of the Yb 4f shell was attributed to YbCp3 , based on experimental magnetic and optical data.10 A mixing of an ionic with a charge-transfer configuration in the wavefunction helped to rationalize the experimental data. This configurational 1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

mixing argument is somewhat reminiscent of the controversy regarding the electronic structure of cerocene and bis(η 8 -pentalene)cerium,11–15 i.e. whether these are closed-shell Ce(IV) systems, whether the ligands donate an electron to form Ce(III) or whether the wavefunction represents a mixture of these configurations. Dolg and Mooßen argued recently, based on data obtained from multi-reference wavefunction theory (WFT) calculations,13, 15 that the electronic structure of these systems is best described by a Ce(III) configuration with a weak Ce – ligand interaction that is not unlike a stretched covalent bond with pronounced multi-reference character. Two main mechanisms are considered to drive metal-ligand covalency, viz. orbital overlap and orbital degeneracy.16 The overlap argument is familiar to chemists and underlies the usual textbook explanations of covalency. The degeneracy argument is based on the interaction of two fragments, such as a metal and a ligand, with their mixing estimated via perturbation theory (PT). In an orbital model framework, the fragment orbital energies then appear as differences in the denominators of the PT expressions for the mixing coefficients, while the interaction matrix elements appear in the numerators. Hence, a near-degeneracy is thought to be capable of inducing a strong mixing of the fragment orbitals and an apparent covalency even if the fragment orbitals are essentially not overlapping / interacting. Kohn-Sham (KS) theory and its generalized versions utilizing – by necessity – approximate functionals remains the computational tool of choice for systems without pronounced multireference character. It offers familiar orbital-based descriptions of the electronic structure and bonding along with an arsenal of analysis tools for dissecting the interactions between fragments in molecules. KS orbital energies depend on the chosen functional, and for hybrid functionals they depend strongly on the fraction of exact exchange (eX). It is then a question whether a calculation with a given functional may accidentally cause degeneracy-based mixing while a calculation with a different functional would not. Regarding overlap-driven mixing, approximate KS calculations – in particular with non-hybrid functionals – generally suffer from a delocalization error (DE).17, 18 This means that the extent of orbital mixing may be severely exaggerated in KS calculations,19 depending on the functional. Hartree-Fock (HF) theory tends to produce an opposite ‘localization error’, which is considered here as a negative delocalization error accompanied by weaker mixing in molecular calculations. We previously noticed a small 4f population in natural bond orbital (NBO) analyses of trigonal planar LaX3 model complexes (X = F, Cl, Br, I) based on non-hybrid KS calculations.9 The 4f populations were used in part to quantify and rationalize calculated trends for the Ln NMR chemical shifts. A systematic trend in the NMR shift dependent on the halide was previously observed experimentally for the solids and qualitatively reproduced by calculations for the same model LaX3 complexes as well as larger cluster models.20 Bogart et al. recently synthesized an 8coordinate cerium(IV) complex and reported a 4f population of 0.8 from hybrid KS calculations 2 ACS Paragon Plus Environment

Page 2 of 29

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with a relativistic pseudopotential for Ce.7 Dolg et al. noted unreasonably large excess 4f populations in KS calculations of Eu(III) complexes, leading to overestimated optimized Eu–ligand distances.21 Clavaguera et al. noted a pronounced mixing of Lu 4f and F 2p orbitals in LuF3 in quasi-relativistic KS calculations with a non-hybrid functional,22 which was attributed to a degeneracy-driven mixing. Ramakrishnan et al. subsequently employed quasi-relativistic KS calculations with a Hubbard U parameter (‘DFT+U’) to remove the degeneracy, which essentially eliminated the orbital mixing in LuF3 .23 The usual interpretation of the Hubbard U parameter for 4f (and 3d) metal shells is that it corrects their unphysically strong delocalization in KS calculations. Therefore, the finding of Reference 23 may also be interpreted in terms of a strong reduction of the DE. Xu et al. studied the LnF3 series by non-hybrid quasi-relativistic KS calculations and concluded that the Ln3+ −F– covalency arises from a dative interaction of F 2p orbitals with Ln 5d orbitals, with a small admixture of Ln 6s. Any Ln 4f contributions were attributed to polarization of the 5d orbitals.24 The notion of degeneracy-driven mixing and whether it should be considered as covalency has generated some debate.24–26 For instance, Ji et al. questioned the concept by pointing out that it does not cause an energetic stabilization of Lu−X bonds.25 Concerning actinide complexes, Kaltsoyannis emphasized that covalency is intuitively associated with orbital overlap and internuclear charge buildup.26 The case for 4f orbital covalency loses traction when discarding orbital mixing without overlap as a measure of covalency and given the small spatial extent of Ln 4f orbitals. However, experimental evidence for both 5d and 4f covalency in octahedral lanthanide hexachlorides was recently obtained by Löble et al. using X-ray absorption spectroscopy (XAS),8 with the observed trends being supported by the results from KS calculations. Questions surrounding 4f covalency in lanthanide complexes remain, and they warrant further study. For instance, would a degeneracy-driven L → M dative bonding lead to a partial electron transfer to the metal and thereby significantly alter the electron density distribution? Or would the DE be able to give rise to a functional-dependent variation of Ln 4f density via overlap-driven L → M dative bonding? For late lanthanides, and for Lu(III) (4 f14 ) in particular, how may the DE or degeneracy impact any mixing between metal 4f and ligand orbitals in KS calculations? Are the 4f populations assigned to complexes of early lanthanides based on KS calculations physically and chemically meaningful, or mainly an artifact of the DE? Are sizable excess 4f populations from L → M dative bonding obtained in KS calculations for which the DE is minimized as best as possible? The main focus of this work is the impact of the DE in relation to lanthanide 4f bonding, via numerical quantification of the DE caused by various approximate XC functionals in correlation with any Ln 4f and 5d contributions in ligand-centered orbitals. The extent of the DE is quantified by the curvature of the energy E(N) between integer N. We also employ non-empirical ‘tuning’ of 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

hybrid functionals with range-separated exchange (RSE) in order to generate molecule-specific functional parametrizations that produce small E(N) curvature. Populations from KS and HF results are compared with those from coupled cluster singles and doubles (CCSD) calculations. Section 2 provides details about the procedure used to quantify the DE in a system as well as other computational aspects. The results are collected and discussed in Section 3. First, the correlation between the DE and the 4f populations is established. NPA results from KS and CC calculations are compared next, followed by a discussion of the impact of the DE on occupied orbitals in comparison with previous works. Population donation into lanthanide 5d versus 4f orbitals from ligand p orbitals is studied in more detail for LaX3 and CeCl62– . Section 3 subsequently considers the possibility of degeneracy-driven orbital mixing in KS calculations of some of the complexes studied and analyzes results from a two-orbital model. Section 4 summarizes our findings.

2 Theoretical and Computational Details Non-empirical tuning of range separated exchange (RSE) functionals has received considerable attention within the computational chemistry community.27 RSE functional tuning has also emerged within the diagnosis of the DE as a prescription to improve calculated molecular groundstate and response properties. This has been shown for properties which typically exhibit a strong functional dependence, such as electric field gradients, excitation energies, optical and fundamental band gaps, optical rotation, hyperfine coupling, and others.19, 28–32 If the Ln−X covalent interactions predicted from KS calculations show a pronounced functional dependence, then it is possible that a tuned RSE functional yields more reliable results than standard functionals. The tuning procedure has been detailed elsewhere19, 27 and is therefore only briefly outlined −1 here. In the range-separation function used in this work,33 the inter-electronic distance r12 is split into a short-range and a long-range component via 1 − [α + β erf(γ r12 )] α + β erf(γ r12 ) 1 = + r12 r12 r12

(1)

γ is the range separation parameter, α determines the fraction of eX in the short-range, and α + β is the fraction of eX at large interelectronic separation. The parameter γ may be varied to minimize the square of the difference between −εH (Ni ), the highest occupied molecular orbital (HOMO) energy of a reference system with Ni electrons, and its vertical ionization potential IP(Ni ) calculated as the difference between the total energies of the Ni−1 system and the Ni system. This non-empirical tuning establishes numerically as close as possible an exact KS condition for any many-electron system, namely that the negative HOMO energy equals the vertical 4 ACS Paragon Plus Environment

Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ionization potential. The minimization may be carried out with different numbers of electrons simultaneously by minimizing the sum of the square differences (the criterion is often referred to as J 2 ). With N0 being the electron number of the target molecule, we optimize the sum of the square differences for i = 0, 1 which has the added benefit that the negative LUMO energy of the N0 -electron system, −εL (N0 ), represents its electron affinity. The tuning of RSE functionals for molecules requires that α + β = 1, i.e. the functional switches to full eX asymptotically. In one-dimensional (1D) tuning, for fixed α = 1 − β the parameter γ minimizing J 2 is found. Originally, the procedure was set up for a functional not affording short-range eX i.e. α = 0, β = 1. Two-dimensional (2D) tuning repeats the procedure for different values of α = 1 − β and selects the parametrization with the least DE, quantified, as mentioned in the Introduction, by the curvature of E(N) calculated between integer Ni .29 The exact electronic energy should afford straight-line segments with slopes of −IP(Ni ) for Ni−1 < N < Ni , i.e. it has no curvature and changes its slope discontinuously at the Ni . A negative curvature of E(N) (concave, like that of a parabola opening downward) reflects too much localization (typical of HF). A large positive curvature (convex, like that of a parabola opening upward) is indicative of unphysical delocalization and particularly pronounced in calculations with non-hybrid functionals. There is currently less experience with the tuning of open-shell systems, and the fact that the HOMO of the α -spin set of orbitals may be above the LUMO of the β -spin set, or vice versa, may cause ambiguities.34 We considered the orbitals and orbital energies used for the tuning procedure in spin-unrestricted calculations to be from the same set of α or β spin orbitals, respectively, as Janak’s theorem is valid for spin-density functionals.35 The E(N) curvatures were obtained in two different ways. In one approach, the energies were calculated explicitly for fractional N varying in steps of 0.2 electrons. The coefficients of the N 2 term in quadratic fits of the data were then used to assess the extent of the DE. Alternatively, employing a cubic interpolation proposed by Johnson et al. E(N) was determined from the energies and frontier orbital energies of the N−1 , N0 , and N+1 systems via36 E(∆N) = (∆E)∆N + [(ε0 − ∆E)(1 − ∆N) + (∆E − ε+1 )∆N]∆N(1 − ∆N)

(2)

Here, ∆N is the deviation in electron number from the N0 reference system, ∆E = EN+1 − EN , ε0 is the LUMO energy of the N0 system and ε+1 is the HOMO energy of the N+1 system. The equation was originally defined for addition of an electron to a system. When applied to the addition of an electron to the N+1 system, a constant of −∆E = −(EN − EN−1 ) was added to retain the total energy difference between N−1 and N+1 systems. The KS calculations for this work were performed with NWChem version 6.3.39, 40 Recent developments have improved the scaling order of CC calculations by discarding cluster ampli5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1: Bond lengths (Å) used for D3h LnX3 model systems. X La a Gd b Lu c F 2.220 2.053 1.968 Cl 2.620 2.489 2.396 Br 2.740 2.640 2.544 I 2.990 2.841 2.779 a References 9, 37. b Reference 38. c Reference 22.

Figure 1: Structures of some experimentally known complexes used in this study. 1 is an octahedral cerium(IV) hexachloride. 8 point charges of +0.25 are positioned at the vertices of the dashed cube. 2 is an 8-coordinate redox active Ce complex. On the right is the famework of 3a-c where M is La, Gd, and Lu, respectively. See Table S4 for details on additional ligands complexed to the metal in 3. Also note that the central metal in 3a-c resides above the plane of the porphyrin, with the deviation from planarity decreasing across the series. 4 is triscyclopentadienyl ytterbium(III).

6 ACS Paragon Plus Environment

Page 6 of 29

Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tudes between orbitals that are well separated spatially.41 The method is called the domain-based local pair natural orbital coupled cluster singles and doubles (DLPNO-CCSD) and calculations were performed using ORCA version 3.0.3.41, 42 The DLPNO method is not currently implemented for open-shell systems, for which the LPNO scheme was used instead. The population analyses for all calculations were carried out with a locally modified version of the Natural Bond Orbital (NBO) 5.043–45 code. Stuttgart small core (28 electron) effective core potentials (ECPs) were used for the lanthanide atoms,46, 47 along with matching segmented contracted valence basis sets.48 One set of molecules is comprised of lanthanide trihalides, LnX3 [Ln = La, Gd, Lu; X = F, Cl, Br, I], with idealized planar D3h geometries and Ln–X distances as given in Table 1. Additional computations were performed for EuX3 and YbX3 . Since the trends for these systems are the same as for the other lanthanide systems they are not discussed in detail in Section 3. The 6311G*49, 50 Pople basis sets were employed for the light halides in the KS calculations. A 46 electron ECP and accompanying Stuttgart RLC basis set was used for iodide.51 The aug-ccpVDZ basis was employed for all non-lanthanide atoms in the CC calculations, along with a cc-pVDZ fitting basis set.52, 53 Changing the auxiliary basis to that accompanying aug-cc-pVDZ had only a minor role in increasing the weight of the reference by around 1 percentage point and was therefore deemed unnecessary. Tests were performed to verify that the correlation-consistent Dunning basis sets and truncation of the coupled cluster expansion are suitable for the purpose of this study. The corresponding population data can be found in Figure S2. Additionally, the results using the Pople basis sets in the KS calculations were validated against the same calculations performed using the Dunning basis sets (Figure S1 in the SI). In the test calculations on LaCl3 , the 4f populations with the Dunning basis sets are too small compared to 6-311G* without the augmentation. The smallest Cl 3p exponent in the Pople basis is 0.109 but in the cc-pVDZ basis it is 0.162, and the more diffuse 3p shell appears to facilitate stronger donation of electron density to the metal. For cc-pVTZ, the smallest Cl 3p exponent is 0.130 and the effect on the 4f population of the LaCl3 test system from augmenting this basis with diffuse functions is much smaller than for the DZ basis. The second set of complexes studied consists of experimentally characterized complexes, referred to in this paper as 1, 2, 3a-c, and 4. See Figure 1 and Table S4 for details regarding the structures. 1 is an octahedral CeCl62– with an experimental bond length of 2.599 Å that was studied by XAS.8 Eight +0.25 point charges were positioned above the normal vectors of the octahedral faces to stabilize the negative charge of 1. Without compensating charges, the system is not stable with respect to the loss of an electron. 2 is an 8-coordinate Ce(IV)-centered redox active compound, noted for its calculated 4f population in the ground state.7 3a-c are expanded porphyrin complexes containing La, Gd, and Lu, respectively.54 4 is a YbCp3 complex with a 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

pseudo trigonal-planar coordination environment of the metal,10 mentioned in the Introduction, to which significant 4f covalency has been attributed. Geometries for 2 and 3a-c were optimized as in previous works at the B3LYP level of theory with 6-31G*55 basis sets.7, 56 The geometry for 4 was taken from crystallographic data and used without further optimization.10 Test calculations were performed to compare the full complexes of 2 and 3a-c including extended alkyl groups R with model structures using truncated substituents. In the case of 2, R groups were replaced by CH3 and for 3a-c R groups were replaced by H. The calculations indicated no significant change in the electronic structure around the central metal or in the properties of interest. Therefore, the truncated structures were used in subsequent calculations. Similarly, test calculations showed insignificant changes in the metal populations upon changing the ligand basis sets from 6-31G* to 6-31G57 basis, and the smaller 6-31G basis was used for complexes 2 and 3a-c. For the smaller system 4 the 6-311G**58 basis sets were used for C and H. A variety of functionals were tested, ranging from HF exact exchange (eX) only (i.e. HartreeFock theory), RSE functionals, global hybrid functionals with varying fractions of eX in the exchange functional, and non-hybrid generalized gradient approximation (GGA) functionals. A fully long-range corrected (LC) RSE hybrid version of the Perdew, Burke, and Ernzerhof (PBE)59 GGA with α = 0.25, β = 0.75, γ = 0.3 dubbed LCPBE0 and the ‘Coulomb attenuating method’ (CAM) version of the B3LYP60–63 functional, CAMB3LYP,33 with α = 0.19, β = 0.46, γ = 0.33 constitute the RSE functionals studied. The tuning procedure determined optimal γ (1D tuning) or γ , α = 1 − β (2D tuning) values for the PBE-based hybrid. The tuned parametrizations are referred to herein as tuned LC (TLC). The numerical values of α and γ that were used for each complex are given in Table S3. The popular global hybrids PBE064 (25% eX) and B3LYP (20% eX) were selected as well, along with the non-hybrid GGA functionals PBE and BP.65, 66 Scalar relativistic effects resulting from high nuclear charges were treated in the calculations indirectly, via the applied ECPs. Spin-orbit effects are neglected for the present study, as it has been in most of the computational studies cited in the Introduction.

3 Results and Discussion 3.1

Delocalization Error

First, we investigate the accuracy of the interpolation function of Equation (2) for the purpose of obtaining a measure of the curvature of E(∆N). Here, ∆N = N −N0 is the electron number relative to N0 of the system of interest. A comparison between calculations of E(∆N) using fractional N explicitly versus interpolated is shown in Figure 2. The agreement is very good. Figures S4 and S5 provide additional data showing that for our samples the interpolation is safe to use. 8 ACS Paragon Plus Environment

Page 8 of 29

Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: Energy of LaCl3 as a function of electron number relative to the energy of the neutral system (∆N = N − N0 = 0). The lines represent the energy determined from the interpolation of Equation (2) and the markers correspond to calculations with fractional electron numbers. The curvature measures of E(∆N) listed next to the plot legend are the second derivative of Equation (2) averaged over ∆N between two consecutive integers, in eV. The numbers listed correspond to the negative and positive ∆N region, respectively. Inset is a magnified view of the slope of the energy around ∆N = 0. Deviations between the energies of the systems calculated with fractional N and the interpolated energies predicted are on the order of 0.1 eV or less for essentially all cases. Accordingly, we proceed with the interpolation. Numerical data indicating the sign and extent of the DE for selected complexes are collected in Figures 3 and 4 in the form of curvature measures for E(∆N). These were determined as the second derivative of Equation (2) averaged over ∆N between two consecutive integers. The trends for the systems not shown are very similar and therefore omitted in order to avoid clutter. Data for the full set of the LnX3 complexes, including EuX3 and YbX3 , can be found in Figures S3 and S6). Complexes which underwent the 2D tuning procedure are 1, 2, LaF3 , LaCl3 , GdCl3 , LuCl3 , and 4. For the remaining systems, TLC indicates a 1D tuned parametrization. A periodic trend across the halides is observed in the ∆N < 0 regime where the DE decreases when the ligand become chemically ‘softer’, i.e. easier to polarize. In the ∆N > 0 regime the trends along the 9 ACS Paragon Plus Environment


HF LCPBE0

PBE0 B3LYP

PBE BP

N0 → N+1

2.5 2.0

3 0.8

1.5

2 0.6

1.0

1 0.5

0.4 0

0.0

0.2 −1

Lu F

Lu I3

3

0.8

G dI 3

3

G dF

La F

3

0.6

La I3

−1.0

0.4

Lu I3

3

Lu F

G dI 3

3

0.2

G dF

La I3

0.0 −2 0.0

3

−0.5

La F

Curvature Coefficient (eV)

TLC CAMB3LYP

N−1 → N0

1.0 4

1.0

Figure 3: Curvature measures of E(∆N) (eV, see caption of Fig. 2) of selected LnX3 systems. Left: ‘Cationic’ range ∆N < 0. Right: ‘Anionic’ range ∆N > 0. The functionals in the bar charts are ordered for each molecule from left to right as HF, LCPBE0, TLC, CAMB3LYP, PBE0, B3LYP, PBE, BP. LCPBE0 TLC

CAMB3LYP PBE0

N−1 → N0

1.0 3.0

B3LYP PBE

BP N0 → N+1

4

2.5

3

0.8

2.0 2

1.5 0.6 1.0

1

0.4 0.5

0

0.0 0.2

4

3c

0.8

3b

1

0.6

3a

−2

0.4

4

3c

0.2

3b

2

0.0

1

0.0 −1.0

2

−1

−0.5

3a

Curvature Coefficient (eV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

1.0

Figure 4: Curvature measures of E(∆N) (eV, see caption of Fig. 2) for larger complexes 1 through 4. Left: ‘Cationic’ range ∆N < 0. Right: ‘Anionic’ range ∆N > 0. There was no convergence of the HF calculations for 1 to 4 for the N−1 and N+1 systems and the corresponding values are omitted in the charts. halide series are less pronounced and not as systematic. Figures 3 and 4 demonstrate that there is a pronounced and systematic dependence of the DE, as measured by E(N) curvature, on the choice of functional. The trends are in line with other numerical studies of the DE, namely that HF and functionals with large fractions of eX produce negative curvatures indicating too localized electronic structures. In contrast, positive 10 ACS Paragon Plus Environment

Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Isosurfaces (±0.03) of selected canonical MOs of LuF3 from calculations with different functionals. Also shown are two natural localized molecular orbitals (NLMOs) from the PBE calculation. Notation for the symmetry species as in Reference 67. See Figure 9 regarding the symmetry classification of different d and f orbitals in the D3h point group. curvature is characteristic of too much delocalization, which is most pronounced with the nonhybrid functionals PBE and BP but still severe with the two global hybrids. CAMB3LYP tends to perform better than the global hybrids but not entirely satisfactory. The LC functionals do well to minimize the DE, in particular upon tuning. The tuning procedure specifically involves the HOMO and LUMO, and one may expect that it mainly reduces the DE as far as it relates to the frontier orbitals. The question is then whether the tuning and the conclusions drawn from the data of Figure 3 has any bearing on the 4f shell of lanthanide complexes, in particular for late lanthanides where the 4f shell is more core-like. Figure 5 graphically demonstrates how changes in the sign and magnitude of the DE correlate with changes to some of the canonical MOs of LuF3 . The PBE set is very similar to that of Clavaguera et al.22 who previously noted the mixing of 4f with ligand orbitals attributed to a degeneracy-driven mechanism. The PBE orbitals calculated for the present study exhibit the same – considerable – extent of mixing. The HF calculation, on the other hand, generates little mixing. The tuned functional produces orbitals between the two extremes, which goes along with the small magnitude of the DE as quantified by E(N) curvature. Compared to standard functionals, TLC reduces the 4f mixing with ligand orbitals in LuF3 but it does not completely suppress it. For comparison, upon application of DFT+U by Ramakrishnan et al. 23 the orbital mixing in 14a′1 / 16a′1 shifted from 42.8% / 54.0% to 99.6% / 1.2% Lu 4f and from 53.7% / 41.3% to 0.2% / 91.0% F 2p based on Mulliken gross atomic populations. 11 ACS Paragon Plus Environment


0.25 HF LCPBE0 TLC CAMB3LYP PBE0

0.20 4f Population

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

0.15

B3LYP PBE BP DLPNO-CCSD

0.10 0.05 0.00

LaF3

LaCl3

LaBr3

LaI3

GdF3

GdCl3

GdBr3

GdI3

Figure 6: Variation of 4f occupations for selected LnX3 from a natural population analysis (NPA) for different KS functionals and CCSD calculations. When a population analysis as discussed in the next section is carried out, the combined 4f natural atomic orbital (NAO) population for all LuX3 systems is 14.00 irrespective of the functional used. Moreover, in the localized MO set produced by this calculation (i.e. the natural localized molecular orbital, or NLMO set), the 4f shell appears as a set of semi-core orbitals without ligand contributions. Two examples from the PBE calculation for LuF3 are shown in Figure 5. The same findings apply to the Lu complex 3c. This means that for these systems the DE correlates with the tendency of the KS calculations to produce occupied canonical orbitals in which 4f mixes with ligand orbitals. However, this orbital mixing can be undone by taking localized linear combinations. In this case, if mixing occurs in the canonical set of orbitals it is no measure of covalency or 4f delocalization no matter if it is being driven by the DE or by a degeneracy. 3.2

Natural Population Analysis (NPA)

As discussed in the previous subsection, despite a varying extent of the DE in the different calculations for the Lu complexes, there is no accompanying dramatic change in the electron density or in the localized MOs. At least for Lu, the 4f shell can be considered as localized. Mixing of 4f with ligand orbitals in the canonical MOs may occur via a degeneracy or because of the DE, but without apparent consequences. The mixing can be undone via orbital localization. The reader is reminded that for a delocalized electronic structure the orbitals, or a subset of them, cannot be localized to the same degree as it is possible in systems with localized bonding. If there were consequences to the 4f mixing with ligand orbitals in LuX3 this should be apparent also after localizing the orbitals. For the complexes of other – in particular early – lanthanides, empty or partially filled metal 4f shells, along with the empty 5d and 6s shells, may potentially act as acceptors of covalent 12 ACS Paragon Plus Environment

Page 13 of 29

1.4 HF LCPBE0 TLC CAMB3LYP PBE0

1.2 4f Population

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.0 0.8 0.6

B3LYP PBE BP DLPNO-CCSD

0.4 0.2 0.0

1

2

3a

3b

4

Figure 7: Variation in 4f occupation of the larger complexes with different KS functionals. Plotted are deviations from formal f occupations of 0 for La, Ce, Gd β -spin and 6 for Yb β -spin. The CC results of the larger systems were deemed unreliable and are omitted here, as explained in the text. dative (donation) bonding from filled ligand orbitals. Dative bonding can be viewed as a specific form of electron delocalization. The lone pair of a ligand atom facilitating dative bonding with a metal can in this case not be localized to a one-center orbital (i.e. a proper lone-pair). Rather, upon localization of the orbitals of the complex the ‘lone pair’ exhibits two- (or multi-) center character, with metal orbital contributions reflecting the donation bonding that goes along with an overall stabilization of the system. The DE (or – potentially – a degeneracy-driven mechanism) may cause the extent of the electron donation to be exaggerated, which would truly impact the electron density and/or spin density. Such changes can be quantified, for instance, via a population analysis. It is worthwhile mentioning that the set of natural atomic orbitals (NAOs) used in the NPA is orthonormal, i.e. metal-centered NAOs are orthogonal to ligand orbitals and vice versa. For the lanthanide 4f orbitals that are of prime concern here, the overlap with the ligands prior to the orthonormalization step is relatively weak and should not affect the interpretation of the results (further information is provided in Section 3.3). Further, the NAOs are created from a system’s density matrix after covalent interactions have taken place, and any electron density build-up or loss relative to the lanthanide ion’s formal 4f and 5d occupation must be from or to the ligands, respectively. The underlying basis sets used in the calculations are non-orthogonal. Figures 6 and 7 display the 4f populations from the NPA for our samples, (excluding the Lu complexes) as deviations from the formal 4f occupations. CC calculations were performed for all systems, but for complexes 2, 3a-c and 4 the results were deemed unreliable due to low weights of the references in the CC wave functions and large values of the T1 diagnostic. Attempts to use Brueckner or DFT orbitals had negligible effects on these criteria. The CC calculations for the lanthanide halides for which populations are given in Figure 6 employed a larger than default 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

set of non-frozen orbitals (orbitals with energies between ±350 eV were correlated) because of a noticeable impact on the 4f populations compared to using default settings. The La and Ce systems exhibit sizable 4f populations due to L → M dative bonding. The question of whether the 4f populations are a consequence of donation into polarized 5d orbitals, or ligand → 4f donation, is addressed in the following subsection. Irrespective of the mechanism by which electron density ends up in 4f orbitals, the populations are consistent with the functional-dependence of the DE: For a given complex, increasing 4f population goes along with an increasingly positive DE. For the Ce(IV) complexes, this trend is very pronounced and the 4f populations are considerable, presumably owing to the higher positive charge of the ion. The situation for high-spin Gd(III) (4 f7 ) is more complicated for the following reason: The calculations produce both a degree of dative L → M bonding with donation into the empty 4f β -spin orbitals, as well as some delocalization of the occupied 4f α -spin orbitals. A plot for the total 4f populations of GdX3 can be found in the SI, Figure S8. When only the β -spin 4f orbitals are considered, as in Figure 6, clear trends in line with the other systems emerge. The Gd 4f populations are overall smaller than those of the early lanthanides, which seems intuitive as the 4f shell contracts along the lanthanide series with increasing atomic number. The data for the Yb 4f13 complex 4 are also for the β -spin orbitals. The KS calculations give large excess 4f populations with a strong dependence on the functional. The data for 4 in Figure 7 are very similar to Figure S8 in the SI giving the total 4f population. We note that the large excess 4f populations calculated for 4 with the non-hybrid functionals are somewhat larger than those of YbX3 model systems with X = Cl, Br, I, shown in the SI, but not dramatically so. The 4f populations for the Ce and La systems obtained from the CC calculations lie in between those of HF and the TLC functional calculations. The tuned LC population is reasonably close to the CC population for LaX3 and GdX3 . The population from the TLC calculation for 1 is much larger than the CC population. As the latter remains rather close to HF, this raises the question whether CCSD sufficiently corrects over HF, and whether TLC sufficiently corrects the DE of the KS calculations as far as the results from the NPA are concerned. We are presently unable to settle this question. 2D-tuning for 1 increases the predicted 4f population by 0.25 relative to LCPBE0 and above those obtained with the global hybrids, but as for all other compounds the non-hybrid functionals produce the largest 4f populations. For complex 2, the populations from tuned KS calculations are similar to those calculated with LCPBE0, while CAMB3LYP, PBE0 and B3LYP give increasingly large 4f populations but not dramatically so. Bogart et al.7 reported a 4f population of ca. 0.8 using B3LYP, which is consistent with our hybrid KS data when accounting for differences in the basis sets used. For the Yb complex 4 specifically, a single-reference CC calculation is problematic even though the T1 diagnostic value did not raise serious concerns (0.009). The HF reference enters 14 ACS Paragon Plus Environment

Page 14 of 29

Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8: Fractional spin calculations E(α − β ) for YbCp3 for selected KS functionals. The value of α − β indicates 2hSz i; for ±1 the KS determinants are MS = ±1/2 eigenfunctions. the CCSD wavefunction with a weight of only 53%, which indicates multi-reference character and likely renders any NPA-derived CC 4f population unreliable. Indeed, as mentioned in the Introduction, 4 has previously10 been assigned a double-reference ground state mixing a LMCT 4f14 configuration with the 4f13 configuration underlying the present HF and KS calculations. According to Reference 10, the weight of the LMCT configuration is on the order of 11% and therefore we may expect a corresponding 4f population in excess of 13 (or in excess of 6 for β -spin only). The HF calculation produces essentially no deviation from the idealized population of the 4f13 configuration, whereas all KS calculations yield excess populations. The non-hybrid calculations appear to be severely deficient in the sense that the population is around 0.55 and therefore strongly overestimated relative to TLC and the global hybrid functionals. However, the non-hybrid functionals are potentially less affected by the static correlation error,18 which is likely to play a role for 4 due to the assigned multiconfigurational ground state. The static correlation error can be probed via calculations with fractional spins: The unpaired electron can be described via fractional occupations of an α versus β spin orbital in a spin-unrestricted calculation such that the total electron number remains the same. Details can be found in Reference 18. The energy from fractional spin calculations for 4 with selected functionals is plotted in Figure 8. The energies for the MS = ±1/2 systems are identical. For exact KS calculations, the energy should also be the same in the fractional spin systems. As Figure 8 clearly demonstrates, this is not the case in calculations with approximate spin-density functionals or eX only (HF). Large deviations of the energies of the fractional spin systems from the energies for the MS = ±1/2 systems indicate a more pronounced static correlation error. Among the tested functionals, the non-hybrid functional performs best. However, as the GGA calcula15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9: Symmetry species of different d and f orbitals in D∞h and D3h symmetry. The calculated a′2 and a′′2 canonical MOs of LaF3 are also shown (TLC functional, isosurfaces ±0.03, contour lines in selected planes ±0.01 to 1 logarithmically spaced). Notation for the symmetry species as in Reference 67. tions suffer the most from the DE the 4f populations are nonetheless likely to be too large in these calculations. To summarize thus far: The Lu systems show no evidence of mixing of the 4f shell with ligand orbitals such that it would cause 4f delocalization. For the La(III) and especially the Ce(IV) systems, all calculations, including CCSD and TLC, produce sizable 4f populations. For highspin Gd(III), donation into 4f β -spin orbitals still amounts to on the order of 5% of an electron, with some calculations also producing a degree of delocalization of the α -spin 4f orbitals. The KS calculations for the Yb(III) complex 4 and YbX3 model systems give surprisingly large excess 4f populations. For 4, the multi-reference character of the ground state complicates the picture, and approximation in the functionals are likely to cope differently with the static correlation error versus the DE. The non-hybrid functionals appear to exaggerate the 4f population when compared with the findings of Reference 10; the TLC calculation produces more reasonable results. Overall, for the set of complexes studied the 4f populations correlate with the magnitude of the DE as it manifests itself in the curvature of E(N). 3.3

Donation into 5d versus 4f

It has been argued that the 4f populations in complexes of early lanthanides are the result of a polarization of 5s,p orbitals and of the polarization of 5d orbitals into which the ligands donate electron density.24 An apparent population of 4f ‘orbitals’ would then be due to fully and partially occupied d orbitals acquiring some 4f character, in order to describe their polarization and to increase the overlap with the ligands. As it has been stated previously,68, 69 ‘a little hybridization goes a long way’ in the formation of covalent interactions. Focusing on LaX3 first, Figure 9 lists the irreducible representations of d and f metal atomic orbitals in D3h and D∞h symmetry. We 16 ACS Paragon Plus Environment

Page 16 of 29

Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 10: 5d and 4f orbital populations (NPA) in LaF3 from calculations with a PBE-based global hybrid functional with varying percentage of exact exchange (eX). For the remaining LaX3 systems see Figure S10 in the SI. The ‘F 2p loss’ data represents the loss of population relative to fully occupied fluorine 2p orbitals. Below is a set of contour plots of a canonical a′2 MO (see Figure 9) in steps of 25% eX illustrating the decreasing metal contribution. Contour lines ±0.01 to ±1, logarithmically spaced. The circles drawn around the La center have a radius of roughly half the covalent radius of lanthanum. posit that the population of 4f arising from ligand-4f mixing in canonical MOs belonging to the a′2 and a′′2 species gives a better idea about ligand → 4f versus 5d donation versus d-polarization because for these symmetries d and f do not mix. Figure 9 shows plots for two canonical MOs of LaF3 where ligand p orbitals mix with 4fσ and one of the 4fφ , respectively, in canonical MOs that exclude d-character because of their symmetry. Figure 10 shows the 5d and 4f populations of LaF3 calculated with a PBE-based global hybrid with varying percentage of eX in the exchange functional, as well as the corresponding loss of fluorine 2p population. The latter is practically equal to the sum of donation into La 5d and 4f orbitals. Donation from F 2p into La 5d is overall more pronounced than donation to 4f; the 5d populations are about double the 4f populations across the range from the least eX (positive DE) to largest eX percentage (negative DE). Polarization of the filled outer core shells by 4f and 5d basis functions should lead to corresponding 4d and 5s,p NAO populations of less than the formal ones. The numerical data (not shown) indicate that this effect is negligible for 4d, and small, around 0.015 for 5s,p. The bulk of the 4f population must be attributed to a different mechanism. We note in passing that the energies of the canonical MOs that are dominant in 4f character vary considerably in energy among the different functionals. For example, using the 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 11: Symmetry species of different d and f orbitals in Oh symmetry. Isosurfaces (±0.03) of selected canonical MOs of CeCl62– (TLC functional). PBE global hybrid with 0 to 100 % eX, the energy of the lowest 4f-dominant canonical MO in LaF3 varies between -2.5 and 6.75 eV. For GdF3 , the range is -10 to -25 eV. For comparison, the HOMO energies in these systems vary considerably as well, between -15.75 and -7.75 and -7.5 to -16.5 eV for the La and Gd complex, respectively. In these two examples, occupied f orbitals are stabilized by increasing eX while unoccupied f orbitals are destabilized by increasing eX. Of the five NBOs calculated for LaF3 with dominant 5d character and assigned as ‘low occupancy lone pairs’, i.e. 5d acceptor orbitals of dative bonding (occupations between 0.025 and 0.04), two have only 5% f character. This may qualify as predominantly a polarization of the 5d orbitals. The other three NBOs have between 18 to 31% f character, which may be more appropriately interpreted as a d-f hybridization where a significant portion of electron density that may initially be donated to 5d ends up in the 4f shell. Finally, the sum of 4f NAO populations corresponding to symmetries other than those of the d orbitals, i.e. a′2 and a′′2 , constitute between 30 and 40% of the total f populations given in Figure 6. The situation is similar for the other LaX3 complexes studied herein and therefore we forego further discussion of these results. Additional relevant NBO data on select complexes can be found in Figure S2. Based on these calculations and the data shown for LaF3 , it would seem appropriate to assign at least some of the 4f populations determined by the NPA to ligand → 4f donation rather than mere polarization of a partially occupied 5d shell. In the case of octahedral symmetry, due to the presence of the inversion center, the gerade d and the ungerade f orbitals belong to separate irreducible representations which can mix with ligand-centered combinations of the same symmetry. Therefore, even if d-f hybrid orbitals are generated in a subsequent localization procedure (such as NBO/NLMO used in this study), the 4f population arises from canonical MOs that do not mix with d orbitals. Figure 11 lists the irreducible representations to which the d and f orbitals belong in the Oh point group. Also shown are two isosurfaces corresponding to one of each of the t1u and t2u canonical MOs of 18 ACS Paragon Plus Environment

Page 18 of 29

Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CeCl62– (1). There are two sets of t1u and one set of t2u occupied ligand orbitals with pronounced Ce 4f character, while the LUMO has a2u symmetry. In contrast to LaF3 , where all 5d and 4f populations arise from low-occupancy non-bonding one-center orbitals at the metal, as classified by the NBO algorithms, the L → M donation in 1 is strong enough such that a threshold is crossed and the algorithms instead construct six two-center ligand-metal bonding NBOs along with additional low-occupancy metal non-bonding orbitals with populations donated by the ligands. Löble et al. deduced 10% Cl 3p character for 1 in unoccupied metal-centered orbitals according to experimental XAS data and 15.5% Cl 3p character according to KS calculations.8 The results are usually interpreted in terms of a corresponding set of occupied ligand-centered orbitals with sizable metal character. The bonding Ce−Cl NBOs from our TLC calculation each afford a 10% contribution from Ce. This Ce contribution is assigned 50% f character and 33% d character (these localized orbitals do not belong to any of the symmetry species and may therefore afford f-d hybridization). The bonding NBOs have essentially double occupation, resulting in about 0.1 4f population per NBO. Accompanying the bonding NBOs are complementary low-occupancy (ca. 0.05) antibonding NBOs with 90% contribution from Ce with the same percentages of d and f character, creating f populations of 0.023 each. The remaining source of f population resides in low-occupancy non-bonding NBOs of pure f character. Figure 6 shows that TLC and B3LYP calculations for 1 create very similar 4f populations; the latter functional was employed for the study of Löble et al. From the CC density matrix, bonding NBOs are created as well, with about 3% Ce f character, and the remaining f population is assigned to low-occupancy one-center metal NBOs. See Figure S2 for a compilation of pertinent NBOs in select lanthanide halides. It is worthwhile reiterating that whether bonding or one-center NBOs are created in the calculation depends on internal thresholds in the NBO algorithms but this should not have any influence on the population analysis. The NPA projects the part of the total electron density (or spin density) associated with f basis functions onto a basis of 4f, 5f, etc. atomic orbitals constructed by a different set of criteria. In the canonical orbitals that parametrize the density in the KS calculations there is no hybridization of f and d in Oh symmetry. If the 4f atomic orbitals created by the NBO code are representative of a Ce atom, e.g. in their spatial compactness, there must be ligand → 4f donation be present in 1. The result varies with the type of KS functional vs. HF vs. CCSD, but all of them produce substantial 4f populations (Figure 6). The radial density maximum of the 4f shell of Ce in a 6s2 5d1 4f1 configuration is at 0.72 a0 (atomic units) according to numerical relativistic HF calculations.70 Lanthanum in a 4f1 configuration would have a 4f shell with a comparable location of the radial density maximum. The La and Ce basis sets afford a contracted f function with large exponents to describe this 4f shell, along with more diffuse 4f functions. The 4f functions with the lowest exponents in the basis sets 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: 6s, 4f and 5d NAO populations, in parts per thousand (ppt), for selected 4f0 systems calculated with different numbers of 4f basis functions. ‘Full 4f’ denotes the populations from a calculation with all 4f basis functions included. ‘Diffuse 4f’ refers to inclusion of only the two most diffuse 4f functions. ‘No 4f’ indicates a basis without 4f functions. Full 4f Diffuse 4fa No 4fb 6s 4f 5d 4f + 5d 6s 4f 5d 4f + 5d 6s 5d LaF3 HF 13 83 185 268 13 61 192 252 14 199 LCPBE0 16 130 284 414 17 63 306 369 19 330 B3LYP 15 142 323 466 16 54 360 414 12 396 PBE 15 199 374 573 16 52 440 491 18 486 LaCl3 HF 50 104 320 424 51 77 332 409 54 336 LCPBE0 56 148 403 550 59 76 429 506 63 445 B3LYP 50 162 460 622 53 63 505 568 55 532 PBE 51 224 500 723 53 60 577 637 55 610 2– CeCl6 HF 240 500 462 961 248 309 514 823 253 555 LCPBE0 257 904 492 1396 277 313 669 982 283 731 B3LYP 232 1080 493 1573 255 285 772 1057 259 847 PBE 235 1302 490 1791 264 279 875 1153 267 959 Only kept La (Ce) exponents α = 0.6098, 0.1973 (α = 0.6364, 0.2164). b The NBO program requires at least one 4f function in the basis in order to complete without errors. A dummy 4f function with exponent 10−6 was supplied in order to be able to complete the NBO step of the calculation.

a

cover the spatial range of the 5d set (Rmax (5d) ∼ 2a0 ) and therefore mainly serve as polarization functions. Table 2 provides NAO populations (in ppt) of two Ln(III) halides and CeCl62– , all formally in 4f0 configurations. The calculations were performed with the full set of 4f functions, with the diffuse set only, and without 4f functions. The impact of the DE, giving large 5d and 4f NAO populations, is readily apparent. Adding the ‘diffuse’ 4f functions mainly polarizes the 5d shell and increases the overall donation from the ligand, presumably due to better overlap. Compared to the pronounced variations in the 5d populations with the functional in these calculations, the 4f polarization contribution stays relatively constant. This is not the case with the full 4f basis including the contracted functions representing the La and Ce 4f atomic orbitals. The functional dependence of the 4f population is now even more pronounced than for 5d. Moreover, the 5d populations are markedly lower than in the calculation with only polarizing 4f functions. This may indicate that the donation to 4f is facilitated by the 5d shell. Table S1 in the SI provides additional population data from calculations for the same La(III) 20 ACS Paragon Plus Environment

Page 20 of 29

Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12: Selected 4f NAOs of LaF3 and CeCl62– obtained with different functionals. The insets show isosurfaces of the NAOs (±0.03). Plots of the NAO from the metal nucleus along a horizontal line to the right with respect to the molecular orientations shown in the insets. The graphs reflect the NAO radial functions R(r) and radial densities r2 R2 . The positions of the inner radial density maxima are compatible with numerical atomic relativistic HF data from Reference 70. and Ce(IV) systems as in Table 2. In these calculations, small-exponent d functions were removed from the basis and the full set of f functions kept. One set of calculations used only one contracted d function representing the occupied 4d shell, and another set used the 4d contracted basis function plus one additional uncontracted d function with a radial density maximum of 0.86 a0 , i.e. slightly larger than the 4d and 4f density maxima but much smaller than that for the 5d shell. In both sets of calculations there are appreciable 4f populations. With only the contracted basis functions representing 4d, the 4f populations are far above those in Table 2 for ‘Diffuse 4f’ and not very far below the ‘Full 4f’ data. Interestingly, with the one additional uncontracted d function, the 4f populations of Ce alone in CeCl62– become comparable to the combined 4f and 5d populations for the full basis in Table 2. Moreover, for all three systems the 4f populations increase sharply with the presence of the added d function and become larger than the ‘Full’ 4f populations in Table 2. These trends suggests that the second d basis function acts as a conduit for donation into the 4f shell if there is no 5d shell available for L → M donation. 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The information collected so far indicates that there is donation into 4f ‘proper’ taking place, and that its magnitude is rather sensitive to the DE. The donation is somewhat less pronounced than the calculated 4f NAO populations would seem to indicate, because the NAOs afford a diffuse part in the spatial 5d polarization region. These features are displayed graphically in Figure 12 for a representative NAO of LaF3 and CeCl62– , respectively. Corresponding plots for the β -spin NAOs of GdF3 and GdCl3 appear very similar to those of LaF3 (Figure S11). The orbitals are plotted on an axis pointing from the metal nucleus out along a maximum of the corresponding angular 4f function and the graphs therefore reflect the shapes of the NAO radial function and radial density. The insets display corresponding isosurfaces. The NAO set is orthonormal and therefore the orbitals have small orthogonalization ‘tails’ on the ligand atoms. In the radial functions, the pronounced peaks positioned less than 1 a0 away from the metal, and the corresponding radial density maxima are compatible with the numerical atomic relativistic HF data of Reference 70 and therefore correspond to a 4f atomic orbital. But there is also density further out than 2a0 that we must attribute primarily to the 5d polarization by 4f basis functions. The integrated density is reflected accordingly in the NAO populations. However, if the TLC calculations can be trusted not to suffer excessively from the DE affecting the metal - ligand interactions and the balance between 5d vs. 4f populations, then a considerable fraction (in the case of Ln(III)) or even most (in the case of Ce(IV)) of the 4f population is not attributed to 5d polarization but to ligand → 4f donation. 3.4

Degenerate orbital mixing

In order to investigate if a degeneracy-driven 4f mixing with ligand orbitals may cause noticeable electron density changes, a situation was constructed where a formally empty 4f shell and filled ligand orbitals may become degenerate. We focused initially on LaF3 and varied the eX fraction in a PBE-based global hybrid functional in steps of 5%, as in Figure 10. Similar calculations were performed for the other LaX3 , CeCl62– , and LuF3 , with the results shown in the SI, Figure S10. Further, La3+ orbital energies were computed in a field of three negative point charges corresponding to the LaF3 geometry, and corresponding orbital energies for a system of three F– were computed in the crystal field of a +3 point charge accounting for the La ion. The metal 4f and ligand 2p fragment orbital energies showed an energetic crossover between 5% and 10% eX. However, across the range of eX percentages in Figures 10 and S10 there is no sudden increase in either 5d or 4f population that may indicate degeneracy-driven covalency causing a shift in the electron density. The SI contains a numerical analysis of the solutions of the generalized 2 × 2 eigenvalue problem with parameters chosen to match approximately the fragment orbital energies of LaX3 22 ACS Paragon Plus Environment

Page 22 of 29

Page 23 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from the point charge models mentioned above. The model indicates that in the regime of weak overlap a near-degeneracy driven mixing would have a rather sudden onset compared to a case where the fragment orbitals differ more. However, without any overlap the interaction matrix element between the fragment orbitals would be zero as well and mixing may occur in this case only for an exact degeneracy. This situation would not correspond to covalency as understood by chemists, just like the canonical MO wavefunction of H2 , which can be expressed as a 5050 mixture of a covalent and an ionic valence-bond structure, would not appear to chemists as representing covalency at internuclear separations large enough that the overlap is negligible. If a near-degeneracy driven as opposed to overlap/DE driven mixing can be established via variations of the parameters in a KS functional, the model suggests that one may perhaps expect a rather sudden increase or decrease of the 4f populations. The gradual decrease of the 4f populations in Figures 10 and S10 in unison with those of 5d with increasing fraction of eX suggests would seem to favor a DE driven mechanism with small – but non-negligible – overlap instead. The L → M donation bonding is in this case understood in terms of the usual ‘chemist’s orbitals’ picture, overlap (even if small), non-vanishing Fock matrix elements between donor and acceptor orbitals, and not too dissimilar fragment orbital energies. The analysis in the previous subsection does not rule out a L → M donation mechanism that takes place via involvement of the lanthanide 5d shell.

4 Summary and Conclusions In KS calculations of lanthanide complexes, the extent of mixing of ligand orbitals with orbitals from the 4f and 5d shells is dependent on the functional. The 4f and 5d populations obtained from a natural population analysis (NPA) follow the same trend as the KS delocalization error (DE). The latter can be quantified numerically via the curvature of the energy E(N) for electron numbers N between integers, either by calculations with explicit fractional electron numbers of by the interpolation of Equation (2). However, for lanthanide complexes with empty or partially filled 4f shells, L → M donation bonding may occur and raise the 4f populations above the formal 4f electron count of the lanthanide ion. This goes along with dative bonding to the metal 5d shell. While some of the 4f populations arise from basis functions polarizing the 5d shell, there also appears to be genuine ligand → 4f donation. Since the extent of the donation bonding is sensitive to the DE it may be important to investigate the DE numerically in applications of KS theory that are targeted at interpreting experimental data in terms of 4f (or 5f) involvement in covalent metal - ligand interactions. For the closed 4f shell of the Lu systems studied herein, there is no evidence of 4f delocalization in the population analysis or in the localized orbital representation of the electronic structure. 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

However, the extent of the DE correlates well with the extent of mixing between 4f and ligand orbitals in the canonical MOs. We should point out that the correlation of the 4f populations in the NPA with the DE-measures derived from the E(N) curvatures is not a direct proof that the DE causes the excess 4f population, although the fact that the correlation is present for a variety of different lanthanide complexes strengthens such an argument. A more traditional ‘orbital diagram’ based rationale would require that, when the KS Fock matrix is represented in a basis of metal and ligand fragment orbitals, the trends in the 4f populations among the different types of functionals go along with a closer match of ligand donor and Ln 4f acceptor orbital Fock matrix elements. However, such a match would likewise correlate with the DE. Attempts were made to create a degeneracy-driven mixing for small model systems of La, Ce, and Lu with – perhaps – a more sudden appearance and disappearance of mixing. However, the population data exhibits a monotonic decrease of the 4f (and 5d) population for La and a constant 4f population of 14.00 for Lu, with increasing fractions of eX in the functional. The results of the present study pertain to metal complexes as they would be encountered in solution or gas phase. For complex ions in a solid, the functional-dependent extent of L → M donation and the localization of the Lu 4f shell are likely to be very similar. However, in general the KS functional requirements for the computation of the electronic structure of periodic solids of f-element compounds, in particular of ionic solids, are notably different from those for molecules,71, 72 and the case of L → 4f donation in solids with lanthanide ions would appear to warrant further study.

Acknowledgments We thank Dr. Eric Schelter for providing us with optimized coordinates of complex 2 and the Center for Computational Research (CCR) at the University at Buffalo for providing computational resources. This research has been supported by the National Science Foundation, grant CHE-1265833.

Supporting Information Available Additional computational data, xyz coordinates for all systems. This material is available free of charge via the Internet at http://pubs.acs.org.

References [1] Hargittai, M. Chem. Rev. 2000, 100, 2233–2302. 24 ACS Paragon Plus Environment

Page 24 of 29

Page 25 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[2] Kovács, A.; Konings, R. J. M. J. Phys. Chem. Ref. Data 2004, 33, 377–404. [3] Choppin, G. R. J. Alloys Compd. 2002, 344, 55–59. [4] Maron, L.; Eisenstein, O. J. Phys. Chem. A 2000, 104, 7140–7143. [5] Eisenstein, O.; Maron, L. J. Organomet. Chem. 2002, 647, 190–197. [6] Desreux, J. F.; Reilley, C. N. J. Am. Chem. Soc. 1976, 98, 2105–2109. [7] Bogart, J. A.; Lewis, A. J.; Medling, S. A.; Piro, N. A.; Carroll, P. J.; Booth, C. H.; Schelter, E. J. Inorg. Chem. 2013, 52, 11600-11607. [8] Loble, M. W.; Keith, J. M.; Altman, A. B.; Stieber, S. C. E.; Batista, E. R.; Boland, K. S.; Conradson, S. D.; Clark, D. L.; Pacheco, J. L.; Kozimor, S. A.; Martin, R. L.; Minasian, S. G.; Olson, A. C.; Scott, B. L.; Shuh, D. K.; Tyliszczak, T.; Wilkerson, M. P.; Zehnder, R. A. J. Am. Chem. Soc. 2015, 137, 2506-2523. [9] Moncho, S.; Autschbach, J. Magn. Reson. Chem. 2010, 48, S76-85. [10] Denning, R. G.; Harmer, J.; Green, J. C.; Irwin, M. J. Am. Chem. Soc. 2011, 133, 2064420660. [11] Dolg, M.; Fulde, P.; Stoll, H.; Preuss, H.; Chang, A.; Pitzer, R. M. Chem. Phys. 1995, 195, 71-82. [12] Kerridge, A.; Coates, R.; Kaltsoyannis, N. J. Phys. Chem. A 2009, 113, 2896-2905. [13] Dolg, M.; Mooßen, O. J. Organomet. Chem. 2015, 794, 17-22. [14] Edelstein, N. M.; Allen, P. G.; Bucher, J. J.; Shuh, D. K.; Sofield, C. D.; Kaltsoyannis, N.; Maunder, G. H.; Russo, M. R.; Sella, A. J. Am. Chem. Soc. 1996, 7863, 13115–13116. [15] Mooßen, O.; Dolg, M. Chem. Phys. Lett. 2014, 594, 47–50. [16] Neidig, M. L.; Clark, D. L.; Martin, R. L. Coord. Chem. Rev. 2013, 257, 394–406. [17] Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Science 2008, 321, 792-794. [18] Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Chem. Rev. 2012, 112, 289-320. [19] Autschbach, J.; Srebro, M. Acc. Chem. Res. 2014, 47, 2592-2602.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[20] Ooms, K. J.; Feindel, K. W.; Willans, M. J.; Wasylishen, R. E.; Hanna, J. V.; Pike, K. J.; Smith, M. E. Solid State Nuc. Mag. Reson. 2005, 28, 125-134. [21] Dolg, M.; Cao, X.; Ciupka, J. J. Electron Spectrosc. Relat. Phenom. 2014, 194, 8–13. [22] Clavaguéra, C.; Dognon, J.-P.; Pyykkö, P. Chem. Phys. Lett. 2006, 429, 8–12. [23] Ramakrishnan, R.; Matveev, A. V.; Rösch, N. Chem. Phys. Lett. 2009, 468, 158–161. [24] Xu, W.; Ji, W.-X.; Qiu, Y.-X.; Schwarz, W. H. E.; Wang, S.-G. Phys. Chem. Chem. Phys. 2013, 15, 7839–7847. [25] Ji, W.-X.; Xu, W.; Schwarz, W. H. E.; Wang, S.-G. J. Comput. Chem. 2015, 36, 449–458. [26] Kaltsoyannis, N. Inorg. Chem. 2013, 52, 3407-3413. [27] Baer, R.; Livshits, E.; Salzner, U. Ann. Rev. Phys. Chem. 2010, 61, 85–109. [28] Pritchard, B.; Autschbach, J. Inorg. Chem. 2012, 51, 8340-8351. [29] Srebro, M.; Autschbach, J. J. Phys. Chem. Lett. 2012, 3, 576-581. [30] Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. J. Chem. Theory Comput. 2012, 8, 1515-1531. [31] Sun, H.; Autschbach, J. J. Chem. Theory Comput. 2014, 10, 1035–1047. [32] Moore II, B.; Srebro, M.; Autschbach, J. J. Chem. Theory Comput. 2012, 8, 4336-4346. [33] Yanai, T.; Tew, D. P.; Handy, N. C. Chem. Phys. Lett. 2004, 393, 51-57. [34] Stepanov, N. F.; Ustenko, A. A.; Dement’ev, A. I. Vestn. Mosk. Univ., Ser. 2: Khim. 1973, 14, 102-104. [35] Parr, R. G.; Yang, W. Density functional theory of atoms and molecules; Oxford University Press: New York, 1989. [36] Johnson, E. R.; Otero-de-la Roza, A.; Dale, S. G. J. Chem. Phys. 2013, 139, 184116. [37] Myers, C. E.; Graves, D. T. J. Chem. Eng. Data 1977, 22, 436–439. [38] Tsukamoto, S.; Mori, H.; Tatewaki, H.; Miyoshi, E. Chem. Phys. Lett. 2009, 474, 28–32.


Page 26 of 29

Page 27 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[39] Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Dam, H. J. J. V.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. Comp. Phys. Comm. 2010, 181, 1477–1489. [40] Bylaska, E. J. et al. “NWChem, A Computational Chemistry Package for Parallel Computers, Version 6.3 (2014 developer’s version)”, Pacific Northwest National Laboratory, Richland, Washington 99352-0999, USA., 2014. [41] Riplinger, C.; Neese, F. J. Chem. Phys. 2013, 138,. [42] Neese, F. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73–78. [43] Reed, A. E.; Weinhold, F. J. Chem. Phys. 1985, 83, 1736-1740. [44] Weinhold, F. Natural bond orbital methods. In Encyclopedia of computational chemistry; von Ragué Schleyer, P., Ed.; John Wiley & Sons: Chichester, 1998 pages 1792–1811. [45] Glendening, E. D.; Badenhoop, J. K.; Reed, A. E.; Carpenter, J. E.; Bohmann, J. A.; Morales, C. M.; Weinhold, F. “NBO 5.0, http://www.chem.wisc.edu/˜nbo5”, Theoretical chemistry institute, University of Wisconsin, Madison, 2001. [46] Dolg, M.; Stoll, H.; Preuss, H. J. Chem. Phys. 1989, 90, 1730-1734. [47] Cao, X.; Dolg, M. J. Chem. Phys. 2001, 115, 7348-7355. [48] Cao, X.; Dolg, M. J. Mol. Struct.: THEOCHEM 2002, 581, 139–147. [49] Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650-654. [50] Curtiss, L. A.; McGrath, M. P.; Blandeau, J.-P.; Davis, N. E.; R. C. Binning, J.; Random, L. J. Chem. Phys. 1995, 103, 6104-6113. [51] Bergner, A.; Dolg, M.; Küchle, W.; Stoll, H.; Preuß, H. Mol. Phys. 1993, 80, 1431–1441. [52] Dunning, T. H. J. Chem. Phys. 1989, 90, 1007-1023. [53] Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1993, 98, 1358-1371. [54] Sessler, J. L.; Mody, T. D.; Hemmi, G. W.; Lynch, V. Inorg. Chem. 1993, 32, 3175–3187. [55] Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213. [56] Cao, X.; Dolg, M.; Stoll, H. J. Chem. Phys. 2003, 118, 487-496. 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[57] Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257-2261. [58] Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650-654. [59] Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865-3868. [60] Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. [61] Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–789. [62] Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200-1211. [63] Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623-11627. [64] Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158-6170. [65] Becke, A. D. Phys. Rev. A 1988, 38, 3098-3100. [66] Perdew, J. P. Phys. Rev. B 1986, 33, 8822-8824. [67] Altmann, S. L.; Herzig, P. Point-group theory tables; Clarendon Press: Vienna, 2nd ed.; 2011. [68] Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. [69] Claxton, T. A. Nature (London) 1965, 208, 891–892. [70] Desclaux, J. P. At. Data Nucl. Data Tab. 1973, 12, 311-406. [71] Hay, P. J.; Martin, R. L.; Uddin, J.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 034712. [72] Wen, X.-D.; Martin, R. L.; Henderson, T. M.; Scuseria, G. E. Chem. Rev. 2013, 113, 1063–96.


Page 28 of 29

Page 29 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents / Abstract Graphics


Impact of the KohnâSham Delocalization Error on ... - ACS Publications

Recommend Documents

Impact of the KohnâSham Delocalization Error on ... - ACS Publications