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

Keywords: lead, adsorption, water purification, UiO-66, metal–organic frameworks, computational chemistry, heavy metals, defects.

As a bivalent cation, lead has diverse and extremely detrimental impacts on the body in all quantities.^[1,2] The high toxicity of Pb²⁺ 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–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²⁺ 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⁻¹ reported Pb²⁺ adsorption capacity for tripolyphosphate-modified kaolinite clay^[15]) and activated carbon (Pb²⁺ 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₆(OH)₄(O)₄ 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²⁺ 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₂ using a microwave-promoted synthesis technique and found a maximum Pb²⁺ adsorption capacity of 116.74 mg g⁻¹.^[22] Saleem et al. conducted post-synthetic modification in order to replace the amine functional group on UiO-66–NH₂ with a range of sulfur- and nitrogen-containing groups such as NCO, NCS, and NHCSNHCH₃.^[23] The latter bulky functional group led to a maximum adsorption capacity of 232 mg g⁻¹. A significant adsorption capacity of 420 mg g⁻¹ by UiO-66–(COOH)₂ (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²⁺ and the deprotonated –COO^– groups, which are likely to dominate at pH 6. Using the more functionalised UiO-66–(COOH)₄ led to less Pb²⁺ 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²⁺ 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²⁺ adsorption.^[25]

A relatively unexplored factor affecting adsorption performance is the presence of defects. A recent study using UiO-66–NH₂ found that an increased concentration of defects due to missing linkers connecting the metal oxide nodes led to enhanced Pb²⁺ 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²⁺ with a maximum adsorption capacity of 245.8 mg g⁻¹.^[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²⁺. 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)₂ and UiO-66–(COOH)₄ 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²⁺ adsorption. This emphasises the need for more research efforts directed at untangling the multitude of factors that influence Pb²⁺ adsorption. This article seeks to address this gap by comparing the viability of linker functional groups and node defect sites as Pb²⁺ 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²⁺ 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²⁺ by considering the energetics of sequential water addition, solvation energies, and Pb–O bond lengths; we compare our findings with previous literature on Pb²⁺ hydration. This is a crucial step to ensure the accuracy of our computational model. We then employ our model of solvated Pb²⁺ 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²⁺ 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 (ΔGs) 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₂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²⁺ 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²⁺ adsorption. From this, we determine general design principles for the fine-tuning of UiO-66 and other MOFs for Pb²⁺ removal.

In order to accurately model Pb²⁺ 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²⁺ has been studied extensively both experimentally and computationally in previous literature with the main aims of determining the coordination number of Pb²⁺ 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²⁺ (Fig. 1). Holodirected geometries involve a symmetric organisation of ligands around the central Pb²⁺ 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²⁺ complex attempting to minimise energetically unfavourable antibonding interactions between Pb²⁺and the coordinating ligands.^[34] There are several factors that influence whether a given Pb²⁺ complex exhibits a hemidirected or holodirected geometry. One factor is the coordination number; it is often reported in the literature that Pb²⁺ 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²⁺, 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²⁺ 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²⁺ 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²⁺ first solvation shell.^[34]

Fig. 1.  A visual depiction of (a) holodirected, and (b) hemidirected geometries for Pb(H₂O)₆²⁺. 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 ωB97X-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.

A flurry of computational work has also aimed to conclusively establish the structure of the first solvation shell of Pb²⁺,^{[35,36,38–40]} but the final geometries obtained are heavily dependent on the computational method. Molecular dynamics simulations indicate that Pb²⁺ 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²⁺ ion optimises to a holodirected or hemidirected geometry. Wander and Clark found holodirected geometries for Pb(H₂O)_n²⁺ (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 ωB97X-D/6–311++G(d,p) and ωB97X-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₂O)₆²⁺ 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 ωB97X^[44] and ωB97X-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²⁺ 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²⁺ hydration (aug-cc-pVDZ on H and O, aug-cc-pVDZ-PP with the Stuttgart–Koeln small-core multiconfiguration Dirac–Hartree–Fock-adjusted (SK MCDHF RSC) ECP on Pb^[51–54] as well as 6-31G(d,p) on H and O, LANL2DZ with the HayWadt ECP on Pb^[55–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 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 ωB97X-D3BJ always produced hemidirected geometries with all tested basis sets, as did the ωB97X 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 ωB97X-D3BJ/def2-SVP (def2-TZVPP with the def2-ECP on Pb)) in excellent agreement with experimental findings.^[34] Using the ωB97X-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 (ωB97X-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²⁺. Additionally, as illustrated in the following, this approach led to optimised hemidirected coordination geometries also for Pb²⁺ coordinated by ligands other than water.

Table 1.  Coordination geometry of the Pb(H₂O)₆²⁺ 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

In addition to the coordination geometry of hydrated Pb²⁺, 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 kcal mol⁻¹ = 4.186 kJ mol⁻¹); Fig. 2). A negative ΔG 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 solvation energies for Pb²⁺ within a 4 kcal mol⁻¹ range and that they are all in good agreement with the experimental value (see Table S1 and related discussion in the Supplementary Material). Among the coordination numbers investigated, we found that the Pb–O bond lengths in Pb(H₂O)₆²⁺ (2.55 Å average) are significantly closer to experimental values (2.54 Å average^[34]) than those associated with Pb(H₂O)₇²⁺ (2.62 Å), Pb(H₂O)₈²⁺ (2.66 Å), and Pb(H₂O)₉²⁺ (2.70 Å). This result is consistent with a previous computational study on Pb²⁺ 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₂O)₆²⁺ as the best model of solvated Pb²⁺, and therefore chose to use it to model solvated Pb²⁺ interaction with UiO-66 in the following.

Fig. 2.  The reaction free energies (ΔGs) for sequential water addition onto the Pb²⁺ cation to create Pb(H₂O)_n²⁺ (n = 1 – 9). The ΔGs plotted are calculated for the following reaction: Pb(H₂O)_n–1²⁺ + H₂O → Pb(H₂O)_n²⁺.

We started our investigation of Pb²⁺ 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–24,26] We selected a variety of functionalised linkers that have either demonstrated success for adsorption of Pb²⁺ (UiO-66–COOH/COO^–[24], UiO-66–NHCSNHCH₃,^[23] UiO-66–NCS,^[23] UiO-66–NCO,^[23] and UiO-66–NH₂^[22,23,26]) or possess attributes that could improve Pb²⁺ adsorption (UiO-66–SO₃H/SO₃^–, UiO-66–CONH₂, 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₆H₆ 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²⁺ and the functional group at the adsorption site, thus allowing us to determine the relative favourability of each functionalised linker for Pb²⁺ absorption. The optimised geometries of the functionalised linkers coordinating solvated Pb²⁺ are illustrated in Fig. 3 and the ΔGs for Pb(H₂O)₆²⁺ adsorption by these linkers are compared in Fig. 4. A coordination number of six for Pb²⁺ is maintained across reactant and product to avoid biasing the ΔGs; the adsorption reactions are hence modelled as a ligand exchange where solvating water molecules are replaced by functionalised linkers.

Fig. 3.  A visual depiction of solvated Pb²⁺ adsorbed onto UiO-66 functionalised with a variety of linkers modelled using small cluster models. The various functionalised linkers are labelled according to the functional group for brevity. Pb²⁺ maintains a coordination number of six and a hemidirected coordination geometry across all adsorption configurations. For functionalised linkers with multiple potential adsorption sites (–SO₃^–, –CONH₂, –NCS, –NCO), the most stable geometry is presented. Pb atoms are in grey, N in pale blue, S in yellow, O in red, C in brown, and H in off-white.

Fig. 4.  Comparison of the reaction free energies (ΔGs, kcal mol⁻¹) for Pb(H₂O)₆²⁺ binding onto different adsorption sites on functionalised linkers (labelled according to the functional group for brevity) while a coordination number of six is maintained. The ΔGs plotted are calculated for the following reaction: Pb(H₂O)₆²⁺ + UiO-66–Lⁿ → UiO-66–L–Pb(H₂O)_aⁿ⁺² + bH₂O where L is the functional group; n = 0 for all functional groups except –COO^– and –SO₃^– where n = –1; a = 5 in all cases except for –COO^– and –NHCSNHCH₃, which coordinate Pb²⁺ in a bidentate fashion leading to a = 4; b = 1 in all cases except for –COO^– and –NHCSNHCH₃ for which it is 2. The geometries of the product complexes are illustrated in Fig. 3.

Out of all studied functionalised linkers, –COO^– (–7.8 kcal mol⁻¹, Fig. 4), –NHCSNHCH₃ (–3.2 kcal mol⁻¹, Fig. 4), and –SO₃^– (+2.0 kcal mol⁻¹, Fig. 4) display the largest thermodynamic driving force for binding to Pb(H₂O)₆²⁺ and displacing solvating water molecules. Factors that lead to enhanced Pb²⁺ affinity for these functional groups are negative charge (–COO^–, –SO₃^–) and the ability to chelate the Pb atom (–COO^– and –NHCSNHCH₃, with the latter coordinating Pb²⁺ with both the S site and the benzene as shown in Fig. 3); the combination of these two factors for –COO^– leads it to be the most favourable functionalised linker for Pb²⁺ adsorption, consistent with experimental reports showing a high Pb²⁺ adsorption capacity by UiO-66–(COOH)₂.^[24] We attempted to stabilise a multidentate binding geometry also for –SO₃^–, 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⁻¹. Similarly, the only feasible adsorption geometry for –CONH₂ involves the Pb²⁺ atom bound to the carbonyl group, as we found that geometries with Pb²⁺ 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⁻¹, Fig. 4), the additional interaction with the benzene ring in –NHCSNHCH₃ 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₂ (+6.4 kcal mol⁻¹), –NCO (+14.5 kcal mol⁻¹), and –NCS (+12.3 kcal mol⁻¹) perform worse than pristine UiO-66 but have been experimentally demonstrated to be more successful for Pb²⁺ removal from water,^[22,23] thus suggesting that factors other than the thermodynamic binding strength influence Pb²⁺ 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₃, –SO₃^–). 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²⁺ is a borderline acid. Our trend shows Pb²⁺ 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₃, and –SO₃^– 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²⁺ prefers to bind to soft groups such as benzene over hard groups such as ammonia and water.^[64]

Overall, while the absolute ΔGs are generally unfavourable (Fig. 4), this does not mean it is impossible for Pb²⁺ 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²⁺, 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²⁺ or by limiting diffusion. Another aspect not considered by our model thus far is the fact that Pb²⁺ 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²⁺ adsorption, likely owing to easier diffusion of Pb²⁺ 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²⁺ 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₂O)₄²⁺ was selected to model partially desolvated Pb²⁺ 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²⁺ is not excessively high (–8.6 kcal mol⁻¹ to remove two water molecules from Pb(H₂O)₆²⁺, 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₂, –CONH₂). Not surprisingly, ΔGs for Pb(H₂O)₄²⁺ adsorbing without losing water molecules are overall ~10 kcal mol⁻¹ more favourable than the ΔGs obtained for fully solvated Pb(H₂O)₆²⁺ undergoing linker/water exchange (Fig. 4). We additionally examined the extreme case of adsorption of fully desolvated Pb²⁺ and found that removing all water molecules predictably leads to significantly more favourable adsorption ΔGs (e.g. ΔG = –64.8 kcal mol⁻¹ for adsorption on –COO^– as 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²⁺ is removed (see Section 4 of the Supplementary Material). In particular, we were able to stabilise adsorption geometries where the bare Pb²⁺ 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.

Fig. 5.  Comparison of the reaction free energies (ΔG, kcal mol⁻¹) for Pb(H₂O)₄²⁺ adsorbing onto different functionalised linkers (labelled according to the functional group for brevity) with no loss of water molecules. The ΔGs plotted are calculated for the following reaction: Pb(H₂O)₄²⁺ + Lⁿ → L–Pb(H₂O)₄ⁿ⁺², where L is the ligand and n = 0 for all functionalised linkers except –COO^– and –SO₃^– where n = –1.

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²⁺ enters the pore. In Fig. 6, we report the energetics of adsorption of partially desolvated and fully solvated Pb²⁺ onto –COO^– with varying final coordination number. The –COO^– adsorption 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 ΔGs become significantly more favourable for species with a lower initial coordination number. A final coordination number of six is preferred (–16.4 kcal mol⁻¹ for a Pb(H₂O)₄²⁺ initial species) compared with lower coordination numbers such as four (–10.8 kcal mol⁻¹) and five (–15.3 kcal mol⁻¹) or a higher final coordination number of seven (–13.6 kcal mol⁻¹); this trend is consistent across all initial coordination numbers investigated. The reason for a final coordination number of six being preferred over seven for –COO^– is the steric strain generated by attempting to maintain a hemidirected geometry at higher coordination numbers, as indicated by the elongation of the Pb–O bond to the fifth water in –COO^––Pb(H₂O)₅²⁺ (2.88 Å compared with the maximum Pb–O bond length of 2.76 Å in –COO^––Pb(H₂O)₄²⁺). A significant result is that it is more favourable for the destabilised Pb(H₂O)₄²⁺ species to adsorb onto –COO^– to a final coordination number of six (–16.4 kcal mol⁻¹, Fig. 6) relative to rehydrating to form Pb(H₂O)₆²⁺ (–8.6 kcal mol⁻¹ in total, Fig. 2). This suggests that the desolvation of Pb²⁺ possibly occurring under experimental conditions could truly enhance the affinity of Pb²⁺ 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²⁺, 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.

Fig. 6.  Comparison of the reaction free energies (ΔG, kcal mol⁻¹) for the solvated Pb²⁺ ion with varying initial coordination number adsorbing onto UiO-66–COO^– with varying final coordination number. The ΔGs plotted are calculated for the following reaction: Pb(H₂O)_a²⁺ + UiO-66–COO^– + bH₂O → UiO-66–COO^––Pb(H₂O)_c²⁺ + dH₂O, where a is the initial coordination number; c is the final coordination number with 2 subtracted owing to UiO-66–COO^– coordinating Pb²⁺ in a bidentate fashion, and b, d are the required number of water molecules for balancing the equation.

To summarise, our analysis of the Pb²⁺ 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–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²⁺ 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²⁺ adsorption has yet to be thoroughly determined. Our results indicate that defects could play a significant role in the adsorption of Pb²⁺ in aqueous solutions by forming thermodynamically favourable Zr–O–Pb bonds. Shown in Fig. 7a is the original hydroxide–water capped defect site (UiO-66–[OH^–/H₂O]), and Fig. 7b–e depicts the adsorption geometries of Pb²⁺ on this defect site investigated in this study with the coordination sites indicated in the caption. These four unique geometries were identified after a thorough search of stable adsorption geometries, using several different initial guesses in which the Pb²⁺ cation was placed close to different potential adsorption sites on UiO-66–[OH^–/H₂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²⁺ 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²⁺ preferred to retain when interacting with the defective UiO-66 nodes. In all of the identified adsorption geometries (Fig. 7b–e), Pb²⁺ has a coordination number of six.

Fig. 7.  Comparison of the reaction free energies (ΔG, kcal mol⁻¹) for solvated Pb²⁺ (Pb[H₂O]₆²⁺) adsorbing onto the hydroxide–water capped defect site ([UiO-66]–OH^–/H₂O) with different adsorption configurations while maintaining a coordination number of six. (a) Cluster model used to simulate the [UiO-66]–OH^–/H₂O site. (b) Adsorption configuration with Pb²⁺ bound to the OH^– group belonging to the defect site and to a carboxylic oxygen belonging to a 1,4-benzenedicarboxylic acid (BDC) linker. (c) Adsorption configuration with Pb²⁺ bound to the OH^– and H₂O groups belonging to the defect site. (d) Adsorption configuration with Pb²⁺ bound to the OH^– and H₂O groups belonging to the defect site as well as the μ3-OH group. (e) Adsorption configuration with Pb²⁺ bound to the OH^– and H₂O groups belonging to the defect site, two carboxylic oxygens belonging to two BDC linkers, and the μ3-O group. The ΔGs plotted are calculated using the following reaction: Pb(H₂O)₆²⁺ + UiO-66–[OH^–/H₂O] → UiO-66–[OH^–/H₂O]–Pb(H₂O)_(6–a)²⁺ + aH₂O, where a is the number of water molecules from the Pb²⁺ hydration shell in Pb(H₂O)₆²⁺ that are lost when Pb²⁺ binds to UiO-66–[OH^–/H₂O] (a = 2 for (b) and (c); a = 3 for (d); and a = 4 for (e)). Pb atoms are in grey, Zr in green, O in red, C in brown, and H in off-white.

In the most energetically favourable structure (Fig. 7b), Pb²⁺ forms two bonds with the MOF: one to the –OH^– group 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 –OH^– group and one to the –H₂O group on the defect site – but is less stable by 4 kcal mol⁻¹. We suggest the latter structure is less stable because of bond constraints on Pb²⁺; more specifically, when Pb²⁺ binds to the –OH^– group on the defect site (which is extremely favourable, given that this group is negatively charged), it is much easier for Pb²⁺ to bind to the carboxylic oxygen on the BDC linker rather than the oxygen on the –H₂O group. In fact, the bond between Pb²⁺ and the carboxylic oxygen in Fig. 7b (2.44 Å) is much closer in length to the bonds in Pb(H₂O)₆²⁺ (see Fig. S1 in the Supplementary Material) than the bond between Pb²⁺ and the –H₂O group on the defect site in Fig. 7c (2.89 Å).

In the structure in Fig. 7d, Pb²⁺ makes three bonds: one to the –OH^– group and another to the –H₂O group on the defect site, and a third bond to μ3-OH group in the MOF. This structure is 8.3 kcal mol⁻¹ less stable than the most stable structure (Fig. 7b). The structure in Fig. 7e is the least stable structure (ΔG 13.7 kcal mol⁻¹ less favourable than the most stable structure) and is characterised by five bonds between Pb²⁺ and the MOF: one to the –OH^–group and one to the –H₂O group on the defect site, two bonds to different carboxylic oxygens from BDC linkers, and a fifth bond to μ3-O in the MOF. From this, a clear negative trend between stability and the number of bonds between Pb²⁺ and UiO-66 emerges. We believe that this trend is explained by the increased steric constraints on Pb²⁺ coordination geometry, since making bonds to UiO-66 would prevent Pb²⁺ 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²⁺ 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–μ3-O bonds in the MOF are lengthened upon Pb²⁺ adsorption, with ~0.05 Å average distortion (Fig. 7a).

Overall, we observe that Pb²⁺ chose to bind to the –OH^– group on the defect site in all cases. This is easily explained by the electrostatic attraction between Pb²⁺ and the negative charge on the –OH^– group and is consistent with our results showing that Pb²⁺ 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²⁺ is the carboxylic oxygens (i.e. the same adsorption site as pristine UiO-66). However, as demonstrated by the fact that all the identified adsorption configurations involve Pb²⁺ bound to the negatively charged hydroxide group on the UiO-66–[OH^–/H₂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₂O] defect rather than other defects generated by monocarboxylic acid modulators and solvents acting as capping ligands. While binding to the OH^–group on the UiO-66 defect site is always preferred, it is also favourable to adsorb onto the H₂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. μ3-O or μ3-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⁻¹) for Pb²⁺ 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²⁺ adsorption, with the defective UiO-66–[OH^–/H₂O] node showing the greatest thermodynamic affinity for Pb²⁺. 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²⁺, as indicated by the thermodynamic preference of Pb²⁺ for the –COO^– and –SO₃^– linkers and the –OH^– group at the defect site. Based on these results, we propose that the deprotonated μ3-OH group may also act as an adsorption site. In fact, in the acidic conditions typically utilised for Pb²⁺ adsorption,^[24,25] some of μ3-OH groups are likely to be deprotonated according to the pK_as determined in previous work (3.15, 4.43, 6.9, and 11).^[70] Furthermore, adsorption modes involving coordination by both the OH^– group at the defect site and the μ3-O^– group may lead to even stronger adsorption. This aspect will be the subject of future investigations. Chelation of Pb²⁺ (i.e. the formation of multiple bonds between Pb²⁺ and the adsorption site) also plays a significant role, as demonstrated by our results for the –COO^– and –NHCSNHCH₃ adsorption sites as well as the UiO-66–[OH^–/H₂O] site. Given these results, we hypothesise that multiple functionalised linkers may coordinate Pb²⁺ at once, leading to an increase of Pb²⁺ binding affinity to the MOF. On the other hand, coordination by multiple linkers may lead to significant constraints on the coordination geometry of Pb²⁺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²⁺ coordination shell.

The comparable ΔGs 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)₂ of 420 mg g^−1 [24] may be due to a synergistic effect with Pb²⁺ adsorption occurring both at –COO^– functionalised 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²⁺ 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²⁺ 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²⁺.

In this work, we present a DFT analysis of Pb²⁺ 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²⁺ 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²⁺ found in experimental reports. The functional ωB97X-D3BJ produced a hemidirected geometry for Pb²⁺ 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 (ωB97X-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²⁺ species giving a solvation energy in best agreement with experiments. Our study on functionalised UiO-66 for Pb²⁺ 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²⁺ affinity are electrostatic attraction and the ability to chelate Pb²⁺, and identified UiO-66–COO^– and UiO-66–NHCSNHCH₃ as the most favourable absorption sites for Pb²⁺ removal, in agreement with experimental reports.^[23,24] We further explored the possibility that solvated Pb²⁺ 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²⁺ 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²⁺ 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²⁺ adsorption using UiO-66–COO^– determined 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²⁺ 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²⁺ 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₂O] site, formed by hydroxide and water molecules capping the undercoordinated Zr atoms created by the missing BDC linker, is able to coordinate Pb²⁺ using the many oxygen atoms available at this site. The capacity of the UiO-66–[OH^–/H₂O] defect site to directly adsorb Pb²⁺ could contribute to explaining high adsorption capacities reported in some experimental reports. Future work could be aimed at modelling the process of Pb²⁺ 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²⁺ 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²⁺ without unduly constraining its hydration shell are key factors that must be considered when designing MOFs for water purification from toxic cations.

All of the calculations in this study were conducted using Kohn–Sham DFT^[71,72] as implemented in Orca 4.2.1.^[73–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 ωB97X-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²⁺ ion found a similar result using the ωB97X 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.

We incorporated the effect of solvation using a mixed implicit/explicit solvation model. Implicit solvation was included using the Conductor-like Polarizable Continuum Model (CPCM) with a scaled van der Waals surface and Gaussian charge scheme.^[79,80] The first solvation shell of the Pb²⁺ cation was accurately modelled by including explicit water molecules as thoroughly discussed in the Results and Discussion section. Similarly to our previous study,^[47] we described water molecules in solution using a nonamer; the reasoning behind this choice is thoroughly discussed in Fig. S3 as well as Table S2 and associated discussion in the Supplementary Material.

Two types of cluster models were used to model UiO-66 depending on the adsorption site of interest. This approach has been used commonly and successfully in the literature to model UiO-66 nodes and linkers, and we direct the interested reader to some recent reviews on this topic.^[81,82] The cluster model used to simulate adsorption on the functionalised linker is constructed using a benzene ring with the relevant functional group attached (see Fig. 3 and accompanying discussion) to model functionalised UiO-66. This is a valid approach for modelling interactions with the adsorption site on the functionalised linkers as demonstrated by previous studies in the literature.^[61,62,83] The cluster model used to study adsorption on the UiO-66 defective node (Fig. 7 and accompanying discussion) has been used widely in the literature in the past,^{[27,30,84–87]} including in recent work by our group where we report the details of how the cluster model was obtained.^[47] For reasons thoroughly discussed in the Introduction section, we chose to model our defective node as UiO-66–[OH^–/H₂O] (i.e. a UiO-66 node with a missing linker defect and the undercoordinated Zr atoms saturated by a OH^– group and H₂O). This cluster model was then used to calculate the energetics of Pb²⁺ adsorption onto the defect site (see Fig. 7 and accompanying discussion).

The lack of imaginary frequencies in the vibrational frequency analysis was used to verify whether the optimised structure was at a local minimum on the potential energy surface. The frequency calculations also provided thermal corrections to the final electronic energy. For all species involving the MOF cluster models, only the vibrational corrections were considered in order to account for the constraints exerted by the ordered crystal structure that prevent translational and rotational degrees of freedom. The translational, rotational, and vibrational corrections were all included for the species free in solution, such as the water cluster or solvated Pb²⁺.

ΔGs were calculated according to the equation

where G_products is the sum of the free energies of all the products and G_reactants is the sum of the free energies of all the reactants. Favourable reaction free energies are associated with a negative ΔG value. G_products and G_reactants were calculated using the standard thermodynamic cycle approach in order to incorporate thermal and solvation effects.^[88,89] Using this cycle, the Gibbs free energy for each solvated species is given by:

where the superscripts B1 and B2 indicate that the energy is determined using the basis set choice for optimisations and frequencies (B1) or the basis set choice for refined single-point energy calculations (B2) respectively; the subscripts gas and solv signify whether the calculation was performed in gas phase or with implicit solvation; E denotes electronic energies while G denotes free energies at room temperature; ΔG^0→* is the correction for the change in standard states from gas phase to liquid phase as detailed in ref. [89].

More benchmark results, calculations details, and Cartesian coordinates of all the geometries used to obtain the data in this work are available on the Journal’s website.

The data that support this study are available in the article and accompanying online supplementary material.

The authors declare no conflicts of interest.

This research was supported by an Australian Government Research Training Program (RTP) Scholarship awarded to C. S. C.