All DFT and AIMD simulations used the Quickstep code that is integrated in the CP2K software. This code uses the Gaussian and plane wave method (GPW), employing an auxiliary plane wave basis set within a Gaussian orbital scheme. All calculations used "molecularly optimised" Gaussian basis sets from the work of Vandevondele and Hutter, employing basis sets of double-zeta quality (DZVP-MOLOPT-SR) for AIMD simulations and triple-zeta basis sets (TZVP-MOLOPT and TZVP-MOLOPT-SR) for DFT optimisations and calculations of vibrational spectra. To model the core electrons, Gaussian-type Goedeker-Teter-Hutter (GTH) pseudopotentials derived by Krack were used. As in the previous study of 5-FU adsorption in H-FAU, the dispersion-corrected rev-vdW-DF2 functional was used, which was found to perform well in a benchmarking of DFT adsorption energies against results from high-level coupled-cluster calculations. Structure optimisations and vibrational calculations used triple-zeta basis sets and a cutoff energy of 900 Ry. Fixing the unit cell parameters, all atomic coordinates were optimised using a BFGS optimiser, enforcing the following convergence criteria: Maximal gradient: 5•10 -6 Ha bohr -1 ; maximal displacement: 10 -5 bohr (1 Ha = 2,625.5 kJ mol -1 , 1 bohr = 0.52918 Å). All structure visualisations were prepared using VESTA. Calculations of infrared (IR) spectra of adsorbed 5-FU used the partial Hessian approach, including only displacements of 5-FU and the closest framework atoms. IR spectra were plotted using MOLDEN, employing a Lorentzian broadening with a full width at half maximum of 10 cm -1 . AIMD simulations used double-zeta basis sets and a cutoff energy of 600 Ry. They were carried out in the canonical ensemble (constant N, V, and T), using a time step of 0.5 fs and a Nosé-Hoover thermostat with a time constant of 50 fs. Different temperatures (298, 373, and 473 K) and different durations (from 12.5 ps up to 25 ps, i.e., between 25,000 and 50,000 steps) were considered in different simulations, as detailed in the results section. From the AIMD trajectories, snapshots were extracted every 2,500 steps (= every 1.25 ps). Each of these snapshots was fully optimised as described above, thereby allowing to probe several local minima. Finally, Al NMR isotropic chemical shifts were computed for selected DFT-optimised structures. These calculations used the gauge-including projector augmented-wave (GIPAW) method as implemented in the CASTEP code. The PBE exchange-correlation functional was used. Isotropic chemical shieldings were converted into chemical shifts using an equation proposed in a recent DFT study of zeolites, which employed a very similar computational approach. Further details are provided in section S1.3.