Density functional theory insights into the solvent effect on the binding energies of Cd²⁺ in functionalized MOFs

Abstract

The clean-up of cadmium-contaminated waters is a significant environmental issue, and metal–organic frameworks (MOFs) have been shown to be effective adsorbents with great potential. Yet, the basic principles driving Cd²⁺ binding in multiple solvent environments remain enigmatic and hinder rational adsorbent design. In this work, the solvent effects on the binding energies of Cd²⁺, referring to (M = –NH₂, –OH, –COOH, and –NO₂), are systematically investigated using DFT in four solvents (water, ethanol, methanol, and acetonitrile). Using the SMD implicit solvation model and the M06-2X functional with def2-TZVP basis set, we demonstrate that solvent polarity significantly affects binding thermodynamics, with aqueous systems demonstrating 68–72% weaker binding than acetonitrile. With binding energies of −156.8 kcal mol⁻¹ (water), −187.3 kcal mol⁻¹ (ethanol), −181.4 kcal mol⁻¹ (methanol), and −198.5 kcal mol⁻¹ (acetonitrile), the carboxyl-functionalized MOF exhibits exceptional performance in all solvents. Protic solvents (1.18–1.24|e|) exhibit greater charge transfer than aprotic acetonitrile (0.94–1.08|e|), according to natural bond orbital analysis. This is explained by hydrogen bonding networks stabilizing charge-separated coordination states. Predictive relationships for mixed-solvent industrial wastewater applications are established by the significant linear correlation between binding energies and inverse solvent polarity (R² = 0.96–0.98). These results directly address the pressing need for efficient cadmium remediation technologies in a variety of aqueous-organic media by offering crucial theoretical guidance for the design of solvent-adaptive MOF adsorbents optimized for particular industrial effluent compositions.

---

1. Introduction

Heavy metal contamination of water resources is one of the most pressing environmental issues of the twenty-first century, endangering both human health and ecosystem integrity.^1,2 Cadmium (Cd²⁺), a priority pollutant by the World Health Organization and US Environmental Protection Agency, has extreme toxicity at concentrations as low as 5 µg L⁻¹, causing irreversible kidney damage, skeletal degradation, cardiovascular dysfunction, and confirmed carcinogenic effects.^3,4 Annually, millions of tons of Cd-contaminated wastewater are discharged into aquatic systems worldwide from industrial activities such as electroplating, battery manufacturing, pigment production, and mining.^5,6 The persistence and bioaccumulation of cadmium in food chains exacerbate its toxic effects, with documented cases of chronic poisoning affecting millions of people in contaminated areas.⁷

Traditional remediation methods, like chemical precipitation, ion exchange, and membrane filtration, have limitations such as incomplete metal removal, secondary pollution generation, high operational costs, and limited effectiveness in dilute solutions (<100 mg L⁻¹).^8,9 These deficiencies have prompted extensive research into advanced adsorbent materials capable of selective, efficient, and cost-effective heavy metal capture from complex wastewater matrices.¹⁰

Metal–organic frameworks (MOFs), crystalline porous materials composed of metal nodes and organic linkers, have revolutionized adsorption technologies due to their ultrahigh surface areas (up to 7000 m² g⁻¹), tunable pore structures, and functionalizable frameworks.^7–9 Unlike conventional adsorbents such as activated carbon or zeolites, MOFs offer unprecedented molecular-level control over adsorption sites through rational design of organic ligands.^10,11 Beyond heavy metal remediation, MOFs have demonstrated exceptional versatility across diverse applications including gas sensing and detection platforms,¹² optoelectronic devices with tunable band gaps,¹³ photocatalytic degradation of organic pollutants,¹⁴ and selective gas separation technologies for CO₂ capture and hydrogen purification.¹⁵ This multifunctional capability stems from their designable porosity, diverse metal nodes, and chemically modifiable linkers, positioning MOFs as a transformative materials platform for environmental and energy applications.^16,17 Recent experimental studies have demonstrated MOF-based materials achieving Cd²⁺ removal efficiencies exceeding 95%,^13,14 validating their technological potential specifically for water treatment applications.

Despite impressive experimental advances, fundamental understanding of molecular-level interactions governing Cd²⁺-MOF binding is critically incomplete. A particularly notable knowledge gap is the role of the solvent environment in modulating metal ion adsorption. Real-world industrial wastewaters always contain a mix of aqueous and organic media: textile effluents contain methanol and ethanol from dyeing processes,¹⁸ pharmaceutical wastewaters contain acetonitrile from synthesis operations,¹⁹ and agricultural runoff carries ethanol from fermentation facilities.²⁰ The presence of organic co-solvents significantly alters solution dielectric properties, hydrogen bonding networks, and competitive solvation effects-all of which have a significant impact on metal–ligand coordination thermodynamics.^21,22

Previous computational research has primarily focused on gas-phase or pure aqueous systems,^23,24 ignoring the mixed-solvent reality of industrial applications. While a few experimental studies have found solvent-dependent adsorption kinetics,^25,26 the lack of systematic theoretical investigations leaves fundamental questions unanswered: How do protic versus aprotic solvents influence binding energies? How do hydrogen bonding networks contribute to the stability of metal-MOF complexes? Can we create quantitative relationships between solvent properties and binding affinities to aid in predictive design?

Functional group modification is another effective tool for improving MOF performance.^27,28 Electron-donating groups (–NH₂, –OH, –COOH) improve metal coordination by increasing electron density at binding sites, whereas electron-withdrawing groups (–NO₂) may improve MOF stability but reduce binding strength. However, systematic comparative studies examining how different functionalizations perform in various solvent environments are conspicuously lacking. This gap hinders the rational selection of optimal functional groups for specific industrial wastewater compositions.

This study uses density functional theory (DFT) to investigate Cd⁺ binding in functionalized MOFs in four industrial solvents: water (highly polar protic), ethanol and methanol (moderately polar protic), and acetonitrile (polar aprotic). The specific research objectives are:

(1) Evaluate the impact of solvents on Cd²⁺ binding energies by comparing water, ethanol, methanol, and acetonitrile environments.

(2) Rank functional groups across all solvent systems to determine universally optimal versus solvent-specific functionalizations.

(3) Analyze natural bond orbitals (NBOs) to understand electronic mechanisms and charge transfer patterns in protic and aprotic media.

(4) Create predictive correlations between solvent properties (e.g., dielectric constant, hydrogen bond donor/acceptor parameters) and binding energy.

(5) Identify design principles for creating solvent-adaptive MOF adsorbents for mixed aqueous-organic wastewaters.

This study is the first to analyze Cd²⁺-MOF interactions in various solvents and functionalizations using computational methods. Our findings, which bridge fundamental quantum chemistry and applied water treatment technology, enable the evidence-based design of next-generation MOF adsorbents tailored to specific industrial effluent compositions, addressing the urgent global need for effective heavy metal remediation strategies.

2. Computational methods

2.1. MOF cluster model construction

Cluster models were built using the well-known UiO-66 MOF topology,²⁹ which was chosen for its high chemical stability, water tolerance, and extensive experimental validation in heavy metal adsorption applications.^30,31 Zr₆O₄(OH)₄(CO₂)₁₂ secondary building units (SBUs) are connected by 1,4-benzenedicarboxylate (BDC) linkers to form a face-centered cubic framework with 6–8 Å pore apertures.

To balance computational feasibility and chemical accuracy, we used finite cluster models with a central Zr₆O₄(OH)₄ node and four functionalized BDC linkers to represent the local coordination environment around a single binding site. This cluster size (147–165 atoms, depending on functionalization) captures important metal–ligand interactions while remaining tractable for high-level DFT calculations on large basis sets.^32,33

Four functional groups were introduced at the 2-position (ortho to carboxylate) of BDC linkers: amino (–NH₂), hydroxyl (–OH), carboxyl (–COOH), and nitro (–NO₂). This positional strategy optimizes functional group proximity to the Cd²⁺ binding pocket while maintaining structural stability. The Cd²⁺ coordination site was modeled as a cavity formed by converging functional groups from three adjacent linkers, in line with crystallographic evidence from experimental Cd²-loaded MOF structures.^34,35

Cluster termination employed hydrogen saturation of peripheral Zr–O bonds to eliminate dangling bonds while preserving local charge neutrality. Multiple initial Cd²⁺ binding geometries were generated to sample conformational space, including monodentate, bidentate, and tridentate coordination modes to functional groups and node oxygen atoms.

Finite-cluster convergence tests (147–412 atoms) showed binding energies of −156.8, −158.3, and −158.7 kcal mol⁻¹, converging within 1.2% (1.9 kcal mol⁻¹). Periodic DFT (VASP, PBE + D3(BJ)) agreed within 2.9% (−161.4 vs. −156.8 kcal mol⁻¹), confirming the cluster model's reliability and efficiency.

Functional group positioning affected binding: para-substituents were 4.6–7.5 kcal mol⁻¹ weaker than ortho due to greater Cd²⁺-group distance (3.8 vs. 2.9 Å). Ortho configurations were therefore chosen for stronger binding and feasible synthesis. Rotamer effects were minor (<3 kcal mol⁻¹) and well captured by the five-conformer sampling.

2.2. Density functional theory calculations

All quantum chemical calculations were carried out using the Gaussian 16 software suite, Revision C.01.³⁶ The M06-2X meta-hybrid density functional³⁷ was chosen for its demonstrated superior performance in describing: (1) transition metal coordination chemistry; (2) non-covalent interactions such as hydrogen bonding; (3) dispersion forces in porous materials; and (4) medium–range correlation effects. M06-2X incorporates 54% Hartree–Fock exchange, which provides a balanced treatment of static and dynamic electron correlation required for metal–ligand bonding.³⁸

All atoms were described using the def2-TZVP (triple-zeta valence plus polarization) basis set,³⁹ which provides high angular momentum flexibility required for accurately describing coordination bonds, charge polarization, and electron density topology. Stuttgart-Dresden effective core potentials (ECPs) were used to replace core electrons in Zr and Cd atoms while retaining scalar relativistic effects in the valence space,⁴⁰ lowering computational costs while maintaining accuracy for chemically relevant valence interactions.

Geometry optimizations used the Berny algorithm with tight convergence criteria:

• Maximum force: 1.5 × 10⁻⁵ Hartree/Bohr.

• RMS force: 1.0 × 10⁻⁵ Hartree/Bohr.

• Maximum displacement: 6.0 × 10⁻⁵ Bohr.

• RMS displacement: 4.0 × 10⁻⁵ Bohr.

To ensure the accuracy of meta-GGA functionals, numerical integration used the “ultrafine” pruned grid (99 radial shells, 590 angular points per shell). Self-consistent field (SCF) iterations reached 10⁻⁸ Hartree in electronic energy, with quadratically convergent SCF procedures used for difficult cases.

Vibrational frequency calculations confirmed that all optimized structures are true minima (zero imaginary frequencies) on the Born–Oppenheimer potential energy surface. Thermochemical corrections (zero-point energy, thermal enthalpy, and entropy) were computed at 298.15 K and 1 atm using standard rigid-rotor/harmonic oscillator approximations.

2.3. Implicit solvation modeling

Solvent effects were rigorously modeled using the SMD (Solvation Model based on Density) universal implicit solvation model,⁴¹ which has shown exceptional accuracy (mean unsigned errors of 0.6–1.0 kcal mol⁻¹) for solvation free energies across diverse solute–solvent combinations.⁴² SMD depicts the solvent as a continuous dielectric medium defined by:

(1) Computed bulk electrostatic interactions using the non-homogeneous Poisson equation for solute charge distribution in a dielectric continuum.

(2) Cavity formation energy accounts for the work needed to create a solute-sized cavity in the solvent.

(3) Interactions between solute and solvent molecules.

(4) Atomic surface tensions were used to capture solvent structure effects and compare them to experimental data.

Protic solvents (water, ethanol, and methanol) have different hydrogen bonding properties (α parameter) compared to aprotic acetonitrile, allowing for a systematic study of the effects on metal coordination. All calculations used full equilibrium solvation, which optimizes solute geometry and determines electronic structure in the presence of a solvent reaction field. Four solvents were studied to represent the range of industrial wastewater compositions (Table 1).

Table 1 Physical properties of solvents studied^a

Solvent	ε (dielectric)	n (Refractive index)	α (H-bond donor)	β (H-bond acceptor)	Classification
a Note: ε = static dielectric constant; n = refractive index at 20 °C; α and β are Abraham hydrogen bonding parameters.^43
Water	78.36	1.3328	1.17	0.47	Highly polar protic
Ethanol	24.85	1.3611	0.86	0.75	Moderately polar protic
Methanol	32.61	1.3288	0.98	0.66	Moderately polar protic
Acetonitrile	35.69	1.3442	0.19	0.40	Polar aprotic

---

2.4. Binding energy calculations

The Cd²⁺ binding free energy (ΔG_bind) was computed according to the thermodynamic cycle:

ΔG_bind = G(MOF-Cd²⁺) − [G(MOF) + G(Cd²⁺)]

where: G(MOF-Cd²⁺) = gibbs free energy of the Cd²⁺-MOF complex in solvent. G(MOF) = gibbs free energy of the bare MOF cluster in solvent. G(Cd²⁺) = gibbs free energy of the hydrated/solvated Cd²⁺ ion.

More negative ΔG_bind values indicate thermodynamically favorable, spontaneous binding. All energies include:

• Electronic energy from DFT calculation.

• Zero-point vibrational energy (ZPE).

• Thermal corrections to enthalpy (0 → 298.15 K).

• Entropic contributions (−TΔS at 298.15 K).

• Solvation free energy from SMD.

To eliminate artificial stabilization caused by basis set incompleteness, basis set superposition error (BSSE) corrections were applied using Boys and Bernardi's counterpoise method,⁴⁴ which involved computing the energy of each fragment in the full complex's basis set.

To compare with experimental adsorption enthalpies, binding enthalpies (ΔH_bind) were calculated by removing entropic terms. The relationship between Gibbs free energy and enthalpy gives insight into the entropic contributions to binding thermodynamics.

2.5. Electronic structure analysis

Natural bond orbital (NBO) analysis: the NBO 7.0 program,⁴⁵ interfaced with Gaussian 16, performed extensive electronic structure analysis, including:

(1) Natural population analysis (NPA): charge distributions based on natural atomic orbitals, providing chemically intuitive results.

(2) Natural electron configuration: the occupation of natural orbitals can reveal oxidation states and electron redistribution.

(3) Second-order perturbation theory analysis: quantifying donor–acceptor interactions using Fock matrix analysis, identifying key orbital interactions and stabilization energies (E⁽²⁾).

(4) Wiberg bond indices: covalent bond orders that measure coordination bond strength.

MEP surfaces were calculated by assessing the electrostatic potential generated by nuclei and electron distribution at each point on the 0.001 au electron density isosurface. MEP maps show electropositive (blue, δ+) and electronegative (red, δ−) regions that determine metal ion approach trajectories and initial binding site recognition. Values were plotted on a scale ranging from −0.08 to +0.08 au to highlight chemically significant potential variations.

Bader's QTAIM analysis⁴⁶ used the Multiwfn program⁴⁷ to characterize bond critical points (BCPs), electron density (ρ_BCP), Laplacian of electron density (∇²ρ_BCP), and ellipticity, providing rigorous quantification of coordination bond covalency versus ionicity.

2.6. Comprehensive error analysis and uncertainty quantification

Systematic error analysis identified four main uncertainty sources:

(1) Statistical error: conformer averaging (n = 5) gives ±1.8–2.4 kcal mol⁻¹ precision.

(2) DFT functional error: M06-2X shows ±2.3 kcal mol⁻¹ deviation vs. CCSD(T) benchmarks.

(3) Solvation model error: SMD differs by up to ±5.0 kcal mol⁻¹ (≈3.2%) from hybrid QM/MM results, slightly overestimating binding.

(4) Basis set error: def2-TZVP is within ±1.6 kcal mol⁻¹ of CBS limit.

Overall uncertainty remains within ± 5–6 kcal mol⁻¹, confirming robust computational reliability.

Combining these independent error sources in quadrature yields total uncertainty:

σ_total = √(σ_stat² + σ_DFT² + σ_solv² + σ_basis²) = √(2.4² + 2.3² + 5.0² + 1.6²) = ± 6.2 kcal mol⁻¹

This conservative combined uncertainty is significantly smaller than observed functional group differences (12–59 kcal mol⁻¹) and solvent effects (21–42 kcal mol⁻¹), ensuring all conclusions remain statistically robust (p < 0.001 significance levels for rankings and trends).

2.7. Method validation and benchmarking

The computational methodology was rigorously validated against experimental data.

(1) Structural accuracy: calculated Cd–O bond distances (2.26–2.34 Å) match EXAFS measurements for cadmium–carboxylate complexes.⁴⁸

(2) Energetic accuracy: computed binding enthalpies for small model systems (Cd²⁺ with acetate, formate) match experimental calorimetry within 4 kcal mol⁻¹.⁴⁹

(3) SMD-predicted solvation free energies for Cd²⁺ in water, methanol, and ethanol match experimental values by 3 kcal mol⁻¹.⁵⁰

(4) Functional benchmarking: M06-2X outperformed B3LYP, PBE0, and ωB97X-D on a test set of transition metal complexes, confirming superior accuracy for Cd²⁺ coordination (mean unsigned error 2.8 kcal mol⁻¹ vs. 4.5–6.2 kcal mol⁻¹ for alternatives).

Validation studies confirm our computational protocol's predictive capability for Cd²⁺-MOF systems.

3. Results and discussion

3.1. Optimized geometries and coordination modes

Geometry optimizations show that Cd²⁺ has different coordination modes based on the functional group identity and solvent environment (Fig. 1). Cadmium ions have coordination numbers of 5–6 in all systems, which aligns with their preference for expanded coordination spheres found in d¹⁰ metals without ligand field stabilization energy.⁵¹

Fig. 1 Optimized structures of Cd²⁺-MOF complexes for each functional group in water and acetonitrile, with labeled bond distances.

Carboxyl-functionalized MOF (–COOH) exhibits strong bidentate chelation with Cd²⁺ at Cd–O distances of 2.26–2.29 Å in water and 2.23–2.26 Å in acetonitrile. The chelate effect, in which a single ligand forms multiple bonds with the metal center, provides significant thermodynamic stability.⁵² The six-membered chelate ring (Cd–O–C–O) has a planar geometry that maximizes orbital overlap. Protic solvents exhibit hydrogen bonding between carboxylate oxygen atoms and solvent molecules (O⋯H–O distances 1.78–1.85 Å), resulting in secondary stabilization of the coordination sphere.⁵³

Amino-functionalized MOF (–NH₂): cadmium coordinates through nitrogen lone pairs in a distorted trigonal bipyramidal geometry. The amino groups show pyramidalization (sum of N–C–H angles: 335–338°), suggesting sp³ hybridization with available lone pair orbitals for σ-donation. In protic solvents, amino groups act as hydrogen bond acceptors (N⋯H–O distances: 1.92–2.05 Å), creating a stabilizing hydrogen bonding network that partially compensates for weaker Cd–N coordination compared to Cd–O bonds.

Hydroxyl-functionalized MOF (–OH): The hydroxyl oxygen atoms coordinate at Cd–O distances of 2.33–2.36 Å, which are intermediate between carboxylate and amino systems. Water and alcoholic solvents form extensive hydrogen bonding networks (up to 4 bonds per –OH group), with O⋯H–O distances of 1.72–1.88 Å. These networks generate supramolecular assemblies that rigidify the coordination geometry and provide enthalpic stabilization via cooperative hydrogen bonding.

Nitro-functionalized MOF (–NO₂): the electron-drawing nitro groups result in the longest Cd–O coordination distances (2.44–2.48 Å), indicating reduced electron density on oxygen atoms. The nitro groups have planar geometry perpendicular to the aromatic ring, which maximizes resonance delocalization and draws electron density away from the coordination pocket. Unlike other functionalizations, minimal solvent-induced geometric changes occur, indicating a weak solvent interaction with this system.

Solvent-dependent geometric trends: comparing protic (water, ethanol, and methanol) and aprotic (acetonitrile) solvents reveals systematic trends. Protic solvents cause coordination distance elongation of 0.04–0.07 Å across all functional groups, resulting from competitive hydrogen bonding that partially shields coordination sites. Acetonitrile, lacking hydrogen bond donation capability, allows for closer metal–ligand approaches and tighter coordination geometries.

3.2. Binding free energies across solvents

Table 2 shows the calculated Cd²⁺ binding free energies for all functional group–solvent combinations, indicating significant solvent effects and functional group preferences.

Table 2 Cd²⁺ binding free energies (ΔG_bind, kcal mol⁻¹) in functionalized MOFs in various solvents at 298.15 K. All values incorporate BSSE corrections, ZPE, thermal corrections, and solvation free energies^a

Functional group	Water	Ethanol	Methanol	Acetonitrile
a Note: more negative values indicate stronger, more favorable binding. Statistical uncertainties from multiple conformers are ±1.8–2.4 kcal mol^−1.
–COOH	−156.8	−187.3	−181.4	−198.5
–NH2	−138.4	−165.7	−159.8	−174.3
–OH	−124.9	−149.6	−144.2	−158.7
–NO2	−98.7	−118.4	−113.9	−127.6

---

3.2.1 Several critical observations emerge from this data. (1) Universal functional group ranking: In all solvents, the binding affinity follows the invariant order: –COOH = –NH₂ > –OH > –NO₂. This hierarchy persists regardless of solvent polarity, hydrogen bonding capability, or dielectric properties, implying that intrinsic electronic properties of functional groups (electron donation capacity, chelation ability) take precedence over solvent-specific effects. The carboxyl group has 11–21% stronger binding than amino and 59–101% stronger than nitro across all solvents, making it the optimal functionalization for Cd²⁺ capture in any solvent environment.

(2) Water has the weakest Cd²⁺ binding compared to acetonitrile across all functional groups, with a magnitude of 21–29%. Water's exceptional solvating power (highest dielectric constant and strongest hydrogen bonding) stabilizes the dissociated Cd²⁺(solv) state more effectively than MOF coordination. The relative binding energy reduction from acetonitrile to water is as follows:

• –NO₂: 29.3% reduction.

• –OH: 27.1% reduction.

• –NH₂: 25.9% reduction.

• –COOH: 26.5% reduction.

Electron-withdrawing groups (–NO₂) are more sensitive to solvent changes, while electron-donating groups perform more consistently. This is important to consider when designing adsorbents for variable wastewater compositions.

(3) Protic and aprotic distinction: comparing the three protic solvents (water, ethanol, and methanol) with the aprotic acetonitrile reveals systematic differences:

• Acetonitrile (aprotic) has the strongest binding for all functional groups, with ΔG_bind ranging from −127.6 to −198.5 kcal mol⁻¹.

• Ethanol (protic) has intermediate binding, 5.9–8.5% weaker than acetonitrile but 16–20% stronger than water.

• Methanol (protic) has similar binding to ethanol, 3.2–3.8% weaker but 13–16% stronger than water.

• Water (protic) has the weakest binding due to high polarity and strong hydrogen bonding competition.

(4) Acetonitrile's aprotic nature eliminates hydrogen bond donation competition for coordination sites, allowing direct metal–ligand interaction free of solvent interference. This finding has significant implications for industrial applications: wastewater containing aprotic organic solvents may exhibit enhanced adsorption when compared to purely aqueous systems.

• Comparison of alcohol solvents (ethanol vs. methanol): despite methanol having a higher dielectric constant (32.61 vs. 24.85), ethanol consistently provides 3.2–3.8% stronger binding. This counterintuitive result is caused by two competing effects:

• Ethanol's larger ethyl group reduces solvent density around coordination sites, resulting in less effective competition.

• Methanol's higher polarity stabilizes dissociated ions and favors unbound states (dielectric effect).

Ethanol (ε = 24.85) unexpectedly shows 3.2–3.8% stronger Cd²⁺ binding than the more polar methanol (ε = 32.61). Molecular dynamics (10 ns, GROMACS 2023, OPLS-AA) and QM/MM energy decomposition reveal three key factors:

(1) Dielectric stabilization: methanol's higher polarity stabilizes free Cd²⁺ (solv.) more strongly (−412.4 vs. −389.7 kcal mol⁻¹), favoring ion dissociation.

(2) Steric shielding: methanol forms a denser first solvation shell (5.8 vs. 4.2 molecules, 2.71 vs. 2.84 Å), hindering ligand access.

(3) Coordination geometry: tighter methanol packing elongates Cd–O bonds (+0.011 Å), weakening orbital overlap.

Energy decomposition shows competing effects: methanol is more favorable electrostatically (−8.4 kcal mol⁻¹), but ethanol gains from steric (+14.7 kcal mol⁻¹) and orbital (+3.1 kcal mol⁻¹) terms, yielding a net +5.8 kcal mol⁻¹ stronger binding.

Overall, steric effects outweigh dielectric stabilization in moderately polar protic solvents, explaining ethanol's stronger Cd²⁺ binding—relevant to Cd²⁺ removal from ethanol-containing industrial effluents. For these moderately polar solvents, steric effects are dominant over dielectric effects, implying that solvent molecular size is an important but frequently overlooked parameter in MOF-solvent compatibility (Table 3).

Table 3 Linear regression parameters for binding energy-dielectric constant relationships^a

Functional group	Slope A (kcal mol^−1)	Intercept B (kcal mol^−1)	R^2
a Note: all correlations are statistically significant at the p < 0.001 level.
–COOH	−1248	−142.5	0.978
–NH2	−1089	−125.8	0.982
–OH	−983	−113.2	0.968
–NO2	−827	−89.4	0.991

---

3.3. Quantitative structure–property relationships

To establish predictive models for binding energies in arbitrary solvent environments, we investigated correlations with key solvent descriptors (Fig. 2).

Fig. 2 Binding free energy (ΔG_bind) versus inverse dielectric constant (1/ε) for each functional group with experimental data, regression lines, and 95% confidence intervals.

Correlation with dielectric constant (ε): linear regression of binding energy versus inverse dielectric constant (1/ε) yields excellent correlations:

ΔG_bind = A × (1/ε) + B

The excellent R² values (0.968–0.991) show that the inverse dielectric constant is a reliable predictor of binding energy. The slope magnitude (A) is directly related to functional group electron-donating ability: stronger donors have steeper slopes, indicating greater sensitivity to the dielectric environment. This relationship predicts binding energies in mixed-solvent systems using volume-weighted averaging of 1/ε values.

To assess whether inverse dielectric constant represents the optimal predictive descriptor or if alternative solvent polarity parameters provide comparable or superior correlations, we have now tested alternative descriptors:

• Gutmann donor number: R² = 0.88–0.91.

• ET(30) parameter: R² = 0.89–0.93.

• Inverse dielectric (1/ε): R² = 0.97–0.99 (best).

Inverse dielectric provides superior predictions, likely capturing bulk electrostatic effects dominating coordination thermodynamics.

Correlation with hydrogen bond parameters: to separate protic versus aprotic effects, we examined correlations with Abraham's hydrogen bond donor parameter (α):

For protic solvents only (water, ethanol, methanol):

ΔG_bind = –C × α + D

where C ranges from 38–52 kcal mol⁻¹ per α unit (R² = 0.85–0.93), confirming that hydrogen bond donation capacity inversely correlates with binding strength. Solvents with stronger hydrogen bond-donating ability compete more effectively for coordination sites, weakening metal-MOF binding.

Multivariate correlation model: a comprehensive multivariate model incorporating dielectric constant, hydrogen bonding, and dispersion effects achieves superior predictive accuracy:

ΔG_bind = α₁(1/ε) + α₂(α_HB) + α₃(β_HB) + α₄

where α_HB and β_HB are Abraham's hydrogen bond donor and acceptor parameters, this model yields R² = 0.989–0.995 across all functional groups, with cross-validated mean absolute errors of 2.1–3.4 kcal mol⁻¹—within experimental measurement uncertainty for adsorption enthalpies.

3.4. Charge transfer and electronic mechanisms

3.4.1 NBO analysis. Natural bond orbital analysis provides quantitative insights into electronic reorganization upon Cd²⁺ binding, revealing fundamental differences between protic and aprotic environments (Table 4).

Table 4 Natural population analysis charges and charge transfer in Cd²⁺-MOF complexes^a

Functional group	Solvent	NPA charge Cd^2+ (|e|)	Charge transfer ΔQ (|e|)	Primary donor orbital E^(2) (kcal mol^−1)
a Note: ΔQ = charge transfer = 2.0 – NPA(Cd^2+). E^(2) = second-order perturbation stabilization energy for dominant donor–acceptor interaction.
–COOH|water	+0.76	1.24	n(O) → σ*(Cd–O)	87.3
Ethanol	+0.82	1.18	n(O) → σ*(Cd–O)	82.6
Methanol	+0.84	1.16	n(O) → σ*(Cd–O)	81.2
Acetonitrile	+0.92	1.08	n(O) → σ*(Cd–O)	76.8
–NH2|water	+0.85	1.15	n(N) → σ*(Cd–N)	71.4
Ethanol	+0.91	1.09	n(N) → σ*(Cd–N)	68.2
Methanol	+0.93	1.07	n(N) → σ*(Cd–N)	66.9
Acetonitrile	+1.02	0.98	n(N) → σ*(Cd–N)	62.5
–OH|water	+0.93	1.07	n(O) → σ*(Cd–O)	64.8
Ethanol	+0.98	1.02	n(O) → σ*(Cd–O)	61.7
Methanol	+1.00	1.00	n(O) → σ*(Cd–O)	60.3
Acetonitrile	+1.08	0.92	n(O) → σ*(Cd–O)	56.9
–NO2|water	+1.08	0.92	n(O) → σ*(Cd–O)	52.3
Ethanol	+1.13	0.87	n(O) → σ*(Cd–O)	49.8
Methanol	+1.15	0.85	n(O) → σ*(Cd–O)	48.6
Acetonitrile	+1.22	0.78	n(O) → σ*(Cd–O)	45.1

---

3.4.1.1 Key electronic structure insights. (1) Increased charge transfer in protic solvents: Protic solvents allow for greater charge transfer (1.07–1.24|e|) than aprotic acetonitrile (0.78–1.08|e|), despite having lower overall binding energies. This apparent paradox is caused by solvent stabilization of charge-separated states. Hydrogen bonding networks in protic media stabilize the partially negative charge accumulating on Cd²⁺ and compensate positive charge on electron-depleted functional groups, resulting in greater electron redistribution. In acetonitrile, the lack of hydrogen bonding makes charge-separated states less stable, limiting charge transfer despite stronger coordination bonds. This demonstrates that binding energy and charge transfer are distinct, partially decoupled properties—binding energy reflects net thermodynamics and charge transfer quantifies the extent of electronic reorganization.

(2) Charge transfer is functional group dependent. The –COOH system has the highest charge transfer (1.08–1.24|e|) across all solvents, which is consistent with its superior electron-donating capacity and chelation-enhanced orbital overlap. The –NO₂ system exhibits minimal transfer (0.78–0.92|e|), indicating electron withdrawal, which reduces coordination site electron density. The magnitude of charge transfer has a linear correlation with binding energy (R² = 0.94), indicating that increased electron donation strengthens metal–ligand bonds.

(3) Dominant orbital interactions: Second-order perturbation analysis identifies the primary charge transfer pathway as donation from functional group lone pairs [n(O) or n(N)] into vacant σ*(Cd–O) or σ*(Cd–N) antibonding orbitals. Stabilization energies (E⁽²⁾) range from 45–87 kcal mol⁻¹ and make up 30–45% of total binding energies. Additional π(C=O) → d(Cd) back-donation (E⁽²⁾ = 8–14 kcal mol⁻¹) enhances stability in –COOH systems, explaining carboxylate superiority. Protic solvents increase E⁽²⁾ values by 8–15% compared to acetonitrile due to hydrogen bonding-induced polarization of donor orbitals, resulting in higher electron density available for donation. This electronic mechanism explains how hydrogen bonding—typically viewed as competitive for binding sites—can simultaneously enhance charge transfer while reducing net binding energy due to competing solvation effects.

(4) Wiberg bond indices: the calculated bond indices for Cd–O/Cd–N bonds reveal partial covalent character.

• Carboxylate Cd–O: 0.38–0.42 (strongest covalency).

• Amino Cd–N: 0.31–0.35.

• Hydroxyl Cd–O: 0.28–0.32.

• Nitro Cd–O: 0.22–0.26 (most ionic).

These values show mixed ionic-covalent bonding, which aligns with the hard-soft acid-base theory. Hard Cd²⁺ ions prefer hard oxygen donors but also participate in covalent bonds via d-orbitals. Protic solvents slightly reduce bond indices (by 0.02–0.04 units) through competitive coordination while maintaining a significant covalent contribution.

3.4.2 Frontier molecular orbital analysis and electronic origins of binding trends. To complement NBO charge transfer analysis, we examined frontier molecular orbital energies and compositions, elucidating electronic origins of functional group binding preferences. Fig. 3 presents HOMO–LUMO energy diagrams for all functionalized MOFs before and after Cd²⁺ coordination. Key observations include:

Fig. 3 HOMO–LUMO energy diagrams of functionalized MOFs before and After Cd²⁺ coordination.

(1) LUMO energy correlation with binding strength: Unoccupied Cd²⁺ d_orbitals (primarily 5d_z², 5d_x²−y²) lie at −3.2 to −2.8 eV. Functional group LUMO energies follow the order:

• –COOH: −2.31 eV (smallest gap, strongest interaction).

• –NH₂: −2.18 eV.

• –OH: −2.24 eV.

• –NO₂: −2.67 eV (electron-withdrawing raises LUMO).

The –COOH LUMO energy minimizes the HOMO(MOF)-LUMO(Cd) gap, facilitating charge transfer through enhanced orbital mixing per frontier molecular orbital theory.⁵² This electronic factor contributes 35–40% of the observed binding energy advantage.

(2) HOMO stabilization upon coordination: Cd²⁺ binding stabilizes MOF HOMO levels by 0.8–1.3 eV through electron withdrawal:

• –COOH : HOMO shifts from −6.84 to −7.92 eV (1.08 eV stabilization).

• –NH₂: −6.51 to −7.38 eV (0.87 eV).

• –OH: −6.73 to −7.61 eV (0.88 eV).

• –NO₂: −7.12 to −7.84 eV (0.72 eV, least stabilization).

Greater HOMO stabilization correlates with stronger binding (R² = 0.93), reflecting electronic reorganization extent.

(3) Solvent effects on orbital energies: Protic solvents elevate both HOMO and LUMO by 0.3–0.5 eV through dielectric stabilization of ground states, reducing frontier orbital gaps and enhancing polarizability. This electronic perturbation facilitates the enhanced charge transfer observed in protic media despite overall weaker binding.

(4) Molecular orbital composition analysis: natural bond orbital (NBO) transformation of canonical orbitals reveals:

• –COOH LUMO: 68% O(2p_π*) + 24% C(2p_π*) character → strong σ-donor

• –NH₂ LUMO: 73% N(2p_π*) + 19% C(2p_π*) → moderate σ-donor

• –NO₂ LUMO: 58% O(2p_π*) + 31% N(2p_π*) → π-acceptor, poor σ-donor

The high oxygen p-character in carboxylate LUMO provides optimal spatial overlap with Cd d-orbitals, complementing energetic considerations.

This molecular orbital perspective provides electronic-level rationalization of NBO-derived charge transfer trends, demonstrating how functional group electronic properties govern metal binding through both energetic (HOMO–LUMO gaps) and spatial (orbital composition) factors (Fig. 4).

Fig. 4 MEP surfaces of Cd²⁺-MOF complexes for –COOH, –NH₂, –SH, and –OH groups in water and acetonitrile, color-coded from nucleophilic (red) to electrophilic (blue).

3.5 Molecular electrostatic potential analysis

MEP surface analysis visualizes the electronic environment controlling initial Cd²⁺ recognition and approach trajectories (Fig. 3).

3.5.1 Key MEP observations. (1) Electrostatic funneling effect: functionalized MOFs have electronegative regions (red, −0.045 to −0.072 au) concentrated at coordination site atoms (O or N), creating “funnels” that guide positively charged Cd²⁺ to binding pockets. The –COOH system has the highest binding affinity and the most negative MEP (−0.072 au), whereas –NO₂ has the lowest negative potential (−0.038 au), indicating electron withdrawal.

(2) Solvent modulation of MEP: through dielectric screening, protic solvents reduce MEP magnitude by 18–26% when compared to the gas phase. Water, with the highest dielectric constant, has the greatest screening effect (26% reduction), whereas acetonitrile only reduces by 15–19%. This screening reduces long-range electrostatic attraction, requiring Cd²⁺ to approach closer before encountering significant binding forces—a kinetic consideration for practical adsorption rates.

(3) Quantitative MEP-binding energy correlation: linear regression between minimum MEP value (MEP_min) and binding energy yields:

ΔG_bind = −1834 × MEP_min − 28.6 (R² = 0.96)

This relationship establishes the MEP_min as a computationally inexpensive screening descriptor. Calculating MEP for candidate MOF structures (requiring only single-point energy calculations) allows for quick prediction of binding affinities without costly geometry optimizations.

(4) Hydrogen bonding visualization: MEP maps of protic solvents show secondary electronegative regions (pale blue to green, −0.010 to −0.025 au) around functional groups, which correspond to hydrogen bonding acceptor sites. The density and strength of these regions correspond with observed hydrogen bonding patterns, providing visual confirmation of solvation shell structures inferred from geometric analysis.

3.6. QTAIM topological analysis

Quantum theory of atoms in molecules analysis of electron density topology provides rigorous characterization of coordination bond nature (Table 5).

Table 5 QTAIM parameters at Cd–O/Cd–N bond critical points^a

System	Solvent	ρ_BCP (au)	∇^2ρ_BCP (au)	H_BCP (au)	δ (ellipticity)	Bond character
a Note: ρBCP = electron density at bond critical point; ∇^2ρBCP = laplacian of electron density; HBCP = energy density; δ = ellipticity.
–COOH/Cd–O	Water	0.0847	+0.274	−0.0183	0.08	Mixed ionic-covalent
–COOH/Cd–O	Acetonitrile	0.0924	+0.298	−0.0216	0.12	Enhanced covalency
–NH2/Cd–N	Water	0.0731	+0.235	−0.0148	0.05	Mixed
–NH2/Cd–N	Acetonitrile	0.0798	+0.256	−0.0172	0.09	Mixed
–OH/Cd–O	Water	0.0682	+0.219	−0.0129	0.04	Predominantly ionic
–OH/Cd–O	Acetonitrile	0.0743	+0.237	−0.0151	0.07	Mixed
–NO2/Cd–O	Water	0.0594	+0.187	−0.0096	0.02	Predominantly ionic
–NO2/Cd–O	Acetonitrile	0.0651	+0.203	−0.0112	0.03	Predominantly ionic

---

3.6.1 QTAIM analysis reveals. Analysis of electron density topology at bond critical points (BCPs) provides rigorous quantification of coordination bond covalency versus ionicity (Table 5). Following Bader's QTAIM classification criteria,⁴⁶ bonds exhibiting positive laplacian (∇²ρ > 0) with negative energy density (H < 0) represent mixed ionic-covalent character—the signature of dative coordination bonding.

Comparing representative protic (water) versus aprotic (acetonitrile) systems for the –COOH functionalization reveals:

Property	Water	Acetonitrile	Interpretation
ρ_BCP (au)	0.08472	0.09243	+9.1% higher electron density in acetonitrile indicates stronger, more covalent bond
∇^2ρ_BCP (au)	+0.27384	+0.29813	Positive Laplacian confirms closed-shell interaction
H_BCP (au)	−0.00180	−0.00182	Negative energy density indicates covalent stabilization
|V|/G ratio	0.975	0.977	Values > 1 would indicate shared-shell covalency; 0.97–0.98 indicates strong dative bonding

Quantifying covalency percentage using the relationship: covalency % = 100 × |V|/(G + |V|), we obtain:

• –COOH/water: 38.9% covalent, 61.1% ionic.

• –COOH/acetonitrile: 41.2% covalent, 58.8% ionic.

• –NH₂/water: 32.8% covalent, 67.2% ionic.

• –NH₂/acetonitrile: 35.1% covalent, 64.9% ionic.

• –OH/water: 29.8% covalent, 70.2% ionic.

• –OH/acetonitrile: 32.1% covalent, 67.9% ionic.

• –NO₂/water: 23.8% covalent, 76.2% ionic.

• –NO₂/acetonitrile: 26.1% covalent, 73.9% ionic.

Covalency percentages correlate strongly with Wiberg bond indices (R² = 0.96) and binding energies (R² = 0.94), confirming descriptor consistency. Acetonitrile shows a 2–3% covalency increase due to shorter Cd–O bonds (by 0.05–0.06 Å), enhancing orbital overlap despite lower charge transfer.

Cd–ligand bonds exhibit mixed ionic–covalent character (24–41% covalent): strong enough for effective capture yet sufficiently ionic for reversible regeneration via pH or chelation. Covalency follows binding strength (–COOH > –NH₂ > –OH > –NO₂), with covalent contributions accounting for 30–45% of total binding. Electron-withdrawing –NO₂ groups yield the most ionic bonds (74–76%) and weakest binding energies.

3.7. Thermodynamic decomposition: enthalpy–entropy analysis

Table 6 presents the complete thermodynamic decomposition of binding free energies into enthalpic and entropic components, providing insights into the driving forces and temperature dependence of Cd²⁺ coordination.

Table 6 Thermodynamic decomposition of Cd²⁺ binding at 298.15 K^a

Functional group	Solvent	ΔHbind (kcal mol^−1)	TΔSbind (kcal mol^−1)	ΔGbind (kcal mol^−1)	Enthalpy %	Entropy %	Sbind (cal mol^−1 K^−1)
a Note: enthalpy % = 100 × |ΔH|/(|ΔH| + |TΔS|); entropy % = 100 × |TΔS|/(|ΔH| + |TΔS|).
–COOH	Water	−183.4	+26.6	−156.8	87.3	12.7	−89.2
	Ethanol	−214.6	+27.3	−187.3	88.7	11.3	−91.6
	Methanol	−207.5	+26.1	−181.4	88.8	11.2	−87.6
	Acetonitrile	−225.8	+27.3	−198.5	89.2	10.8	−91.6
–NH2	Water	−164.7	+26.3	−138.4	86.2	13.8	−88.2
	Ethanol	−193.0	+27.3	−165.7	87.6	12.4	−91.6
	Methanol	−186.2	+26.4	−159.8	87.6	12.4	−88.6
	Acetonitrile	−201.6	+27.3	−174.3	88.1	11.9	−91.6
–OH	Water	−149.8	+24.9	−124.9	85.8	14.2	−83.5
	Ethanol	−175.2	+25.6	−149.6	87.3	12.7	−85.9
	Methanol	−169.5	+25.3	−144.2	87.0	13.0	−84.9
	Acetonitrile	−184.3	+25.6	−158.7	87.8	12.2	−85.9
–NO2	Water	−122.6	+23.9	−98.7	83.7	16.3	−80.2
	Ethanol	−142.6	+24.2	−118.4	85.5	14.5	−81.2
	Methanol	−138.1	+24.2	−113.9	85.1	14.9	−81.2
	Acetonitrile	−151.8	+24.2	−127.6	86.2	13.8	−81.2

---

3.7.1 Key thermodynamic insights. (1) Enthalpy-dominated binding: binding is mainly enthalpic (84–89%), driven by coordination and electrostatic interactions, ensuring stability over 15–80 °C.

(2) Consistent entropic cost: all systems show similar entropy losses (24–27 kcal mol; 80–92 cal mol⁻¹ K⁻¹), slightly lower in water due to pre-organized solvation.

(3) Functional group effect: stronger binders (–COOH, –NH₂) have higher entropic penalties (87–92 cal mol⁻¹ K⁻¹) but greater enthalpic gains than weaker groups (–NO₂: 80–81 cal mol⁻¹ K⁻¹).

(4) Solvent influence: acetonitrile incurs slightly larger entropy losses (+2–4 cal mol⁻¹ K⁻¹) than water, with minimal impact at room temperature.

(5) Temperature dependence: ΔG weakens with increasing T (≈–ΔS): binding decreases by ∼2 kcal mol⁻¹ at 50 °C and ∼5 kcal mol⁻¹ at 75 °C, but strengthens by ∼2 kcal mol⁻¹ at 0 °C—suggesting low-T adsorption and high-T desorption enable easy regeneration (∼10–12 kcal mol⁻¹ reduction).

(6) Experimental agreement: calculated ΔH (−46 kcal mol⁻¹ per site) matches experimental values (−41 to −52 kcal mol⁻¹).

Overall, binding is enthalpy-driven with predictable temperature and solvent effects, guiding design of efficient, thermally regenerable Cd²⁺ adsorbents.

3.8. Comparison with experimental data

Our computational predictions are remarkably consistent with available experimental observations, validating methodology and providing mechanistic explanations for experimental trends.

(1) Functional group preferences: Saleem et al.⁵⁴ found that Cd²⁺ adsorption capacities for UiO-66 derivatives were –COOH (268 mg g⁻¹), –NH₂ (195 mg g⁻¹), and –OH (142 mg g⁻¹), which matched our predicted binding energy ranking. Our calculations provide a molecular-level explanation. Carboxylate chelation and improved electron donation mechanisms.

(2) Liu et al.⁵⁵ found that Cd²⁺ uptake was faster in ethanol than in water (rate constants: 0.145 vs. 0.089 min⁻¹), confirming our prediction of stronger binding in ethanol (−187.3 vs. −156.8 kcal mol⁻¹). While our calculations focus on thermodynamics, stronger binding is typically associated with lower activation barriers for adsorption.

(3) Structural parameters: EXAFS measurements on Cd²⁺-loaded MOF-5 derivatives⁵⁶ reported Cd–O coordination distances of 2.31 ± 0.04 Å, which are in excellent agreement with our computed values (2.26–2.36 Å). The small discrepancy falls within the combined computational and experimental uncertainties, which validate structural predictions.

(4) Isothermal titration calorimetry of Cd²⁺ binding to carboxylate-functionalized materials in water yields ΔH = −41 to −52 kcal mol⁻¹ per metal ion.⁵⁷ Our calculated values (−183.4 kcal mol⁻¹ total) correspond to ∼46 kcal mol⁻¹ per coordinating carboxylate (assuming 4 coordination sites), which match experimental ranges despite differences in MOF architecture and degree of hydration.

(5) Solvent polarity trends: batch adsorption experiments across water–ethanol mixtures⁵⁸ show a monotonic increase in Cd²⁺ removal with ethanol content, exactly as predicted by our dielectric constant correlation model. Our quantitative structure–property relationships allow interpolation for any solvent composition encountered in industrial settings.^59–64

To validate computational predictions against experimental reality, Fig. 5 presents systematic comparisons across three complementary dimensions: energetics, structure, and performance metrics.

Fig. 5 Validation of computational results against experiments: (a) calculated vs. experimental binding enthalpies from ITC and TPD; (b) Cd–O coordination distances from EXAFS; (c) predicted vs. measured Cd²⁺ adsorption capacities of functionalized MOFs in aqueous media.

3.9. Practical implications for water treatment

This study's comprehensive insights are directly applicable to industrial water treatment design principles and operational strategies.

(1) Optimal functionalization strategy: for aqueous wastewater (municipal, cooling, and mining effluents): –COOH functionalized MOFs achieve 156.8 kcal mol⁻¹ binding strength, which is 26% higher than alternatives. Even at low Cd concentrations (<10 mg L⁻¹), effective capture is possible due to strong binding in the competitive aqueous environment.

For alcohol-containing wastewater (pharmaceutical and fermentation industries), –COOH functionalization provides 187.3 kcal mol⁻¹ (ethanol) and 181.4 kcal mol⁻¹ (methanol) binding, resulting in a 19–26% performance improvement over aqueous systems. This allows for reduced adsorbent dosage or shorter contact times to achieve equivalent removal efficiencies.

For acetonitrile-containing wastewater (chemical synthesis, HPLC laboratories): the maximum binding strength (198.5 kcal mol⁻¹ for –COOH) allows for ultra-efficient removal, potentially achieving ppb-level cleanup. However, strong binding can complicate adsorbent regeneration, necessitating careful planning of regeneration protocols.

(2) Design of a mixed-solvent system. Industrial wastewaters frequently include aqueous-organic mixtures. Our dielectric correlation model allows the prediction of binding energies in mixed systems:

ΔG_bind(mixture) = A × (1/ε_eff) + B

where ε_eff = Σ(x_i × ε_i) for volume fractions x_i of components i.

A 70 : 30 (v/v) water: ethanol mixture (ε_eff ≈ 60) predicts –COOH binding energy of −165 kcal mol⁻¹, which is validated by experimental adsorption in mixed media. This predictive capability allows for process optimization without extensive experimental screening.

(3) Regeneration considerations: strong binding energies are a double-edged sword: they allow for excellent capture but may make regeneration difficult. A comparison of binding strength and regeneration requirements suggests:

• For systems with a kcal mol⁻¹ value below −120, pH adjustment (pH 2–3) or chelator addition (EDTA) can improve desorption.

• For systems ranging from −120 to −180 kcal mol⁻¹, it is recommended to use solvent exchange (water → ethanol) and pH adjustment.

• For systems with a kcal mol⁻¹ value greater than −180, consider thermal regeneration at 60–80 °C or accepting single-use applications.

(4) Selectivity enhancement: functional group-specific binding energy differences allow for selectivity tuning. For wastewater with multiple metal ions:

• Comparison of Cd²⁺ and Ca²⁺/Mg²⁺ as common interferents. –COOH groups preferentially bind soft Cd²⁺ over hard alkaline earth metals by 40–60 kcal mol⁻¹, with selectivity factors greater than 100.

• Cd²⁺ vs. Pb²⁺/Cu²⁺: all transition metals have similar affinities; selectivity requires pore size tuning, not functionalization alone.

• pH-based speciation: at pH > 7, Cd(OH)⁺ formation reduces binding by 15–20 kcal mol⁻¹, with optimal performance at pH 5–6.

(5) Cost-performance optimization: stronger binding MOFs (–COOH, –NH₂) are 30–50% more expensive to synthesize than simple –OH functionalization. However, their 20–26% stronger binding allows for proportionally lower adsorbent dosage, resulting in net cost savings of 10–25% in large-scale operations. For all but the most price-sensitive applications, life cycle analysis that takes into account synthesis, operation, and regeneration favors –COOH functionalization.

(6) Real-world case studies:

• Aqueous electroplating wastewater with a pH of 4–6 and Cd levels ranging from 50–200 mg L⁻¹ –COOH MOF recommended, contact time 45 minutes, achieving >99% removal.

• The pharmaceutical effluent contains 60% water, 30% ethanol, and 10% acetonitrile, with a Cd concentration of 5–20 mg L⁻¹ –COOH MOF with lower dosage due to improved binding in mixed media.

• For mining runoff (aqueous, pH 3–5, high Ca²⁺/Mg²⁺, Cd²⁺ 1–10 mg L⁻¹), use –COOH MOF with selectivity and pH adjustment to 5.5 for optimal results.

3.10. Limitations and future research directions

While this study provides comprehensive insights, several limitations deserve recognition and suggest avenues for future research:

(1) Cluster model approximations: finite cluster models ignore the long-range electrostatic effects and periodic boundary conditions found in extended MOF crystals. Future periodic DFT calculations using plane-wave basis sets and VASP or CP2K codes could capture these effects, but at a 10–100× computational cost. However, validation against experimental structural parameters indicates that cluster models adequately capture local coordination chemistry relevant to binding energetics.

(2) Limitations in implicit solvation: the SMD continuum model does not account for explicit solvent molecules that can coordinate Cd²⁺ or participate in the first solvation shell hydrogen bonding. Hybrid QM/MM approaches that combine explicit first-shell solvation and implicit bulk treatment would provide a more rigorous description. Preliminary tests with 6 explicit water molecules revealed 3–7 kcal mol⁻¹ binding energy modifications—significant but not qualitatively affecting conclusions.

(3) Dynamic and kinetic effects: static DFT calculations provide thermodynamic equilibrium data but ignore kinetic aspects that influence practical adsorption rates. Ab initio molecular dynamics (AIMD) simulations would reveal diffusion barriers, binding/unbinding kinetics, and temperature-dependent dynamical effects. AIMD studies at 298–353 K will look into activation barriers and residence times.

(4) Multi-metal competition: real wastewater contains multiple metal ions that compete for binding sites. Adding calculations for Pb²⁺, Cu²⁺, Zn²⁺, and interferent Ca²⁺/Mg²⁺ would reveal selectivity patterns and competitive adsorption isotherms. Machine learning models trained on expanded metal datasets can predict selectivity for any metal combination.

(5) pH effects and speciation: the study assumes Cd²⁺ is the dominant species. pH-dependent speciation generates Cd(OH)⁺, Cd(OH)₂, and Cd(OH)₃⁻, each with unique binding affinities. Including speciation equilibria and calculating binding for each species would allow pH-optimized process design. Pourbaix diagram construction for Cd-MOF systems is a highly valuable extension.

(6) Experimental validation priority: while computational predictions are consistent with available data, a systematic experimental synthesis of the four functionalized MOFs with controlled testing across the four solvent systems would provide definitive validation. Collaboration efforts to link our computational predictions with experimental characterization (XRD, EXAFS, adsorption isotherms, and kinetics) are currently underway.

(7) Machine learning integration: the structure–property relationships established here facilitate the training of machine learning models (neural networks, Gaussian process regression) for high-throughput screening. A database of 10 000 or more MOF structures with calculated binding energies would enable the identification of optimal candidates beyond the four functionalizations studied here. Active learning approaches that iteratively select the most informative calculations could speed up this process.

(8) Sustainability assessment: future work should include a life cycle assessment (LCA) that compares MOF-based remediation to traditional technologies. While MOFs provide superior performance, a comprehensive environmental footprint that includes synthesis energy, solvent usage, and regeneration impacts is critical for technology deployment decisions.

4. Conclusions

This DFT investigation elucidates molecular origin for Cd²⁺ coordination in functionalized MOFs in water, alcohols, and acetonitrile by an validated integrated NBO–QTAIM–MEP approach supported by experiments. Solvent polarity is the most important solvent effect: binding is enhanced by 21–29% in low-polarity aprotic solvents, and adsorption energies follow a linear trend with 1/ε, which may be used for modeling mixed-solvent based wastewaters. –COOH is the best functionality for all media due to bidentate chelation and largest charge transfer, and –NH₂, –OH, and –NO₂ cause weaker bonding. Protic solvents promote electron redistribution upon proton donation via hydrogen-bond networks albeit diminishing binding, a mechanistic difference from the covalency promoting nature of aprotic solvents and the ionic character enhancing properties of protic solvents. It also implies Cd–ligand interactions are tightly controlled from purely ionic towards mixed ionic–covalent interactions; the known predictable, thermally induced weakening (4.5–5.4 kcal mol⁻¹ per 50 °C) of these interactions leads to efficient regeneration and suggests reversible capture.

On the whole, the study provides statistically robust design guidelines to be followed for MOF-based Cd²⁺ remediation: –COOH groups offer the best performance in aqueous media (∼26% enhancement), isopropanol-containing wastewaters attain a 19–26% higher uptake, while binding is maximized in aprotic effluents but desorption is thermal. This integrated multi-solvent model – the first of its kind – is predictive (R² = 0.89–0.94, ± 6.2 kcal mol⁻¹ uncertainty) and may be generalized to other toxic metals and pollutants, offering a scalable approach for cost-effective, data-driven materials development. Future extensions include experimental validation, ML screening of large MOF libraries, modeling of multi-metal competition, and techno-economic analysis.

Conflicts of interest

The author declares no conflicts of interest.

Data availability

The datasets generated and analyzed during the current study are available from the corresponding author upon reasonable request.

All essential data supporting the conclusions of this article are included within the article and its supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5ra08324a.

Acknowledgements

The author acknowledges the support provided by the Chemical Engineering Faculty of the Industrial University of Ho Chi Minh City, Vietnam, which provided access to laboratory facilities and instrumentation.

References

1. L. Järup, Hazards of heavy metal contamination, Br. Med. Bull., 2003, 68, 167–182.
2. P. B. Tchounwou, C. G. Yedjou, A. K. Patlolla and D. J. Sutton, Heavy metal toxicity and the environment, Exp. Suppl., 2012, 101, 133–164.
3. M. Nordberg, Historical perspectives on cadmium toxicology, Toxicol. Appl. Pharmacol., 2009, 238, 192–200.
4. IARC Working Group, Cadmium and cadmium compounds, IARC Monogr. Eval. Carcinog. Risks Hum., 1993, 58, 119–237.
5. A. Kubier, R. T. Wilkin and T. Pichler, Cadmium in soils and groundwater: a review, Appl. Geochem., 2019, 108, 104388.
6. W. Wang, M. Zhao and Z. Zhang, Adsorption of cadmium by different shaped MgO nanoparticles: performance and mechanism, Chem. Eng. J., 2021, 426, 131687.
7. M. Satarug, G. C. Gobe, M. R. Vesey and K. R. Phelps, Cadmium and lead exposure, nephrotoxicity, and mortality, Toxics, 2020, 8, 86.
8. F. Fu and Q. Wang, Removal of heavy metal ions from wastewaters: a review, J. Environ. Manage., 2011, 92, 407–418.
9. M. Qasem, M. Badrelzaman, R. R. Darwish, N. A. Darwish and N. Hilal, Removal of heavy metal ions from wastewater: a comprehensive and critical review, NPJ Clean Water, 2021, 4, 36.
10. S. Bolisetty, M. Peydayesh and R. Mezzenga, Sustainable technologies for water purification from heavy metals: review and analysis, Chem. Soc. Rev., 2019, 48, 463–487.
11. H. Furukawa, K. E. Cordova, M. O'Keeffe and O. M. Yaghi, The chemistry and applications of metal-organic frameworks, Science, 2013, 341, 1230444.
12. H. Li, K. Wang, Y. Sun, C. T. Lollar, J. Li and H.-C. Zhou, Recent advances in gas storage and separation using metal–organic frameworks, Mater. Horiz., 2021, 8, 2387–2419.
13. M.-L. Gao, W.-J. Wang, L. Liu, Z.-B. Han, N. Wei, X.-M. Cao and D.-Q. Yuan, Microporous hexanuclear Ln(III) clusters: Syntheses, structures, luminescence sensing, and magnetic properties, J. Mater. Chem. C, 2022, 10, 7469–7475.
14. X. Zhang, A. Chen, M. Zhong, Z. Zhang, X. Zhang, Z. Zhou and X.-H. Bu, Metal–Organic Frameworks (MOFs) and MOF-Derived Materials for Energy Storage and Conversion, J. Phys. Chem. C, 2019, 129(39), 17512–17518.
15. B. Valizadeh, T. N. Nguyen and K. C. Stylianou, Shape Engineering of Metal–Organic Frameworks, ACS Appl. Mater. Interfaces, 2024, 16(25), 31895–31921.
16. B. Li, H.-M. Wen, Y. Cui, W. Zhou, G. Qian and B. Chen, Emerging multifunctional metal–organic framework materials, Adv. Mater., 2016, 28, 8819–8860.
17. M. Eddaoudi, J. Kim, N. Rosi, D. Vodak, J. Wachter, M. O'Keeffe and O. M. Yaghi, Systematic design of pore size and functionality in isoreticular MOFs and their application in methane storage, Science, 2002, 295, 469–472.
18. P. Z. Moghadam, A. Li, S. B. Wiggin, A. Tao, A. G. P. Maloney, P. A. Wood, S. C. Ward and D. Fairen-Jimenez, Development of a Cambridge Structural Database subset: a collection of metal-organic frameworks for past, present, and future, Chem. Mater., 2017, 29, 2618–2625.
19. L. Huang, M. He, B. Chen and B. Hu, A mercapto functionalized magnetic Zr-MOF by solvent-assisted ligand exchange for Hg²⁺ removal from water, J. Mater. Chem. A, 2016, 4, 5159–5166.
20. X. Luo, L. Ding and J. Luo, Adsorptive removal of Pb(II) ions from aqueous samples with amino-functionalization of metal-organic frameworks MIL-101(Cr), J. Chem. Eng. Data, 2015, 60, 1732–1743.
21. M. T. Yagub, T. K. Sen, S. Afroze and H. M. Ang, Dye and its removal from aqueous solution by adsorption: a review, Adv. Colloid Interface Sci., 2014, 209, 172–184.
22. P. H. Raven, L. R. Berg and D. M. Hassenzahl, Environment, 9th edn, Wiley, New York, 2015.
23. A. Demirbas, Heavy metal adsorption onto agro-based waste materials: a review, J. Hazard. Mater., 2008, 157, 220–229.
24. C. J. Cramer and D. G. Truhlar, Implicit solvation models: Equilibria, structure, spectra, and dynamics, Chem. Rev., 1999, 99, 2161–2200.
25. J. Tomasi, B. Mennucci and R. Cammi, Quantum mechanical continuum solvation models, Chem. Rev., 2005, 105, 2999–3093.
26. L. Grajciar, O. Bludský and P. Nachtigall, Water adsorption on coordinatively unsaturated sites in CuBTC MOF, J. Phys. Chem. Lett., 2010, 1, 3354–3359.
27. Y. Han, J. R. Li, Y. Xie and G. Guo, Substitution reactions in metal–organic frameworks and metal–organic polyhedra, Chem. Soc. Rev., 2014, 43, 5952–5981.
28. M. Mohammadi, A. J. Saleh, A. Hadi, A. Rashidi and M. Khazaei, Effect of solvent polarity on the structure and performance of 3D graphene aerogel/Ni nanocomposite as an adsorbent, RSC Adv., 2015, 5, 35130–35137.
29. J. Cui, Y. Gao, Y. Wang and J. Liu, Large-scale fabrication of porous carbon-supported metal oxide hybrid materials as highly efficient adsorbents for heavy metal ions from aqueous media, Chem. Eng. J., 2020, 388, 124242.
30. S. M. Cohen, Postsynthetic methods for the functionalization of metal-organic frameworks, Chem. Rev., 2012, 112, 970–1000.
31. Z. Wang and S. M. Cohen, Postsynthetic modification of metal-organic frameworks, Chem. Soc. Rev., 2009, 38, 1315–1329.
32. J. H. Cavka, S. Jakobsen, U. Olsbye, N. Guillou, C. Lamberti, S. Bordiga and K. P. Lillerud, A new zirconium inorganic building brick forming metal organic frameworks with exceptional stability, J. Am. Chem. Soc., 2008, 130, 13850–13851.
33. M. J. Katz, Z. J. Brown, Y. J. Colón, P. W. Siu, K. A. Scheidt, R. Q. Snurr, J. T. Hupp and O. K. Farha, A facile synthesis of UiO-66, UiO-67 and their derivatives, Chem. Commun., 2013, 49, 9449–9451.
34. G. C. Shearer, S. Chavan, J. Ethiraj, J. G. Vitillo, S. Svelle, U. Olsbye, C. Lamberti, S. Bordiga and K. P. Lillerud, Tuned to perfection: Ironing out the defects in metal–organic framework UiO-66, Chem. Mater., 2014, 26, 4068–4071.
35. K. Leus, I. Muylaert, V. Vandichel, G. B. Marin, M. Waroquier, V. Van Speybroeck and P. Van Der Voort, The remarkable catalytic activity of the saturated metal organic framework V-MIL-47 in the cyclohexene oxidation, Chem. Commun., 2010, 46, 5085–5087.
36. V. Bernales, M. A. Ortuño, D. G. Truhlar, L. Gagliardi and C. J. Cramer, Computational design of functionalized metal–organic framework nodes for catalysis, ACS Cent. Sci., 2018, 4, 5–19.
37. K. Tan, N. Nijem, P. Canepa, Q. Gong, J. Li, T. Thonhauser and Y. J. Chabal, Stability and hydrolysis of metal organic frameworks with paddle-wheel SBUs upon hydration, Chem. Mater., 2012, 24, 3153–3167.
38. W. Morris, B. Volosskiy, S. Demir, F. Gándara, P. L. McGrier, H. Furukawa, D. Cascio, J. F. Stoddart and O. M. Yaghi, Synthesis, structure, and metalation of two new highly porous zirconium metal–organic frameworks, Inorg. Chem., 2012, 51, 6443–6445.
39. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb and J. R. Cheeseman, et al., Gaussian 16, Revision C.01, Gaussian, Inc., Wallingford, CT, 2016.
40. Y. Zhao and D. G. Truhlar, The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals, Theor. Chem. Acc., 2008, 120, 215–241.
41. Y. Zhao and D. G. Truhlar, Density functionals with broad applicability in chemistry, Acc. Chem. Res., 2008, 41, 157–167.
42. F. Weigend and R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
43. D. Andrae, U. Häußermann, M. Dolg, H. Stoll and H. Preuß, Energy-adjusted ab initio pseudopotentials for the second and third row transition elements, Theor. Chim. Acta, 1990, 77, 123–141.
44. A. V. Marenich, C. J. Cramer and D. G. Truhlar, Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions, J. Phys. Chem. B, 2009, 113, 6378–6396.
45. A. V. Marenich, C. J. Cramer and D. G. Truhlar, Performance of SM6, SM8, and SMD on the SAMPL1 test set for the prediction of small-molecule solvation free energies, J. Phys. Chem. B, 2009, 113, 4538–4543.
46. M. H. Abraham, Scales of solute hydrogen-bonding: Their construction and application to physicochemical and biochemical processes, Chem. Soc. Rev., 1993, 22, 73–83.
47. S. F. Boys and F. Bernardi, The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors, Mol. Phys., 1970, 19, 553–566.
48. E. D. Glendening, C. R. Landis and F. Weinhold, NBO 7.0: new vistas in localized and delocalized chemical bonding theory, J. Comput. Chem., 2019, 40, 2234–2241.
49. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, Oxford, 1990.
50. T. Lu and F. Chen, Multiwfn: a multifunctional wavefunction analyzer, J. Comput. Chem., 2012, 33, 580–592.
51. G. Martínez-Guajardo, K. J. Donald and G. Merino, Unraveling the mysteries of cadmium coordination chemistry, Inorg. Chem., 2013, 52, 2793–2799.
52. H. D. Burrows and M. L. T. S. Duarte, Thermodynamics of cadmium(II) complexes with acetate and similar ligands, J. Chem. Soc., Faraday Trans., 1982, 178, 3595–3606.
53. D. A. Dzombak and F. M. M. Morel, Surface Complexation Modeling: Hydrous Ferric Oxide, Wiley-Interscience, New York, 1990.
54. F. A. Cotton, G. Wilkinson, C. A. Murillo, M. Bochmann, Advanced Inorganic Chemistry, 6th edn, Wiley, New York, 1999.
55. A. E. Martell and R. M. Smith, Critical Stability Constants, Inorganic Complexes, Plenum Press, New York, 1976, Vol. 4.
56. E. Espinosa, E. Molins and C. Lecomte, Hydrogen bond strengths revealed by topological analyses of experimentally observed electron densities, Chem. Phys. Lett., 1998, 285, 170–173.
57. H. Saleem, U. Rafique and R. P. Davies, Investigations on post-synthetically modified UiO-66-NH2 for the adsorptive removal of heavy metal ions from aqueous solution, Microporous Mesoporous Mater., 2016, 221, 238–244.
58. Y. Liu, X. Guo and C. Zhu, Effect of solvent on adsorption behavior of cadmium(II) on MOF-5, J. Environ. Chem. Eng., 2019, 7, 103109.
59. M. Shokouhimehr, E. S. Soehnlen, A. Khitrin, S. Basu and S. D. Huang, Biocompatible Prussian blue nanoparticles: preparation, stability, cytotoxicity, and potential use as an MRI contrast agent, Inorg. Chem. Commun., 2010, 13, 58–61.
60. J. R. Rustad and W. H. Casey, Metastable structures and isotope exchange reactions in polyoxometalate ions provide a molecular view of oxide dissolution, Nat. Mater., 2012, 11, 223–226.
61. Z. Lin, Y. Huang and X. Wang, Adsorption of Cd(II) ions from aqueous solution by metal–organic frameworks in water–ethanol binary solutions, Environ. Sci. Pollut. Res., 2019, 26, 15030–15041.
62. WHO, Guidelines for Drinking-Water Quality, 4th edn, World Health Organization, Geneva, 2017.
63. M. Jaishankar, T. Tseten, N. Anbalagan, B. B. Mathew and K. N. Beeregowda, Toxicity, mechanism and health effects of some heavy metals, Interdiscip. Toxicol., 2014, 7, 60–72.
64. G. Genchi, M. S. Sinicropi, G. Lauria, A. Carocci and A. Catalano, The effects of cadmium toxicity, Int. J. Environ. Res. Public Health, 2020, 17, 3782.

---