Pore-filling contamination in metal–organic frameworks †

Tuneable pore sizes, ordered crystal structures, and large surface areas are some of the main attractive features of metal–organic frameworks (MOFs). To fully understand the structure–property relationships of these materials, accurate characterisation of their structural features is essential. The surface areas of MOFs are routinely estimated from the physical adsorption of gases. By applying the Brunauer, Emmett & Teller (BET) theory to an adsorption isotherm, the surface area is calculated from the amount of gas that forms a monolayer on the pore surface. While this technique is used ubiquitously within the porous solid community, its accuracy can be greatly affected by pore-filling contamination. This process causes an overestimation of the BET surface area from the overlap of surface and pore-filling adsorption as molecules that are not in contact with the surface are erroneously included into the surface area calculation. Experimentally, it is rather challenging to examine the effects of pore-filling contamination, which typically rely on accurate atomistic simulations to provide insight. In this work, we employ grand canonical Monte Carlo simulations and other theoretical approaches to assess the impact of pore-filling contamination on MOF surface areas. With a focus on the rht and nbo topologies, we show how experimental studies that suggest MOF surface areas can be increased by replacing phenyl rings for alkynes are largely influenced by the pore-filling contamination effect.


Introduction
Metal-organic frameworks (MOFs) represent a large group of porous coordination polymers that have received much attention for their potential applications in gas separation, storage and sensing as well as in catalysis and electronic devices. [1][2][3][4][5][6][7][8] These hybrid materials contain metal ions or clusters that are bridged together by polydentate organic ligands, 9 forming a series of three dimensional crystalline frameworks which typically possess ultra-high surface areas and porosities. 10 Unlike traditional nanoporous materials and minerals, such as zeolites, there is a wider variety of building blocks to choose from when developing new MOFs. 11 Therefore, a greater range of pore architectures and sizes can be achieved through judicious selection of the inorganic and organic building blocks, providing a more rational approach to design. [12][13][14] The wide variety of inorganic and organic moieties makes it hypothetically possible to design a near-infinite number of metal-organic frameworks. Great flexibility in design of MOFs structure is highly desirable for selective functionality of porous hosts in adsorption based applications, such as energy storage 15 or carbon sequestration. 16 For these applications, it is crucial that MOFs materials provide high and efficient uptake of adsorbate, and the use of open metal sites [17][18][19] has been shown to provide strong adsorption sites at low pressures, where uptake is governed by the strength of hostguest interactions. In addition, low pressure uptake can be improved through the inclusion of functionalised ligands 20,21 or through post-synthetic modification of the framework. 22 At higher pressures, uptake is correlated with the amount of the surface area and pore volume provided by the adsorbent, and the accurate estimation of these structural properties becomes increasingly important in understanding the mechanisms of adsorption.
The surface area of MOFs is routinely determined by applying the Brunauer, Emmett and Teller (BET) theory 23 to nitrogen adsorption isotherms. The BET theory aims to explain the physical adsorption of guest molecules through the formation of adsorption layers. The surface area of a MOF is calculated by determining the amount of guest molecules, which form a single adsorbed layer across the pore surfaces. This quantity, referred to as the monolayer capacity, is typically extracted from nitrogen adsorption experiments performed on the activated MOF sample. Alternatively, a nitrogen isotherm can be simulated using the grand-canonical Monte Carlo method (GCMC), which explores adsorption thermodynamics of classical systems at equilibrium. Using either method, the monolayer capacity can be extracted from the adsorption isotherms using the linearized BET equation: where, p 0 and p are the saturation and equilibrium pressures, respectively, v is the amount of gas adsorbed at the equilibrium pressure p, v m is the monolayer capacity, and c is a constant related to the heat of adsorption of adsorbed monolayer. Using eqn (1) (3); a tolerance of 20% is suggested for this latter criterion. From the appropriate pressure range, the value of v m can be readily determined and used to calculate the specific surface area (SSA), i.e. the surface area per gram of MOF, using the following equation: where A BET is the SSA, s is the nitrogen adsorption crosssectional area, N A is Avogadro's constant and m s is the mass of the MOF sample. In most calculations, the monolayer is assumed to form a close-packed state, making s = 0.162 nm 2 . However, some surfaces can influence the packing of nitrogen more strongly than others, creating a range of s values for different matierals. 24 As eqn (2) shows, the SSA is directly influenced by the nitrogen adsorption cross section, and any change to its value will, therefore, affect the SSA. Although the BET method remains a popular characterisation tool, the underlying theory relies on numerous assumptions which, for some adsorbents, oversimplify or incorrectly describe the adsorption process. 24,26,27 For example, MOFs containing multiple pore sizes violate a core assumption in BET theory that states that all adsorption sites are homogenous. This is because at low pressure an adsorbate will preferentially adsorb in a smaller pore where there are a greater number of MOF-adsorbate interactions. Therefore, alternative methods that can efficiently and accurately determine surface areas of porous materials are highly desirable. In 2004, Düren et al. used Monte Carlo (MC) methods to determine the geometric surface areas of various porous adsorbents. 28 In each MC simulation, an atom-sized probe is randomly inserted across the surface of an adsorbent atom numerous times. If an insertion does not overlap with any other adsorbent atoms, it can be used to determine the accessible surface area (ASA) of that atom. By repeating this process for each adsorbent atom, the ASA of the entire adsorbent can be determined based on the number of accessible insertions and the size of the probe. 29 As the ASA is dependent on the probe size, it is appropriate to use a nitrogen-sized probe when compared with the BET surface area.
As organic ligands make up most of a MOFs structure, their shape and size contribute significantly to the amount of surface area formed. By expanding the size of the organic ligand, a framework with larger surface area, pore volume and porosity can be created. In many examples, this expansion is achieved by bonding multiple phenyl rings together, forming part of an isoreticular series of MOFs. 14,30,31 While this method is effective, previous work has shown that using multiple, ''area efficient'' alkyne groups is a better approach to increase MOF surface areas. 32 By employing the latter strategy, two new MOFs, NU-109 and NU-110, have been produced with BET surface areas exceeding 7000 m 2 g À1 ; two of the largest values recorded in the literature. Since this publication, 32 a number of MOFs and other porous materials have been designed with alkyne units forming a core part of their structures. [33][34][35][36] One of these structures, NU-111, developed by Hupp and co-workers, 37 possesses large dendritic ligands that contain six alkyne linkers coordinated through copper paddlewheel nodes to form a cubic framework with rht topology (Fig. 1). Hupp and co-workers find that NU-111 has a BET surface area of about 5000 m 2 g À1 , almost 20% larger than that of NOTT-119, 38 an isoreticular MOF that replaces each alkyne group with a phenyl ring (Fig. 1a).
This indicates that replacing phenyl rings with alkyne groups may be more beneficial for increasing the SSA of MOFs. In a later study, these conclusions were supported by Bichoutskaia et al. 39 who computed the ASAs of a hypothetical series of MOFs with rht topology. For seven MOFs structures considered in this study, the results also indicate that substituting phenyl rings with alkynes groups leads to a general increase in specific surface area.
In earlier studies, Zhou and co-workers synthesised PCN-16 40 and PCN-46 41 MOF structures with nbo topology, which is formed by coordinating tetra-topic carboxylate ligands through copper paddlewheel nodes. The resulting framework possesses a rhombic unit cell formed from a combination of short  (Fig. 2b). The ligands used to form PCN-16 and PCN-46 were designed with one and two alkyne moieties, respectively (Fig. 2a). By replacing these alkynes moieties with a single phenyl ring, the isoreticular NOTT-101 MOF is formed as synthesised by Schröder and co-workers. 14 The BET surface areas of PCN-16, PCN-46 and NOTT-101 are 2273, 2500 and 2316 m 2 g À1 , respectively, further demonstrating examples where MOF SSAs may be increased through the replacement of phenyl rings with alkyne functional groups.
While, the above observations support the effectiveness of the adopted strategy, recent work 42,43 has demonstrated that many MOF BET surface areas can be largely overestimated, despite fulfilling the Rouquerol consistency criteria. The effect responsible for this inaccuracy is known as pore-filling contamination (PFC) that typically affects MOFs containing multiple pores, which adsorb nitrogen at different relative pressures. PFC can lead to an overlap of monolayer formation with porefilling, causing large amounts of nitrogen to be included into the monolayer capacity calculation, despite being located in the centre of the pore cavities. The overall effect of PFC is an overestimated value of the BET surface area. Experimentally, it is extremely challenging to establish and quantify how PFC affects the BET calculation, as the exact locations of each adsorbed nitrogen molecule are not known. In recent years, GCMC simulations have been extensively used as a tool to study adsorption thermodynamics and simulate isotherms for MOFs. 44 The quality and accuracy of these simulations are determined by the choice of force-field used to describe the chemical interactions within the system. If a good match is observed between the experimental and simulated isotherms, then it is assumed that the force-field provides an accurate description of the underlying chemical interactions within the system, 45 allowing exact positions of each adsorbate molecule to be extracted from simulation. For this reason, GCMC simulations provide a reliable method for assessment of PFC effects on calculations of the BET surface area of MOFs.
In this study, we determine how PFC affects the SSA calculation of nine different MOFs within the nbo and rht topologies. 14,37,38,40,41,[46][47][48] These MOFs ( Fig. 1a and 2a) were chosen to also study how SSA changes when replacing phenyl rings with alkynes. We employ GCMC simulations to simulate the nitrogen isotherms and use the BET method to determine the surface areas of each investigated MOF. In addition, the ASA is calculated for comparison with the BET surface areas.
With the obtained results, we discuss what impact this may have on the observed experimental trends.

Results and discussion
The experimental and simulated surface areas of the nine investigated MOFs are recorded in Table 1. A list of computational details used for each simulation may be found in the ESI. † In general, the BET method provides good agreement between the experimental and simulated surface areas, with simulations frequently predicting a larger surface area. This occurs because simulations use perfect crystal structures that are free of defects, adsorbed solvent molecules or collapsed regions of framework, all of which may be found to some extent in experimental structures. In the case of NOTT-119, we find a larger than anticipated discrepancy between the experimental and simulated BET surface areas (Fig. S19, ESI †). We observe that for NOTT-119 the experimental isotherm shows typical type IV behaviour with a pore filling step at p p 0 ¼ 0:25. However, our simulations predict a sharp increase in nitrogen adsorption at p p 0 o 0:08 followed by a plateau in loading. This shape is more consistent with a type Ib isotherm, which is usually recorded for materials with wide micropores. 49 The nitrogen monolayer is however predicted to form at p p 0 % 0:11, which correlates to a 68% or 560 cc(STP) g À1 difference in nitrogen loading. Simulations using generic force-fields such as OPLS 50 and UFF 51 also fail to capture the correct experimental isotherm shape, however the simulation snapshots that we collected remained useful for the trend analysis demonstrated later in this paper (Fig. S30, ESI †). The data on surface areas recorded in Table 1 show that for the rht topology, the BET method tends to predict larger values than the geometric method, whilst the opposite trend is observed for the nbo topology. These trends arise as a consequence of pore-filling contamination, which leads to overestimation of BET surface areas in some MOFs. To verify this, we generate snapshots from the GCMC simulations at the pressure of monolayer formation, as determined from the consistency criteria. By measuring the distance of each nitrogen molecule from the framework atoms, we are able to distinguish between monolayer and pore-filling molecules (Fig. 3), using a distance cut-off technique described previously. 43 The number of nitrogen molecules contributing to PFC can then be calculated.
For MOFs with the nbo topology, we observe small amounts of PFC which can be attributed to the narrow widths of micropore, shown in the pore size distributions (Fig. S1, ESI †). The surface adsorption of nitrogen in these MOFs saturates the majority of the pore space, leaving space for only a few additional molecules to adsorb. Of the four MOFs, PCN-16 possesses the smallest amount of PFC, in which 3.00 nitrogen molecules per unit cell (2.6%) are erroneously included into the calculation of the monolayer capacity. As the linker size is increased, we observe a small increase in the amount of PFC affecting the BET calculation. NOTT-101, PCN-46 and NOTT-102 contain 5.00, 8.00 and 9.25 nitrogen molecules per unit cell, respectively, which contribute to PFC; positions of these molecules are typically found to reside at the centre of each pore as illustrated in Fig. 4.   In contrast to these trends, the snapshots of the rht topology MOFs show much larger proportions of PFC. These MOFs contain three different pore types, whose diameters range between 1.0 and 2.4 nm (Fig. S2, ESI †). Consequently, the beginning of monolayer formation in the larger pores can overlap quite strongly with pore-filling or saturation in the smaller pores, leading to high levels of PFC. Although the consistency criteria were used for all calculations, the snapshots of the nitrogen monolayer show that pore-filling has already begun in the largest pore, before the completion of an adsorbed monolayer (Fig. S31-S35 Although each monolayer capacity has been derived from linear regions that fulfil or minimise the error within the consistency criteria, we have shown so far that these values of BET surface area can still be susceptible to error from PFC. This evidence from our simulations prompted a more accurate determination of the BET surface area by removing the porefilling contaminant nitrogen molecules from the monolayer capacity calculation (labelled as method 1 in Table 1). For the nbo topology MOFs, the corrections to the original BET surface area are quite subtle, ranging from 15-180 m 2 g À1 . However, for the rht MOFs, large reductions were observed, with the SSAs of NOTT-119, NOTT-116 and NU-111 decreasing by more than 1000 m 2 g À1 . Whilst the PFC correction leads to a general improvement, in the nbo MOFs there are still some noticeable differences between the values of the ASAs and these corrected BET surface areas. In addition, we now observe that the corrected BET surface areas of all investigated MOFs are smaller than their corresponding ASAs. As discussed previously, most BET calculations assume that the adsorbed monolayer occupies the surface in a closed packed state. However, this packing may be changed or influenced across different surfaces leading to a different value of the nitrogen adsorption cross-section area, s. In previous work, 43 it has been shown that the accuracy of this assumption may be improved by using snapshots at the saturation pressure, p 0 . At p 0 , it is expected that the nitrogen adsorbed in the pores packs like a liquid, i.e. s value is more accurate, and can be used to estimate the maximum number of nitrogen molecules that form a close packed monolayer. This method has also been shown to give very similar surface areas to the geometric ASAs, over a range of different topologies. The monolayer capacity from this method can be calculated by counting the total number of nitrogen molecules in contact with the pore walls at p 0 . As such, the method relies very little on the assumptions of the BET theory and the calculation becomes very similar to that of a geometric calculation. The surface areas determined using these monolayer capacities correspond to a more accurate upper limit of surface area, that can be achieved from a perfect crystal structure. We chose to use this method in our work as it also allows us to observe the difference that is made to the PFC corrected surface areas when the packing of the monolayer is also considered (labelled as method 2 in Table 1).
Using method 2, an excellent agreement between the ASAs and the PFC corrected surface areas is observed for all nine MOFs (Fig. 5). This suggests that these surface areas provide a more accurate characterisation than those affected by PFC. For a few MOFs, we observe that the PFC corrected surface areas are larger than their geometric ASAs. This occurs because the ASA calculation uses a hard-sphere probe, which for some MOFs cannot fully access certain areas of the framework. However, these areas may be fully accessible to the softer nitrogen molecules used in GCMC simulations, leading to an overall larger surface area (Fig. S36, ESI †). While the agreement between the geometric method and method 2 is very good, we note that computational expense of method 2 is much greater than it is for calculating an ASA. However, the use of a hard-sphere probe for the ASA calculation does not always provide accurate surface areas. As a result, method 2 can be used to verify the ASA results and it provides an alternative method for calculating accurate surface areas when the geometric method fails to do so.
With the excellent agreement between the ASA and the PFC corrected surface areas from method 2, we can now compare the surface areas of the investigated MOFs. Table 2 contains both experimental and simulated results, which measure the changes in surface area between MOFs containing phenyl moieties and MOFs that replace them for alkyne groups. A  Table 2 implies that the SSA is increased upon phenyl substitution for alkynes, whilst a negative change implies a decrease in the SSA. We note that the comparison of two experimental surface areas can be influenced by both PFC and the quality of the sample. The presence of defects, adsorbed solvent or collapse regions of framework may be present in one sample more than others, leading to a biased comparison. However, by using perfect crystals and accounting for pore-filling contamination, a more balanced comparison can be made, which in this work can be used to assess the efficiency of replacing phenyl moieties with alkyne groups. Due to the small amounts of pore-filling contamination present in each nbo topology MOF, a good general agreement between the experimental and simulated BET surface areas is already observed, regardless of any PFC corrections to the surface area (Table 2). In contrast, the higher levels of pore-filling contamination present in each rht topology MOF strongly influences how efficient alkynes are at boosting SSAs. In these MOFs, we expect that the comparisons generated from the ASAs and PFC corrected surface areas provide a more reasonable assessment of this strategy, which appears to be less efficient than previously suggested through experimental analysis. For both topologies, the changes to the SSA are strongly dependent on the number of substituted alkynes. To explain these trends, the SSA is represented as a product of the number of unit cells per gram and the total surface area of a MOFs unit cell. MOFs are periodic structures represented in a simulation by a unit cell infinitely repeated in three dimensions. The SSA of these crystals can be thus calculated as the product of the surface area of a unit cell and the number of unit cells per gram of MOF material. The latter quantities can be directly extracted from both the crystal structure data and the computation of the surface area, making it simple to compare how mass and total surface area change when substituting phenyl rings with alkyne groups (Table 3).
A single phenyl ring has more than three times the mass of an internal alkyne. Therefore, by replacing any phenyl with three or less alkynes, the mass of the unit cell is reduced and hence, an increase in the number of unit cells per gram is observed. However, the SSA is also dependent on the total surface area of the unit cell. In general, this is affected by the size of the ligands and their accessibility for adsorption. As each alkyne is smaller than each phenyl ring, the coupling of multiple alkynes is required to increase the overall unit cell surface area. For example, PCN-46 has two alkyne groups and NOTT-101 has one phenyl ring per ligand. By coupling two alkynes together, PCN-46 has a lower mass and longer ligand     (Table 3). When a phenyl ring is substituted for a single alkyne, the decrease in mass is generally too small to offset the decrease in the unit cell surface area, resulting in an overall smaller SSA. However, NU-111 shows exception to this trend when compared with NOTT-119. From experiments, the BET surface area of NU-111 is suggested to be approximately 20% larger than NOTT-119. However, our simulations show that the BET SSAs of NU-111 and NOTT-119 can be overestimated by PFC, suggesting that these MOFs possess very similar SSAs. This result is also confirmed by the ASAs, which show a small 83 m 2 g À1 difference between the two MOFs. While, the total unit cell surface areas differ by B7000 Å 2 , the reduced unit cell mass of NU-111 provides an additional 5.4 Â 10 18 unit cells per gram, offsetting the reduction to the unit cell surface area and providing a minor increase in the overall SSA.
From our results shown in Tables 2 and 3, it is clear that the substitution of phenyl rings for alkynes is a valid approach to increasing SSA in two scenarios. The first scenario is when a single phenyl ring is substituted for two or more alkyne groups. This provides an increase in unit cell surface area and a decrease in mass, resulting in a large SSA. However, if the reduction in mass is large enough to compensate the reduction in unit cell surface area, it is also possible to observe increases in the SSA from substituting phenyls and alkynes in a 1 : 1 ratio. Although, this remains a valid approach to increasing SSA, our simulations show that this process is not as efficient as previously suggested from experimental BET analysis.

Conclusions
In summary, GCMC simulations were employed to determine how pore-filling contamination affects the specific surface area calculations of nine MOFs. Whilst, the experimental and simulated BET surface areas generally agree, a poorer agreement between the simulated BET surface areas and accessible surface areas is observed. Assuming this to be the effect of PFC, GCMC simulation snapshots were used to measure the amounts of PFC in each MOF by calculating the distance of each nitrogen molecule from the pore walls. The molecules which do not adsorb onto the pore walls were classified as pore-filling contaminants and deemed responsible for the over-estimation of the BET surface area. As each MOF with nbo topology considered in this work possesses narrow micropores, PFC provided a negligible difference to their BET surface areas. However, a much larger overestimation of the surface area was found for rht topology MOF due to the multiple pore sizes which are susceptible to large amounts of PFC. We correct the surface area for the effects of PFC by using two methods. In the first method, a simple removal of each pore-filling contaminant was made, and the monolayer that remained in the snapshots at the predicted monolayer formation pressure was used to recalculate the surface area. Due to the poor assumption of the monolayer packing in this method, a second method which used snapshots from the saturation pressure was employed. By removing the molecules that were classed as pore-filling contaminants in these snapshots, a denser monolayer was observed which provided a more accurate computation of the surface area, using a similar calculation to that of the geometric method. As a result, the PFC corrected surface areas that were obtained from the second method were found to be in excellent agreement with the ASAs, validating the results we obtained from both of these methods. In addition, the short amount of time required to assess PFC in a GCMC snapshot suggests that this analysis could be implemented to be routinely calculated after any nitrogen adsorption simulation. While, the experimental BET surface area of a MOF may be affected by its crystal quality and the frameworks susceptibility to pore-filling contamination, these simulated surface areas provide a fair even comparison between each investigated MOF, eliminating any deceptive trends that may arise from experimental analysis. Overall, our simulation results show that replacing phenyls with multiple smaller alkynes is a valid approach to increasing the SSA when multiple alkynes are used. However, the relative gain in SSA is much smaller than predicted from experimental results. This is rationalised by expressing the SSA as a product of the number of unit cells per gram and the total unit cell surface area, showing how alkynes provide much lighter unit cells, which is the main driving force for the gains in SSA.

Conflicts of interest
There are no conflicts to declare.