Free energy of ligand removal in the metal-organic framework UiO-66

We report an investigation of the"missing-linker phenomenon"in the Zr-based metal-organic framework UiO-66 using atomistic forcefield and quantum chemical methods. For a vacant benzene dicarboxylate ligand, the lowest energy charge capping mechanism involves acetic acid or Cl-/H2O. The calculated defect free energy of formation is remarkably low, consistent with the high defect concentrations reported experimentally. A dynamic structural instability is identified for certain higher defect concentrations. In addition to the changes in material properties upon defect formation, we assess the formation of molecular aggregates, which provide an additional driving force for ligand loss. These results are expected to be of relevance to a wide range of metal-organic frameworks.


Substitution by monodentate ligands
When substituting a ligand for 2 monodentate molecules (one charge capping and one neutral) there is a choice on whether to do this so that the different groups are trans or cis. We identify trans to be favourable as a higher symmetry is maintained (Table S1).

DMF dielectric constant
The temperature dependence of the dielectric constant as reported by Bass et al. was used to obtain the free energies of solvation of molecules using the COSMO program in NWChem. S1 Table S2 gives the values extracted from reported values as published in the referenced paper. Figure S12 gives the published data points.  Figure S1: Dielectric constant of DMF with temperature as published by Bass et al.

Forcefield description of bulk UiO-66
Approaches to forcefield parametrisation can differ. S2-S4 In this instance forcefield parameters are taken from CHARMM S5,S6 with the total energy (E CHARM M ) calculated under the harmonic approximation considering both non-bonding and bonding parameters. We approached the interaction between metal node and organic linker by considering these components are separate entities bonded by van der Waals and charge interactions. Organic linkers including those used for charge capping were fully described by a bonding forcefield with CHARMM parameters with Lennard-Jones non-bonding interactions. The explicit parameters will be provided in the SI. Water molecules were described with the TIP3P (original S7 and modified S8 ) model.
Formal charges were used for the metal ions, the head of the carboxylic acid used charges for similar carboxylic acid systems. S9 The remaining charges were calculated in GULP S10 using the Gasteiger charge equilibrium model considering bond increments and nearest neighbour electronegativities. S11 Buckingham non-bonding terms were then used to describe the interaction between metal node and organic linker to reproduce bond lengths from density functional theory (DFT): All Buckingham terms were fitted to recreate both structural and materials properties as previously described.

Defect formation energies
The thermodynamic cycle considered is given in Figure S2. The description of each Gibbs free energy value used to calculate the final defect energy of formation are provided in Table   S4.
gas phase double protonation of BDC NIST database ∆G 5 solvation free-energy of BDC in DMF NWChem S5 Figure S2: Thermodynamic cycle for the removal of 1 BDC linker from UiO-66, where X is the anion used for charge compensation.      Figure S3: Configurations for 2 vacant BDC linkers. Red lines represent the missing linkers with the names of the configurations corresponding to those plotted in Figure 2 for the acetate/Cl − /H 2 O capping mechanisms, respectively. Configurations are numbered in order of defect energy magnitude and therefore differ for either capping mechanism. In brackets is the degeneracy of configuration used to calculate the Boltzmann distributions. There are 4 metal nodes in the cubic unit-cell, individual metal nodes are highlighted as green, blue and pink with all 8 corners of the plotted cubic cell making the fourth metal node. Figure S4: Configurations for 3 vacant BDC linkers. Red lines represent the missing linkers with the names of the configurations corresponding to those plotted in Figure 2 for the acetate/Cl − /H 2 O capping mechanisms respectively. Configurations are numbered in order of defect energy magnitude and therefore differ for either capping mechanism. In brackets is the degeneracy of configuration used to calculate the Boltzmann distributions. There are 4 metal nodes in the cubic unit-cell, individual metal nodes are highlighted as green, blue and pink with all 8 corners of the plotted cubic cell making the fourth metal node. S10 Table S11: Helmholtz free energies (kJmol −1 ) and overall defect formation energy (E f (defect)) for removal of 1 linker in UiO-66 for the Cl −1 charge capping models. ∆A was calculated to include contributions from enthalpy and entropy for the specified temperatures (∆A = ∆U -T∆S). Values given in brackets follow full volume relaxation at the stated temperature and therefore are the Gibbs free-energy (∆G = ∆H -T∆S).  Table S12: Helmholtz free energies (A) and overall defect formation energy (E f (defect)) for the removal of 1 linker from UiO-66. Given in this table are results considering OH − and CH 3 COO − charge capping models. ∆A was calculated to include contributions from enthalpy and entropy for the specified temperatures (∆A = ∆U -T∆S). Values given in brackets follow full volume relaxation at the stated temperature and therefore are the Gibbs free-energy (∆G = ∆H -T∆S).  4 193.4 193.3 193.6 194.1 410 178.7 193.0 193.2 192.6 196.7 196.5 196.9 197.4 S13  5 191.4 192.6 192.5 192.8 193.11 191.8 196.9 410 192.5 194.6 195.8 195.6 196.0 196.28 194.9 200.0 S15

XRD for Cl − /H 2 O charge compensation
Simulated XRD patterns (λ = 1.5418 nm) for the cubic and monoclinic phases of UiO-66 calculated with 1 and 7 linkers missing are given ( Figures S5 and S6).

Boltzmann distributions
The equilibrium distributions of defects in the structure can be estimated from where g n is the degeneracy of each configuration, ∆G n is the defect free energy (Tables S12, S13, and S14), and Z is the vibration partition function   Figure S12: Atom types used for UiO-66 framework and charge compensating molecules.  Table S21: Non-bonding parameters for the UiO-66 framework. Note the absence of Lennard-Jones parameters for Zr and O(oxide) to avoid double counting of interactions associated with the metal node. Lennard-Jones total energy is calculated using epsilon/sigma rather than A and B typical constants. E L−J = (σ/r) 12 -2(σ/r) 6 . Cut-off was set to 12.5Å and 12.0Å for Lennard-Jones and Buckingham parameters, respectively.   Table S24: Non-bonding parameters for DMF. Lennard-Jones total energy is calculated using epsilon/sigma rather than A and B typical constants. E L−J = (σ/r) 12 -2(σ/r) 6 . Cut-off was set to 12.5Å and 12.0Å for Lennard-Jones and Buckingham parameters, respectively. Table S27: Non-bonding parameters for Cl − . Lennard-Jones total energy is calculated using epsilon/sigma rather than A and B typical constants. E L−J = (σ/r) 12 -2(σ/r) 6 . Cut-off was set to 12.5Å and 12.0Å for Lennard-Jones and Buckingham parameters, respectively.  Table S28: Charges used for acetic acid. Atom types are given in Figure S9. Note that C(methyl) and H(methyl) are the same atom types as were used for DMF.

Lennard-Jones
Atom type charge O(aa) -1.13 H(methyl) 0.01 C(methyl) 0.06 Ccarb 1.17 Table S29: Non-bonding parameters for acetic acid. Lennard-Jones total energy is calculated using epsilon/sigma rather than A and B typical constants. E L−J = (σ/r) 12 -2(σ/r) 6 . Cut-off was set to 12.5Å and 12.0Å for Lennard-Jones and Buckingham parameters, respectively.   Table S32: Non-bonding TIP3P and new parameters for water. Lennard-Jones total energy is calculated using epsilon/sigma rather than A and B typical constants. E L−J = (σ/r) 12 -2(σ/r) 6 . Cut-off was set to 12.5Å and 12.0Å for Lennard-Jones and Buckingham parameters, respectively.  Figure S13: Measured intermolecular distances between components of clusters between molecular precursors.