To achieve a relatively stable adsorption site for each adsorption, we utilized several geometries sampled by Monte Carlo (MC) as our initial guess. After optimization of each geometry, the most stable geometries were further used for adsorption energy calculations. Geometry optimizations and single-point energy calculations were performed using density functional theory (DFT) at the ωB97X-D 66-68 /Def2-SVP and ωB97X-D/Def2-TZVPPD 69-71 levels of theory, respectively, unless noted otherwise. The zero-point energy and thermal corrections were evaluated by the rigid rotor harmonic oscillator (RRHO) method. All calculations were performed using a development version of the Q-Chem software package. Monte Carlo simulations. Using the open-source RASPA package, the Widom particle insertion method was employed to calculate the heats of adsorption for a diverse set of gas molecules in UiO-66 and MOF-808. MC simulations in a canonical ensemble (i.e., NVT) at a low temperature of 100 K were carried out to probe the favorable adsorption configurations as well as the associated energies. As discussed in the previous subsection, these configurations were subsequently used as the starting point in the QM/MM calculations. In the MC calculations, the structures of UIO-66 and MOF-808 were adopted from the Cambridge Crystallographic Data Centre (CCDC: 1018045 and 1002672, respectively) with some modifications (i.e., adding missing hydrogen atoms on the metal node ). The framework structure was considered rigid. Intermolecular interactions were described using the pair-wise sum of 12-6 Lennard-Jones potential and Coulombic energies with static point charges. The former was truncated and shifted to zero at a cutoff radius of 12 Å, while the long-range nature of the latter was computed using the Ewald summation method. The Lennard-Jones parameters for the framework atoms were taken from the UFF, while their partial charges were assigned with the values reported in the study of Yang et al. To describe guest molecules considered herein, various models were adopted. Sulfur dioxide (SO2) and carbon monoxide (CO) were both modeled using the ESP-MM model. Cycloheptane, methyl cyclopentane, methylcyclohexane, cyclohexene, and cyclohexanol were described using the OPLS-AA model. All of the rest of the molecules were modeled using the TraPPE model. Lennard-Jones parameters between dissimilar atoms were estimated with the Lorentz-Berthelot mixing rule. The dimensions of the simulation domain were at least twice the cutoff radius in each perpendicular direction.