Determination of structures and binding energies of europium(III) complexes bound to biologically relevant anions

Abstract

Anion sensing is a crucial area of research with applications in medical diagnostics, environmental monitoring and industrial analysis. Lanthanide complexes are actively studied for these purposes, particularly in biological media, as their positively charged central cations can provide strong electrostatic interactions with the respective analytes and because their unique photophysical properties allow for efficient detection. Within this work, we investigate the suitability of various computational methods to determine the geometries of europium(III) complexes and their binding energies to anions of biological relevance, putting special emphasis on the influence of conformational flexibility. We find that the composite r²SCAN-3c method produces accurate geometries at a low computational cost and that numerical stability as well as computational efficiency are increased when using a Y atom as proxy for Eu. We highlight the importance of conformational flexibility, necessitating an appropriate conformer search algorithm. Finally, when studying binding energies in host–guest complexes, we find that, despite investigating a range of methodological influences, we are unable to find a computational protocol that achieves even semi-quantitative accuracy. In summary, this work provides a solid starting point for future investigations on luminescent lanthanide(III) complexes and their anion sensing properties while also illustrating the open challenges.

---

1. Introduction

The properties of lanthanide complexes have been investigated extensively in the last two decades because of their numerous useful applications.¹ Lanthanide complexes are employed in several important fields such as bioimaging,^2–4 optoelectronic devices,^5–7 memory devices,^8–10 ionic liquids^11–13 and selective anion sensing.^14–17 In the field of anion sensing, lanthanide-based complexes are advantageous because they feature a positively charged central cation that can form strong electrostatic interactions with anions.¹⁸ Anions, in turn, play a crucial role in biological, industrial and agricultural processes, and hence, a significant effort is made by researchers for the development of novel and improved anion sensors.¹⁹ Europium(III) complexes are particularly well suited for anion sensing in biological media, and hence, these have been used as sensors for detecting kinase protein,²⁰ biomarker proteins in human serum,^21,22 adenosine²³ and bicarbonates.^17,24

The electronic structure of lanthanides is largely responsible for their unique photophysical properties,¹¹ such as sharp emission bands²⁵ and large Stokes shift of about 200 nm,²⁶ making them useful to serve as sensors. Sensors based on these complexes have distinctive selectivities, and their detection limits are similar to those of commonly employed organic fluorophores like coumarin, fluorescein and rhodamine.²⁷ However, they are considered superior probes to conventional organic fluorophores due to their very narrow emission bands having a full-width half maximum of less than 10 nm which reduces spectral interference and enhances the precision with which emission changes can be monitored during sensing.²⁸ In addition, these complexes have long-lived luminescence, which allows sensitive time-resolved detection of emission and prevents interference from background excitation and scattered radiation.^29–32

The structural design of ligands to make efficient lanthanide complexes has been an important goal of several researchers,^3,33 as the properties of these complexes can be improved for more sensitive and selective detection of specific ions and molecules by carefully designing the structure of ligands around the central ion. For this purpose, quantum mechanical modelling can help by providing molecular insights and energetics of the relevant processes.

Researchers have used quantum chemical methods to identify important structural features of these complexes, for instance, the binding modes of guests to the complexes.³⁴ The quantum chemical studies have also been proven to be useful in verifying the presence of important chemical interactions such as π–π stacking and hydrogen bonding in them.^35,36 They have assisted in understanding the increased selectivity of one guest molecule over others by providing a deeper understanding of the binding geometries.³⁷ In addition to structural aspects, theoretical methods have provided an understanding of the photophysical behaviour of the complexes, such as sensitisation of the lanthanide ion by ligands.^3,38,39

Having a computational method which can accurately describe the structure of discrete lanthanide complexes is crucial for advancing their applications. Earlier work involved the development of various variants of the Sparkle model as extensions to different semi-empirical methods.^40–42 Subsequently, Grimme and co-workers have shown that their general purpose GFN-xTB method outperforms Sparkle/PM6 for geometries of lanthanide complexes.⁴³ The next step by this group amounts to the development of a force-field (GFN-FF) for lanthanides that can describe the lanthanide-to-ligand bond lengths with a mean absolute deviation of 0.14 Å.⁴⁴ As an alternative option, Okayasu and Yuasa revealed that DFT based predictions, in combination with experimental techniques, can be used to determine the configuration of europium(III) complexes in solution.⁴⁵

The modelling of lanthanide complexes is challenging due to their open-shell electronic structure, large number of possible conformers,²⁰ flexible coordination geometries spanning coordination numbers usually from 7 to 10 as well as the presence of multiply charged species and the influence of solvation. Considering the importance of these complexes and the difficulties in their design process, our aim is to present a workflow that closely reflects the objectives and challenges identified in our previous studies of related complexes,^34–36,46 thereby establishing an approach optimised for anion-recognition applications.

Our work begins by comparing geometries obtained using various computational methods and evaluating them against experimental crystal structures. For this purpose we choose the four experimentally synthesised europium(III) complexes [Eu.1],⁴⁷ [Eu.L2],⁴⁸ [Eu.L3],⁴⁹ and [Eu.L5],⁴⁹ based on either cyclen (1,4,7,10-tetraazacyclododecane) or TACN (1,4,7-triazacyclononane) macrocyclic scaffolds, as shown in Fig. 1. [Eu.1] has been developed as a probe for peroxynitrite and features a coordinatively saturated nonadentate ligand.⁴⁷ [Eu.L2] was synthesised by Parker and coworkers from the tris-(4-bromopyridyl) europium(III) complex in a two-step, one-pot sequence and has three-fold symmetry in solution; [Eu.L3] and [Eu.L5] were also synthesised from the same precursor complex using Sonogashira coupling reactions. [Eu.1] has an approximate square antiprismatic geometry, while [Eu.L2], [Eu.L3] and [Eu.L5] have tricapped trigonal prismatic coordination geometry differing in the nature of the donor atoms (phosphinate or carboxylate donors) and peripheral substituents.^48,49 After testing how well geometries are reproduced when starting from the crystal structure, we move to the more realistic setting of evaluating how well they can be determined in a completely blind approach. In this context, we emphasise the importance of appropriate conformer search algorithms. As a final step, we move to the even more challenging question of predicting binding energies and geometries in host–guest complexes. For this purpose, we choose two derivatives of [Eu.1], the complexes [Eu.pB(OH)₂]⁺ and [Eu.mB(OH)₂]⁺, also shown in Fig. 1.^36,46 These complexes are experimentally known to bind to phosphate, adenosine monophosphate (AMP), as well as adenosine diphosphate (ADP) and we wanted to investigate their binding properties using computation.⁴⁶

Fig. 1 Chemical structures of the six europium(III) complexes investigated within this work.

2. Computational details

2.1. Electronic structure

Geometries of the lanthanide(III) complexes were optimised with the GFN-FF^44,50 force field, the semi-empirical GFN2-xTB⁵¹ method, and density functional theory (DFT) using various functionals (ωB97M-V,⁵² B3LYP⁵³ with D4 dispersion correction^54,55 and r²SCAN-3c⁵⁶) and basis sets (def2-SVP and def2-TZVPP;⁵⁷ 6-31G* and 6-311G**⁵⁸). Here, r²SCAN-3c utilises mTZVPP basis set, the D4 London dispersion corrections,^54,55 and a geometric counter poise (gCP)⁵⁹ correction to account for basis set superposition error (BSSE).⁶⁰ Noting that gCP corrections are already present in r²SCAN-3c, we do not include an explicit counterpoise correction. The basis sets def2-SVP and def2-TZVPP are paired with the small-core ECP28MWB (=def2-ECP)⁶¹ replacing 28 core electrons. In this instance, we attempt computations on Eu(III) in heptet multiplicity and, as a simpler analogue with no f-electrons, also Y(III) as a singlet. Y(III) is considered here as a computational model for Eu(III) and all experimental data considered are based on Eu(III). Indeed, Y(III) is often used instead of Eu(III) in DFT computations because it is diamagnetic, has a 4f⁰ configuration, while its charge and ionic radius are similar to Eu(III), giving comparable geometries without f-electron complications.^39,62–65 It is also worth noting that the f-electrons in lanthanides are usually assumed to not take part in the bonding and have negligible angular effects.⁶⁶ Thus, geometries should not be notably affected by replacing the central ion as long as the ionic radii are similar.

The basis sets 6-31G* and 6-311G** are paired with the large-core ECP52MWB,⁶⁷ which also represents the f-electrons and replaces in total 52 electrons of Eu(III) and computations are performed in a pseudo-singlet configuration. The notations Y/28, Eu/28 and Eu/52 are used in the discussion to refer to the use of small and large-core potentials, respectively. We, thus, have three different levels of representing the europium atom and it is not a priori clear which one is the best. While the small core computation with Eu(III) seems to be the highest level, we note a severe limitation. The true ⁷F₀ ground state of Eu(III) is roughly spherically symmetric due to static correlation within the f-orbitals. As a consequence, the seven f-orbitals should have an overage occupation of 6/7 electrons each. But within the present computations, six of the f-orbitals are singly occupied, and one is left empty, thus breaking the spherical symmetry. Even more, the system now contains seven different quasi-degenerate ⁶f-configurations with the potential of causing severe problems in terms of numerics and reproducibility. By contrast, the ECP52MWB and Y(III) computations preserve the spherical symmetry of the core and semi-core electrons per construction and are, thus, expected to be numerically more stable (as also seen in our preliminary computations).

To perform DFT calculation with ωB97M-V functional (6-31G* and 6-311G** basis sets and ECP52MWB) we used Q-Chem 6.1⁶⁸ while, for the other calculations, Orca 5.0^69–71 interfaced to libxc⁷² was used. In these calculations, water was included as solvent via the conductor-like polarizable continuum model (CPCM).⁷³ Semi-empirical calculations were carried out using xTB 6.5.1.⁷⁴ Water solvation in the xTB calculations is included via analytical linearized Poisson–Boltzmann (ALPB) implicit solvation model.⁷⁵

To explore the conformational space of the complexes, the conformer-rotamer ensemble sampling tool (CREST version 2.12), based on a direct sampling scheme for the generation of conformational ensembles at the xTB level using a Y(III) central atom was employed.⁷⁶ Settings used in CREST are: non-covalent interaction (NCI) mode, scaling factor of 0.6 for outer wall potential, metadynamics length of 50 picoseconds and energy threshold of 10 kcal mol⁻¹ (=42 kJ mol⁻¹). It is worth mentioning that it is important to optimise the conformations generated by CREST, at DFT level because their relative energies at semiempirical and DFT levels can be very different.^77,78 If more than 50 structures were generated, only the 50 lowest-energy ones were considered for further optimisation with r²SCAN-3c to keep computational effort at a reasonable level.

Subsequently, for the anion binding energies, geometries were optimised with r²SCAN-3c and single-point energy calculations were performed with the previously mentioned methods on these optimised geometries. In addition to these methods, ωB97M-V/def2-TZVPP(D) and ωB97M-V/6-311(+)G** computations were also performed; here we used diffuse basis functions on all oxygen and phosphorus atoms as studies have indicated that diffuse functions are useful to describe interactions with anions.^79,80 Computations of binding energies used the SMD⁸¹ model to include the solvent effects of water. Note that SMD is formally a more advanced model than CPCM as it does not only include a PCM contribution but also models the effect of hydrogen bonds,⁸¹ but we did not notice a notable difference in energies. Corrections for free energy contributions were included in electronic energies, from single-point Hessian with the application of biasing potential in Cartesian space to retain the initial geometry of the complex,⁸² with the GFN2-xTB method.

Analysis of geometries was performed with the MATILDA package as described below.⁸³

2.2. Complex geometries

Lanthanide(III) complexes often feature eight or nine coordination sites out of which two sets of four coordinated atoms may arrange as squares in parallel planes in cyclen based complexes or three sets of three atoms as triangles in TACN based complexes. The case of two squares is shown schematically in Fig. 2 (left). The twisting angle between the squares is a crucial structural feature commonly used to explain spectroscopic and magnetic features.^84–87 The structures can adopt a staggered conformation, which is denoted square antiprism (SAP) and corresponds to a twisting angle of 45°. Structural constraints can inhibit the complex from adopting this idealised twist angle leading to a twisted square antiprism (TSAP) structure. Note that in the limit of a twisting angle of 0° one would obtain a fully eclipsed square prism structure; but such structures are not usually observed in practice. In Gd(III) complexes, increasing the twist toward the TSAP geometry weakens the axial coordination environment, leading to faster inner-sphere water exchange and corresponding changes in relaxivity.⁸⁸ In Eu(III) complexes, the SAP–TSAP twist alters the crystal-field symmetry, producing characteristic modifications in Stark splitting and the overall spectral shape.⁸⁹ We have selected four Eu(III) complexes for analysing the twist in their geometries which are [Eu.1],⁴⁷ [Eu.L2],⁴⁸ [Eu.L3],⁴⁹ and [Eu.L5].⁴⁹ Geometries were obtained from the Cambridge Crystallographic Data Centre (CCDC), deposition numbers: 1965719, 982687, 1849022 and 1849026, respectively.

Fig. 2 Schematic depiction of twisting angles in lanthanide complexes considering square (left) and triangular (right) coordination geometries. The pink sphere is the auxiliary atom, and the green sphere is the central lanthanide atom. The dashed line represents the angle being measured.

In practice, there is not a unique definition of the twisting angle as it depends on the relative positions of 9 individual atoms. For SAP/TSAP complexes, we use the construction as outlined in Fig. 2 (left). We first place an auxiliary atom at the center of the bottom square (which is a cyclen macrocycle in the present cases). Subsequently, we compute four individual dihedral angles connecting bottom square - auxiliary - Eu - top square, iterating over the four atoms of bottom and top square. The twisting angle is now defined as the average of those for dihedrals. If the absolute angle is greater than 35° the complex is assigned SAP and if it is lower than 25°, TSAP is assigned. A borderline status is assigned to any value between 25° and 35°.

Whereas the SAP/TSAP discussions, based on square coordination geometries, are more common in the literature, we can also define a similar parameter for complexes with a trigonal structure, such as [Eu.L2], [Eu.L3], and [Eu.L5] as shown in Fig. 1. In this case, we perform the analogous operations using the bottom/middle and bottom/top triangles (Fig. 2, right). Here, we always follow along individual coordination arms of the ligand to match the corners of the different triangles.

To characterise the structures, we also use the rotational constant, following ref. 76. For the present purpose, rotational constants are used as these provide a unique and easily computable descriptor that provides some information about how similar two structures are. The rotational constants are computed in MATILDA using

where h is Planck's constant and I_k is k-th eigenvalue of the moment of inertia tensor. Below, we report the maximal B_k value, denoted B.

Root-mean-square deviations (RMSD) between two molecular structures are computed in MATILDA in connection with a quaternion algorithm for superimposing the structures.⁹⁰ RMSDs are computed as mass-weighted averages of the individual atomic displacements.

2.3. Binding energies

Binding energies were computed according to the following general reaction scheme where an H₂O molecule is displaced by the guest anion G⁻

Eu·L–H₂O + G⁻ ⇄ Eu·L–G⁻ + H₂O (2)

where Eu·L refers to the europium atom along with its polydentate ligand.

We first consider a “raw” binding constant K_b that is derived directly from eqn (2).

here, K_b depends on the concentration of host–guest complex [Eu·L–G⁻], of water [H₂O], of the water-bound complex [Eu·L–H₂O] and of the free guest [G⁻] anion. The standard free binding energy is subsequently given as

Directly measuring these quantities from experiment is not convenient and one usually determines an apparent binding constant. Here, we use the relation

for the apparent binding constant. Compared to eqn (3), the concentration of water [H₂O] is integrated into the binding constant. Furthermore, in the denominator, we include the sum of [G⁻] and [HG], that is, the total amount of analyte present in solution independently of its protonation state. Combining eqn (4) and (5), we can now derive the following relationconnecting the standard free binding energy with the experimentally determined apparent binding constant . Here, the concentration of water is inserted as 55.6 mol L⁻¹, the [H⁺] value is chosen according to experimental pH, and the K_a of the guest is taken from literature. Whereas the right-hand-side of eqn (6) is determined by experiment, we compute the binding free energy on the left-hand-side using quantum chemistry computations in a standard rigid-rotor/harmonic-oscillator model (with rotor cut-off value of 50 cm⁻¹) for the original reaction described in eqn (2).

We noted that the binding energies tended to be strongly overestimated and assigned this fact to the problem that the electrostatic interaction between the positively charged Eu.L complex and the guest G⁻ may be overestimated in our implicit solvent models as solvation or desolvation of the interacting species can have great influence on the binding free energies.^91,92 We attempted to correct for this effect in a way that did not significantly enhance the computational effort by including an ad-hoc correction using an Na⁺ counterion in the equation

Eu·L–H₂O + G–Na ⇄ Eu·L–G⁻ + H₂O–Na⁺ (7)

here, it is assumed that the guest is not actually solvated freely but forms an ion pair with a cation in solution. We believe that such a process would be particularly relevant for multiply charged anions as discussed below. However, we also want to point out that this is only a very rough approximation. We computed the standard free energies from eqn (7) and again compared this to the right-hand-side of eqn (6), as determined experimentally.

3. Results and discussion

3.1. Binding geometry

Obtaining accurate ground state geometry of the lanthanide complexes is the first and highly relevant step for complex design.⁹³ Particularly, the accuracy of the distances between the lanthanide ion and the ligand atoms directly coordinated to it are of high interest. This knowledge is crucial for an accurate assessment of the influence of the chemical environment on the properties of the complex.⁹⁴ In addition, the twisting angle and particularly the differentiation between SAP and TSAP is of high interest.⁹⁵ The structure determination of these paramagnetic complexes is challenging by commonly employed experimental techniques, i.e., single-crystal X-ray and nuclear magnetic resonance spectroscopy.⁴⁵ Therefore, we were interested in finding a computational method that can give accurate geometries of europium(III) complexes with a reasonable computational cost.

We investigate the four complexes [Eu.1],⁴⁷ [Eu.L2],⁴⁸ [Eu.L3],⁴⁹ and [Eu.L5],⁴⁹ as presented in Fig. 1, using various computational methods, as discussed in Section 2.1. In all cases, we start the geometry optimisation from the experimental crystal structure and investigate how much the final converged geometry deviates from this. An analysis of the resulting bond lengths is presented in Fig. 3. For this comparison, one should keep in mind that finite-temperature effects can affect the measured equilibrium bond-lengths r_e compared to the 0 K value r₀ predicted by DFT. Typical deviations between r₀ and r_e are on the order of 0.01 Å.^96,97 Therefore, experimental r_e values should be considered as an upper bound to r₀, meaning that the values predicted by DFT should be slightly smaller.

Fig. 3 Analysis of computed bond lengths between europium(III) and ligand donor atoms for [Eu.1], [Eu.L2], [Eu.L3], and [Eu.L5] for various computational methods: (a) individual errors; (b) mean absolute errors. The x-axis label specifies the, density functional, basis set, central atom, and number of electrons in the ECP.

In the following, we investigate DFT in conjunction with the B3LYP and ωB97M-V functionals, considering that these were used for previous lanthanide studies.^34,37,39,46 In addition, as computationally cheaper options, we apply the composite r²SCAN-3c method, the semi-empirical GFN2-xTB method, and the force field GFN-FF.

Fig. 3a shows the individual errors in bond lengths versus the crystal structure grouped by the type of bond (O–Eu or N–Eu); Fig. 3b shows the mean absolute error (MAE) averaged over all values in panel a. Starting with GFN-FF, the computationally cheapest method investigated here, we find an MAE of 0.101 Å, and as Fig. 3a shows there is a clear trend for overbinding with most of the bond lengths underestimated.

Proceeding to the semi-empirical GFN2-xTB method, we attempted geometry optimisations using either an Y(III) or an Eu(III) atom. The results obtained between both methods are similar but, interestingly, we obtain a slightly better MAE (0.060 Å vs. 0.068 Å) when using Y in the computations. Continuing with r²SCAN-3c, we find very accurate bond lengths with an error of about 0.028 Å (Y) and 0.024 Å (Eu). While the formal replacement of europium(III) with yttrium(III) in the computation does not have a notable effect on the MAE (always with respect to Eu(III) experiments), there is a clear difference in bond lengths as shown in Fig. 3a. Using Y the bond lengths are predominantly underestimated, using Eu they are predominantly overestimated. In line with the comments made above, a slight underestimation might actually correctly reflect finite temperature effects, thus, favouring the Y results.

Proceeding to the widely used B3LYP functional, in conjunction with the D4 dispersion correction, we find that it provides an MAE of 0.040 Å (Y) and 0.037 Å (Eu), which is somewhat increased compared to r²SCAN-3c. D4-B3LYP underestimates O–Eu distances while overestimating N–Eu distances with similar performance for using Y and Eu atoms.

Finally, we investigated the ωB97M-V functional using five different basis/ECP/central atom combinations. The use of ωB97M-V/def2-SVP and def2-ECP (28 electrons in core) with yttrium, again, resulted in an MAE of about 0.038 Å. A slight improvement is obtained when using Eu as the central atom. Increasing the basis set to def2-TZVPP has only a negligible effect. Lastly, ωB97M-V with ECP52MWB (52 electrons in core) gives values of MAE of 0.030 and 0.029 with 6-31G* and 6-311G**, respectively. Interestingly, as opposed to before, the bond lengths are now mostly overestimated.

In addition to bond lengths, we have also determined the twist angles for the same four complexes. The coordination sphere of [Eu.1] is determined by two squares and we investigate the twisting between the top square formed by the three acetate oxygens and one ether oxygen and the bottom square formed by the four macrocyclic nitrogen atoms. In the tricapped trigonal prismatic [Eu.L2], [Eu.L3] and [Eu.L5] complexes, we investigate the three triangles formed by the three oxygens of the acetate groups (top), the three nitrogens of the quinoline (middle) and the three nitrogens of the macrocycle (bottom). The values of the twisting angles are presented in Fig. 4 with the experimental reference values as solid lines and the computed results as solid shapes (circles, triangles, diamonds and squares) connected by dashed lines. The experimental twist for the squares of [Eu.1] is 23° while for [Eu.L2], [Eu.L3] and [Eu.L5] the twist between bottom and middle (ϕ^tri_bm) is slightly above or below 41° and between bottom and top (ϕ^tri_bt) is in a range of 95–101°.

Fig. 4 Twisting angles between (a): bottom/top (bt) triangle (b): bottom/middle (bm) triangle (c): bottom/top (bt) square of the complexes determined using MATILDA (smooth lines are for crystal structures).

Viewing first GFN-FF, we find that its predicted torsion angles are quite far away from the reference values. For [Eu.1] the torsion angle is too high, meaning that GFN-FF does no longer characterise this as clear TSAP. For the other three complexes, both ϕ^tri_bm and ϕ^tri_bt are again significantly higher than the references. GFN2-xTB also shows some clear deviations, especially concerning the ϕ^tri_bt in tricapped trigonal complexes. By contrast, all DFT methods predict ϕ^tri_bt very accurately for [Eu.L2], and only show a slight overestimation for [Eu.L3] and [Eu.L5]. For [Eu.1], the DFT results are always very close to the reference value. Interestingly, for ϕ^tri_bt in the tricapped complexes [Eu.L2/3/5] we find a zigzag pattern where twisting angles using Y as central atom are always higher than when using Eu. This trend is very pronounced for GFN2-xTB but also clearly visible for the DFT methods. Finally, we note that one should not overinterpret the results from Fig. 4. The experimental twisting angles were measured in the solid state at finite temperature whereas our computations represent the solution phase and formally 0 K. As such one cannot expect perfect agreement. However, the clear conclusion remains that the GFN-FF and GFN2-xTB results are of a somewhat lower quality whereas all DFT methods perform satisfactorily.

Summarising this section, we note that the GFN-FF and GFN2-xTB methods provide computationally inexpensive means for getting approximate geometries and that GFN2-xTB performs slightly better when using a Y atom compared to an Eu atom. Improved results are obtained using DFT and within our dataset the most accurate results are obtained for r²SCAN-3c with almost equal performance for using Y and and Eu. We will, therefore, use r²SCAN-3c as the workhorse for the remaining part of this study. Noting the almost equal performance but significant numerical advantages, we will use a central Y atom rather than an Eu atom in the geometry optimisations.

3.2. Conformational analysis

Whereas the previous section used a somewhat artificial starting point, i.e., testing how well structures are reproduced given the crystal structure as input, we now move to a more realistic setting. Within the following, we use no experimental input for our primary investigation. In a first step, we create a hand-drawn 3D structure for the four complexes of interest ([Eu.1], [Eu.L2], [Eu.L3] and [Eu.L5]) and optimise it at the r²SCAN-3c level. Subsequently, we perform a CREST/GFN2-xTB conformer search based on this structure. Finally, these structures are re-optimised using r²SCAN-3c. We are now interested in evaluating three interrelated questions: (i) is a hand-drawn structure generally good enough to represent the crystal structure, (ii) are improved results obtained using a conformer search algorithm, and (iii) is it reasonable at all to represent such complexes using a single structure?

The results of this investigation are presented in Fig. 5. Here we plot the relative energy against the rotational constant (used as a generic parameter to indicate by how much two structures differ). Within Fig. 5 the hand-drawn and crystal structures are shown as a red cross and blue plus sign, respectively. Note that both of these structures are optimised using r²SCAN-3c (using hand-drawn and crystal structure as a starting point). The CREST ensembles using r²SCAN-3c and GFN2-xTB optimisations are shown using green circles and grey diamonds respectively.

Fig. 5 Analysis of conformational flexibility for four selected europium(III) complexes: left – energy (relative to crystal structure) plotted against rotational constant; right – 3D crystal structures for [Eu.1] and [Eu.L3].

Starting with [Eu.1] (shown on the upper left in Fig. 5), we find that that this complex possesses quite a dense distribution of conformers. Indeed, within an energy window of 42 kJ mol⁻¹, a total of 357 structures were found using CREST/GFN2-xTB (shown as grey diamonds). 50 of these structures were further optimised using r²SCAN-3c, as shown using green circles. For [Eu.1] we find that the hand-drawn structure lies very closely with the crystal structure and that both are embedded within the ensemble of CREST structures. CREST produces some low-energy structures with the lowest one being 2.9 kJ mol⁻¹ below the crystal structure when computed at the r²SCAN-3c level. This finding should be understood in the sense that in solution a range of structures are energetically accessible and that the crystal structure is similar in energy to one of the lowest energy conformers but structurally quite different.

Moving to [Eu.L2], we again find a very dense distribution of conformers with 382 structures generated below 42 kJ mol⁻¹. Reoptimisation using r²SCAN-3c leads to a dense set of structures clustered around the crystal structure with only one structure higher up near 40 kJ mol⁻¹. In this case, it is also noteworthy that the lowest energy conformer using r²SCAN-3c corresponds to the 19th conformer from GFN2-xTB highlighting the importance of investigating larger numbers of conformers. Strikingly, our hand-drawn structure lying at 58.4 kJ mol⁻¹ is completely removed from any of the other structures. This finding illustrates the danger of directly using hand-drawn structures without further scrutiny which would, in this case, yield completely unphysical energies. A more detailed analysis of [Eu.L2] is presented in Fig. 6 illustrating the central part of the complex as displayed for the crystal structure, the lowest energy conformer and the hand-drawn structure. The figure illustrates that in terms of, both, energy and RMSD the lowest energy CREST conformer is quite similar to the crystal structure while the hand-drawn structure deviates strongly. It is difficult to spot the underlying structural dissimilarity and only on close inspection one finds a different arrangement of the methyl groups attached to the phosphorus atoms. For the crystal structure and CREST structure, one finds that all three methyl groups are pointing upwards yielding a roughly C₃ symmetric arrangement. By contrast, for the hand-drawn structure, only one methyl group points up and the other two point sideways. The conclusion from this analysis is that subtle changes in conformation and stereochemistry can have a massive impact on the energetics. It is not a priori clear what arrangement is favourable and an automated conformer search algorithm can be greatly beneficial.

Fig. 6 Central part of the [Eu.L2] complex displayed for its optimised geometry from crystal structure, lowest energy CREST conformer and hand-drawn structure. (Energy is in kJ mol⁻¹ and RMSD in Å).

Resuming with [Eu.L3] in Fig. 5, we find that now a significantly smaller set of low energy conformers is generated using CREST/GFN2-xTB (only 33 below 42 kJ mol⁻¹). The crystal structure coincides with the lowest energy structure from the CREST/r²SCAN-3c optimisation. The hand-drawn structure is slightly higher in energy at 3.4 kJ mol⁻¹. Finally, the simplest case is found for [Eu.L5]. In this case only four low-energy conformers were generated. Out of these conformers, one exactly coincided with both the crystal and hand-drawn structures.

In summary, Fig. 5 highlights the challenges in working with flexible systems like europium(III) complexes. Whereas in some cases, such as [Eu.L5], it is simple to identify the lowest energy conformer, there are other cases such as [Eu.1] and [Eu.L2] where a large number of low-energy conformers are involved. Even more, subtle differences in the way structures are generated can have dramatic effects on the outcomes producing even qualitatively incorrect results as exemplified for [Eu.L2]. Specifically, we have highlighted the challenges of directly using a hand-drawn structure. To obtain more robust results, we suggest to always start with a crystal structure if available. If no crystal structure is available we would recommend a systematic conformer search.

3.3. Case study: anion binding

Having established a protocol for obtaining reliable geometries of europium(III) complexes, we now want to present a more extended case study of biological anions binding to [Eu.pB(OH)₂]⁺ and [Eu.mB(OH)₂]⁺. [Eu.pB(OH)₂]⁺ is derived from [Eu.1] by removing one acetate donor to open a binding site at the Eu(III) center for phosphoanions, with a para-positioned boronic acid designed to engage the diol of adenosine monophopshate (AMP). Its isomer, [Eu.mB(OH)₂]⁺, has the boronic acid at the meta position to probe the optimal geometry for simultaneous phosphate and diol binding. Target anions include HPO₄²⁻, AMP, and ADP, which are biologically important, allowing assessment of how ligand geometry modulates recognition of these species. We considered monodentate binding for all these anions as in previous studies the monodentate binding of inorganic phosphate is assumed to be the preferred mode based on the shape of the emission spectra,⁴⁶ which are sensitive to coordination geometry. More broadly, there is no crystallographic evidence supporting bidentate binding of inorganic phosphate or monophosphate esters to Eu(III) centres. Consistent with this, a recent study⁹⁸ using EXAFS and EPR on a cyclen-based Eu(III) complex showed that adenosine monophosphate binds in a monodentate fashion. Whereas ADP could bind in a bidentate fashion in principle, we have previously argued that monodentate binding is dominant for the complexes studied here.⁴⁶

Theoretical determination of accurate binding free energies for such systems is often difficult due to a number of challenges, such as the charged state of molecules, solvation, pH and basis set superposition error (BSSE).^59,99–101 To study these complexes, we first performed conformer analysis on the water-bound europium(III) complexes, host–guest complexes (as binding of guests to the macrocyclic hosts is known to cause conformational variations¹⁰²) and anions. Then, the lowest energy structure (as determined at the r²SCAN-3c level using a Y(III) central ion) was selected to compute the binding free energy. To perform computations on complexes, we have replaced Eu with Y for reasons mentioned before, and geometries of all structures were optimised with the r²SCAN-3c method. A comparison to experimental binding constants was performed following Section 2.3.

The optimised molecular geometries for the lowest energy conformer from CREST/r²SCAN-3c are displayed in Fig. 7. These show the binding of host europium(III) complexes to water, hydrogen phosphate, AMP and ADP. Here, the binding of AMP and ADP is monodentate (binding through only one phosphate oxygen). The guest anions are shown on the right side of the complex, binding to the central cation via one of their negatively charged oxygen atoms. Selected bond lengths and the SAP twist angles are listed below.

Fig. 7 Structure of the [Eu.pB(OH)₂]⁺ complex bound with (a) a water molecule, (b) HPO₄²⁻, (c) AMP, (d) ADP and [Eu.mB(OH)₂]⁺ with (e) a water molecule, (f) HPO₄²⁻, (g) AMP, (h) ADP. Bond distances of central metal to ligand atoms are represented as N(m) macrocycle, N(q) quinoline, O(a) acetate, O(e) ether, O(w) water, O(p) phosphate.

Starting with the water-bound [Eu.pB(OH)₂]⁺ complex in Fig. 7a, we note that all Eu–N bonds are in the range of 2.5–2.6 Å. There is more variation with the bonds involving oxygen. The Eu–O(ac) bonds involving the charged acetate groups are the shortest (2.3 Å), followed by the bonds to water (2.4 Å). Conversely, the Eu–O(e) bond involving the ether is significantly longer highlighting the weaker interaction with this group.

Upon binding of HPO₄²⁻ the binding geometry around the europium atom stays largely intact but some crucial adjustments to bond lengths occur, as shown in Fig. 7b. Most importantly the Eu–O(p) bond (2.24 Å) involving phosphate is significantly shorter than the analogous Eu–O(w) bond involving water was in the previous example (2.39 Å) reflecting the enhanced charge on this atom. As a consequence the Eu–O(a) bonds are slightly elongated. We also find that one of the macrocyclic nitrogens is now farther out at 2.71 Å while the distance to the ether oxygen is now only 2.81 Å reflecting some further changes in the coordination sphere. Viewing the three-dimensional structure in Fig. 7b, we also note that the phenylboronic acid group is now pointing in a different direction as compared to before highlighting the influence of the conformer search procedure.

Moving on to the AMP and ADP complexes (Fig. 7c and d), we again note that the overall coordination geometry remains intact, but also observe some changes in bond lengths. Especially, we note an increase in the Eu–N(q) and Eu–O(e) bond lengths, indicative of steric crowding caused by the larger size of AMP and ADP compared to water and HPO₄²⁻. Conversely, the distance of the phosphate oxygen is similar to the HPO₄²⁻-bound [Eu.pB(OH)₂]⁺ complex. Interestingly, there is also a change in the twist angles for the lowest energy conformers from −37° to 25° changing the formal classification from SAP to TSAP. Reviewing finally the phenylboronic acid group in all four complexes in Fig. 7, we note that it is arranged in quite divergent directions in all four cases highlighting again the comformational flexibility of this complex.

Next, we proceed to the [Eu.mB(OH)₂]⁺ complex, as shown in the same Fig. 7. We generally observe similar trends as in the case of the para isomer but some notable differences are present. Viewing the water bound complex, we note that the phenylboronic acid group is now pointing away, apparently to form a hydrogen bond with the associated acetate group. In addition, we find that the Eu–O(e) bond is significantly shortened for meta compared to para (2.83 Å vs. 3.08 Å). Comparing the phosphate bound [Eu.mB(OH)₂]⁺ complex in Fig. 7f with its para counterpart, we find a different arrangement of the phenylboronic acid group, which is now pointing toward the phosphate group to form a hydrogen bond. This goes along with a quite dramatic elongation of the formal Eu–O(e) bond, which is now 3.15 Å. Proceeding to the AMP bound [Eu.mB(OH)₂]⁺ complex in Fig. 7g, we note that this stands out due to its particularly short Eu–O(p) bond of only 2.217 Å. Viewing finally the ADP bound [Eu.mB(OH)₂]⁺ complex in Fig. 7h, we find that the coordination geometry is quite similar to the para analogue but that the conformer adopted by ADP is quite different with a more folded arrangement for meta and a more linear arrangement for para.

To summarise the results from Fig. 7, we note that binding of the analytes leaves the coordination geometry around the europium atom largely intact, with only some of the bond lengths being altered. At the same time, the conformational flexibility of the ligands is illustrated, particularly through the phenylboronic acid group, which is moved in various directions possibly to optimise hydrogen bond interactions. The variability in the structures illustrates the importance of using a conformer search algorithm when computing the geometries of host–guest complexes. Finally, we note that the lowest energy conformers as shown here, cannot be assumed to be the only relevant conformers considering the existence of a number of similar energy conformers. Nonetheless, we believe that the presented results provide a reasonable representative binding geometry.

Having discussed the geometries, we were next interested whether it is possible to obtain accurate binding free energies. Experimentally determined values [following eqn (6)] are presented in the first column of Table 1. There is not much variation in these values, all lie around 35 kJ mol⁻¹, with a slightly higher preference of [Eu.pB(OH)₂]⁺ for ADP and somewhat less for HPO₄²⁻.

Table 1 Evaluation of various protocols for computing binding free energies for [Eu.pB(OH)₂]⁺ (para) and [Eu.mB(OH)₂]⁺ (meta) bound to phosphate, AMP, and ADP. Results are always based on r²SCAN-3c geometries and energies

Complex	Anion	Exp.^a	Drawn^b	CREST^c	Na^+^d
a Experimental binding energy computed via eqn (6) using binding constants from ref. 46.b Hand-drawn and optimised without conformer search.c CREST conformer search carried out for each component.d Estimation of the effect of Na^+ counterions via eqn (7) (including CREST conformer search)
Para	[HPO4]^−2	−31.1	−79.4	−79.2	−34.3
	AMP	−35.1	−81.0	−67.3	−26.7
	ADP	−37.1	−81.9	−67.3	−15.0
Meta	[HPO4]^−2	−31.2	−78.2	−91.7	−46.9
	AMP	−34.3	−84.2	−79.6	−38.9
	ADP	−34.2	−49.5	−71.8	−18.7

---

We evaluated three different strategies for the computation of the binding free energies. In the first approach, we optimised the hand-drawn europium(III) complexes directly using r²SCAN-3c. The second method includes performing conformer search on complexes using CREST and then optimising the resulting conformers at the r²SCAN-3c level. Subsequently the lowest energy conformer (at r²SCAN-3c level) is used for calculating the binding free energy. The results from these two approaches, as presented in Table 1 in the columns “Drawn” and “CREST”, are clearly not satisfactory severely overestimating binding energies producing values about twice as large as they should be. This overbinding could be attributed to the fact that charged ions are screened in solution, due to the effect of the solvent molecules, and interact with counterions. To account for this effect, we apply an ad hoc correction by considering the binding of an Na⁺ counterion to the analyte anion in the computation following eqn (7). Interestingly, this correction does indeed lead to the values being close to experimental results. However, the qualitative trends and relative numbers are still not satisfactory.

Noting that none of the above methods produce satisfactory results, we now investigate the effect of different computational methods. Fig. 8 presents energies with the inclusion of the Na⁺ counter ion, i.e., via eqn (7). We investigate the same set of methods as above but do not include the small-ECP Eu computations as these produced quite erratic results due to numerical instabilities related to the ⁶f electron configuration. Conversely, we included two more methods, denoted ωB97M-V/def2-TZVPP(D) and ωB97M-V/6-311(+)G** including diffuse basis functions on all the O- and P-atoms. In Fig. 8, the r²SCAN-3c values represent the results also shown in Table 1. The other results are based on the same geometries using single-point computations with the respective method. We first note that GFN2-xTB produces extremely high binding energies (around 300 kJ mol⁻¹, not shown here). Continuing with D4-B3LYP/def2-SVP and ωB97M-V/def2-SVP, we find that some of the energies are severely overestimated reaching above 70 kJ mol⁻¹ in three cases, each. We can rationalise this overestimation in terms of basis-set superposition error (BSSE) and, indeed, when using the larger def2-TZVPP basis set in conjunction with ωB97M-V, binding energies are strongly reduced, more in line with r²SCAN-3c. Increasing the basis set further and adding diffuse functions (def2-TZVPP(D)) has only a slight further influence. Interestingly, when using ωB97M-V/6-31G* in conjunction with ECP52MWB, we do not observe strong overbinding despite the fact that 6-31G* is of polarised double-ζ quality analogous to def2-SVP. Increasing the basis set to 6-311G** now somewhat increases the binding, which is the opposite trend to what would be expected if the BSSE between Eu on phosphate is the main influence and what was observed with the def2 basis sets. The addition of diffuse functions in the basis set (6-311(+)G**) somewhat lowers binding energies putting them on the right order of magnitude.

Fig. 8 Binding free energies of anions with europium(III) complexes (using Na⁺ counterions) determined with different computational methods.

Finally, we briefly discuss the results obtained when not using a counterion, as shown in Fig. S2. In this case, the binding energies are even more severely overestimated, as also seen in Table 1. Interestingly, GFN2-xTB underestimates the binding energies while all other methods overestimate them.

Reviewing the results presented here, we are left with the conclusion that while geometries can be meaningfully reproduced, none of the presented computational approaches yields satisfactory results for binding energies. We discussed two potential sources of error, (i) use of the wrong conformer and (ii) BSSE. We showed that these can be avoided via (i) an appropriate conformer search algorithm and (ii) larger basis sets and/or appropriate BSSE corrections (such as the gCP correction in r²SCAN-3c). However, even after including these, the results are neither of the correct order of magnitude nor reflect the experimentally observed trends. We hypothesise that the missing contribution lies in the solvation of the multiply charged ions investigated here. Apparently, the relevant solvation effects are not fully described by the SMD continuum model employed. The ad hoc addition of an Na⁺ ion moves the results in the right direction, but this approximation is too coarse to produce quantitative results. A more realistic solvation environment might be obtained via dynamics with explicit solvent or free energy perturbation approaches. However, such approaches are out of the scope of this work.

4. Conclusions

Within this work, we investigated various computational approaches for studying europium(III) complexes with a focus on their potential as anion sensors in aqueous media. Aside from investigating density functionals and basis sets, we also put an emphasis on the most practical description of the f-electrons within the europium atom, described either explicitly when using a small-core ECP, implicitly within a large-core ECP, or avoided altogether when substituting the Eu atom for Y. To analyse geometries of the complexes, we have developed tools for computing bond lengths, twist angles and rotational constants and made them available via the MATILDA package. Within this work three interrelated questions were investigated: (i) how well are binding geometries reproduced when starting from the crystal structure, (ii) what is the influence of conformational flexibility in generating appropriate structures, and (iii) how well binding energies in host–guest complexes can be predicted.

Starting with geometries, we find that the composite DFT-based r²SCAN-3c method comes out as the optimal approach being both accurate and computationally efficient. We also find that substituting Eu for Y does not noticeably affect the results while providing significant advantages in terms of numerical stability and computational efficiency. To reduce computational cost, the GFN2-xTB and GFN-FF methods were shown to be effective albeit with the downside of notably increased errors in bond lengths and twisting angles.

Progressing to the influence of conformational flexibility, we illustrated that some of the investigated complexes possess a large number of low-energy conformers highlighting the importance of sampling these appropriately to get quantitative results. Even more, we noticed in one instance that simply using a hand-drawn structure led to qualitatively wrong results and an energy overestimated by more than 50 kJ mol⁻¹. We found CREST to be a suitable conformer search tool noting, however, that the generated structures had to be reoptimised at a higher level for an accurate estimation of the energetic ordering and we used r²SCAN-3c for this purpose.

We concluded this work with a case study on anion binding, considering phosphate, AMP, and ADP, using two previously synthesised probe complexes. First, we studied the host–guest structures obtained with the different anions illustrating changes in bond distances, twist angles and more wide ranging changes in conformer structure upon binding, underscoring the structural variability of these complexes. Subsequently, we investigated the possibility of computing binding energies. Despite investigating a variety of protocols and computational methods, we were not able to pinpoint any approach that yielded satisfactory results. We have highlighted the importance of basis set superposition error and conformational flexibility both of which, when neglected introduce strong errors. However, even after taking both into account appropriately, the binding energies are still strongly overestimated. We attribute the remaining error to challenges of appropriately describing charged anions in solution, especially as far as charge screening and the effect of counterions is concerned. Describing such effects would require dynamics simulations or free energy perturbation approaches, both of which are beyond the scope of this work. We are, thus, left with the conclusion that within standard quantum chemistry approaches the description of binding energies is out of reach in these systems.

To summarise, we conclude that r²SCAN-3c proves to be an effective method for obtaining accurate geometries of europium(III) complexes. In connection with CREST, also conformationally flexible systems can be treated well. Accurate binding energies are unfortunately out of reach with the presented methodologies. As such we believe that the presented methods are beneficial for rationalising and illustrating experimental results. However, at this point it is not possible to predict binding affinities, e.g., to develop anion selective probes. Nonetheless, we hope that the results presented here serve as a good representation of the state of the art and will serve as an inspiration to develop more accurate methods.

Author contributions

Conceptualization – F. P., U. R., S. J. B.; investigation – U. R.; methodology – F. P., J-M. M., U. R.; software – F. P., U. R., J-M. M.; supervision – F. P., S. J. B.; writing – original draft – U. R., F. P.; writing – review & editing – F. P.; J-M. M.; S. J. B.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

Data for this article (input, output and coordinate files for all computations performed) are available at Loughborough University's institutional repository at https://doi.org/10.17028/rd.lboro.31281982.

Supplementary information (SI): detailed derivation for binding energy (Section S1), structures of anions (Fig. S1), complexation of sodium cation with guests (Fig. S2), and binding energy without sodium cations (Fig. S3). See DOI: https://doi.org/10.1039/d6cp00457a.

Acknowledgements

The authors acknowledge the 'lovelace' high performance computing service at Loughborough University.

References

1. J.-C. G. Bünzli, J. Lumin., 2016, 170, 866–878.
2. S. V. Eliseeva and J.-C. G. Bünzli, Chem. Soc. Rev., 2010, 39, 189–227.
3. D. Parker, J. D. Fradgley and K.-L. Wong, Chem. Soc. Rev., 2021, 50, 8193–8213.
4. E. G. Moore, A. P. Samuel and K. N. Raymond, Acc. Chem. Res., 2009, 42, 542–552.
5. A. K. Singh, Coord. Chem. Rev., 2022, 455, 214365.
6. A. Dalal, K. Nehra, A. Hooda, D. Singh, P. Kumar, S. Kumar, R. S. Malik and B. Rathi, Inorg. Chim. Acta, 2023, 550, 121406.
7. J. Kido and Y. Okamoto, Chem. Rev., 2002, 102, 2357–2368.
8. J. Fang, H. You, J. Chen, J. Lin and D. Ma, Inorg. Chem., 2006, 45, 3701–3704.
9. R. Jankowski, M. Wyczesany and S. Chorazy, Chem. Commun., 2023, 59, 5961–5986.
10. S. Roy, A. Chakraborty and T. K. Maji, Coord. Chem. Rev., 2014, 273, 139–164.
11. K. Binnemans, Chem. Rev., 2007, 107, 2592–2614.
12. G. Tao, Y. Huang, J. A. Boatz and J. M. Shreeve, Chem. – Eur. J., 2008, 14, 11167–11173.
13. S.-F. Tang and A.-V. Mudring, Inorg. Chem., 2019, 58, 11569–11578.
14. S. E. Bodman and S. J. Butler, Chem. Sci., 2021, 12, 2716–2734.
15. S. J. Butler and D. Parker, Chem. Soc. Rev., 2013, 42, 1652–1666.
16. N. B. Lima, S. M. Goncalves, S. A. Junior and A. M. Simas, Sci. Rep., 2013, 3, 2395.
17. T. L. Martinon, A. Sinha, A. W. Schlimgen, A. D. Persell, T. E. Garner, K. L. Peterson, J. I. Siepmann and V. C. Pierre, JACS Au, 2025, 5, 5512–5523.
18. Y. Zheng, C. Tan, G. P. Drummen and Q. Wang, Spectrochim. Acta, Part A, 2012, 96, 387–394.
19. C. M. G. Dos Santos and T. Gunnlaugsson, Dalton Trans., 2009, 4712–4721.
20. Z. Liu, J. Bai, C. Liang, Y. Wang, H. Nie, X. Liu and H. Yan, Biosens. Bioelectron., 2022, 203, 114055.
21. E. Al-Enezi, A. Vakurov, A. Eades, M. Ding, G. Jose, S. Saha and P. Millner, Sensors, 2021, 21, 831.
22. H. Li, S. L.-W. Ng, D. J. Black, W. Han, R. Pal and D. Parker, Chem. Sci., 2025, 16, 14544–14552.
23. D.-W. Huang, C.-G. Niu, G.-M. Zeng and M. Ruan, Biosens. Bioelectron., 2011, 29, 178–183.
24. L. Bridou, L. G. Nielsen and T. J. Sorensen, J. Rare Earths, 2020, 38, 498–505.
25. D. Wang, J. Zhang, Q. Lin, L. Fu, H. Zhang and B. Yang, J. Mater. Chem., 2003, 13, 2279–2284.
26. M. R. Felício, P. D. Vaz, C. D. Nunes and M. M. Nolasco, New J. Chem., 2022, 46, 21311–21323.
27. H. A. Azab, G. M. Khairy, N. Abd El-Ghany and M. A. Ahmed, J. Photochem. Photobiol., A, 2019, 374, 1–9.
28. S. Pandya, J. Yu and D. Parker, Dalton Trans., 2006, 2757–2766.
29. P. Coppo, M. Duati, V. N. Kozhevnikov, J. W. Hofstraat and L. De Cola, Angew. Chem., Int. Ed., 2005, 44, 1806–1810.
30. R. Zhang and J. Yuan, Acc. Chem. Res., 2020, 53, 1316–1329.
31. M. Starck, J. D. Fradgley, R. Pal, J. M. Zwier, L. Lamarque and D. Parker, Chem. – Eur. J., 2021, 27, 766–777.
32. H. Ma, K. Chen, B. Song, Z. Tang, Y. Huang, T. Zhang, H. Wang, W. Sun and J. Yuan, Biosens. Bioelectron., 2020, 168, 112560.
33. G. De Sa, O. Malta, C. de Mello Donegá, A. Simas, R. Longo, P. Santa-Cruz and E. da Silva Jr, Coord. Chem. Rev., 2000, 196, 165–195.
34. M. L. Shipton, F. A. Jamion, S. Wheeler, A. M. Riley, F. Plasser, B. V. L. Potter and S. J. Butler, Chem. Sci., 2023, 14, 4979–4985.
35. S. E. Bodman, P. Stachelek, U. Rehman, F. Plasser, R. Pal and S. J. Butler, Chem. Sci., 2025, 16, 5602–5612.
36. S. E. Bodman, C. Breen, S. Kirkland, S. Wheeler, E. Robertson, F. Plasser and S. J. Butler, Chem. Sci., 2022, 13, 3386–3394.
37. S. H. Hewitt, G. MacEy, R. Mailhot, M. R. Elsegood, F. Duarte, A. M. Kenwright and S. J. Butler, Chem. Sci., 2020, 11, 3619–3628.
38. I. Georgieva, T. Zahariev, A. J. Aquino, N. Trendafilova and H. Lischka, Spectrochim. Acta, Part A, 2020, 240, 118591.
39. S. Shuvaev, M. A. Fox and D. Parker, Angew. Chem., Int. Ed., 2018, 130, 7610–7614.
40. A. V. de Andrade, N. B. da Costa Jr, A. M. Simas and G. F. de Sá, Chem. Phys. Lett., 1994, 227, 349–353.
41. R. O. Freire, G. B. Rocha and A. M. Simas, Inorg. Chem., 2005, 44, 3299–3310.
42. M. A. Filho, J. D. L. Dutra, H. L. Cavalcanti, G. B. Rocha, A. M. Simas and R. O. Freire, J. Chem. Theory Comput., 2014, 10, 3031–3037.
43. M. Bursch, A. Hansen and S. Grimme, Inorg. Chem., 2017, 56, 12485–12491.
44. T. Rose, M. Bursch, J. M. Mewes and S. Grimme, Inorg. Chem., 2024, 63, 19364–19374.
45. Y. Okayasu and J. Yuasa, J. Phys. Chem. Lett., 2021, 12, 6867–6874.
46. S. E. Bodman, C. Breen, F. Plasser and S. J. Butler, Org. Chem. Front., 2022, 9, 5494–5504.
47. C. Breen, R. Pal, M. R. Elsegood, S. J. Teat, F. Iza, K. Wende, B. R. Buckley and S. J. Butler, Chem. Sci., 2020, 11, 3164–3170.
48. S. J. Butler, M. Delbianco, N. H. Evans, A. T. Frawley, R. Pal, D. Parker, R. S. Puckrin and D. S. Yufit, Dalton Trans., 2014, 43, 5721–5730.
49. K. Mason, A. C. Harnden, C. W. Patrick, A. W. Poh, A. S. Batsanov, E. A. Suturina, M. Vonci, E. J. McInnes, N. F. Chilton and D. Parker, Chem. Commun., 2018, 54, 8486–8489.
50. S. Spicher and S. Grimme, Angew. Chem., Int. Ed., 2020, 59, 15665–15673.
51. C. Bannwarth, S. Ehlert and S. Grimme, J. Chem. Theory Comput., 2019, 15, 1652–1671.
52. N. Mardirossian and M. Head-Gordon, Phys. Chem. Chem. Phys., 2014, 16, 9904–9924.
53. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
54. E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017, 147, 034112.
55. E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2019, 150, 154122.
56. S. Grimme, A. Hansen, S. Ehlert and J.-M. Mewes, J. Chem. Phys., 2021, 154, 064103.
57. R. Gulde, P. Pollak and F. Weigend, J. Chem. Theory Comput., 2012, 8, 4062–4068.
58. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222.
59. H. Kruse and S. Grimme, J. Chem. Phys., 2012, 136, 154101.
60. M. Bursch, J.-M. Mewes, A. Hansen and S. Grimme, Angew. Chem., Int. Ed., 2022, 61, e202205735.
61. M. Dolg, H. Stoll and H. Preuss, J. Chem. Phys., 1989, 90, 1730–1734.
62. M. Soulie, F. Latzko, E. Bourrier, V. Placide, S. J. Butler, R. Pal, J. W. Walton, P. L. Baldeck, B. Le Guennic and C. Andraud, et al., Chem. – Eur. J., 2014, 20, 8636–8646.
63. K. Sénéchal-David, A. Hemeryck, N. Tancrez, L. Toupet, J. G. Williams, I. Ledoux, J. Zyss, A. Boucekkine, J.-P. Guégan and H. Le Bozec, et al., J. Am. Chem. Soc., 2006, 128, 12243–12255.
64. T. T. Da Cunha, J. Jung, M.-E. Boulon, G. Campo, F. Pointillart, C. L. Pereira, B. Le Guennic, O. Cador, K. Bernot and F. Pineider, et al., J. Am. Chem. Soc., 2013, 135, 16332–16335.
65. F. Pointillart, B. L. Guennic, O. Maury, S. Golhen, O. Cador and L. Ouahab, Inorg. Chem., 2013, 52, 1398–1408.
66. R. O. Freire and A. M. Simas, J. Chem. Theory Comput., 2010, 6, 2019–2023.
67. M. Dolg, H. Stoll, A. Savin and H. Preuss, Theor. Chim. Acta, 1989, 75, 173–194.
68. E. Epifanovsky, A. T. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons and A. L. Dempwolff, et al., J. Chem. Phys., 2021, 155, 084801.
69. F. Neese, J. Comput. Chem., 2003, 24, 1740–1747.
70. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12, e1606.
71. F. Neese, J. Comput. Chem., 2023, 44, 381–396.
72. S. Lehtola, C. Steigemann, M. J. Oliveira and M. A. Marques, SoftwareX, 2018, 7, 1–5.
73. V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995–2001.
74. S. Grimme, C. Bannwarth and P. Shushkov, J. Chem. Theory Comput., 2017, 13, 1989–2009.
75. S. Ehlert, M. Stahn, S. Spicher and S. Grimme, J. Chem. Theory Comput., 2021, 17, 4250–4261.
76. P. Pracht, F. Bohle and S. Grimme, Phys. Chem. Chem. Phys., 2020, 22, 7169–7192.
77. S. Finta, A. V. Kalikadien and E. A. Pidko, J. Chem. Theory Comput., 2025, 21, 5334–5345.
78. A. A. Otlyotov, T. P. Rozov, A. D. Moshchenkov and Y. Minenkov, Organometallics, 2024, 43, 2232–2242.
79. J. Simons, J. Phys. Chem. A, 2008, 112, 6401–6511.
80. J. Simons, J. Phys. Chem. A, 2023, 127, 3940–3957.
81. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
82. S. Spicher and S. Grimme, J. Chem. Theory Comput., 2021, 17, 1701–1714.
83. S. Ghosh, U. Rehman and F. Plasser, MATILDA: scripts for manipulation, comparison and conversion of molecular structures, 2025; available from: https://github.com/felixplasser/matilda.
84. D. Parker, E. A. Suturina, I. Kuprov and N. F. Chilton, Acc. Chem. Res., 2020, 53, 1520–1534.
85. A. C. Harnden, E. A. Suturina, A. S. Batsanov, P. K. Senanayake, M. A. Fox, K. Mason, M. Vonci, E. J. McInnes, N. F. Chilton and D. Parker, Angew. Chem., 2019, 131, 10396–10400.
86. M. Woods, S. Aime, M. Botta, J. A. Howard, J. M. Moloney, M. Navet, D. Parker, M. Port and O. Rousseaux, J. Am. Chem. Soc., 2000, 122, 9781–9792.
87. L. G. Nielsen and T. J. Sørensen, Inorg. Chem., 2019, 59, 94–105.
88. S. Aime, A. Barge, J. I. Bruce, M. Botta, J. A. Howard, J. M. Moloney, D. Parker, A. S. De Sousa and M. Woods, J. Am. Chem. Soc., 1999, 121, 5762–5771.
89. O. A. Blackburn, R. M. Edkins, S. Faulkner, A. M. Kenwright, D. Parker, N. J. Rogers and S. Shuvaev, Dalton Trans., 2016, 45, 6782–6800.
90. C. F. Karney, J. Mol. Graph. Model., 2007, 25, 595–604.
91. J. H. Jensen, Phys. Chem. Chem. Phys., 2015, 17, 12441–12451.
92. S. C. Patrick, P. D. Beer and J. J. Davis, Nat. Rev. Chem., 2024, 8, 256–276.
93. J. D. L. Dutra, I. F. Gimenez, N. B. D. C. Junior and R. O. Freire, J. Photochem. Photobiol., A, 2011, 217, 389–394.
94. R. O. Freire, G. B. Rocha and A. M. Simas, J. Braz. Chem. Soc., 2009, 20, 1638–1645.
95. S. E. Bodman and S. J. Butler, Chem. – Eur. J., 2025, e02251.
96. J.-M. Mewes, A. Hansen and S. Grimme, Angew. Chem., Int. Ed., 2021, 60, 13144–13149.
97. H. Ohmagari, N. Marets, J. Kamata, M. Yoneyama, T. Miyauchi, Y. Takahashi, Y. Yamamoto, Y. Ogihara, D. Saito and K. Goto, et al., Front. Chem., 2022, 10, 1047960.
98. H. K. Pyle, M. Judd, A. Barancewicz, A. J. Mayer, N. Cox, S. A. Kondrat and S. J. Butler, Inorg. Chem., 2024, 63, 20726–20736.
99. F. Jensen, J. Chem. Theory Comput., 2024, 20, 767–774.
100. R. C. Lochan and M. Head-Gordon, Phys. Chem. Chem. Phys., 2006, 8, 1357–1370.
101. Y. Inada and H. Orita, J. Comput. Chem., 2008, 29, 225–232.
102. J. W. Steed, Coord. Chem. Rev., 2001, 215, 171–221.

---