Computational Investigation of Adsorptive Removal of Pb 2 1 from Water by the UiO-66 Metal–Organic Framework: Comparison of Adsorption Sites on Defects and Functionalised Linkers

Adsorption using metal–organic frameworks (MOFs) such as UiO-66 has shown great promise in remediating water sources contaminated with toxic heavy metals such as Pb 2 þ , but detailed information about the adsorption process remains limited. In this article, we gained mechanistic insights into Pb 2 þ adsorption using both functionalised and defective UiO-66 by performing density functional theory calculations using cluster models. Our benchmarked approach led to a computational model of solvated Pb 2 þ (a hemidirected Pb(H 2 O) 6 2 þ complex) fully consistent with experimental reports. The analysis of Pb 2 þ adsorption using functionalised UiO-66 determined that factors such as electrostatic attraction, chelation, and limited constraints on the Pb 2 þ coordination geometry lead to enhanced binding affinity. For these reasons, UiO-66-COO – was identified as the most promising functionalised MOF, consistent with experimental literature. We additionally explored a novel aspect of Pb 2 þ adsorption by UiO-66: the role of missing linker defects that often characterise this MOF. We found that the defects expected to form in an aqueous environment can act as excellent adsorption sites for Pb 2 þ and the preferred adsorption geometry is again determined by electrostatic attraction, chelation, and constraints on the Pb 2 þ coordination geometry. Overall, we conclude that functional groups and defect sites can both contribute to Pb 2 þ adsorption and our study provides crucial design principles for improving the UiO-66 MOF performance in toxic Pb 2 þ removal from water. defects.


Introduction
As a bivalent cation, lead has diverse and extremely detrimental impacts on the body in all quantities. [1,2] The high toxicity of Pb 2þ is attributed to this metal's mimicry of calcium in the body and subsequent ability to disrupt enzymes, resulting in protein malfunction that interferes significantly with the processes of the central nervous system. Water poisoning with lead remains a modern and global issue. [3] Dangerous elevated blood lead levels that lead to developmental delays in children are found across the world, including countries such as Australia, Canada, and the US. [4][5][6][7] Identifying an efficient way to remove this toxic heavy metal from water is imperative. This urgent need is heightened by the fact that water pollution is rising with increasing wastewater discharge from the agriculture industry, urbanisation, and other practices. [8] Water purification techniques include chemical, physical, and biological means such as filtration, electrodialysis, ion exchange, adsorption, and membrane bioreactors. [8] Adsorption techniques represent a particularly promising option as they are easy to operate, typically inexpensive, and do not produce as much sludge or other toxic pollutants as other methods. [9] As a consequence, adsorption today is an effective, simple process that has become one of the most common ways to purify water sources. [8] The efficacy and selectivity of the absorbent materials can be manipulated by changing pore size, surface area, and functional groups. [10] As an example, adsorption of Pb 2þ ions from aqueous solutions has been successfully achieved using carbon nanotubes as well as other carbonaceous materials functionalised with oxygen and nitrogen. [11,12] Metal-organic frameworks (MOFs) are one class of materials that has proven success in water purification applications, often demonstrating adsorption capacities for toxic heavy metals greater than those of conventional adsorbents [13,14] such as modified clays (78.74 mg g À1 reported Pb 2þ adsorption capacity for tripolyphosphate-modified kaolinite clay [15] ) and activated carbon (Pb 2þ adsorption capacity of 17.5 mg g À1 [16] ). MOFs are characterised by unprecedented porosity and surface area, explaining their excellent performance in separations and water purification applications. Significant research effort has been dedicated to a subclass of MOFs -UiO-66 and its functionalised derivatives -for the removal of a variety of heavy metals from water; we direct the interested reader to a recent review on this topic. [17] UiO-66 is a highly porous and water-stable MOF composed of Zr 6 (OH) 4 (O) 4 nodes connected by 1,4-benzenedicarboxylic acid (BDC) linkers. [18,19] In addition to the pristine form of UiO-66, stable structures with functionalised BDC linkers as well as defects consisting of missing linkers and nodes have been synthesised. UiO-66 with defects produced by missing linkers requires other groups to saturate the undercoordinated Zr sites and balance charge; these capping ligands can be monocarboxylic acids used in some synthesis procedures [20] or hydroxide anions and water molecules ubiquitously present in aqueous solution. [21] Interestingly, these functionalised and defective MOFs often have adsorption capacities, catalytic performance, and gas separation capabilities superior to that of pristine UiO-66. [20] UiO-66 with functionalised linkers has shown significant aptitude for the adsorption of Pb 2þ from water sources and recent studies suggest that functional groups, pore size, and defects may play a significant role in influencing its adsorption capacity. Wang et al. obtained UiO-66-NH 2 using a microwavepromoted synthesis technique and found a maximum Pb 2þ adsorption capacity of 116.74 mg g À1 . [22] Saleem et al. conducted post-synthetic modification in order to replace the amine functional group on UiO-66-NH 2 with a range of sulfur-and nitrogen-containing groups such as NCO, NCS, and NHCSNHCH 3 . [23] The latter bulky functional group led to a maximum adsorption capacity of 232 mg g À1 . A significant adsorption capacity of 420 mg g À1 by UiO-66-(COOH) 2 (i.e. two -COOH groups present on each BDC linker) at pH 6 was reported by Zhao et al. [24] The authors propose that this outstanding adsorption capacity is due to electrostatic attraction between Pb 2þ and the deprotonated -COOgroups, which are likely to dominate at pH 6. Using the more functionalised UiO-66-(COOH) 4 led to less Pb 2þ adsorption, [24] potentially owing to steric effects and reduced pore size, which then limits the diffusion of ions through the material. This indicates that there is a balance that needs to be reached for maximum Pb 2þ adsorption: increasing the number of functional groups increases the number of binding sites but decreases the pore size, thus limiting access to more binding sites. A separate study by Morcos et al. found UiO-67 (an analogue of UiO-66 with larger apertures and pores) to have increased adsorption capacity relative to UiO-66, highlighting the importance of pore size as a factor for Pb 2þ adsorption. [25] A relatively unexplored factor affecting adsorption performance is the presence of defects. A recent study using UiO-66-NH 2 found that an increased concentration of defects due to missing linkers connecting the metal oxide nodes led to enhanced Pb 2þ adsorption. [26] Defects might enhance adsorption performance simply by facilitating diffusion through the MOF but it is also possible that they might play an active role in adsorption. For instance, in another recent study, the defect sites created at the metal oxide nodes by absent linkers were directly functionalised to be capped by thiourea and amidinothiourea groups; these functional groups adsorbed Pb 2þ with a maximum adsorption capacity of 245.8 mg g À1 . [25] Defects spontaneously formed under experimental conditions may also act as adsorption sites. Previous reports suggest that in aqueous conditions, the most common type of defect likely to form consists of a hydroxide and water group capping the undercoordinated Zr atoms. [21,27] These capping groups could potentially coordinate Pb 2þ . However, it remains difficult to establish the effect of defects conclusively simply based on experimental adsorption capacities of samples reported in the literature owing to the many confounding factors that influence defect concentration. For example, many synthesis steps that are not directly targeted towards defect engineering can lead to increased defect concentration, such as washing the MOF with solvent [28] or synthesising UiO-66 with functionalised and bulky linkers [29,30] (e.g. the precursors to UiO-66-(COOH) 2 and UiO-66-(COOH) 4 MOFs previously mentioned [24] ), consequently influencing adsorption capacity.
Overall, the recent experimental work in the field indicates that both defective metal oxide nodes with unsaturated metal sites capped by nucleophilic groups as well as functionalised linkers could play a complementary role in Pb 2þ adsorption. This emphasises the need for more research efforts directed at untangling the multitude of factors that influence Pb 2þ adsorption. This article seeks to address this gap by comparing the viability of linker functional groups and node defect sites as Pb 2þ adsorption sites on UiO-66 using a robust computational model centred on density functional theory (DFT) simulations with cluster models. We first examine the coordination geometry of solvated Pb 2þ generated by a range of computational methods and benchmark our computational method against relevant experimental results. Additionally, we explore the preferred coordination number of solvated Pb 2þ by considering the energetics of sequential water addition, solvation energies, and Pb-O bond lengths; we compare our findings with previous literature on Pb 2þ hydration. This is a crucial step to ensure the accuracy of our computational model. We then employ our model of solvated Pb 2þ to determine the energetics of adsorption onto a range of functionalised linkers and compare the results with non-functionalised UiO-66. To complement this study, we also investigate how the adsorption thermodynamics is affected by partial desolvation, which may be required for solvated Pb 2þ to pass through the UiO-66 apertures and diffuse through the MOF. Finally, we examine adsorption onto UiO-66 defective nodes and compare the reaction free energies (DGs) of adsorption with those obtained for the adsorption sites on the functionalised linkers. Out of a wide range of potential defect sites present in UiO-66 depending on synthesis conditions, [19,31] we chose to focus our analysis on the defect site with a hydroxide and water molecule bound to the open Zr metal sites created by the missing BDC linker (i.e. UiO-66-[OH -/H 2 O] adsorption site). We selected this site for several reasons, including its ubiquity in aqueous conditions [21] and superior stability compared with other hydroxide and water containing defect sites, [27] as well as the presence of several oxygen atoms available to coordinate the Pb 2þ cation. We conclude by discussing the key aspects that lead to more favourable adsorption energies based on our analysis to help disentangle the complex interplay of factors that could lead to enhanced Pb 2þ adsorption. From this, we determine general design principles for the fine-tuning of UiO-66 and other MOFs for Pb 2þ removal.

Results and Discussion
In order to accurately model Pb 2þ interactions with UiO-66, it is important to first determine how the cation prefers to be solvated in aqueous solution. The solvation shell of Pb 2þ has been studied extensively both experimentally and computationally in previous literature with the main aims of determining the coordination number of Pb 2þ and determining the orientation of the solvating water molecules around the cation. In the following, we first give an overview of these results from the literature and then compare them with our own computational findings to benchmark our model.
There are two relevant solvation geometries for Pb 2þ (Fig. 1). Holodirected geometries involve a symmetric organisation of ligands around the central Pb 2þ atom (Fig. 1a), due to the stereochemically inactive 6s lone pair. However, it is possible for there to be a visible void in the structure, leading to a hemidirected geometry (Fig. 1b). [32] It has been suggested in the literature that this is due to the hybridisation of the inert 6s orbital with the 6p orbitals to form a stereochemically active lone pair. [33] Alternatively, it has been proposed that the structural distortions visible in the hemidirected geometries are due to the Pb 2þ complex attempting to minimise energetically unfavourable antibonding interactions between Pb 2þ and the coordinating ligands. [34] There are several factors that influence whether a given Pb 2þ complex exhibits a hemidirected or holodirected geometry. One factor is the coordination number; it is often reported in the literature that Pb 2þ complexes with coordination numbers less than six exhibit hemidirected geometries while coordination numbers greater than nine always produce holodirected geometries. [32,33,35] A hemidirected geometry is also more likely if there is higher charge transfer between the ligands and Pb 2þ , there are attractive interactions between the ligands, or if the Pb-ligand bond is more ionic in nature. [32,36] It is possible for more unusual geometries such as bisdirectional structures to form where the lone pair of Pb 2þ is split over the equatorial plane [37] but these structures are formed owing to bulky ligands such as protein fragments and hence are not relevant to the present study. A study by Persson et al. used extended X-ray absorption fine structure (EXAFS) spectroscopy to determine the Pb-O bond lengths of Pb 2þ in aqueous solution. The authors found a large distribution of Pb-O bond lengths and this result along with the absence of multiple scattering contributions in the spectra indicated asymmetry in the structure, suggesting a hemidirected arrangement. [34] Furthermore, the derived ionic radii and average Pb-O bond lengths compared with the spectral data led the authors to conclude the most likely coordination geometry involves six water molecules in the Pb 2þ first solvation shell. [34] A flurry of computational work has also aimed to conclusively establish the structure of the first solvation shell of Pb 2þ , [35,36,[38][39][40] but the final geometries obtained are heavily dependent on the computational method. Molecular dynamics simulations indicate that Pb 2þ has a coordination number of nine and a holodirected geometry. [38,40,41] On the other hand, static quantum mechanical calculations find that coordination numbers between six and nine are all isoenergetic and thermodynamically accessible at room temperature. [35,36,39] The chosen DFT functional, basis set, and dispersion correction scheme greatly affect whether the solvated Pb 2þ ion optimises to a holodirected or hemidirected geometry. Wander and Clark found holodirected geometries for Pb(H 2 O) n 2þ (n ¼ 6-8) structures using B3LYP/aug-cc-pvdz-PP geometry optimisations conducted in the gas phase. [36] Kuznetsov et al. attempted to replicate their geometries using a range of model chemistries and found that incorporating implicit solvation in the geometry optimisation transformed initial holodirected structures into hemidirected structures; [39] they additionally found that the vB97X-D/6-311þþG(d,p) and vB97X-D/aug-cc-pVDZ model chemistries led to an average Pb-O bond length of 2.56 Å , extremely close to the experimental value of 2.54 Å found by Persson et al. [34] In order to disentangle the effect of model chemistry on coordination geometry, we optimised the Pb(H 2 O) 6 2þ structure using a range of functionals and basis sets as well as testing the effect of including the D3BJ dispersion correction [42,43] ( Table 1). The four functionals chosen were vB97X [44] and vB97X-D3BJ [45,46] (as both of these functionals have been used successfully to model heavy metals [39,47] ), B3LYP [48] (a common functional used in the literature for Pb 2þ hydration studies [36,39] ), and PBE0 [49] (a functional recommended specifically for optimising geometries of late transition metal structures [50] ). We studied a range of basis sets and effective core potentials (ECPs) chosen because they have been used in the literature for modelling either Pb 2þ hydration (aug-cc-pVDZ on H and O, aug-cc-pVDZ-PP with the Stuttgart-Koeln smallcore multiconfiguration Dirac-Hartree-Fock-adjusted (SK MCDHF RSC) ECP on Pb [51][52][53][54] as well as 6-31G(d,p) on H and O, LANL2DZ with the HayWadt ECP on Pb [55][56][57][58] ) or other heavy metals (def2-SVP on H and O, def2-TZVPP with the def2-ECP on Pb [53,59,60] ). We also tested removing the larger basis set on Pb and using def2-SVP on all atoms with either the def2-ECP or the HayWadt ECP for Pb. The initial geometry for each calculation was an octahedral holodirected configuration. The geometries generated by PBE0 are dependent on basis set choice: the Ahlrichs basis sets lead to hemidirected geometries while the Dunning basis sets lead to holodirected ones. B3LYP produces holodirected geometries with all tested basis sets and also leads to bond Pb-O lengths significantly smaller than the experimental ones (2.44 Å average with B3LYP/def2-SVP (def2-TZVPP with def2-ECP on Pb) v. 2.54 Å experimental value [34] ). Adding dispersion corrections (D3BJ) to B3LYP . Note the octahedral symmetry in the holodirected geometry (obtained using the B3LYP/def2-SVP (def2-TZVPP on Pb) model chemistry) and conversely the void in between the ligands for the hemidirected geometry (obtained using the vB97X-D3BJ/def2-SVP (def2-TZVPP on Pb) model chemistry). Table 1 provides a full summary of the effect of the computational methods tested in this study on the solvation geometry. Pb atoms are in grey, O in red, and H in off-white. significantly improves the model, creating a hemidirected geometry with an average Pb-O bond length of 2.57 Å (with the def2-SVP basis set (def2-TZVPP with def2-ECP on Pb)). This finding is consistent with a previous benchmarking paper on transition metal geometries that stressed the importance of dispersion corrections. [50] We found that vB97X-D3BJ always produced hemidirected geometries with all tested basis sets, as did the vB97X functional with the def2-SVP (def2-TZVPP with the def2-ECP on Pb) basis set. These hemidirected geometries have average Pb-O bond lengths (2.55 Å with vB97X-D3BJ/def2-SVP (def2-TZVPP with the def2-ECP on Pb)) in excellent agreement with experimental findings. [34] Using the vB97X-D3BJ/aug-cc-pVDZ (aug-cc-pVDZ-PP with the SK MCDHF RSC ECP on Pb) model chemistry slightly worsens the average Pb-O bond length (2.58 Å ), indicating that the Ahlrichs basis sets are the most suitable for these calculations. Overall, we concluded that our chosen model chemistry (vB97X-D3BJ/def2-SVP (def2-TZVPP with def2-ECP on Pb)), which is further discussed in the Computational Details section, provides the best description of the coordination geometry of hydrated Pb 2þ . Additionally, as illustrated in the following, this approach led to optimised hemidirected coordination geometries also for Pb 2þ coordinated by ligands other than water.
In addition to the coordination geometry of hydrated Pb 2þ , we also investigated the preferred coordination number. We found that all coordination numbers ranging from six to nine are thermodynamically accessible at room temperature, owing to the small energy cost or gain for transitioning between these coordination numbers (maximum of 3.8 kcal mol À1 (1 kcal mol À1 ¼ 4.186 kJ mol À1 ); Fig. 2). A negative DG value indicates exergonicity. This finding is fully consistent with previous literature [35,36,39] and reflects the extremely dynamic nature of the solvation shell. Additionally, we also found that all of these coordination numbers lead to similar computed  6 21 complex generated with a variety of functionals and basis set combinations The effect of dispersion corrections in the form of the D3BJ correction is also considered. All geometries were optimised from an octahedral holodirected geometry. Additionally, all geometry optimisations were conducted in the gas phase. The acronyms ECP and SK-MCDHF-RSC ECP stand for effective core potential and Stuttgart-Koeln small-core multiconfiguration-Dirac-Hartree-Fock-adjusted effective core potential, respectively Functional Basis set on O, H Basis set on Pb ECP on Pb Coordination geometry of optimised structure . This result is consistent with a previous computational study on Pb 2þ hydration. [39] Additionally, we found that coordination numbers below seven led to hemidirected geometries consistent with experimental findings, [34] while coordination numbers eight and nine had holodirected geometries due to the steric crowding of water molecules forcing a symmetric arrangement (see Fig. S1 and accompanying discussion in the Supplementary Material). Overall, based on these findings and the comparison with the literature, we identified hemidirected Pb(H 2 O) 6 2þ as the best model of solvated Pb 2þ , and therefore chose to use it to model solvated Pb 2þ interaction with UiO-66 in the following.
We started our investigation of Pb 2þ adsorption by UiO-66 by focusing on the functionalised linkers as the adsorption sites, as the bulk of experimental work has been devoted to these systems. [22][23][24]26] We selected a variety of functionalised linkers that have either demonstrated success for adsorption of Pb 2þ (UiO-66-COOH/COO - [24] , UiO-66-NHCSNHCH 3 , [23] UiO-66-NCS, [23] UiO-66-NCO, [23] and UiO-66-NH 2 [22,23,26] ) or possess attributes that could improve Pb 2þ adsorption (UiO-66-SO 3 H/SO 3 -, UiO-66-CONH 2 , UiO-66-SH). For conciseness, throughout the manuscript we refer to the functionalised linkers by their functional groups (e.g. UiO-66-SH is referred to as -SH). While -NHCSNHC 6 H 6 was also tested by Saleem et al., [23] we excluded it from our study owing to its significant steric bulk and relatively poor experimental adsorption performance. As discussed in the Computational Details section, the functionalised linkers adsorption sites were modelled using small cluster models consisting of benzene rings with the functional group attached (see Fig. 3); this approach is well documented in the literature. [61,62] This model captures the bonding interaction between Pb 2þ and the functional group at the adsorption site, thus allowing us to determine the relative favourability of each functionalised linker for Pb 2þ absorption. The optimised geometries of the functionalised linkers coordinating solvated Pb 2þ are illustrated in Fig. 3 and the DGs for Pb(H 2 O) 6 2þ adsorption by these linkers are compared in Fig. 4. A coordination number of six for Pb 2þ is maintained across reactant and product to avoid biasing the DGs; the adsorption reactions are hence modelled as a ligand exchange where solvating water molecules are replaced by functionalised linkers.
Out of all studied functionalised linkers, -COO -(-7.8 kcal mol À1 , Fig. 4), -NHCSNHCH 3 (-3.2 kcal mol À1 , Fig. 4), and -SO 3 -(þ2.0 kcal mol À1 , Fig. 4) display the largest thermodynamic driving force for binding to Pb(H 2 O) 6 2þ and displacing solvating water molecules. Factors that lead to enhanced Pb 2þ affinity for these functional groups are negative charge (-COO -, -SO 3 -) and the ability to chelate the Pb atom (-COOand -NHCSNHCH 3 , with the latter coordinating Pb 2þ with both the S site and the benzene as shown in Fig. 3); the combination of these two factors for -COOleads it to be the most favourable functionalised linker for Pb 2þ adsorption, consistent with experimental reports showing a high Pb 2þ adsorption capacity by UiO-66-(COOH) 2 . [24] We attempted to stabilise a multidentate binding geometry also for -SO 3 -, but the significant steric strain imposed by this planar tridentate binding mode caused this binding configuration to be less favourable than a monodentate binding mode by 2 kcal mol À1 . Similarly, the only feasible adsorption geometry for -CONH 2 involves the Pb 2þ atom bound to the carbonyl group, as we found that geometries with Pb 2þ bound to the amine group were not stable upon geometry optimisation. While binding solely to the sulfur atom is unfavourable (-SH, 8.1 kcal mol À1 , Fig. 4), the additional interaction with the benzene ring in -NHCSNHCH 3 stabilises the complex significantly; this interaction is possible due to the size and flexibility of the latter functional group. Owing to these steric considerations, there is no other functionalised linker that allows this interaction with the benzene ring. Several linkers such as -NH 2 (þ6.4 kcal mol À1 ), -NCO (þ14.5 kcal mol À1 ), and -NCS (þ12.3 kcal mol À1 ) perform worse than pristine UiO-66 but have been experimentally demonstrated to be more successful for Pb 2þ removal from water, [22,23] thus suggesting that factors other than the thermodynamic binding strength influence Pb 2þ adsorption. This aspect is further discussed later in the manuscript.
The overall trend in Fig. 4 is mostly consistent with hard-soft acid-base (HSAB) theory, [63] with three exceptions (-COO -, -NHCSNHCH 3 , -SO 3 -). HSAB theory classifies electron acceptors (acids) and electron donors (bases) as 'hard' or 'soft' based on criteria such as atomic radius, effective nuclear charge, and polarisability. Soft acids and bases typically possess large atomic radii and high polarisability, while hard acid and bases have the opposite characteristics. Acids form strong bonds with bases of the same class. Pb 2þ is a borderline acid. Our trend shows Pb 2þ base preference in the order N (borderline when bound to benzene) . benzene (soft) . S (soft) . O (hard), hence mostly as expected from HSAB theory except for the three linkers noted earlier. The preference for the -COO -, -NHCSNHCH 3 , and -SO 3 groups due to electrostatic and chelating effects has already been thoroughly explained earlier.
The overall trend is also fully consistent with previous work that found that Pb 2þ prefers to bind to soft groups such as benzene over hard groups such as ammonia and water. [64] Overall, while the absolute DGs are generally unfavourable (Fig. 4), this does not mean it is impossible for Pb 2þ to interact with these functionalised linkers. It is important to note that the relative ordering of functionalised linkers is a more significant result than the absolute thermodynamic binding strength. While extensive effort has been made on our part to accurately model the coordination geometry and energetics of solvated Pb 2þ , it is impossible to thoroughly account for the full chemical environment. This could include a second solvation shell, counterions, or solvents such as water or dimethylformamide (DMF) captured in the pores of the MOF during synthesis, which could all hamper adsorption either by direct competition with Pb 2þ or by limiting diffusion. Another aspect not considered by our model thus far is the fact that Pb 2þ may need to partially desolvate in order to pass through the small triangular apertures of pristine UiO-66 (6 Å side length) [18,65,66] despite the dynamic rotation of the BDC linkers. [67] The need for desolvation may also be exacerbated in the presence of functional groups on the BDC linkers that would lead to even smaller or more crowded apertures. The partial desolvation hypothesis is consistent with the fact that experimental studies indicate using UiO-67 [25] and the presence of defects [26] can lead to increased Pb 2þ adsorption, likely owing to easier diffusion of Pb 2þ through the MOF. Furthermore, partial desolvation in order to pass through apertures in porous materials is a well-documented process, [68,69] and moderate dehydration only has a minor energy cost (see Fig. 2). Given the plausibility of the partial desolvation hypothesis both in terms of sterics and energetics, it is important to analyse the possible fate of partially desolvated Pb 2þ once it has passed through the narrow UiO-66 apertures. We therefore decided to investigate the energetics of direct adsorption when the solvation shell was partially removed.
Pb(H 2 O) 4 2þ was selected to model partially desolvated Pb 2þ and study the effect of partial desolvation on the energetics of adsorbing onto the full range of studied functionalised linkers with no loss of water molecules (Fig. 5). We chose this initial species as the energetic cost of desolvation relative to fully solvated Pb 2þ is not excessively high (-8.6 kcal mol À1 to  remove two water molecules from Pb(H 2 O) 6 2þ , Fig. 2). The trends of favourability across linkers are maintained from Fig. 4, with some minor changes in ordering for nearly isoenergetic species (non-functionalised UiO-66, -NH 2 , -CONH 2 ). Not surprisingly, DGs for Pb(H 2 O) 4 2þ adsorbing without losing water molecules are overall ,10 kcal mol À1 more favourable than the DGs obtained for fully solvated Pb(H 2 O) 6 2þ undergoing linker/water exchange (Fig. 4). We additionally examined the extreme case of adsorption of fully desolvated Pb 2þ and found that removing all water molecules predictably leads to significantly more favourable adsorption DGs (e.g. DG ¼ -64.8 kcal mol À1 for adsorption on -COOas shown in Fig. S2 and accompanying discussion in the Supplementary Material). The latter study also showed that desolvation may lead to increased accessibility of additional binding sites on the functionalised linkers with multiple adsorption sites once the steric bulk caused by the solvation shell of Pb 2þ is removed (see Section 4 of the Supplementary Material). In particular, we were able to stabilise adsorption geometries where the bare Pb 2þ ion is bound to the nitrogen atom for the -NCO and -NCS functionalised linkers in addition to the adsorption geometries on the terminal oxygen and sulfur atoms.
We further explored the adsorption thermodynamics by studying the different extent of desolvation that may be induced by the small apertures as well as different degrees of rehydration that may occur once Pb 2þ enters the pore. In Fig. 6, we report the energetics of adsorption of partially desolvated and fully solvated Pb 2þ onto -COOwith varying final coordination number. The -COOadsorption site was chosen for this case study as it led to the most thermodynamically favourable adsorption (Fig. 4). As already noted in Fig. 5, the DGs become significantly more favourable for species with a lower initial coordination number. A final coordination number of six is preferred (-16.4 kcal mol À1 for a Pb(H 2 O) 4 2þ initial species) compared with lower coordination numbers such as four (-10.8 kcal mol À1 ) and five (-15.3 kcal mol À1 ) or a higher final coordination number of seven (-13.6 kcal mol À1 ); this trend is consistent across all initial coordination numbers investigated. The reason for a final coordination number of six being preferred over seven for -COOis the steric strain generated by attempting to maintain a hemidirected geometry at higher coordination numbers, as indicated by the elongation of the (-8.6 kcal mol À1 in total, Fig. 2). This suggests that the desolvation of Pb 2þ possibly occurring under experimental conditions could truly enhance the affinity of Pb 2þ to functional groups, making adsorption competitive, if not more favourable, than rehydration. Overall, our results suggest that desolvation induced by small apertures may favour adsorption and removal of Pb 2þ , but there is clearly a fine balance between driving desolvation and impeding diffusion. Molecular dynamics simulations may be used to further explore this aspect in future work.
To summarise, our analysis of the Pb 2þ cation interacting with different functionalised linkers of UiO-66 shows that changing the functional group may lead to significantly different adsorption outcomes in agreement with the literature on UiO-66 [22][23][24] as well as other materials. [11,12] However, our observations do not fully account for all reported experimental results in terms of which functional groups lead to the greatest  adsorption, [23,24,26] even when the effects of potential Pb 2þ desolvation are taken into account. This indicates other factors such as the presence of defects should also be critically analysed. This aspect is investigated next. Defects in UiO-66 are ubiquitous, [21] but their impact on Pb 2þ adsorption has yet to be thoroughly determined. Our results indicate that defects could play a significant role in the  adsorption geometries, using several different initial guesses in which the Pb 2þ cation was placed close to different potential adsorption sites on UiO-66-[OH -/H 2 O]. To improve the accuracy of our geometry search, we used both the solvated and the fully desolvated cation when building initial guesses. Working with the fully desolvated ion allowed us to easily probe all the possible adsorption sites, as Pb 2þ would not be restricted in making new bonds by its hydration shell. On the other hand, the use of the solvated cation allowed us to determine how much of its first solvation shell Pb 2þ preferred to retain when interacting with the defective UiO-66 nodes. In all of the identified adsorption geometries (Fig. 7b-e), Pb 2þ has a coordination number of six.
In the most energetically favourable structure (Fig. 7b), Pb 2þ forms two bonds with the MOF: one to the -OHgroup in the defect site and another bond to the carboxylic oxygen belonging to a BDC linker. The structure in Fig. 7c is also characterised by two bonds -one to the -OHgroup and one to the -H 2 O group on the defect site -but is less stable by 4 kcal mol À1 . We suggest the latter structure is less stable because of bond constraints on Pb 2þ ; more specifically, when Pb 2þ binds to the -OHgroup on the defect site (which is extremely favourable, given that this group is negatively charged), it is much easier for Pb 2þ to bind to the carboxylic oxygen on the BDC linker rather than the oxygen on the -H 2 O group. In fact, the bond between Pb 2þ and the carboxylic oxygen in Fig. 7b (2.44 Å ) is much closer in length to the bonds in Pb(H 2 O) 6 2þ (see Fig. S1 in the Supplementary Material) than the bond between Pb 2þ and the -H 2 O group on the defect site in Fig. 7c (2.89 Å ).
In the structure in Fig. 7d, Pb 2þ makes three bonds: one to the -OHgroup and another to the -H 2 O group on the defect site, and a third bond to m3-OH group in the MOF. This structure is 8.3 kcal mol À1 less stable than the most stable structure (Fig. 7b). The structure in Fig. 7e is the least stable structure (DG 13.7 kcal mol À1 less favourable than the most stable structure) and is characterised by five bonds between Pb 2þ and the MOF: one to the -OHgroup and one to the -H 2 O group on the defect site, two bonds to different carboxylic oxygens from BDC linkers, and a fifth bond to m3-O in the MOF. From this, a clear negative trend between stability and the number of bonds between Pb 2þ and UiO-66 emerges. We believe that this trend is explained by the increased steric constraints on Pb 2þ coordination geometry, since making bonds to UiO-66 would prevent Pb 2þ from freely orientating the coordination shell as if it was an ideal solvation shell. This is especially evident for the structure in Fig. 7e, where the Pb 2þ ion sits deep within the node. This results in one out of the six Pb-O bonds (the bond to the carboxylic O; 2.93 Å ) being significantly elongated relative to others (2.49 Å average). Additionally, this structure is characterised by the greatest distortion of the UiO-66 node as the Zr-m3-O bonds in the MOF are lengthened upon Pb 2þ adsorption, with ,0.05 Å average distortion (Fig. 7a).
Overall, we observe that Pb 2þ chose to bind to the -OHgroup on the defect site in all cases. This is easily explained by the electrostatic attraction between Pb 2þ and the negative charge on the -OHgroup and is consistent with our results showing that Pb 2þ prefers to bind to negatively charged functional groups. While it is possible for monocarboxylic acids such as formate, acetate or trifluoroacetate to act as capping ligands for defects, the only available adsorption site for Pb 2þ is the carboxylic oxygens (i.e. the same adsorption site as pristine . However, as demonstrated by the fact that all the identified adsorption configurations involve Pb 2þ bound to the negatively charged hydroxide group on the UiO-66-[OH -/H 2 O] site, the hydroxide group belonging to the defect site is a preferred coordination site relative to the carboxylic oxygens of BDC and monocarboxylic acids. This result offers further justification for our analysis being focussed solely on the UiO-66-[OH -/H 2 O] defect rather than other defects generated by monocarboxylic acid modulators and solvents acting as capping ligands. While binding to the OHgroup on the UiO-66 defect site is always preferred, it is also favourable to adsorb onto the H 2 O site (Fig. 7c) as well as to an arm of the UiO-66 (i.e. carboxylic O of the BDC linkers, Fig. 7b). Binding to deeper parts of UiO-66 (i.e. m3-O or m3-OH, Fig. 7d, e) is less favourable, likely owing to increased bond strain as previously discussed. We conclude that there is a significant thermodynamic driving force (up to ,12 kcal mol À1 ) for Pb 2þ to bind to the MOF defect site, as long as a large portion of its hydration shell is retained and Pb-O bond strain is limited.
A comparison of all the energetics of adsorption reported in this work suggests that both defective nodes and functionalised linkers could cooperatively contribute to Pb 2þ adsorption, with the defective UiO-66-[OH -/H 2 O] node showing the greatest thermodynamic affinity for Pb 2þ . The inclusion of negatively charged sites was found to enhance the adsorption strength at both defective nodes and adsorption sites owing to electrostatic attraction between the MOF and Pb 2þ , as indicated by the thermodynamic preference of Pb 2þ for the -COOand -SO 3 linkers and the -OHgroup at the defect site. Based on these results, we propose that the deprotonated m3-OH group may also act as an adsorption site. In fact, in the acidic conditions typically utilised for Pb 2þ adsorption, [24,25] some of m3-OH groups are likely to be deprotonated according to the pK a s determined in previous work (3.15, 4.43, 6.9, and 11). [70] Furthermore, adsorption modes involving coordination by both the OHgroup at the defect site and the m3-Ogroup may lead to even stronger adsorption. This aspect will be the subject of future investigations. Chelation of Pb 2þ (i.e. the formation of multiple bonds between Pb 2þ and the adsorption site) also plays a significant role, as demonstrated by our results for the -COOand -NHCSNHCH 3 adsorption sites as well as the UiO-66-[OH -/H 2 O] site. Given these results, we hypothesise that multiple functionalised linkers may coordinate Pb 2þ at once, leading to an increase of Pb 2þ binding affinity to the MOF. On the other hand, coordination by multiple linkers may lead to significant constraints on the coordination geometry of Pb 2þ and result in less favourable adsorption thermodynamics. This is demonstrated by our results showing that too many bonds to the UiO-66 defective node leads to distortion of either the MOF itself or the Pb 2þ coordination shell.
The comparable DGs between functionalised linkers and nodes hint at the complex interconnectedness of factors that contribute to a measured adsorption capacity. Functionalised linkers could have many subtle effects on the chemistry and structure of the MOF beyond acting as direct adsorption sites. As an example, a study by Wei et al. found that using functionalised linkers during synthesis led to more linker vacancies (i.e. defect sites). [30] This suggests that the extremely high adsorption capacity reported by Zhao et al. with UiO-66-(COOH) 2 of 420 mg g À1 [24] may be due to a synergistic effect with Pb 2þ adsorption occurring both at -COOfunctionalised linkers (as proposed by the authors of the study) and at defective nodes created by the synthesis conditions. The functional groups on linkers could also influence the nucleophilicity of the defect sites, in addition to the quantity. Wei et al. also found that functionalised linkers alter the electron density, and hence the chemical behaviour, of Zr atoms on the nodes; it is hence possible that some functionalised linkers could subtly affect the chemical behaviour of capping hydroxide and water molecules on nearby defect sites. [30] These effects could be investigated in future work by using either larger cluster models that incorporate both functionalised linkers and defective nodes in one structure or periodic calculations. Another indirect effect of functionalised linkers previously discussed is the reduction of the UiO-66 aperture size. This could have a complex effect on Pb 2þ adsorption, as we demonstrated that desolvation potentially induced by these functional groups leads to more thermodynamic driving force for adsorption. At the same time, numerous and bulky functional groups may hamper Pb 2þ diffusion through the MOF, [24] thus limiting adsorption. Overall, our investigation and analysis of all the factors discussed in this work provide critical design principles for the synthesis of UiO-66 MOFs to be employed in water purification from toxic Pb 2þ .

Conclusions
In this work, we present a DFT analysis of Pb 2þ adsorption using functionalised and defective UiO-66. We began our study by thoroughly exploring the effect of the model chemistry on the coordination geometry and coordination number of solvated Pb 2þ and validate our computational model by comparing our results with the literature. We found that B3LYP, while a popular functional in the literature, could not replicate the hemidirected coordination geometry of Pb 2þ found in experimental reports. The functional vB97X-D3BJ produced a hemidirected geometry for Pb 2þ solvated by six water molecules; additionally, the average Pb-O bond length of 2.55 Å in this optimised complex is almost identical to the experimental value of 2.54 Å . [34] Overall, our chosen model chemistry (vB97X-D3BJ/def2-SVP (def2-TZVPP with def2-ECP on Pb)) led to results in agreement with previous computational as well as experimental work. All coordination numbers between six and nine were found to be thermodynamically accessible at room temperature, with the six-coordinated Pb 2þ species giving a solvation energy in best agreement with experiments. Our study on functionalised UiO-66 for Pb 2þ capture identified several functional groups that may lead to more favourable adsorption energetics relative to non-functionalised UiO-66, in which the main adsorption site on the BDC linker is the benzene ring. We determined that the factors that lead to linkers with greater Pb 2þ affinity are electrostatic attraction and the ability to chelate Pb 2þ , and identified UiO-66-COOand UiO-66-NHCSNHCH 3 as the most favourable absorption sites for Pb 2þ removal, in agreement with experimental reports. [23,24] We further explored the possibility that solvated Pb 2þ may not be able to pass through the small UiO-66 triangular window while maintaining its full solvation shell; this effect may be worsened by the further narrowing of the aperture by bulky functional groups, hindering diffusion. The solvated Pb 2þ ion may thus have to partially or fully desolvate in order to diffuse through the MOF and reach adsorption sites. We therefore modelled adsorption of partially desolvated Pb 2þ onto the functionalised linkers without loss of water molecules and found that, not surprisingly, desolvation led to energetics of adsorption that are significantly more exergonic. Our case study on partially desolvated Pb 2þ adsorption using UiO-66-COOdetermined that there is a strong preference for adsorption onto a linker compared with rehydration within the pore. The thorough screening of the functionalised linkers carried out in this work also revealed that some adsorption sites on the functionalised linkers such as UiO-66-NCO were less thermodynamically favourable than non-functionalised UiO-66. However, these functional groups have demonstrated moderate experimental efficacy for Pb 2þ adsorption, indicating that factors other than the thermodynamic binding strength may be at play. There is likely a complex interplay of factors that leads to enhanced adsorption, and one that has not yet been fully explored by the literature is the presence of defects, which are ubiquitous in UiO-66. Defects can likely play a significant role in adsorption by creating larger apertures and allowing Pb 2þ to diffuse more easily through the MOF pores. A novel finding presented in this work is that defect sites created by missing linkers could boost adsorption capacity by directly acting as an adsorption site. The UiO-66-[OH -/H 2 O] site, formed by hydroxide and water molecules capping the undercoordinated Zr atoms created by the missing BDC linker, is able to coordinate Pb 2þ using the many oxygen atoms available at this site. The capacity of the UiO-66-[OH -/H 2 O] defect site to directly adsorb Pb 2þ could contribute to explaining high adsorption capacities reported in some experimental reports. Future work could be aimed at modelling the process of Pb 2þ desolvating, rehydrating, and diffusing through the UiO-66 pores in the presence of defects and functionalised linkers using molecular dynamics simulations. Overall, our work provides several new insights into the removal of toxic Pb 2þ from water using functionalised UiO-66 and suggests that defects may be playing a critical role, although this area of research is still mostly unexplored. We also determined that electrostatic attraction and the ability to chelate Pb 2þ without unduly constraining its hydration shell are key factors that must be considered when designing MOFs for water purification from toxic cations.

Computational Details
All of the calculations in this study were conducted using Kohn-Sham DFT [71,72] as implemented in Orca 4.2.1. [73][74][75] The protocol for these calculations is based on previous work in our group where we successfully modelled As V adsorption using UiO-66. [47] The vB97X-D3BJ range-separated hybrid functional was used as it produced geometries consistent with experimental findings as discussed in the Results and Discussion section; a previous computational study on the hydration shell of the Pb 2þ ion found a similar result using the vB97X functional. [39] The included parameterised D3BJ dispersion correction [42,43] is recommended in the literature for optimising geometries of transition metal complexes. [50] The RIJCOSX approximation [76] was also applied in order to reduce computational cost. Final Gibbs free energies in solution were calculated using a multi-step approach. Geometry optimisations and vibrational frequency calculations were conducted using the def2-SVP basis set on all atoms except for Pb and Zr, which were described with the def2-TZVPP basis set [59,60] in conjunction with the def2-ECP effective core potential. [53,77] The def2-ECP describes the core 60 electrons in the Pb atom and the core 28 electrons for Zr. Single-point refinements were conducted on the optimised geometries using the ma-def2-TZVPP basis set for Pb and Zr and the ma-def2-TZVP for all other atoms. [59,78] The use of def2-TZVPP for transition metals is recommended in the literature; [50] the polarisation and diffuse functions present in these basis sets are ideal for modelling transition metals as well as non-covalent interactions such as the weak interactions between the solvating water molecules.