Analysis of SO2 Physisorption by Edge-Functionalized Nanoporous Carbons Using Grand Canonical Monte Carlo Methods and Density Functional Theory: Implications for SO2 Removal.

Nanoporous carbons (NPCs) are ideal materials for the dry process of flue gas desulfurization (FGD) due to their rich pore structure and high specific surface area. To study the effect of edge-functionalized NPCs on the physisorption mechanism of sulfur dioxide, different functional groups were embedded at the edge of NPCs, and the physisorption behavior was simulated using the grand canonical Monte Carlo method (GCMC) combined with density functional theory (DFT). The results indicated that the insertion of acidic oxygenous groups or basic nitrogenous groups into NPCs could enhance the physisorption of SO2. The influence of edge functionalization on the pore structure of NPCs is also analyzed. To further explore the interaction in the adsorption process, the van der Waals (vdW) interaction and electrostatic interaction between the SO2 molecule and the basic structural unit (BSU) were investigated. Simulated results showed that edge functionalization had limited influence on vdW interaction and did not significantly change the distribution characteristics of vdW interaction. According to the study on electrostatic interaction, edge functionalization was found to promote inhomogeneity of the surface charge of the adsorbent, enhance the polarity of the adsorbent, and thus enhance the physisorption capacity of SO2. More importantly, we provide an idea for studying the difference in adsorption capacity caused by different functional groups connected to carbon adsorbents.

Introduction

SO2 gas is one of the major gaseous pollutants in China and other developing countries.[1] It mainly comes from fossil fuel emissions in thermal power plants and industrial processes. Without the flue gas desulfurization (FGD) process, it may cause acid rain and other natural disasters. According to a study of 74 cities in China, SO2 concentrations declined from 39.9 to 17.0 μg/m3 in 2013–2017 with FGD and coal utilization restrictions.[2] Considering the serious consequences of SO2 emissions, how to effectively carry out the FGD is still an important research direction.[3−7]

Conventional FGD materials are mostly based on the ability of alkaline adsorbents to capture acidic gases, such as sodium hydroxide aqueous solution, limestone–gypsum, etc. However, the above conventional materials have some drawbacks, mainly the difficulty to recycle, probability of causing corrosion of the pipeline, and easy production of liquid or solid wastes. Economical and environmentally friendly FGD materials are of great significance for the energy conservation policy, environmental protection, and resource reuse, so the development of a new adsorbent has attracted wide interest. Some new adsorption materials, such as activated carbon, waste concrete powder (WCP), metal–organic frameworks (MOFs), carbon nanotubes (CNTs), porous aromatic frameworks (PAFs), etc., are considered for the adsorption of SO2.[4,8−15] Compared with conventional adsorption materials, nanoporous carbons (NPCs) have attracted wide attention due to their higher specific surface area, more developed nanopores, lower preparation cost, and more stable adsorption efficiency. By comparing the correlation between saturated sulfur capacity of carbon materials and pore volume of different pore sizes, Sun et al.[16] believed that micropores with the high specific surface area were beneficial to physisorption of SO2. Liu et al.[17] measured the specific surface area of carbon materials by the Brunauer–Emmett–Teller (BET) method and found that the adsorption performance of carbon materials with a specific surface area less than 700 m2/g was mainly affected by their physical structures. On the other hand, the effect of embedding different functional groups in carbon materials on the removal of SO2 in flue gas was studied.

Notably, many researchers have found that adding proper functional groups to NPCs improves the ability to capture gases.[18−21] For example, Raymundo-Piñero et al.[19] found that the embedding of nitrogenous functional groups, especially pyridine nitrogen functional groups, could enhance the adsorption of acidic SO2 gas through experiments. Sun et al.[20] found that the role of nitrogen doping is to change the local electron density of carbon materials, the polarity of carbon atoms, and the charge distribution of the carbon surface, which resulted in enhanced adsorption of SO2. Maurya et al.[21] determined the effect of functional groups on the selective adsorption of SO2 in different compositions of SO2/N2. Note that Kandagal[22] studied the effect of functional groups on the binding of methane by combining ab initio calculations and classical Monte Carlo simulations; the results showed that specific site functionalization could have a significant impact on the local adsorption characteristics. Additionally, some literature studies point out the edge-functional groups can significantly enhance gas binding, making the edges have more potential binding sites in materials containing high concentrations of edge carbon.[23−25] Zhou et al.[25] pointed out that edge functionalization could effectively improve the pore structure and morphological characteristics of NPCs in the investigation of the competitive adsorption of nanoporous carbon, indicating that the embedding position of functional groups had a significant impact on gas adsorption. However, the effect of edge-functional groups on NPCs for SO2 adsorption has not been studied to the best of our knowledge.

Accurate numerical simulation on a small scale helped the abovementioned researchers to investigate the microscopic nature that was difficult to explore through experimental methods. The GCMC method can give the adsorption isotherm of the material on the macroscopic scale, but it is still necessary to explore the deeper adsorption mechanism by employing quantum chemistry on the microscopic scale. Combining the GCMC method with DFT to study the structure and properties of materials will provide a way to explore the mechanism that affects physisorption.[26]

The main content of this paper revolves around the influence of edge-functionalized NPCs on the physisorption of SO2. We construct the structure of this study in the following steps: (1) construct a BSU of edge-functionalized NPCs, optimize its structure by DFT, and then carry out the GCMC simulation to predict the physisorption capacity of edge-functionalized NPCs for SO2; (2) calculate the pore parameters of NPCs and then study the effects of electrostatic interaction and vdW interaction on SO2 adsorption. The effects of electrostatic interaction and vdW interaction between the adsorbate and adsorbent on physisorption capacity are then evaluated by the energy decomposition method.

Computational Methods

As shown in Figure , the BSU of nonfunctionalized NPCs is constructed using coronene-shaped graphene. Four functional groups, hydrogen (−H), hydroxyl (−OH), amine (−NH2), and carboxyl (−COOH), are embedded at the edge of the BSU, respectively. Lu et al.[27] indicated that these four functional groups could effectively enhance the gas adsorption capacity of NPCs. The results of geometric optimization analysis of BSUs show a low sensitivity for the selection of the basis set. The B3LYP hybrid function was selected for geometric optimization of carbon materials in previous studies.[27,28] Therefore, the B3LYP hybrid function and 6-31G(d,p) basis set are selected to implement the geometric optimization of edge-functionalized BSUs.

Figure 1

Hydrogen (−H), hydroxyl (−OH), amine (−NH2), and carboxyl (−COOH) are, respectively, modified on 12 carbon atoms at the edge of coronene-shaped graphene to form BSUs with different functional groups.

Electrostatic potential (ESP) analysis is a very effective method to reveal the electrostatic interaction among molecules. In the molecular system, ESP is defined aswhere the two terms Vnuc(r) and Velec(r) represent the contributions of the nucleus and electrons, respectively; R and r denote coordinates of nucleus A and the point charge, respectively; Z is the charge of nucleus A and ρ(r′) is the electron density of the molecular system at the point r′.[29]

The van der Waals interaction and electrostatic interaction can be considered as complementary relations in the analysis of weak interactions. Following ESP, we defined vdW potential aswhere Vrepul(r) and Vdisp(r) denote exchange-repulsion potential and dispersion-attraction potential in vdW potential, respectively, and R and r denote nuclear coordinates of atoms A and B.[30] Employing the geometric combination rule, the ε (potential well) and R0 (equilibrium distance) terms are evaluated asThe above atomic nonbonded parameters are taken from the universal force field (UFF).[31]

By plotting ESP and vdW potential distributions, we can intuitively describe the electrostatic and vdW interaction characteristics of BSUs, which is conducive for us to get a general understanding of the weak interaction characteristics among study objects. The ESP and vdW potential analyses are carried out in Multiwfn.[32] VMD software is employed for visualization of the obtained results.[33]

The stable configuration of a single SO2 molecule adsorbed on the BSU is set to the gas-framework interaction model. Symmetry adapted perturbation theory (SAPT) and the DFT method are widely adopted to analyze weak interactions.[34−36] We divide the possible adsorption sites of BSU into edge sites and basal planes. The edge sites are located near the embedded functional groups of the BSU, and the molecules adsorbed on edge sites are approximately at the same level as the BSU. The basal planes are located in an area parallel to the level of the BSU.

Geometric optimization of adsorption configurations is at the B3LYP-D3/def2-SVP level. The conventional B3LYP hybrid function is not accurate in calculating dispersion in the research system, so it is necessary to increase the accuracy of the calculation using the B3LYP-D3 hybrid function with dispersion correction.[37] In addition, the use of def2-SVP improved by adding polarization function based on def2-SV ensures the reliability of the calculation results.[38]

The weak interaction energy between the SO2 molecule and the BSU is decomposed into electrostatic interaction and vdW interaction in PSI4 code.[39] The induction also belongs to the weak interaction energy, but its contribution in the physisorption process is generally less than the dispersion effect. The electrostatic and vdW interactions are mainly analyzed in this study without consideration of the effect of induction.[40]

For the construction of NPCs, the commonly used models are the slit-like pore model and the coronene-shaped graphitic BSU random structure model.[41] The latter has successfully predicted the adsorption performance of NPCs for gas, and the calculated adsorption amount is in good agreement with the experimental results.[42] Although this model is random, realizations of this structure can be constructed by different placements of the BSU, and the results (for the same density) appear to have nothing to do with the placement details.[41] As far as this study is concerned, there is another advantage in choosing the random structure model; that is, the functional groups can be easily embedded in the BSU to form various modified NPC frameworks. NPCs are assumed as simulation boxes with periodic boundary conditions. By setting the density of all simulation boxes to 0.542 g/cm3, the number of BSUs added to the simulation boxes can be obtained bywhere M represents the quality of a single BSU. The length of simulation boxes is between 30 and 45 Å to ensure the accuracy and efficiency of the calculation.[43] After all of the parameters of the structure are determined, the Monte Carlo algorithm is used to continuously optimize the structure, and the simulation box with the lowest energy and the most stable state is obtained finally. All of the NPC modeling was carried out by Packmol software.[44]

To optimize the packing of NPC structures, we adopted the Monte Carlo method. First, we randomly placed the BSUs in the simulation box at 298 K. Then, BSUs will perform a Monte Carlo-type movement in the box to produce a new structure. We can thus calculate the acceptance probability of this new structureEnew is the potential energy of the new structure and Eold is the potential energy of the old structure. k is the Boltzmann constant and T is the temperature. If Enew is greater than Eold, the new structure is accepted, otherwise rejected. Repeating these operations 10 000 times, we can finally get the optimized simulation box.

Every single simulation box consists of a set of BSUs as shown in Figure . The physisorption process in every single simulation box is simulated by the GCMC method. The Peng–Robinson equation of the state is chosen to calculate the gas-phase density and experimental fugacity.[45]

Figure 2

(a) A series of BSUs. (b) NPCs-None, (c) NPCs-H, (d) NPCs-OH, (e) NPCs-NH2, and (f) NPCs-COOH are simulation boxes composed of corresponding BSUs.

Table shows the parameter settings in the GCMC simulation. In addition, adsorbate molecules are inserted, deleted, reinserted, rotated, and moved at the same probability in each step.

Table 1

Parameter Settings for GCMC

projects	parameters	units
initialization cycles	10 000
production cycles	10 000
electrostatic	Ewald[46]
potential energy function	Lennard-Jones (LJ)[47]
cutoff distance of vdW interaction	12.5	Å
temperature	298/313/373	K
pressure	0–220	kPa

SO2 is modeled as a rigid molecule with three charged LJ interaction sites. The LJ potential parameters and atomic partial charge were reported by Potoff.[48] For NPC frameworks, we use the rigid model. Information of LJ potential parameters was taken from the UFF.[31] For the atomic partial charge of NPCs, generally speaking, the RESP charge is the most reliable method to calculate the charge. Multiwfn software is used to calculate the RESP charge of NPCs (see Figure S1 and Table S1).[32] In summary, the atomic LJ potential parameters and atomic partial charges are shown in Table .

Table 2

Atomic LJ Potential Parameters and Partial Atomic Charges of Adsorbents and Adsorbates

	gas molecule models		NPCs
atom	S(SO2)	O(SO2)	C	H	O	N
ε (K)	73.80	79.00	29.13	22.12	34.72	34.75
σ (Å)	3.39	3.05	3.40	2.57	3.12	3.26
Q (e)	+0.59	–0.29

All of the interaction parameters conform to the Lorentz–Berthelot (LB) mixing rule, which is defined asIt is worth noting that some researchers believed that the involvement of O2 and H2O would affect the adsorption of SO2 and found that the adsorption type can be approximately considered to be dominated by physisorption under anaerobic and anhydrous conditions.[49,50] In this study, O2 and H2O are not added in the GCMC simulation to maintain that physisorption dominates the adsorption process. Ensuring the accuracy of the GCMC simulation, which is only applicable to systems with physisorption, the simulation of GCMC is carried out by RASPA2 software, and the simulation results are plotted as absolute adsorption isotherms.[51]

Parameters such as accessible surface area and available pore volume are obtained using the Düren and Sarkisov methods.[52,53] The calculation method of the accessible surface area can be simply described as follows: the probe molecule rolls on the adsorbent framework, and the center position of the probe molecule is marked to obtain the trajectory around the entire adsorbent framework and then the accessible surface area is calculated aswhere σ is the sum of the diameter of the adsorbent atom and the diameter of the probe atom, f is the fraction of the probe atom without overlapping points with other atoms in the structure, and the total accessible surface area Aabs is given by the sum of a single area a related to each atom of the adsorbent.

The available pore volume is calculated by the Monte Carlo method. The probe molecule is randomly inserted into the porous material at first, and then the effectiveness of successful insertion is judged. If the interaction energy between the probe molecule and the adsorbent is negative, the insertion is effective. After repeating the insertions, the available pore volume is calculated aswhere N is the total number of insertions, Nsuccess is the number of successful insertions, Nadsorbent is the total volume of adsorbent, V is the available pore volume of porous materials, and the Poreblazer software is employed to calculate the above pore parameters.[54]

Results and Discussion

Absolute Adsorption Isotherms

Maurya et al. studied the adsorption of pure SO2 on a bilayer graphene nanoribbon (GNRs) with different functional molar percentages at 303 K.[21] As shown in Figure , we also calculated the excess adsorption isotherms of SO2 at 303 K and compared the results with those reported by Maurya et al.[21]

Figure 3

Simulation results (this work and the report of Maurya et al.) of excess adsorption isotherms of SO2 at 303 K in carbon-based materials modified by (a) H, (b) OH, (c) NH2, and (d) COOH functional group, respectively.[21]

There is a minor difference in excess adsorption isotherms as shown in Figure . In addition to the differences in the model, we believe that the difference in the functional mole percentage can also cause a slight difference in the results. Obviously, for NPCs-OH, NPCs-NH2, and NPCs-COOH, the functional molar percentage is much larger than that of the materials reported by Maurya et al., which is one of the reasons why the adsorption is greater in the present work.[21] In addition, Maurya et al. ignored the atomic charges of graphene.[21] Compared with our work, the electrostatic interaction between adsorbate and adsorbent is smaller, which may lead to a smaller adsorption capacity. Of course, the force field parameters will also affect the simulation results, but in general, our results are roughly consistent with the results of previous reports of Maurya et al., which has proved the reliability of our method.[21]

By independently analyzing the absolute adsorption isotherms using the GCMC simulation at different temperatures, it can be found that the adsorption amount of SO2 in edge-functionalized NPCs decreases with the increase of temperature in Figure . It is generally believed that chemisorption is caused by the formation of chemical bonds between molecules, and the adsorption rate increases with the increase of temperature. Physisorption is caused by the weak intermolecular force, and the adsorption amount decreases with the increase of temperature. In summary, under the parameters of GCMC set in this study, physisorption is the main adsorption type. In addition, we can estimate the desorption temperature roughly according to the simulation results. When the temperature reaches 373 K, the adsorption capacity of SO2 is very small, and high desorption efficiency can be obtained theoretically. Determination of more accurate desorption temperatures needs repeated experiments.

Figure 4

Absolute adsorption isotherms of SO2 in edge-functionalized NPCs at different temperatures (293, 313, and 373 K): (a) NPCs-None, (b) NPCs-H, (c) NPCs-OH, (d) NPCs-NH2, and (e) NPCs-COOH.

From the comparison of absolute adsorption isotherms of different edge-functionalized NPCs, the SO2 adsorptions of NPCs-OH, NPCs-NH2, and NPCs-COOH are remarkably higher than those of NPCs-H and NPCs-None, indicating that the physisorption capacity of NPCs for SO2 has been significantly enhanced by the modification of acidic oxygenous groups and basic nitrogenous groups at the edge of NPCs. Furthermore, it is found that the SO2 physisorption capacity of edge-functionalized NPCs follows the sequence: NPCs-NH2 > NPCs-COOH > NPCs-OH > NPCs-H. Figure S2 shows the number of adsorbed molecules versus cycles to illustrate that the systems have been sufficiently equilibrated. To ensure that all of the temperature and pressure conditions in the simulation are sufficiently equilibrated, we plotted the number of adsorbed molecules versus cycles as shown in Figure S2 (T = 298 K and P = 220 kPa).

The equilibrium configuration snapshots of SO2 gas at 10 kPa were provided in Figure to further understand the gas distribution adsorbed on edge-functionalized NPCs. First, we find that most of the initial adsorption of SO2 occurs at the edge of the BSUs. In addition, edge functionalization enables some SO2 to be adsorbed near BSUs. In other words, adsorption active sites are added to the region where SO2 molecules cannot be absorbed in the past.

Figure 5

Snapshots of (a) SO2 in (b) NPCs-None, (c) NPCs-H, (d) NPCs-OH, (e) NPCs-NH2, and (f) NPCs-COOH at T = 298 K and P = 10 kPa.

Pore-Size Distributions

Pore parameters have a great influence on the adsorption capacity of materials. Structure parameters used to evaluate pores are (1) available pore volume, (2) accessible surface area, (3) pore limiting diameter (DL) and maximum pore diameter (DM), (4) porosity, Φ (%), (5) pore-size distributions (PSDs).[55] It is worth noting that the reported results are single realization of the optimized simulation box (Table ).

Table 3

Pore Parameters with Different Functional Groups (Gas Probe Molecule = SO2)

	NPCs-None	NPCs-H	NPCs-OH	NPCs-NH2	NPCs-COOH
number of BSUs	50	50	50	50	30
dimensions (Å3)	35.303	35.803	42.303	41.903	42.403
available pore volume (cm3/g)	1.25	1.28	1.34	1.40	1.33
accessible surface area (m2/g)	3516.24	3226.12	3312.96	3426.73	3462.18
DL (Å)	15.60	19.04	19.17	18.37	17.25
DM (Å)	20.31	22.36	24.30	27.06	24.48
porosity, Φ (%)	70.98	70.18	76.44	71.99	76.57

Compared with NPCs-None (1.25 cm3/g), the presence of −H (1.28 cm3/g), −OH (1.34 cm3/g), −NH2 (1.40 cm3/g), and −COOH (1.33 cm3/g) functional groups leads to the enlargement of the available pore volume. The changing trend of the available pore volume is the same as that of DL, DM, and porosity (Φ). In fact, by ignoring some other factors, the density of each model is the same, so the larger the functional groups connected, the less number of filled molecules (with the same volume), and as the distance between each molecule becomes longer, the pore size becomes larger.

The accessible surface area of all of the frameworks is within the range of activated and hypothetical high-surface-area carbons (2000–4600 m2/g), which is larger than that of carbon nanotubes (153–1315 m2/g).[56,57] Compared with NPCs-None (3516.24 m2/g), the accessible surface area (3226.12–3462.18 m2/g) of NPCs with functional groups is slightly reduced. In general, the introduction of functional groups makes the available pore volume and pore size of NPCs expand, but it reduces the accessible surface area of NPCs. By analyzing the PSDs shown in Figure , the change of pore size can be further verified.

Figure 6

PSDs of NPCs with different functional groups: (a) PSDs of NPCs-None, (b) PSDs of NPCs-H, (c) PSDs of NPCs-OH, (d) PSDs of NPCs-NH2, and (e) PSDs of NPCs-COOH.

Figure shows that both NPCs-None and NPCs-H fall within the range of pore size of 15–23 Å, exhibiting a fairly dense pore structure. The PSDs of NPCs-OH (19–24.50 Å) and NPCs-COOH (17–24.50 Å), although estimated with different functional groups, exhibit similar pore sizes, indicating that they all have isolated large pores with broad PSDs. Particularly, NPCs-NH2 is within a pore-size range of 18–27.50 Å, with the broadest PSDs and the largest isolated pore.

Van der Waals Potential

The isosurface map and the color-filled plane map of the vdW potential of this system are shown in Figure . The region close to the BSU is fully enclosed by a green isosurface, indicating that the exchange-repulsion potential dominates the vdW potential in this area. The blue isosurface above and below the BSU indicates that the dispersion-attraction potential dominates the vdW potential in this area.

Figure 7

vdW potential distribution of BSUs; the different types of figures on the left and right are the isosurface map and the color-filled plane map, respectively. In the isosurface map, green represents positive value and blue represents negative value. The color scale is in kcal/mol. Ar is chosen as the probe atom. (a) BSU-None, (b) BSU-H, (c) BSU-OH, (d) BSU-NH2, and (e) BSU-COOH.

As shown in Figure , the distribution of vdW potential for different BSUs is compared, and it is found that the region of dispersion attractive potential dominating vdW potential generally distributes in the basal planes of the BSU, and the strength of the vdW potential does not change significantly. In summary, the effect of edge functionalization on vdW potential is trivial, and the embedding of functional groups has a limited effect on the distribution and strength of the vdW potential.

It is generally believed that vdW interaction represents the weak interaction between two atoms of electrical neutrality, and the most direct relationship with the vdW potential is the geometric structure of the system, while the embedding of functional groups does not significantly change the geometric structure of the BSU. Consequently, the change of vdW potential distribution and strength is not obvious. In summary, edge functionalization has limited influence on the distribution and strength of the vdW potential of the BSU. In other words, the influence of vdW potential is not the main reason why edge-functionalized NPCs can improve the physisorption capacity of SO2.

SAPT Calculations

Figure shows the decomposition of adsorption energy between BSUs and SO2 molecules. EvdW represents vdW interactions and Eele represents the electrostatic interactions. It is found that the EvdW (basal planes: 0.1000–9.1568 kcal/mol, edge sites: 1.4588–48.9754 kcal/mol) of the BSUs is positive in Figure a, which indicates the exchange-repulsion is the dominant factor in vdW interactions. Comparing the EvdW between the basal planes and edge sites, we find that the EvdW of edge sites is larger than that of the basal planes except for the NPCs without edge functionalization. In other words, SO2 molecules adsorbed on edge sites are affected by the larger vdW potential after edge functionalization, which is not conducive to adsorption.

Figure 8

Components of the interaction energy between BSU and the SO2 molecule obtained from SAPT calculations: (a) vdW interaction between BSU and the SO2 molecule and (b) electrostatic interaction between BSU and the SO2 molecule.

On the other hand, the Eele (basal plane: −0.2124 to −13.8384 kcal/mol, edge sites: −2.7505 to −42.1327 kcal/mol) of the BSUs is negative in Figure b. Similarly, we further compare the Eele between the basal planes and edge sites, and it is found that the Eele of edge sites is larger than that of the basal planes, except for the NPCs without edge functionalization. That is to say, the SO2 molecules adsorbed on edge sites are affected by the larger electrostatic potential, which is favorable to adsorption.

Figure a,b shows similar trends, indicating that there is a certain correlation between van der Waals interaction and electrostatic interaction. Combining with the vdW potential distribution of the BSU, it is found that due to the enhancement of electrostatic interaction, the adsorption position of SO2 molecules may be closer to the BSU, causing it to go deeper into the vdW positive potential region.

Compared with the Eele (basal planes: −7.7899 to −13.8384 kcal/mol, edge sites: −8.4679 to −42.1327 kcal/mol) of other types of edge-functionalized BSUs, the Eele (basal planes: −4.0614 kcal/mol, edge sites: −2.7505 kcal/mol) of NPCs-None BSU and the Eele (basal planes: −0.2124 kcal/mol, edge sites: −5.5599 kcal/mol) of NPCs-H BSU are smaller, indicating that NPCs-None and NPCs-H restricted adsorption sites for SO2 molecules and their adsorption capacities are weaker. This conclusion is consistent with the results of the low adsorption amount of NPCs and NPCs-H in GCMC simulations.

By comparing the electrostatic interactions at all positions of the five adsorbents, it is found that the electrostatic interaction at the edge of NPCs-NH2 BSU is the strongest, which is about 15 times stronger than that of the NPCs-None BSU. This finding also adequately explains the results of the maximum adsorption amount of NPCs-NH2 in the GCMC simulation. By comparing the electrostatic interaction between NPCs-COOH and NPCs-OH, it is found that the former is stronger than the latter in both the basal planes and edge sites, which also explains that the adsorption capacity follows the sequence NPCs-COOH > NPCs-OH in the GCMC simulation.

In summary, edge-functionalized NPCs enhance the electrostatic effect on SO2 molecules and supplement the adsorption active sites, and the enhancement of the electrostatic effect is positively correlated with the adsorption amount of SO2.

In addition, we used Multiwfn software to calculate the total adsorption energy under the UFF force field. Figure a shows the total adsorption energy obtained from SAPT calculations, Figure b shows the total adsorption energy calculated based on the universal force field (UFF), and some differences cannot be ignored in the quantitative analysis of Figure a,b. The main reason is that the calculation based on the force field is subject to the simple function form of the force field, and its accuracy is still different from that of the DFT-based calculation. Importantly, the rankings of various BSUs’ adsorption capacities based on force field calculations are generally consistent with the DFT-based calculations, demonstrating the effectiveness of the force field and atomic charge we choose.

Figure 9

Total adsorption energy between BSUs and SO2 molecule: (a) total adsorption energy obtained from SAPT calculations and (b) total adsorption energy based on the universal force field (UFF).

Electrostatic Potential (ESP)

According to the above study, we find it necessary to investigate how edge functionalization causes changes in electrostatic interaction. We perform ESP analysis on the vdW surface with an electron density of 0.001 au (0.63 kcal/mol).[58]

As shown in Figure , the ESP of the whole surface of the BSU without edge functionalization is close to zero, and the charge distribution is relatively uniform. However, the ESP of the whole surface of the edge-functionalized BSU significantly deviates from zero, and the charge distribution is relatively nonuniform. We also find that the ESP on the surface of the H atom is generally positive and that on the surface of the O and N atoms is generally negative. The ESP on the surface of the SO2 molecule shows that the electrostatic potential on the surface of the O atom is negative and that on the surface of the S atom is positive. The electron distribution on the surface of the molecule is nonuniform, which could be regarded as a polar molecule.

Figure 10

Distribution of ESP on the vdW surface with an electron density (ρ) of 0.001 au (0.63 kcal/mol). The color scale is given in kcal/mol. (a) BSU-None, (b) BSU-H, (c) BSU-OH, (d) BSU-NH2, (e) BSU-COOH, and (f) SO2.

To quantify the influence of edge functionalization on the polarity of BSUs, we calculated the molecular polarity index (MPI) of different BSUs, which is defined aswhere A refers to the whole area of the vdW surface of the molecular system and V(r) is the value of the electrostatic potential at a point position r in space. The integration is carried out on the whole molecular vdW surface (S).[59]

Figure a shows the MPI of different BSUs. It is found that the MPI (3.5431 kcal/mol) of the BSU without edge functionalization is relatively small. The MPI (8.3623 kcal/mol) after embedding hydrogen and the MPI (8.0821 kcal/mol) after embedding hydroxyl increased by nearly four times, and the MPI (15.5773 kcal/mol) after embedding amino and the MPI (15.4510 kcal/mol) after embedding carboxyl increased by nearly six times, indicating that edge functionalization is very significant for the enhancement of the polarity of BSUs.

Figure 11

(a) Molecular polarity index (MPI) of BSUs with different functional groups, (b) average strength and size of the positive electrostatic potential regions, and (c) average strength and size of the negative electrostatic potential regions.

It is worth noting that the MPI after embedding hydrogen and the MPI after embedding hydroxyl are approximately equal, indicating that both polarities are of little difference. However, the adsorption capacity for SO2 follows the sequence of NPCs-OH > NPCs-H from the results of the GCMC simulation, indicating that the adsorption capacities of the materials with similar polarity are still different. Similarly, the same conclusion can be obtained by analyzing the BSUs embedded in amino or carboxyl. To study the reasons for such difference, we analyzed the average strength and size of the positive and negative electrostatic potential regions of the above molecules.

Figure b shows the percentage of positive areas and the positive average values for different molecules. Figure c shows the percentage of negative areas and the negative average values for different molecules. It is found that the percentage of the positive area (51.22%) of the SO2 molecule is slightly larger than that of the negative area (48.78%), and the positive average value (28.9407 kcal/mol) is larger than the absolute average negative value (23.3097 kcal/mol), indicating that the SO2 molecule is easier to combine with negatively charged Lewis base.

By observing the BSUs with similar polarities, it is found that the negative average values of NPCs-H BSU and NPCs-OH BSU are similar, and the percentage of the latter (68.93%) is much larger than that of the former (55.76%) in negative areas; thus, it makes SO2 molecules easier to combine with the NPCs-OH BSU. When the percentages of negative areas of NPCs-NH2 BSU and NPCs-COOH BSU are similar, the average strength of the former negative electrostatic potential (20.4579 kcal/mol) is much larger than that of the latter (15.6548 kcal/mol), so NPCs-NH2 BSU is easier to combine with SO2 molecules.

Conclusions

The adsorption behavior of edge-functionalized NPCs at different temperatures and pressures is simulated by the GCMC method. The results show that the adsorption capacity follows the sequence: NPCs-NH2 > NPCs-COOH > NPCs-OH > NPCs-H > NPCs-None. Simulated results about the pore structure of NPCs show that the insertion of acidic oxygenous groups or basic nitrogenous groups will enlarge the pore size of NPCs, but it will also cause a slight reduction in the accessible surface area of NPCs. The SAPT calculation results show that the adsorption position of SO2 molecules is closer to the BSU due to the enhancement of electrostatic interaction, which causes the molecules to penetrate the vdW positive potential region and results in more adsorption active sites. It is also found that the enhancement of the electrostatic effect is positively correlated with the amount of SO2 adsorbed. By analyzing the distribution of electrostatic potential, the surface charge distribution of the edge-functionalized BSU is more nonuniform than that of the non-edge-functionalized BSU.