Surface Modification Strategy for Enhanced NO2 Capture in Metal–Organic Frameworks

The interaction strength of nitrogen dioxide (NO2) with a set of 43 functionalized benzene molecules was investigated by performing density functional theory (DFT) calculations. The functional groups under study were strategically selected as potential modifications of the organic linker of existing metal–organic frameworks (MOFs) in order to enhance their uptake of NO2 molecules. Among the functional groups considered, the highest interaction energy with NO2 (5.4 kcal/mol) was found for phenyl hydrogen sulfate (-OSO3H) at the RI-DSD-BLYP/def2-TZVPP level of theory—an interaction almost three times larger than the corresponding binding energy for non-functionalized benzene (2.0 kcal/mol). The groups with the strongest NO2 interactions (-OSO3H, -PO3H2, -OPO3H2) were selected for functionalizing the linker of IRMOF-8 and investigating the trend in their NO2 uptake capacities with grand canonical Monte Carlo (GCMC) simulations at ambient temperature for a wide pressure range. The predicted isotherms show a profound enhancement of the NO2 uptake with the introduction of the strongly-binding functional groups in the framework, rendering them promising modification candidates for improving the NO2 uptake performance not only in MOFs but also in various other porous materials.


Introduction
Nitrogen dioxide (NO 2 ) belongs to a group of highly reactive gases known as nitrogen oxides (NOx) and primarily gets in the air from the burning of fuels [1]. The NO 2 generated by the exhaust gases of the industries as well as by our private cars become an important air pollutant, the toxicity of which has an impact on the environment and human health [2]. NO 2 is formed in the combustion processes of heating systems. The main pollutant emitted directly from hydrocarbon combustion is nitric oxide (NO), along with a small proportion of nitrogen dioxide (NO 2 ). Nitrogen oxide is oxidized by ozone (O 3 ) in the atmosphere, on a 10 min time scale, to give NO 2 [3]. According to the EPA (Environmental Protection Agency), when NO 2 is present in the air, it can be harmful to human health due to the irritation created in the airways of the human respiratory system [4]. In addition, NO 2 and other NO x molecules can interact with water, oxygen, and other chemicals in the atmosphere to form acid rain [5,6], whose harmful ecological effects are detrimental and most prominent in aquatic environments. Due to the serious impacts of nitrogen pollution [2], there are several separation methods of NO 2 from industrial gas mixtures [7] and various adsorbents which have been used for the removal of NO 2 ; some of them including zeolites [8], wood-based activated carbon [9], graphite oxides and iron composites [10]. Zhu et al., in an attempt to adsorb NO 2 in N-doping activated carbon, observed that N-doping leads to easier adsorption of NO 2 molecules, thus increasing NO 2 physisorption energies [11]. After theoretically investigating NO 2 adsorption on the surface of a silicon carbide (SiC) nanotube, Iranimanesh et al. proposed it for NO 2 gas pollutant sensing and removal [12]. Lately, the separation of NO 2 from industrial gas mixtures through its adsorption has been clearly identified as a possible application of metal-organic frameworks (MOFs) [13,14].
In the last two decades, there has been a rapid growth of metal-organic frameworks (MOFs) in the field of porous materials, and their applications vary from adsorption and separation of gases [15] to catalysis [16] and drug delivery [17]. MOFs are crystalline materials with extremely high porosity (up to 90% free volume) and a huge internal surface area extending beyond 6000 m 2 /g [18]. Due to these properties and considering the plethora of both inorganic and organic building blocks, MOFs have attracted great research interest as high-capacity adsorbents to meet various separation demands [18][19][20]. One of the important structure-to-property flexibilities of the MOF structures is the potential of their organic linkers to be modified with the incorporation of various functional groups in order to tune their interaction with selected molecules [18].
Several previous studies have shown that the introduction of functional groups into MOFs leads to enhanced gas uptake performance. In the work of Frysali et al. [21], the introduction of a sulfate anion in the phenyl ring has the highest interaction energy (−5.4 kcal/mol) with CO 2 , a value almost two times larger than the corresponding binding energy for benzene (−2.9 kcal/mol). Klontzas et al. [22] showed that the gravimetric capacity of the Li modified IRMOF-8 was calculated to 10 wt% at 77 K and 100 bar, while the corresponding values show great promise also at room temperature with an uptake of 4.5 wt%, performances significantly enhanced with respect to the unmodified counterpart (up to three times stronger). For NO 2 , the theoretical studies made by Fioretos et al. [23] have shown that NO 2 interacts stronger with functionalized benzenes such as aniline, phenol, and toluene (with binding energies of −2.26 kcal/mol, −1.72 kcal/mol and −2.02 kcal/mol respectively) than with benzene (−1.67 kcal/mol). The amino (-NH 2 ) substituent can be particularly beneficial as the interaction of strongly polar molecules, such as NO 2 , with the amino-substituted aromatic rings is characterized by the contribution of electrostatic dipole-dipole forces resulting in enhanced adsorption.
Taking into account that many MOF frameworks have a phenyl group in their organic linker, together with the effectiveness of the organic linker functionalization strategy for tuning their interaction with guest molecules, in this work, we investigate the interaction of NO 2 molecules with a series of 43 carefully selected benzene molecules by means of density functional theory (DFT) calculations. The selection of the functional groups was based on chemical intuition and findings of previous similar studies [22,[24][25][26][27]. Subsequently, in order to verify the effectiveness of the functionalization on the enhancement of NO 2 uptake in MOFs, grand canonical Monte Carlo (GCMC) simulations were performed for the IRMOF-8, modified by the three functional groups that showed the strongest interaction with NO 2 .

Density Functional Theory
To investigate the interaction of NO 2 molecules with the organic linkers of MOFs, we start with the simplified model of a benzene ring. A large set of 43 functional groups was examined for their binding strength towards the NO 2 molecule. The conformations of the functionalized benzene molecules were optimized using Gershom Martin's double-hybrid density functional DSD BLYP in the resolution of identity (RI) approximation [28] along with the def2-TZVPP basis set and with the corresponding auxiliary basis set for the RI approximation [29][30][31], including the D3BJ (Becke-Johnson damping version) empirical correction term for the dispersion interactions as proposed by Grimme [32][33][34][35]. All geometries were optimized without any symmetry constraints, and the optimized minimum-energy structures were verified as stationary points on the potential energy surface by performing  [36,37]. The dimer energies were corrected for the basis set superposition error (BSSE) using the counterpoise (CP) method as proposed by Boys and Bernardi [38].
The electron density redistribution was calculated as the difference between the electron density of the functionalized benzene-NO 2 (∆D) complex and the electron densities of the functionalized benzene (D(functionalized benzene)) and NO 2 (D(NO 2 )) molecule according to the formula: All electron densities were calculated at the RI-DSD-BLYP/def2-TZVPP level of theory. Mathematical operations on the electron densities along with the visualization of the electron density difference, were done using gOpenMol graphics program [39,40].

Grand Canonical Monte Carlo
To verify the effectiveness of the strongest interacting functional group candidates obtained from the DFT calculations to enhance the NO 2 uptake in MOFs, we employed Monte Carlo simulations in the grand canonical ensemble. IRMOF-8 was selected as the platform for the uptake calculation and was functionalized with the best performing functional groups, as shown in Figure 1.
the RI approximation [29][30][31], including the D3BJ (Becke-Johnson damping version) empirical correction term for the dispersion interactions as proposed by Grimme [32][33][34][35]. All geometries were optimized without any symmetry constraints, and the optimized minimum-energy structures were verified as stationary points on the potential energy surface by performing numerical harmonic vibrational frequency calculations. All calculations were performed with the Orca 4.2 program package [36,37]. The dimer energies were corrected for the basis set superposition error (BSSE) using the counterpoise (CP) method as proposed by Boys and Bernardi [38].
The electron density redistribution was calculated as the difference between the electron density of the functionalized benzene-NO2 (ΔD) complex and the electron densities of the functionalized benzene (D(functionalized benzene)) and NO2 (D(NO2)) molecule according to the formula: All electron densities were calculated at the RI-DSD-BLYP/def2-TZVPP level of theory. Mathematical operations on the electron densities along with the visualization of the electron density difference, were done using gOpenMol graphics program [39,40].

Grand Canonical Monte Carlo
To verify the effectiveness of the strongest interacting functional group candidates obtained from the DFT calculations to enhance the NO2 uptake in MOFs, we employed Monte Carlo simulations in the grand canonical ensemble. IRMOF-8 was selected as the platform for the uptake calculation and was functionalized with the best performing functional groups, as shown in Figure 1. The ΝΟ2 adsorption was calculated for a pressure range up to 1.2 bar at 298 K. Fugacity coefficients for the different thermodynamic states were defined using the Penge-Robinson equation of state [41]. The framework coordinates were taken from the crystallographic data [42] and a cubic periodic box of size 30.1 × 30.1 × 30.1 Å 3 was used for all the frameworks, functionalized or parent. Simulations were performed in supercells incorporating enough repeat units such that all edge lengths were greater than 25.6 Å, i.e., twice the Lennard-Jones (LJ) cut-off radius. For each simulation point, 50,000 cycles were The NO 2 adsorption was calculated for a pressure range up to 1.2 bar at 298 K. Fugacity coefficients for the different thermodynamic states were defined using the Penge-Robinson equation of state [41]. The framework coordinates were taken from the crystallographic data [42] and a cubic periodic box of size 30.1 × 30.1 × 30.1 Å 3 was used for all the frameworks, functionalized or parent. Simulations were performed in supercells incorporating enough repeat units such that all edge lengths were greater than 25.6 Å, i.e., twice the Lennard-Jones (LJ) cut-off radius. For each simulation point, 50,000 cycles were performed for system equilibration, followed by additional 100,000 cycles for sampling over the ensemble averages. For the description of the interactions between the IRMOF-8 and the NO 2 atoms, LJ + Coulomb potentials were used and each atom of the host or the guest was treated explicitly [43]. The framework of IRMOF-8 was kept rigid during the simulations, while NO 2 molecules were allowed to translate and rotate. Nitrogen dioxides were treated as three rigid center molecules with bond lengths between nitrogen and oxygen atoms held fixed at 1.19 Å. For the electrostatic interactions between NO 2 molecules and the host material, point charges equal to qO = −0.073 and qN = 0.146 were placed at oxygen and nitrogen sites of NO 2 , respectively. The framework atoms charges were defined by employing the CHELPG method [44]. For the van der Waals interactions, potential parameters according to COMPASS Force Field [45,46] model were used, with ε = 50.36 and σ = 3.24 Å for Nitrogen atom and ε = 62.51 and σ = 2.93 Å for the oxygen center.
For each MOF framework, the potential parameters were taken from the UFF force field [47], except for organic linker atoms that were treated separately. Lorenz-Berthelot mixing rules were used to describe the NO 2 −IRMOF-8 interactions. For the functionalized organic linkers, the UFF parameters were found to be inconsistent with DFT interactions. To correct this, the parameters of the classical potential were fitted to reproduce the quantum chemical data (SI). All GCMC calculations were carried out with the RASPA software package [48].

Results and Discussion
The interaction of the NO 2 molecule with the full set of the 43 strategically functionalized benzenes was calculated at the RI-DSD-BLYP/def2-TZVPP level of theory and can be seen in Table S1 of the SI section. In Table 1, we present the nine functional group candidates with the highest binding energy together with the non-functionalized benzene for comparison. In Figure 2, the DFT optimized dimer geometries of the NO 2 . . . C 6 H 5 -X systems are shown. Table 1. Binding energies in kcal/mol of the NO 2 . . . C 6 H 5 -X systems, calculated at the RI-DSD-BLYP D3(BJ)/def2-TZVPP level of theory and percentage of binding energy enhancement with respect to the introduction of the unfunctionalized benzene. All interaction energies are corrected for the BSSE by the full counterpoise method [39].

System
Binding Energy (kcal/mol) The highest interaction energy (5.4 kcal/mol) with NO 2 was found for phenyl hydrogen sulfate (-OSO 3 H), which is almost three times stronger than the corresponding binding energy of the unfunctionalized benzene (2.0 kcal/mol).
The energetically most favorable structures are characterized by a weak interaction between NO 2 and organic linkers. More specifically, the acidic protons of the substituents, especially (-OSO 3 H, -OPO 3 H 2 , -PO 3 H 2 , -SOOH, -C(OH) 3 ), tend to interact with NO 2 s oxygen site with binding distances between 1.89 Å and 2.22 Å. The trend in the best binding energies is also confirmed by the electron density redistribution plots in Figure 3. Due to the electrostatic nature of NO 2 s interactions, these plots can serve as a rule of thumb for predicting the most stable NO 2 complexes with organic molecules. The red regions of the electron redistribution plots correspond to rich electron areas located around oxygen atoms, where the blue regions correspond to poor electron areas located mainly around hydrogen atoms. The electron-rich region of the nitrogen atom of NO 2 interacts with the electron-poor regions of the functionalized benzene, and the electron-rich regions around oxygen atoms of NO 2 interact with the electron-poor regions (around hydrogen atoms) of the functionalized benzenes.  Hydrogen, carbon, oxygen, nitrogen, lithium, sulfur, and phosphorus atoms are depicted as white, gray, red, blue, purple, yellow, and orange spheres, respectively. electron redistribution plots correspond to rich electron areas located around oxygen atoms, where the blue regions correspond to poor electron areas located mainly around hydrogen atoms. The electron-rich region of the nitrogen atom of NO2 interacts with the electron-poor regions of the functionalized benzene, and the electron-rich regions around oxygen atoms of NO2 interact with the electron-poor regions (around hydrogen atoms) of the functionalized benzenes.  To test the effect of the proposed surface modification strategy on the enhancement of the NO 2 capture in MOF structures, we selected IRMOF-08 as the platform and modified its framework by adding one functional group per linker, as shown in Figure 1. We calculated the NO 2 adsorption isotherms by performing GCMC simulations for the three functional groups (-OSO 3 H, -PO 3 H 2 , -OPO 3 H 2 ) with the stronger interactions. From the excess volumetric and gravimetric isotherms at 298 K and pressures up to 1.2 bar (Figure 4a,b), there is a significant enhancement of the uptake due to functionalization. At 1.0 bar, the corresponding volumetric and gravimetric uptake for the IRMOF-8 with the OSO 3 H group was found at 193 cm 3 (STP)/cm 3 and 15 mmol/g, respectively, a value 12.8% larger than that of the unmodified IRMOF-8. The corresponding volumetric and gravimetric uptake for the unmodified IRMOF-8 was found at 7.5 cm 3 (STP)/cm 3 and 0.8 mmol/g. To test the effect of the proposed surface modification strategy on the enhancement of the NO2 capture in MOF structures, we selected IRMOF-08 as the platform and modified its framework by adding one functional group per linker, as shown in Figure 1. We calculated the NO2 adsorption isotherms by performing GCMC simulations for the three functional groups (-OSO3H, -PO3H2, -OPO3H2) with the stronger interactions. From the excess volumetric and gravimetric isotherms at 298 K and pressures up to 1.2 bar ( Figure  4a, b), there is a significant enhancement of the uptake due to functionalization. At 1.0 bar, the corresponding volumetric and gravimetric uptake for the IRMOF-8 with the OSO3H group was found at 193 cm 3 (STP)/cm 3 and 15 mmol/g, respectively, a value 12.8% larger than that of the unmodified IRMOF-8. The corresponding volumetric and gravimetric uptake for the unmodified IRMOF-8 was found at 7.5 cm 3 (STP)/cm 3 and 0.8 mmol/g.  Figure 5 shows representative snapshots taken at 0.1 and 1.0 bar for the parent and the three functionalized frameworks. In both pressures, the modified material hosts considerably more NO2 molecules than the parent structure due to the stronger binding sites introduced to the structure by functionalization. This is also verified by the fact that the NO2 molecules are located closer to the functional groups.  Figure 5 shows representative snapshots taken at 0.1 and 1.0 bar for the parent and the three functionalized frameworks. In both pressures, the modified material hosts considerably more NO 2 molecules than the parent structure due to the stronger binding sites introduced to the structure by functionalization. This is also verified by the fact that the NO 2 molecules are located closer to the functional groups.

Conclusions
In this work, we studied by means of density functional theory the interaction of a chemically diverse set of 43 functionalized benzenes with the NO2 molecule.

Conclusions
In this work, we studied by means of density functional theory the interaction of a chemically diverse set of 43 functionalized benzenes with the NO 2 molecule.
The highest interaction energy with NO 2 (5.4 kcal/mol at the RI-DSD-BLYP/def2-TZVPP level of theory) was found for phenyl hydrogen sulfate (-OSO 3 H)-an interaction almost three times larger than the corresponding binding energy for non-functionalized benzene (2.0 kcal/mol).
The groups with the highest NO 2 binding (-OSO 3 H, -PO 3 H 2 , -OPO 3 H 2 ) were selected for functionalizing the linker of IRMOF-8 and investigating the trend in their NO 2 uptake capacities with grand canonical Monte Carlo (GCMC) simulations. GCMC simulations showed a clear enhancement of the NO 2 uptake both gravimetrically and volumetrically at 298 K and pressures up to 1.2 bar for the functionalized MOF, an enhancement even more pronounced at low pressures; at 0.1 bar, the volumetric uptake becomes 40 or 60 or 110 times larger than the unmodified IRMOF-08 by introducing -OPO 3 H 2 , -OSO 3 H, -PO 3 H 2 functional groups, respectively. Based on this significant enhancement, we propose our surface functionalization as a general strategy for improving the NO 2 adsorption uptake not only in MOFs, but also in various other porous materials. Our theoretical results can serve as high accuracy reference calculations and guide synthesis towards materials with high NO 2 adsorption capacity.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/molecules27113448/s1, Table S1: Sorted binding energies (kcal/mol) of the NO 2 . . . C 6 H 5 -X systems under study, calculated at the RI-DSD-BLYP D3(BJ)/def2-TZVPP level of theory. All interaction energy values have been corrected for the basis set superposition error (BSSE) by the full counterpoise method [38]. Percentage of binding energy enhancement with the introduction of the FG compared to benzene; Figure S1: Global minima geometries and binding energy values (in kcal/mol) of all the systems in this study; Figure S2: Electron density redistribution plots of the optimized geometries of the NO 2 . . . C 6 H 5 -X complexes. With red and blue being the regions that gain and lose electron density upon the formation of the complex, respectively; Figure S3: Fitting of the (ε, σ) parameters of the UFF [47] potential on the QM data obtained from the ab-initio scan of NO 2 over benzene; Figure S4: Fitting of the (ε, σ) parameters for the NO 2 . . . C 6 H 5 -OPO 3 H 2 interaction; Figure S5: Fitting of the (ε, σ) parameters for the NO 2 . . . C 6 H 5 -OSO 3 H interaction; Figure S6: Fitting of the (ε, σ) parameters for the NO 2 . . . C 6 H 5 -PO 3 H 2 interaction.
Funding: This research has been co-financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship, and Innovation, under the call RESEARCH-CREATE-INNOVATE (project code: T2E∆K-01976).
Institutional Review Board Statement: Not applicable.