A litmus test for the balanced description of dispersion interactions and coordination chemistry of lanthanoids †

The influence of long-range interactions on the structure of complexes of Eu( III ) with four 9-hydroxy-phenalen-1-one ligands (HPLN) and one alkaline earth metal dication [Eu(PLN) 4 AE] + (AE: Mg, Ca, Sr, and Ba) is analyzed. Through the [Eu(PLN) 4 Ca] + complex, which is a charged complex with two metals—one of them a lanthanoid—and with four relatively fluxional p -ligands, the diﬃculties of describing such systems are identified. The inclusion of the D3(BJ) or D4 corrections to diﬀerent density functionals introduces significant changes in the structure, which are shown to stem from the interaction between pairs of PLN ligands. This interaction is studied further with a variety of density functionals, wave-function based methods, and by means of the random phase approximation. By comparing the computed results with those from experimental evidence of gas-phase photoluminescence and ion mobility measurements it is concluded that the inclusion of dispersion corrections does not always yield structures that are in agreement with the experimental findings.


Introduction
Due to its efficiency in comparison with correlated wave-function theories, Kohn-Sham density functional theory (KS-DFT) has become the most widely used quantum chemical method in chemistry and physics. Even though convergence in accuracy and predictability are not hierarchically defined as clearly as with correlated wave-function methods, 1 efforts in this direction have for many years been graphically represented by what is called Jacob's ladder. 2 Despite the endeavors made in the design and development of hierarchical enhancements of the KS-DFT methods, there are still certain intrinsic limitations in the exchange-correlation functionals such as the self-interaction error (SIE), deficiencies in describing strong correlation, and the lack of long-range dynamic correlation. 3,4 This latter aspect becomes visible in the difficulties of KS-DFT to describe noncovalent interactions.
Fortunately, during the past two decades, different approaches have been developed in order to overcome these drawbacks. 5 Two of the most extended modi operandi have been the following: First, the attainment of dispersion interactions directly from the electron density by defining non-local correlation (NLC) functionals. This approach has brought about the so-called van der Waals functionals (vdW-DF) and represents a new conception within KS-DFT since it adds a description of the dispersion interactions directly to the DFT functionals by combining all kinds of correlations. Examples in this direction are the functionals vdW-DF1, 6 and vdW-DF2 7,8 and the simpler and more efficient approaches of Vydrov and van Voorhis (VV9 and VV10). [9][10][11][12][13] The second most used approximation, especially within the field of chemistry, is the so-called DFT-D methodology [14][15][16][17][18][19] consisting in adjoining a semi-classical correction [see eqn (1)] to the known local functionals, which accounts for the long-range van der Waals interactions between particles in the gas phase described through the attractive correction À1/R n (n = 6, 8,. . .), that is: Here, R AB is the distance between atoms A and B, C AB n denotes the nth-order dispersion coefficient for the atom pair A and B, s n is a scaling factor that depends on the functional used and ensures the exact asymptotic behavior if C AB n is exact, and f damp,n (R AB ) is a damping function that determines the range of the dispersion correction avoiding near singularities for small R AB and double-counting effects of correlation at mid-range distances. An important part of the work developed within this last approach has been carried out by Grimme and collaborators through their models D1, 14 D2, 15 D3 16 -with the original damping function D3(0) or the Becke-Johnson damping function D3(BJ) 17 -and D4. 18,19 While the DFT-D approach is by far the computationally most efficient, and therefore perhaps the most popular dispersion correction, the approaches D1, D2, and D3 have the limitation of considering only the molecular geometry and the environment of each atom pair so that changes in electron density are not reflected in the dispersion coefficients C AB n [eqn (1)]. A step further has been achieved with the most recent D-model, D4, which includes classical atomic partial charges in addition to the geometrical information. The D4 model shows, to date, a promising performance in case of conformational energies of d-block transition metals, 20 interaction energies of some charged molecules with p-systems, 21 and polarizabilites of periodic systems formed by inorganic salts of elements of the groups 1-5. 22 As mentioned above, the inclusion of dispersion interactions has enabled the use of DFT to describe non-covalent inter-and intra-molecular interactions 23 which have been intensely studied primarily in unsaturated organic molecules, in so-called p-p-stacked structures, in biologically relevant molecules such as nucleic acids and proteins, 24,25 and more scarcely also in complexes where metal-metal and metalligand interactions were analyzed (see, e.g., ref. [26][27][28]. Challenges to be explored with all aforementioned approaches are those in which lanthanoids and actinoids are present and especially when different forms of long-range interactions concur simultaneously. Examples are charged systems in which partially flexible p-groups and diverse metals with a wide coordination sphere get together. In these cases, the dilemma between accuracy and computational feasibility becomes a decisive issue when it comes to finding a qualitatively correct answer. 29 Some years ago, some of us worked on the characterization of Eu(III) complexes with different numbers of the 9-hydroxyphenalen-1-one (HPNL) ligand. [30][31][32][33][34] This ligand was chosen for its high absorption cross section, its long-lived phosphorescence, and its ability to coordinate and transfer energy to Eu atoms. [35][36][37] In these studies, the structures of the complexes were always theoretically predicted making use of DFT by using the BP86 38 functional in combination with the def2-SVP 39,40 and def2-TZVPP 39,40 basis sets. The computed values could confirm the experimental evidence obtained through mass selected gas-phase photoluminescence (PL) emission spectra as well as, in some cases, through ion mobility measurements. In the course of this research we came across a class of heterodimetallic complexes with four ligands, that is, Eu(III) 9-oxo-phenalen-1-one coordination complexes [Eu(PLN) 4 AE] + with AE = Mg, Ca, Sr, and Ba (see Fig. 1), 34 for which the BP86/def2-TZVPP approach predicted four different structural motifs.
The two motifs lowest in energy, hereinafter called motif 1 and motif 2, were consistent with a six-fold coordinated Eu(III) ion and a five-fold coordinated alkaline earth ion. Motif 1 (M1) can be described as a complex where one PLN unit coordinates only with the Eu(III) ion while the other three are coordinated to both metal ions. In motif 2 (M2), on the contrary, one PLN unit coordinates only with the alkaline earth dication, whereas the other three form a propeller-like configuration in which the ligands are coordinated to the Eu(III) ion similar to leaf blades to a stem. 31,41,42 For [Eu(PLN) 4 Ca] + and [Eu(PLN) 4 Sr] + motif 2 is only 7-9 kJ mol À1 more stable while for [Eu(PLN) 4 Mg] + motif 1 is by 8 kJ mol À1 the lower conformation. The experimental evidence arising from the ion mobility analysis and the gas-phase photoluminescence (PL) measurements could not rule out at the time any of these energetically almost degenerate structures. 43 The presence of the PLN ligands, typical p-systems, made us question if long-range correlation interactions in these systems might be relevant in distinguishing between motifs 1 and 2. Thus, the present work focuses on the investigation of dispersion effects in motif 1 (M1) and on the analysis of their influence on the molecular structure. For this purpose, the D3(BJ) and D4 models were primarily used and contrasted, when possible, with reference values obtained employing wave-function based methods, as well as with other approaches for evaluating dispersion effects within the framework of KS-DFT.
The remaining part of this manuscript is organized as follows: in Section 2 the theoretical and experimental methodology followed in this work is described. Section 3 is mainly focused on the structure of the [Eu(PLN) 4 Ca] + complex. For this purpose, the effect of using different functionals, the inclusion of non-covalent interactions, and the influence of the starting structures in the geometry optimization process are first analyzed in Section 3.1. A comparison with the experimental data is carried out in Section 3.2. Section 4 elaborates on the results of Section 3 first by briefly discussing many-body dispersion effects (Section 4.1). A benchmark study of the p-p interactions between the PLN ligands is carried out in Section 4.2, the results are contrasted with values obtained using the random phase approximation (RPA) and with reference values obtained from wave-function-based CCSD(T) 44 and CCSD(F12)(T) 44,46,47 computations. In Section 4.3 the role of the alkaline earth and lanthanoid metals coordination sphere in the evaluation of the dispersion interactions by means of the D3(BJ) and D4 models is discussed. Finally, Section 5 summarizes the present work with the observation that the system studied here provides a litmus test for measuring the delicate balance between dispersion interactions and coordination chemistry of europium complexes.

Computations
All quantum chemical computations were performed using the Turbomole program package. 48 DFT calculations were carried out using a spherical grid 5 for numerical integration as defined in the Turbomole program. 49 Unless stated otherwise, the convergence thresholds were 10 À8 E h and 10 À5 E h /a 0 for the self-consistent field (SCF) energy and for the Cartesian gradient, respectively. The resolution-of-the-identity approximation (also called density fitting) was applied in all cases. [50][51][52] Using a Fermi-type expression dynamically determined, floating occupation numbers in the SCF computation were obtained in order to determine the number of unpaired electrons.
In the present work, some of the most popular DFT functionals were employed. In increasing order of sophistication when describing the exchange-correlation contribution (Perdew's metaphorical Jacob's Ladder), 2 the used functionals can be listed as follows: BP86, 38 BLYP, 53 and PBE 54 functionals within the generalized gradient approximation (GGA), the TPSS, 55 SCAN 56 and M06-L 57 functionals belonging to the metageneralized gradient approximation (MGGA) functionals, the non-empirical PBE0 58 functional, the semi-empirical B3LYP 59 and BHLYP 53,59 functionals, and the meta-hybrid functionals TPSSh, 60 M06, 61 and M06-2X. 61 The latter five belong to the class of hybrid functionals (HYB) with different percentages of exact Hartree-Fock exchange. Finally, two range-separated hybrid functionals including dispersion interactions, that is, oB97M-V 62 and oB97X-D, 63 were also tested. For those functionals that do not include long-range dynamic correlation (dispersion), the pairwise dispersion corrections D3 14-16 with Becke-Johnson damping function [D3(BJ)], 17 and D4 18,19 as well as the nonlocal VV10 approach 13,64 were incorporated to the already established local functionals. Moreover three van der Waals functionals, that is, optB88-vdW, 65 BEEF-vdW 66 and mBEEF-vdW 67 based on a non-local electron density dependent correlation (NLC) were also added to the test set. Additional improvements beyond the pairwise additive description of the long-range electronic correlation were carried out by means of the use of the random phase approximation (RPA) within the framework of the adiabatic-connection fluctuation-dissipation theorem in the DFT formulation [68][69][70][71][72] and as defined for ground-state energies and analytic first-order properties by Furche et al. [73][74][75][76] The RPA approach is directly related to the post-Hartree-Fock methods sharing with them also a high computational demand. 73,[77][78][79] Combining DFT with a post-Hartree-Fock method, the double-hybrid functional B2PLYP 80 with second order perturbation theory correlation energy was also used for energies. In the context of the wave-function theory, second-order Møller-Plesset perturbation theory (MP2) 81 as well as explicitly correlated MP2 (MP2-F12) 82 were considered, too, accounting for dynamic electron correlation effects. As reference values for the energies calculated in this work, the coupledcluster singles and doubles with perturbational treatment of triple excitations [CCSD(T)] method 44 and the CCSD(F12) method [45][46][47] were employed. For all wave-function post-Hartree-Fock methods, the core electrons were not correlated (fc). Unless otherwise stated, the def2-TZVPP basis set 39 was used in all computations.
As in previous works, 31,33,34 the energy levels of the complexes were analyzed and simulated using an effective ligandfield Hamiltonian 83-89 based on the theory of many-electron spectra first described by Racah. [89][90][91][92][93] The computations with this effective Hamiltonian were carried out with the McPHASE package. 94,95 In particular the Wybourne normalized ligandfield parameters, L kp , 94 were obtained with the pointc module employing the natural population analysis 96 point charges calculated with different DFT functionals and the def2-TZVPP basis set. The F k and z parameters were used as optimized for the nonanuclear [Eu 9 (PLN) 16 (OH) 10 + ] complex. 33

Experiment
In the gas-phase photoluminescence (PL) measurements, the desolvated ions were mass-selected, stored in an quadrupole ion trap, thermalized via collisions with a helium buffer gas at about 0.2 mbar and 85 K, and photoexcited by the 458 nm line of a continuous wave (cw) Ar + -laser with its beam orthogonal to the ion trap axis. 31,32 Description of the experimental setup has been reported in detail elsewhere. 34,97 Structural information about gas-phase ions at room temperature was obtained using gas-phase ion mobility. Ion mobility mass spectrometry is a technique that spatially separates ions under the influence of an electric field by affecting their motion via collisions with a buffer gas. 98,99 Via the recording of the ions' arrival time at the detector, and calibrated against known references, their cross-sections can be determined accurately. Furthermore, distinguishing isomers with an identical mass-over-charge ratio was achieved at an ion mobility resolution of 150. As described in ref. 34, the N 2 collision cross-sections inferred from trapped ion mobility spectrometry (TIMS) measurements on a timsTOF (Bruker, Bremen, Germany) at 300 K and a pressure of 0.923 mbar N 2 (exit funnel) was compared to scaled He-collision cross-sections computed using both MOBCAL 100-102 and IMOS. [103][104][105] It is to be noted that the He-collision cross-sections computed using MOBCAL and IMOS are in perfect agreement. The scaling relationship, O N2 = 53.9 + 1.110O He , was inferred from the comparison of the He and N 2 cross sections reported for denatured polyalanine by Bush and coworkers. 106 For the sake of simplicity, only He-collision cross-sections will be reported in this work.
[Eu(PLN) 4 AE] + complexes (Fig. 1). As in previous studies with Eu(III) compounds, the structures of these complexes were optimized at the BP86/def2-TZVPP level and in order to investigate more in detail motif 1, dispersion interactions were added. For this purpose, the D3(BJ) model was applied as a first approach to the problem. Fig. 2 shows the different optimized geometries obtained with and without dispersion corrections for motif 1 (M1-D3 and M1, respectively). [For Motif 2, for example, the geometrical changes are considerably less pronounced by including D3(BJ) dispersion interactions (see Section I of the ESI †).] With D3(BJ) corrections, the two-by-two p-p stacking of the PLN ligands apparently represents a key driving force resulting in a significant molecular reorganization for this motif (M1-D3).
Energetically, for motif 1 the M1-D3 structures are considerably more stable than their M1 and their motif 2 [M2-D3] counterparts [see Table 1 of Section I of the ESI †] suggesting that the formation of the stacked PLN pair is associated with a significant attraction energy gain.
As a major structural change occurs upon the inclusion of dispersion corrections resulting in the stacking of the PLN ligands in the M1-D3 structure, a more detailed analysis of the structure was accomplished for one of the complexes of the series, namely [Eu(PLN) 4 Ca] + . In this case the two structures obtained for this complex, that is, M1 (BP86/def2-TZVPP) and M1-D3 (BP86-D3(BJ)/def2-TZVPP) with open-open and stacked configurations, respectively (see Fig. 3), were taken as starting points for the optimization of the structure of this complex using a selection of DFT-functionals. Cartesian coordinates of M1 and M1-D3 are provided in Section II of the ESI. † In particular, one or two functionals were chosen within each rung of Jacob's ladder, that is, GGA, MGGA, HYB, MHYB, and range-separated hybrid. 107 For the double hybrid B2PLYP functional, gradient calculations and therefore geometry optimizations are currently not implemented in Turbomole. The inclusion of dispersion corrections by means of the DFT-D models D3(BJ) and D4 as well as by incorporation of the nonlocal functional VV10 was also investigated. Additionally, the van der Waals functionals optB88-vdW, BEEF-vdW and mBEEF-vdW, albeit especially suitable for solid-state calculations, were also tested. 108 The results including energy differences and qualitative descriptions of the structures are given in Table 1.
For all tested functionals without extra dispersion corrections, and regardless of the starting point geometry, open-open   4 Ca] + optimized with different density functionals as well as relative energies (DE = E (starting geometry = M1-D3(BJ)) À E(starting geometry = M1)) in kJ mol À1 including harmonic zero-point energy effects (DE 0 ), and the corresponding Gibbs free energies under the conditions of the photoluminiscence (TLIF: P = 0.2 mbar, T = 85 K) (DG TLIF ) and ion mobility (IM: P = 0.923 mbar, T = 300 K) (DG IM ) measurements. All computations were performed with the def2-TZVPP basis set (oo = open-open; ss = stacked; hs = half-stacked) Not only open-open structures but also a stacked configuration as well as a new type of structure (half-stacked) are obtained (see Fig. 4). With methods like TPSS-D4 and TPSSh-D4 the stacked and half-stacked structures are even practically isoenergetic (see Table 1). Although the electronic energy differences between the three structural configurations are of the order of a few kJ mol À1 (at most 24 kJ mol À1 ), the stability trend confirmed in all cases is, starting from the most stable, as follows: With the aim of identifying the nature of the critical points found in the geometry optimization using different methods, second derivatives of the potential energy were computed for all of the optimized geometries, verifying that all structures correspond to minima on the potential energy surface. 109 Apart from the almost isoelectronic character of the three types of structures found, the zero-point vibrational energy (ZPE) of these three configurations did not differ by more than 2 cm À1 and the rotational constants, above all the B constant, depart by no more than 50 MHz within a particular functional and with and without dispersion corrections. The latter can clearly be observed in Table 1 where the different ground state energies are presented.
It is expected that Gibbs free energy contributions (DG) to the electronic energy using the conditions of the photoluminescence experiments (T = 85 K and P = 0.2 mbar) and the ion mobility measurements (T = 300 K and P = 0.923 mbar) will lead to relative stability results closer to what it would be expected experimentally. An estimate of DG was carried out via the module freeh of Turbomole. In a first approximation the harmonic approximation for vibrations and the rigid rotor model for the rotations were used in the calculation of the corresponding partition functions. The results obtained under the latter approach are presented in the column denoted as RR in Table 1. Additionally and considering that among the 264 vibrational frequencies of [Eu(PLN) 4 Ca] + , 12-17 have values below 100 cm À1 and 34-39 lie below 200 cm À1 , a treatment of low-lying vibrational modes was carried out using the free-rotor approximation 110 whose results are denoted as FR in Table 1.
In general, with the inclusion of the experimental conditions, the stability trend already observed continues to be confirmed, that is: stacked 4 half-stacked 4 open-open. The results obtained using the rigid-rotor model (RR) in comparison with those obtained using the free-rotor (FR) model seem to be slightly overestimated. Exceptions to the aforementioned bias are the results obtained with the TPSS-D4, and TPSSh-D4 methods where under the pressure and temperature conditions of both experiments the half-stacked structure is slightly more stable than the stacked structure, although only in the order of 5 kJ mol À1 . Another exception is observed with the B3LYP-D3(BJ) method where under the experimental conditions of the ion mobility experiment the open-open structure is slightly more stable (only 2-3 kJ mol À1 ) than the stacked structure.
Regarding the uncertainties of the very small DFT energy differences presented, as well as the quality of the (harmonic) frequencies, the results from the the analysis of the Gibbs free energies at all conditions under consideration should be considered with caution. We therefore turn our attention in the next section to the comparison with the experimental evidence.

Comparison with the experimental information
The different structures resulting from most of the functionals (with and without dispersion corrections) employed in this work were compared with the experimental values arising from the ion mobility measurements (Fig. 5) and the photoluminescence spectrum (Fig. 6).
The comparison of the He-converted cross-section of all structures obtained computationally with those obtained from the ion mobility measurement, suggest again (cf. Fig. 5) that only the open-open configuration is in agreement with the experimental information and both stacked and half-stacked structures can be clearly ruled out. Ligand field theory is currently the only practical model to analyze and simulate the level structure of lanthanoid ions at an accuracy of about 50 cm À1 . In the present work, the Stark splitting of the 7 F 2 level as well as the related transition energies ( 5 D 0 -7 F 2 band) were computed to accomplish the comparison with experiment. The experimental broadening was simulated using Gaussian functions, whose widths are consistent with the experimental resolution and whose superposition provide the maxima displayed. As it can be observed in the few examples shown in Fig. 6, the open-open structure is consistently the one that is most in line with the experimental information both in energy range and number of splits, whereas the half-stacked and stacked structures progressively separate from the experimental profile regarding both features.
Thus, although the M1-D3 structures are energetically the most stable, they poorly agree with experiment, showing that at least in the present type of experimental conditions and resolution they are not present in the gas-phase measurements.

Many-body dispersion effects
Being [Eu(PLN) 4 Ca] + an open shell system with 458 electrons and requiring 2390 basis functions for calculations using the def2-TZVPP basis set, geometry optimizations using the RPA approach or wave-function-based methods such as MP2 or even CCSD(T) were not possible at this time. Therefore, in order to make an assessment of the performance of these two methods in this complex, the lanthanoid ion metal, Eu 3+ , was replaced by indium [In 3+ ]. In 3+ has a comparable size to Eu 3+ with the advantage of having a closed-shell electronic structure with 444 electrons for the [In(PLN) 4 Ca] + complex. In addition to the above mentioned computationally more demanding methods, a few DFT functionals were also tested for evaluating the comparability of the original complex with Eu 3+ with that with In 3+ (see Table 2 of ESI † and Section SIII). Regarding the results obtained with the RPA approach as well as with the MP2 wavefunction-based method, the BP86/RPA calculations led to an open-open structure when the starting structure is M1. However, they result in a stacked structure when M1-D3 is the starting point in the optimization. Stacked or half-stacked configurations are obtained at the MP2 level.

Benchmark study on the stacked PLN ligands
As the major structural change observed upon inclusion of dispersion corrections was the stacking of the PLN ligands in the M1-D3 structure, a more detailed analysis of the energy contributions in this configuration was accomplished (see Section IV of the ESI †). To investigate the problem more in depth, the first step was to select a model where dispersive interactions play a major role. Therefore, the Ca 2+ and Eu 3+ ions in the stacked M1-D3 structure obtained at the BP86-D3(BJ)/ def2-TZVPP level were removed and a proton was added to each of the four ligand anions (PLN) resulting in two neutral dimers (PLNH) 2 , for which the positions of the added H atoms were optimized using the BP86-B3(BJ)/def2-TZVPP level of theory while leaving the rest of the structure fixed. Thus, two models were designed, which are denoted as Ligandpair12 and Ligand-pair34 (see Fig. 7). The interaction energy of these model structures was based on the reaction 2 HPLN Ð PLNH ð Þ 2 , and computed according to the following expressions: As is well known, the magnitude of the corrections due to dispersion effects depends on the functional used. Therefore, the next step to be considered was to contrast the dispersion interaction contributions obtained with different functionals corresponding to various rungs of Jacob's ladder. 2 In addition, all results were compared with reference values obtained from wave-function based methods that cover electronic correlation in a more complete way. The most accurate method that could be used for this (still fairly large) molecular system is the CCSD(T) method. Moreover, the latter was corrected for the  45,46 For consistency, in all cases the def2-TZVPP basis set was employed. However, in order to assess the adequacy of this basis set for the explicitly correlated calculations (F12), a comparative study between the def2-TZVPP and cc-pVTZ-F12 111,112 basis sets employing various methods with and without F12 corrections was carried out on three different configurations of the benzene dimer [see Section V of the ESI]. In the case of Ligandpair12 and Ligand-pair34, the tests carried out using second-order perturbation theory (MP2) with and without explicit F12 correlation show that the differences obtained with both basis sets amount to only 1 kJ mol À1 , which is consistent with the results obtained for the benzene dimer. In addition, we compared MP2-F12 results obtained with the cc-pVTZ-F12 and cc-pVQZ-F12 basis sets 111,112 to those obtained with the def2-TZVPP and found only a minor difference of 4-5 kJ mol À1 . In view of these results, the use of the def2-TZVPP basis set was justified for the benchmark study on Ligandpair12 and Ligandpair34.
In Fig. 8, the interaction energies of the Ligandpair12 obtained with different functionals in comparison with the CCSD(T) and CCSD(F12)(T) reference values are displayed. (Similar results were obtained for Ligandpar34, see Fig. 6 of the ESI. †) The CCSD(F12)(T) reference values were obtained by adding to the CCSD(F12) energies the triples energy contribution from the CCSD(T) calculation. Without dispersion corrections, many interaction energies are fairly repulsive. Among all functionals without dispersion corrections, the highly parametrized Minnesota functionals, M06, M06-L and M06-2X, 113 which show a certain ability to embrace medium-range electron correlation, [114][115][116] as well as the SCAN functional and the double-hybrid B2PLYP achieve the best performance. As already mentioned in the literature, 19,117 the inclusion of dispersion interactions using the atom-pairwise semi-classical models D3(BJ) and D4 remarkably improves the results for all of the functionals, yielding attractive overall interaction energies. It can be observed that precisely, and only, for the functional used as a starting point in the present work, that is, BP86, the D3(BJ) and D4 models overestimate (in an absolute sense) the dispersion interactions by about À14 kJ mol À1 and À8 kJ mol À1 , respectively, in comparison with the CCSD(F12)(T) reference value. With the other functionals, the interaction energies for Ligandpair12 range between À41 and À57 kJ mol À1 while the CCSD(F12)(T) reference value is À51 kJ mol À1 . As already mentioned, the atom-pairwise D3(BJ) model 16 is based on the molecular geometry to estimate the intra-and intermolecular dispersion energies and therefore the changes in electron density are not considered in the dispersion coefficients. With the model termed D4, 18,19,21,22 this last limitation is corrected with the incorporation of charge dependent dispersion coefficients. This improvement of the D4 model yields for the two cases studied here-two typical p-p-interaction systems-changes of only +2 to À7 kJ mol À1 for the dispersion energy, with the D4 energies slightly more attractive than those with D3(BJ). The good performance of the range-separated hybrid functional oB97X-D 63 that includes dispersion corrections following the DFT-D scheme is particular noteworthy, see Fig. 8.
As mentioned in the introduction, the inclusion of longrange correlation effects can also be carried out through the incorporation of a fully nonlocal functional to the already established (exchange and correlation) functionals. The simpler and computationally cheaper nonlocal functional VV10 by Vydrov and van Voorhis was employed here, the results for Ligandpair12 using some of the functionals are: BLYP (À52 kJ mol À1 ), PBE (À43 kJ mol À1 ), and TPSS (À49 kJ mol À1 ). (see Fig. 9). These results are in agreement with previous studies in which it was found that, in general, nonlocal correlation functionals tend to overestimate the long-range dispersion interactions with errors ranging between 24% and 49%. 118 Here we also stress the good performance of the range-separated hybrid  functionals oB97M-V and oB97X-V of Head-Gordon et al. 119 that include the non-local electron-density dependent correlation functional VV10. oB97M-V and oB97X-V show deviations of less than 5 kJ mol À1 from the CCSD(F12)(T) reference values. As it can be observed in Fig. 8 and 9, the wave-function based methods MP2 and MP2-F12 overestimate the interaction energy. Moreover, more complete accounts of electronic correlation in the framework of KS-DFT such as the random phase approximation (RPA) considerably improve, as expected, the predictions of the interaction energies in all cases.

Lanthanoid and alkaline earth metals coordination
Once the ligand-ligand interactions on the molecular configuration of [Eu(PLN) 4 Ca] + had been analyzed, and in view of the structural results reported in the previous section, the next aspect to be considered is the coordination structure of the lanthanoid (Eu 3+ ) and the alkaline earth metal (Ca 2+ ) ions and their influence on the quantification of the London interaction energies. With its large coordination sphere, europium takes usually part in complexes with coordination numbers (CN) of 8, 9, or even higher, prompting coordination geometries of the dodecahedron, octahedron, square antiprism, cube, and biand tricapped trigonal prism types, 120-127 and molecular point-group symmetries such as O h , T d , D 3d , D 3h , and others. In the case of calcium the coordination number can significantly vary, from 9 to other much smaller, resulting in different coordination geometries. 128 As mentioned in the previous sections, in case of the [Eu(PLN) 4 Ca] + complex, the coordination number of the europium and calcium ions are 6 and 5, respectively, presenting consequently a certain distortion with respect to the ideal 6-coordination polyhedra 129,130 and being close to an octahedral configuration. 34 Having said this, in the framework of DFT-D, the concept of coordination numbers plays a significant role. In fact, the term ''fractional coordination numbers'' (fCN) was already introduced in the D3 model 16 trying to differentiate between the diverse hybridization states of the atoms in a molecule and thus to account for different atom pair dispersion coefficients [C AB 6 , see eqn (1)], making them dependent on the ''chemical'' environment while remaining a purely geometric approximation. Just as this concept is defined, 16 it leads to fCN values that are very similar to the conceptually defined CNs, as they are understood in coordination chemistry, in particular for molecules with atoms of the first and second rows. Typical values are 2, 3, and 4 for the carbon atom with sp, sp 2 and sp 3 hybridizations, respectively. However, the fCN tended to be very large and different from integer values in the case of metallic systems. 19 With the D4 model, the calculation of the fCNs is slightly modified by including in their evaluation the atomic electronegativities of the two atoms involved in a bonding pair and by recasting the fCNs such that they roughly correspond to the Wiberg bond orders 131 for single bonds. In addition to this, some adjustments are made in the Gaussian function weighting procedure. All this leads to a better description of the bonding environment and the hybridization type. 18,19 Using the DFT-D3 16 and DFT-D4 18,19 programs, values for the fCN of europium and calcium in the [Eu(PLN) 4 Ca] + complex were calculated for different DFT functionals and compared with the chemical coordination described in the first paragraph of this section (Table 4 of the ESI †). Values for the total dispersion energy of the complex as well as that between the Ca-Eu metal pair are also included in Table 4 of the ESI †. To point out that in the case of the DFT-D4 model, the electric charges are independent of the functional used and calculated following a classical model based on the electronegativity equilibration of Gaussian-type charge densities. 19,132 In general, in the case of atoms of the first and second row and in comparison with the chemical coordination numbers, the D3(BJ) model tends to slightly overestimate the fCN values, whereas the D4 model slightly underestimates them. In the case of the two metals, the differences are quantitatively more relevant. For all the considered geometries, the D3(BJ) model estimates the fCN for europium close to 8 while that for calcium close to 7. In the case of the D4 model, the fCN is around 5 for the two metals, which is in good agreement with the chemical coordination number for calcium but somewhat small for europium. For the half-stacked and stacked structures, the fCN tends to be slightly larger than in the case of the open-open structures, especially in the case of calcium. As expected, with the model D3(BJ) the dispersion energy between the two metals is between 3 and 4 times larger than with the model D4, but in any case this energy does not represent more than about 1% of the total dispersion energy. The type of structure, open-open, half-stacked, or stacked, does not establish a difference in the dispersion interaction between metals. The total dispersion energy is in general greater in absolute terms with the D4 model except for the results obtained with the BP86 functional where |E disp-D3(BJ) (Total)|4|E disp-D4 (Total)| in the case of stacked structures.
Summarizing, it can be said that the D4 model quantifies the fCN values much more closely to the chemical CNs than the D3 model. However, this does not mean necessarily that the molecular structures resulting from the optimization, in view of the results, are more plausible with the D4 model than with the D3(BJ) model.

Conclusions and summary
The influence of non-bonded interactions on the structure of Eu 3+ heterodimetallic complexes with four ligands, i.e., europium(III) 9-oxo-phenalen-1-one coordination complex [Eu(PLN) 4 Ca] + was studied. The BP86-D3(BJ) and BP86-D4 methods were found to overestimate the energy of the p-p interaction between the PLN ligands in absolute terms. Although the D4 model takes into account not only the molecular geometry but also the atomic partial charges and provides values for the ''fractional coordination numbers'' closer to the conceptional coordination numbers of the metals, the optimized structures obtained with BP86-D4, BLYP-D4, TPSS-D4, and TPSSh-D4 are not in better agreement with the experimental evidence. For the cases analyzed, the nonlocal VV10 approximation leads to structural results (half-stacked or stacked configuration) that do not concur with the photoluminescence and ion mobility data. The three tested van der Waals functionals conclude with structures that are in agreement with the experimental findings. Explorations with the In(III)-model complex, that is, [In(PLN 4 Ca] + indicate that the results with the RPA and MP2 can lead to qualitatively wrong results. In summary, in this work a molecular example was found where the evaluation of non-covalent interactions was not straightforward. Thus, a litmus test for the balanced description of dispersion interactions and coordination chemistry was revealed that may become a helpful benchmark system for future method developments in the framework of DFT as well as improvements in the context of computationally efficient wavefunction based methods. Based on the findings of the comparison with the two experiments in Section 3.2, the following simple litmus test is proposed: for a given method two geometry optimizations are carried out starting from BP86/def2-TZVPP and BP86-D3(BJ)/def2-TZVPP optimized geometries (see ESI †). A method passes the test if at least one structure optimization yields an open-open structure. Of course the test cannot provide any information if a method passes or fails for the right reasons. Of the methods under study BP86-D3(BJ), BP86-D4, BLYP-D4, BLYP-VV10, PBE-VV10, TPSS-D4, TPSS-VV10, and TPSSh-D4 fail the test, because the dominant species found in two experiments cannot be reproduced even starting from a structure that is consistent with both experimental findings. The proposed litmus test is passed by all other methods employed in this work, of these BP86, BLYP, PBE, TPSS, B3LYP, BHLYP, TPSSh, PBE0, optB88-vdW, BEEF-vdW, and mBEEF-vdW yield only openopen structures.
The system under study represents a difficult balancing act between dispersion interactions and metal coordination in which pairwise semi-classical approaches like the D-models can lead to misleading results as previously observed, for example, in the study of cryogenic ion vibrational predissociation (CIVP), 133 bond dissociation energies 134 as well as energies, 135 prediction of vibrational spectra 136 and interconversion barriers. 137 We finally note that as shown in this work the inclusion of dispersion corrections, at least via the D3(BJ), D4, and VV10 approaches, although decisive for the description of the interaction energy between the stacked ligandpairs, represents only a necessary but not sufficient condition for an unambiguous prediction of the molecular structures. Therefore, including dispersion corrections of this kind is not always the preferred option when optimizing geometries with DFT. Special attention must be paid to cases where the dispersion corrections influence results significantly.

Conflicts of interest
There are no conflicts to declare.