Computational Evaluation of the Impact of Incorporated Nitrogen and Oxygen Heteroatoms on the Affinity of Polyaromatic Ligands for Carbon Dioxide and Methane in Metal-Organic Frameworks

Density functional theory is employed to explore the binding of carbon dioxide and methane in a series of isoreticular metal-organic frameworks, with particular emphasis on understanding the impact of directly incorporated nitrogen and oxygen heteroatoms on the affinity of the ligand for CO 2 and CH 4 . While the strongest binding sites for both CO 2 and CH 4 were found to be directly above the aromatic rings of the core of the ligand, the introduction of heteroatoms to the core systems were shown to significantly alter both the binding strength and preferred binding locations of CH 4 and CO 2 . The presence of pyrazine rings within the ligand was observed to create new binding sites for both CO 2 and CH 4 and, in the case of CO 2 , severely reduce the binding strength or entirely eliminate binding sites which were prominent in the analogous carbocyclic ligands. These results suggest that while the presence of framework nitrogen and oxygen heteroatoms provides a route to ligands with enhanced affinity for methane, a similar increase in affinity for CO 2 is not guaranteed.

5][26][27] In the present work, the binding of CO2 and CH4 with a number of real and hypothetical organic ligands inspired by the recently synthesized experimental series of metal-organic frameworks, MFM-18X, 28 is explored in vacuo using dispersion-corrected density functional theory in order to shed light on the influence of directly incorporated heteroatomic species on ligand-guest interactions.
The MFM-18X series of frameworks (where X = 0, 1, 3 or 5) employs octacarboxylate ligands with central polycyclic cores of increasing length (Figure 1) and have been shown experimentally to exhibit promise for both CO2 and CH4 adsorption applications.The structures all contain three types of cavities, of which two are particularly relevant to the present work, being sufficiently large (diameter > 12.4 Å) for the adsorption mechanism at low loading to be dominated by interactions between the gas and the ligand cores which form the pore walls.
While the ligands present in MFM-180 and 181 are formed primarily from purely carboncontaining aromatic hydrocarbon, the ligands of MFM-183 and MFM-185 possess heteroacene cores containing additional nitrogen and, in the case of MFM-185, oxygen atoms as well.This arises from synthetic necessity; the longer ligand structures are significantly more synthetically accessible in heterocyclic form than as pure carbocycles.The pure hydrocarbon aromatic analogues of MFM-183 and MFM-185 have proven to be synthetically unobtainable, as have both heterocyclic and carbocyclic versions of MFM-182 and MFM-184.The simulation-based approach taken in this work enables the binding mechanisms of CH4 and CO2 to be systematically investigated not only for the synthesized linkers used in the MFM-18X series, but also for the hypothetical members of the series shown in Figure 2.This provides significant insight into the influence of linker length and, in particular, composition on the interaction between CH4/CO2 and the pore walls of the framework, as well as valuable information for further experimental efforts.

RESULTS AND DISCUSSION
The interaction of single molecules of CO2 and CH4 with a series of organic linker fragments, based on those used in the synthesis of the MFM-18X series of MOFs was investigated via density functional theory with Grimme 3 dispersion correction (DFT-D3) 29 implemented in the Q-Chem software package. 30In addition to the experimentally synthesized ligand systems, the hypothetical naphthalene, anthracene-, tetracene-and pentacene-based ligands were also investigated.
In order to reduce the computational cost of the calculations, the isophthalic acid portions of the linkers were omitted and replaced with single hydrogen atoms (i.e. the R groups shown in Figures 1 and 2 were replaced with C6H5 in all simulations).] As the aromatic core of the ligand -the focus of the present study -is at least two aromatic units removed from the COOH group, however, the influence of the electron-withdrawing group on guest binding in the systems studied herein is likely to be minimal.In referring to the linker fragments studied in this work, we adopt the form Lnc and Lnh, in which n indicates the length of the central aromatic core (n = 0 to 5), c indicates a carbocyclic molecule and h indicates a heterocyclic molecule.
The interaction of the guest molecule with the linker was evaluated in two steps, both using the B3LYP functional. 33An initial geometry optimization of the guest-linker dimer was undertaken using the 6-31+G* basis set, followed by single-point energy calculations using the larger 6-311+G* basis set, from which the binding energy of the guest molecule was calculated following the counterpoise method for the correction of the basis set superposition error. 34Partial charges were determined using the CHELPG method. 35For each guest-linker system, several initial geometries were evaluated and in all calculations the atoms of the linker fragment were kept fixed, while the guest molecule and its constituent atoms were allowed to adjust position upon optimisation.
In order to validate the adopted computational methodology, a range of benchmarking calculations were performed for a simple weakly interacting system -the CO2 -C6H6 dimer.These results are compared with the literature values for the binding energy and equilibrium separation obtained by coupled-cluster theory with iterative singles-and-doubles excitations and a perturbative treatment of triplet excitations, CCSDT(T), in the limit of the complete basis set (CBS) 36 as well as other available data.Table 1 summarizes a comparison between the performance of B3LYP-D3 method, used in this study, with CCSD(T) method and less accurate Moller-Plesset Perturbation theory, MP2.For both configurations of CO2 -C6H6 dimer, 'stacked' (configuration A in Figure 3) and 'end on' (configuration B in Figure 3), B3LYP-D3 results perform consistently better than MP2 and, in the case of configuration A, are in excellent agreement with the 'golden standard' CCSD(T)/CBS results.
A very good agreement for the binding energies calculated using CCSD(T) and DFT-D3 methods has been previously reported 37 for non-covalent interactions of CO2 molecule with a large set of organic linker molecules including nitrogen-containing heterocycles.The authors 37 concluded that both methods are suitable tools to study such weakly bound systems.The performance of B3LYP-D3 hybrid functional augmented with a dispersion term has been also shown to be reliable for description of the interaction between small molecular adsorbates (CH4, H2, N2, CO2, CO, H2O and NH3) with Cu-and Fe-containing coordinatively unsaturated sites (CUSs), with the root-mean-square deviation (RMSD) values being typically below 3 kJ/mol. 20,38 gure 3. Two configurations of CO2 -C6H6 dimer used in benchmarking calculations presented in Table 1.

Table 1. Benchmarking results for the binding energy and equilibrium distance in CO2 -C6H6
dimer comparing the performance of B3LYP-D3, MP2 and CCSD(T) methods.

Method
Basis In the case of the non-heterocyclic linkers (L0c, L1c, L2c, L3c, L4c and L5c), a steady increase in binding strength is observed as the linker length is decreased, reaching a maximum for L1c (Figure 4; top).For both CH4 and CO2, the highest strength binding locations on the carbocyclic linkers were found to be adjacent to the phenyl rings on each end of the linker, either directly over or near the first aromatic ring of the carbocycle (Figure 4; bottom).The primary interaction in this location is between either a positively charged hydrogen (in the case of CH4) or carbon atom (for CO2) with the electron-rich centre of the aromatic ring.Additional, weaker interactions exist between the guest molecule and the edges of nearby terminal phenyl groups.As the separation between the terminal phenyl groups is reduced by removing aromatic units from the centre of the linker, guest molecules are able to interact not only with the nearest pair of phenyl rings but also with those on the opposite end of the linker, resulting in the observed increase in binding strength.Further reduction in linker length from L1c to L0c eliminates the central aromatic component of the linker entirely and results in increased steric hindrance due to the short separation between terminal phenyl groups, forcing the guest molecule further away from the linker (Figure 4; bottom).As a result of these two factors, both CH4 and, in particular, CO2 were found to be much more weakly bound in L0c than in L1c (Figure 4; top).While guest binding on the preferred binding site within the carbocyclic linker series -the aromatic ring nearest the phenyl groups -is weakened as the linker is extended, the introduction of additional aromatic rings results in an increase in the total number of viable binding sites.In the case of CH4, additional binding sites were found both directly above any additional aromatic rings and above the C-C aromatic bonds shared by adjacent rings (Figure 5).The secondary binding sites were found to exhibit very similar binding energies (within 0.2 kJ/mol of each other) but were all found to be approximately 1 kJ/mol (~10%) weaker than the primary binding location for each linker.It should be noted that while these sites are viable in terms of binding energy, many are likely to be mutually exclusive and, particularly in the shorter L0c-L2c systems, the presence of more than one CH4 molecule on each face of the ligand may be sterically impossible.Given the similarity in binding energies observed along the core of the ligand, it is unlikely that CH4 will remain localized on one particular site under ambient conditions and the interaction energy in the real system is therefore likely to be an average of both primary and secondary sites.Little difference was observed in methane binding energy with the secondary binding sites on extension of the linker system from L3c to L4c and L5c, indicating that while an increase in linker length results in a more weakly bound methane molecule at the primary adsorption site, it is not detrimental to the adsorption of methane on secondary sites.The inclusion of a longer ligand is therefore likely to be advantageous, allowing more than one methane molecule to bind with the same ligand without the steric clashes between primary and secondary sites that exist in L0c-L2c, although a full exploration of the many-molecule system would require a classical statistical mechanics approach.
In the case of CO2, secondary binding sites were also observed either directly over any additional aromatic rings present in the linker or over a shared aromatic C-C bond.CO2 was always found to align itself with the O-C-O axis parallel to the plane of the aromatic core of the linker, and generally adopted a position with the O-C-O axis either aligned with or perpendicular to the long axis of the aromatic core (Figure 6).Unlike methane, however, CO2 binds much more weakly to the secondary sites when compared to the primary binding locations and a steady decrease in binding strength is observed for secondary sites as the linker is extended (Figure 6), a result of the more significant electrostatic interactions between CO2 and the phenyl rings at the end of the linkers, which remains non-negligible at a longer separation distance than the primarily dispersive CH4linker interactions.(bottom).Similar binding modes were observed for L3c and L5c (see SI). L0c and L1c do not contain enough aromatic rings to allow for secondary biding sites.
While both CH4 and CO2 behave similarly in the non-heterocyclic linkers with respect to the available binding sites and dependence of binding energy on linker length, the same is not true upon introduction of N-and O-heteroatoms to the system, such as those present in L3h and L5h.
The presence of nitrogen and oxygen within the aromatic core of the linker alters the distribution of charge in the system, introducing regions of significant positive charge (Figure 7).This has the effect of both creating new viable binding sites and significantly altering -and even eliminatingpreviously identified sites.In the case of the heterocyclic L3h linker, the preferred methane binding site was found to be directly over the pyrazine rings of the linker, near the outer phenyl groups, with a binding energy of 8.50 kJ/mol (Figure 8).Further viable binding geometries were observed along the sides of the linker, in which the hydrogen atoms of the methane molecule are able to form weak hydrogen While the location of the strongest methane binding sites on the carbocyclic L3c and heterocyclic L3h linkers remained very similar, the introduction of additional O-heteroatoms in the longer L5h linker significantly alters the preferred CH4 binding locations.For the CH4 -L5h system, the strongest interactions were observed when methane was located directly above the oxygencontaining, electron-rich central ring of the linker (Figure 9).Although a similar configuration was observed in the carbocyclic L5c, methane exhibits a noticeably stronger interaction with the heterocycle (by between 0.5 to 0.9 kJ/mol).In contrast to the L3h system, no dimers were observed in which methane is bound above the plane of the pyrazine rings of L5h (cf.Fig. 8C for L3h), suggesting that the outer phenyl rings play a role in stabilizing the bound methane in the shorter L3h ligand.As in L3h, the interaction between H of CH4 and the heteroatoms of the linker creates a number of new methane binding sites in comparison to the carbocyclic analogue, L5c (Figure 9).
The strongest of these new sites are located near the N-heteroatoms (BE = 6.34 kJ/mol; Fig. 9C), a result of the interaction between C-H and the nitrogen lone pair.Further sites are located near the O-heteroatoms, either via a direct C-H⋯O interaction (BE = 4.47 kJ/mol; Fig. 9D) or a bridging mode, in which methane is seen to interact weakly with both nitrogen and oxygen (BE = 4.76 kJ/mol; SI).As in L3h, the binding energies for these newly created sites exceed kT at room temperature, indicating that they can be expected to play a non-negligible role in methane adsorption at low pressure.Whereas the strongest interactions in the carbocyclic L3c and L5c were found to be directly over the ring systems (BE = 16.06 kJ/mol and 15.17 kJ/mol, respectively), only comparatively weak interactions (about 7-12 kJ/mol) were found for similar conformations in the heterocyclic systems.
As was the case for methane, several new strong CO2 binding sites were identified near the Nheteroatoms (Figure 10), exhibiting binding energies of around 17.5 kJ/mol (L3h) and 16.7 kJ/mol (L5c), approximately 1-2 kJ/mol stronger than any sites observed in their carbocyclic analogues.
The CO2 molecule is able to bind strongly with the linker via a combination of the interaction of the positively charged C of CO2 with the negatively charged N-heteroatoms and interactions between the negative O of CO2 and the regions of slight positive charge above and below the ring.
The overall effect of heteroatoms on the adsorption of CO2 in these systems is likely to be minimal, however, as the introduction of new binding sites along the sides of the linker is, to a certain extent, counterbalanced by the elimination of binding sites above and below the aromatic core.From the point of view of the selective adsorption of CO2 over CH4 at low loading -in which the selectivity is dominated by the interaction of the two species with the framework -the ratio of the binding energies of the two components (BECO2/BECH4) is of greater interest than the binding energies themselves.In the case of primary adsorption sites on non-heterocyclic linkers, CO2 always binds more strongly with the linker than CH4 and the ratio from the calculations was found to vary from 1.21 (in L0c) to a maximum of 1.91 (in L4h).Little dependence of the ratios of binding energies of the five carbocyclic linkers on linker length was observed, with the lowest ratios being observed in L1c and L3c (1.71 and 1.79, respectively), while the ratios in the remaining three linkers ranged from 1.88 to 1.91.When both primary and secondary binding sites are considered, the selectivity of linkers L1c through L5c was found to be almost identical, with binding energy ratios ranging from 1.63 (L5c) to 1.72 (L4c).The presence of heteroatoms in the linkers appears to be beneficial to CO2 selectivity when sites along the sides of the linkers are considered (BECO2/BECH4 = 2.67-2.91).Although CO2 is expected to be selectively bound over CH4 along the edges of the heterocyclic linkers, the re-distribution of charge above and below the aromatic core of the linker is detrimental to CO2 and methane is preferred in sites above the ring system of L3h as a result (BECO2/BECH4 = 0.82).

CONCLUSIONS
In this work, we have presented an extensive computational study of the binding of CH4 and CO2 with a series of real and hypothetical MOF ligands, exploring the effect of linker length and the presence of N-and O-heteroatoms on the affinity of the linker for the two species.It has been shown that the strongest binding sites for both CH4 and CO2 in the carbocylic linkers are found directly above the aromatic rings the linker core, primarily the ring nearest to the outer phenyl rings of the system.The strongest binding on the carbocyclic linkers were observed for the benzene and naphthalene-based systems, in which strong interactions with the π system of the aromatic rings are combined with additional interactions with the nearby outer phenyl rings.The shortest linker, L0c, was found to introduce significant steric restrictions on the binding locations of CH4 and CO2, and both species were only weakly bound as a result.
Although the introduction of heteroatoms to the linkers was found to be beneficial in the case of CH4, the case was not as clear for CO2.In the heterocyclic linkers, CH4 was able to bind in the regions above and below the heterocyclic core of the linker with only a small reduction in binding strength compared to the carbocyclic analogues.][42] The presence of heteroatoms in the L3h and L5h linkers was found to significantly enhance CO2 binding compared to the carbocyclic L3c and L5c.The strongest binding sites, however, were found to be along the sides of the heterocycle and CO2 was no longer able to bind strongly with sites above and below the aromatic core.The effect of incorporation of heteroatoms into MOFs on CO2 adsorption is therefore likely to be highly system dependent, with the largest benefit being in systems where the edges of the heterocycles are accessible to CO2, for example, MFM-18X, 28 NOTT-122 43 and the quinoxaline-based MOFs of Zhu and 44 It is notable that while three quinoxaline MOFs were reported, considerably enhanced CO2 uptake was observed in the only MOF in which the edges of the pyrazine system were accessible.
Comparison of the relative binding strengths of CH4 and CO2 indicates that any of the carbocyclic linkers studied in this work are suitable for selective adsorption of CO2 over CH4 at low pressure but that the presence of at least one central aromatic ring improves selectivity considerably over the single unsaturated C-C bond present in L0c.Comparison of the heterocyclic and carbocyclic linkers suggests that while the presence of heteroatoms can greatly enhance the selectivity towards CO2 for systems in which the edges of the heterocyclic rings are accessible, the presence of heteroatoms tends to favour methane in regions directly above and below the heterocyclic rings and the incorporation of heteroatoms into ligands in which the heteroatom is not able to interact directly with CO2 may, in fact, reduce the selectivity of the MOF for carbon dioxide over methane at low pressure.
ASSOCIATED CONTENT

Supporting Information.
Calculation of pore size distribution; additional CH4 binding sites on L4c, L5c and L5h; additional CO2 binding sites on L3c and L5c; CHELPG charges for L3h and L3c.

S1: CALCULATION OF PORE SIZE DISTRIBUTION
The pore size distributions (PSD) of the MFM-18X series of MOFs were calculated following the method of Gelb and Gubbins 45 , in which the largest spherical probe which is capable of being inserted into the structure is determined via Monte Carlo simulation.In these calculations, the frameworks were assumed to be rigid, with atoms kept fixed at their crystallographic positions.
The framework atoms were described using the OPLS force field. 46For all frameworks, three peaks were identified in the PSD, the largest two being at least 12.4 Å and 14.5 Å respectively.
The central aromatic cores of the linker fragments explored in this work are only accessible to methane and CO2 in these larger cavities.

Figure 1 .
Figure 1.The MFM-18X ligands as previously described (top) and the resulting three-dimensional

Figure 2 .
Figure 2. The hypothetical linker fragments explored in this work.

Figure 5 .
Figure 5. Representative secondary CH4 binding sites on L3c viewed from above and the side of

Figure 6 .
Figure 6.Dependence of CO2 -linker binding energy on linker length for primary (blue) and

Figure 7 .
Figure 7.Comparison of CHELPG partial charge maps of the carbocyclic (top) and heterocyclic

Figure 8 .
Figure 8. Binding sites identified for the CH4 -L3h dimer: A, B) secondary binding sites created

Figure 9 .
Figure 9. Binding sites identified for the CH4-L5h dimer: A, B) preferred binding sites above or

Figure 10 .
Figure 10.Representative binding sites for CO2 on the heterocyclic L3h and L5h linkers: A) Figure S1.Secondary over-bond (top) and over-ring (bottom) CH4 binding sites on ligand L4c.