Capture of heavy hydrogen isotopes in a metal-organic framework with active Cu(I) sites

The production of pure deuterium and the removal of tritium from nuclear waste are the key challenges in separation of light isotopes. Presently, the technological methods are extremely energy- and cost-intensive. Here we report the capture of heavy hydrogen isotopes from hydrogen gas by selective adsorption at Cu(I) sites in a metal-organic framework. At the strongly binding Cu(I) sites (32 kJ mol−1) nuclear quantum effects result in higher adsorption enthalpies of heavier isotopes. The capture mechanism takes place most efficiently at temperatures above 80 K, when an isotope exchange allows the preferential adsorption of heavy isotopologues from the gas phase. Large difference in adsorption enthalpy of 2.5 kJ mol−1 between D2 and H2 results in D2-over-H2 selectivity of 11 at 100 K, to the best of our knowledge the largest value known to date. Combination of thermal desorption spectroscopy, Raman measurements, inelastic neutron scattering and first principles calculations for H2/D2 mixtures allows the prediction of selectivities for tritium-containing isotopologues.

The data has then be transferred into Van't Hoff plots for different coverages (Supplementary Figure 3). The isosteric heat of adsorption is defined as the negative of the isosteric enthalpy of adsorption by the relation ( . 5 For a given surface coverage, it can be derived from the so-called van´t Hoff form (Eq. 1), which is analogous to the Clausius-Clapeyron equation. A detailed derivation is given elsewhere. 6 Using a linear interpolation in the adsorption isotherms for different temperatures, we calculate the pressure of equilibrium (P eq ) for a fixed amount of adsorbed gas (n a ). The variation of the pressure of equilibrium with the temperature for a fixed amount of gas is the so-called Isoster curve and a linear relation is found between the ln(P eq ) and the reciprocal temperature 1/T (see Supplementary Figure 3). A linear regression model was employed for the isosters at different gas concentrations. Using the Equation 1 with the obtained slopes we calculate the isosteric heat of adsorption in function of the coverage (see Supplementary Figure 4).
An alternative analytical approach is based on the Langmuir model. In the low-pressure range we can restrict our considerations to the strong adsorption sites. The data can be fitted to the linear form of the Langmuir equation (Eq. 2). (2) Supplementary Figure 2 shows the H 2 and D 2 adsorption isotherms for different temperatures and their corresponding Langmuir fitting using the range (0-0.1bar). Supplementary Table 1 shows the equilibrium constants and saturation volumes for the H 2 and D 2 adsorption isotherms using the Langmuir fitting. The adsorption enthalpies and entropies are obtained by fitting temperature-dependent equilibrium constants (k) to the van't Hoff equation  Supplementary Note 2: Thermal desorption spectroscopy (TDS). TDS is commonly used in surface science to measure the activation energy of desorption. The activation energies can be determined by two methods (see below), utilizing measurements with different heating rates. Both methods are equivalent, do not require any assumption about pre-exponential factor, reaction order or specific mechanism, and the values for the desorption energies are typically in agreement with values for the isosteric heat of adsorption. 7 After exposure to an equimolar mixture of 10mbar H 2 /D 2 at 100K the thermal desorption spectra have been recorded applying different heating rates of 0.01K.s -1 , 0.05K.s -1 and 0.1K.s -1 . The desorption maxima for both H 2 and D 2 show a clear shift to lower temperatures for slower heating rates indicating a thermally activated desorption process (Supplementary Figure 6 left). Owing to the exchange of H 2 by D 2 at an exposure temperature of 100K, the D 2 desorption rate is by an order of magnitude higher than for H 2 .
For an adsorbed gas molecule, the desorption process can be described by the Polanyi-Wigner equation for first order with the assumption that the pre-exponential factor and the desorption energy are independent of coverage: with instantaneous coverage θ, frequency factor , desorption energy E des , gas constant R, and temperature T. In the case of a linear heating rate β, the temperature can be described by (4) where is the initial temperature. The Polanyi-Wigner equation can be rewritten to A plot of ( ) versus ( ) yields a straight line and the desorption energy can be extracted from the slope (see Supplementary Figure 6 right). 8 Alternatively, Falconer proposed that for the maximum desorption rate at the peak temperature the following relation holds according to equation (3): 7 Again, the logarithm of the maximum desorption rate versus ( ) yields a straight line, and the desorption energy can be extracted from the slope (see Supplementary Figure 6 right).
Supplementary Table 3 gives the desorption energies determined by the two independent methods together with the error of the linear regression. For TDS the error including all experimental uncertainties is difficult to quantify, however, can be approximated to about ±2 kJ.mol -1 . Owing to the smaller signal of H 2 the uncertainty of the evaluation maybe even higher for H 2 .
The heats of adsorption and the activation energies of desorption obtained by three different methods are compared in Supplementary Table 5. Here possible experimental uncertainties, which are higher than the standard deviation, are included. All methods are in good agreement within the expected error range. Furthermore, the value for D 2 is about 2.5 kJ·mol -1 higher than for H 2 which can be related to the different zero-point energies of H 2 and D 2 .
Supplementary Figure 6: TDS spectra of different heating rates. Comparison of different heating rates. Left: The equimolar mixture of 10mbar has been applied for 10min at 100K. The heating procedure has been performed with different heating rates (blue: 0.01K.s -1 , red: 0.05K.s -1 , black: 0.1K.s -1 ). Right: Procedure used for the determination of the desorption energy.  Figure 7: Higher-energy signals in the INS spectra recorded at different temperatures. Development of D 2 replacement at the strong adsorption site, during heating and cooling. The peaks at 4.9 and 6.5 meV, which belong to the strong adsorption sites, show no change between the first 5 K (black) and higher temperatures 40 K, 77 K, indicating that the same H 2 amount is still bound to this site. At the weaker adsorption site (13.4 and 14.7 meV) there is already a decrease in intensity for the 40 K measurement. At 77 K, in addition to the peaks of the strong adsorption site, a third peak emerges around 7.9 meV. This 7.9meV peak was hidden in the broader 6.5 meV signal (Figure 3a, high loading) and appears clearly after the exchange of H 2 by D 2 at the strong adsorption site. A comparison of the spectra measured at 5 K before (black) and after (red) heating to 200 K shows that the intensity does not change significantly for the 7.9 meV transition. The 7.9 meV signal can thus be assigned to an orthohydrogen vibration, which is bound to the weaker adsorption site. The peak at 23.0 meV is most likely the combination of a phonon mode and the rotational transition at 14.7 meV. Any transition from the ground state to an excited state of parahydrogen has very low cross section for neutrons, 9 however, the transitions between the ground state to orthohydrogen states have the larger cross section (82 barn) and the combination of a rotational transition and a phonon also has a high cross section. The signal at 28.7 meV belongs to the transition (J=1 to J=2) in weakly adsorbed H 2 due to the presence of a small amount of orthohydrogen. The excess signal that is observed at the higher loading corresponds to the recoil of the molecular hydrogen adsorbed on the second site. 1 The lack of recoil observed in the initial low loading reflects that all hydrogen molecules are bound at the very strong first adsorption site.

Supplementary
Supplementary Note 3: Raman measurements. The 514.5 nm line of an Argon ion laser was used as excitation source with an incident power of less than 5 mW to avoid sample degradation and heating effects. The scattered light was analyzed by a T64000 Jobin-Yvon Raman spectrometer equipped with a liquid-nitrogen-cooled charge coupled detector (CCD). The measurements were performed in the single grating configuration yielding a resolution of 1.9 cm -1 in the spectral region of interest. A calibration spectrum of an Ar lamp was recorded after each measurement.
For the application of high gas pressures and low temperatures we used a sample cell specially designed for optical measurements, which is described in reference. 10 After activating at 453 K for 3 h under high vacuum, the sample was loaded in the pressure cell under argon atmosphere in a glove box. Then the cell was connected to a turbo molecular pump, providing high vacuum prior to the gas loading. Only half of the cell was filled with Cu(I)-MFU-4l enabling to record reference spectra of the free gas occupying the empty part of the cell under the same pressure/temperature conditions as the sample. The gas pressure inside the cell was stepwise increased up to 22 bar. For the low temperature measurements (down to 40 K), the cell was immersed to a helium bath cryostat with a silicon-diode temperature sensor.  2)) a redshifted signal is observed, which is almost independent of pressure at 45K, whereas at 100K the shifted signal is much weaker and only visible at higher pressures. The attractive interaction between the adsorption sites of Cu(I)-MFU-4l and the H 2 (D 2 ) molecules reduces the degrees of freedom of these physisorbed molecules, which leads to a modified vibrational transitions of the Q(J)-branch. 11

Supplementary Note 4: Benchmark of DFT calculation protocol.
A small molecule (CuH), which represents the Cu(I) site, has been used as a benchmark for the various methods (PBE0-D3, MP2 and CCSD(T)) and basis sets (def2-TZVP, cc-pVTZ, cc-pVQZ) used. We want also to examine the effect of an all-electron (def2-TZVP and cc-pVQZ) versus a basis set with ECP (cc-pVQZ-PP). Moreover, the effect of anharmonicities on frequencies is explored. H 2 is interacting with CuH with a T-shaped orientation, exactly in the manner that is adsorbed on the copper site of the Cu(I)-MFU-4l. Results are summarized in Supplementary  Table 6. According to Supplementary Table 6 Supplementary Table 7. The binding energies and enthalpies change less by 1kJ mol -1 between the models I-Form and I-Cl and the frequencies by ~10 cm -1 . This validates our approach to use Cl anions instead of formate as capping ligands to the Zn sites. The second benchmark is on the effect of the model size. We compare results for models I-Cl and II-Cl from the PBE-D3/def2-TZVP and PBE0-D3/def2-TZVP methods in Supplementary  Table 8. The results change by -2 and -1 kJ mol -1 respectively for H 2 and N 2 , when the model size is increased from I to II.

Supplementary Note 5: Error estimate for using molecular instead of periodic model.
We performed further calculations to assess the error in the computed binding energies of H 2 , from using a finite cluster model instead of the periodic cell. Due to the big cell size, geometry optimization is computationally time and resource demanding, thus an approximate scheme has been used by taking the crystal structure of the original MFU-4l and replacing one of the units in the corner with the atomic positions of the Model-II-Cl (shown in Supplementary Figure 9). We used the PBE0-D3 optimized coordinates of the Model-II-Cl with H 2 . Then three single point energy calculations are performed: dimer of H 2 @Cu-MFU-4l, Cu-MFU-4l and H 2 at the positions of the dimer. In this case, the interaction energies are calculated rather than the binding energies, i.e. the deformation energies of the individual monomers are not included. That is the reason, why the interaction energies are much higher compared to binding energies. These periodic calculations have been performed with the VASP 5.2 program and using the PBE functional with the -D3 dispersion corrections and 400eV kinetic energy cutoff. We have also used the dispersion parameters derived for the PBE0 functional. The same procedure was employed for the molecular Model-II-Cl inserted in a cubic cell of 25.0Å. In this way, we can estimate the long-range effects of the H 2 binding in the crystal. We further decompose the interaction energies into contributions from PBE and -D3 dispersions only. The results are shown in Supplementary  The energy decomposition analysis (EDA), 12 as implemented in the ADF 13 program, was used to analyze the interactions between H 2 and the Cu(I) site. In this scheme, the binding energy can be written as the sum of the interaction energies (ΔE int ) and the deformation energies (ΔE def ) of the monomers. ΔE int is further decomposed into terms that contain a) the electrostatic interactions (ΔV elst ) between the unperturbed charge distributions of the deformed fragment with the field of the other, b) the Pauli repulsion (ΔE Pauli ), which is associated with going from the unperturbed individual fragments to a symmetrized and orthonormalized wave function of their product Ψ 0 =NA[Ψ A Ψ B ] that obeys the Pauli exclusion principle (N is a normalization constant and A the anti-symmetry operator), and c) orbital interactions (ΔE OI ) that account for electron pair bonding, charge transfer, and polarization effects when going from Ψ 0 to the converged wave function of the complex. The first two terms can be combined and the sum of the electrostatic interactions and Pauli repulsion is called steric interaction. The orbital interactions can be further decomposed by employing the natural orbitals for chemical valence (NOCV) theory in combination with the extended transition state (ETS) method. 14,15 The ETS-NOCV decomposes the orbital interactions into different components (σ, π, δ) of the chemical bond. These calculations have been performed with the ADF program using similar computational details as before, that is, the PBE0-D3 functional and the TZP basis set. The results are presented in Supplementary Table10. According to the EDA, the electrostatic interactions are 62% of the total attractive interactions, which is the sum of the electrostatic, the orbital interactions and the dispersions.
Dispersions are very small, only 2%, and orbital interactions reach 36%. The orbital interactions are further decomposed with the ETS-NOCV scheme. After analyzing the NOCVs, we identified the pair of NOCVs that contribute to forward-and back-donation. Contour plots of the most important NOCVs are shown in Supplementary Figure 11. According to the ETS-NOCV, forward-donation terms are 63.5 kJ mol -1 or 62% of the total orbital interactions. Back-donation terms contribute a significant amount of 34.9 kJ mol -1 or 34%. Supplementary ; zero-point energy of adsorbed H 2 and D 2 ( ); Entropy of adsorbed H 2 and D 2 at 77K ( ); Enthalpy of adsorbed H 2 and D 2 at 77K ( ); Gibbs energy of adsorbed H 2 and D 2 at 77K ( ); frequencies of adsorbed H 2 (ν). The calculations have been done with PBE0-D3/def2-TZVPP (H 2 )/aug-cc-pVTZ-PP (Cu)/def2-TZVP (other atoms) on the model-I-Cl of Cu(I)-MFU-4l (Supplementary Figure 9b) (free D 2 )=18.6 kJ mol -1 ; ν(free H 2 )=4408 cm -1 . ν R,T corresponds to modes associated to rotations and translations of the free molecules, respectively.