Interface-enhanced CO₂ capture via the synthetic effects of a nanomaterial-supported ionic liquid thin film

Abstract

Ionic liquids (ILs) are effective CO₂ capture media and recent experimental evidence has demonstrated that the addition of two-dimensional (2D) nanomaterials into ILs can effectively improve their CO₂ capturing capability. However, an in-depth mechanism on how 2D nanomaterials enhance CO₂ absorption is poorly documented. In this study, the adsorption of CO₂ by a representative IL, namely 1-ethyl-3-methyl-imidazole-tetrafluoroborate ([EMIM][BF₄]), coated on graphene (GRA, the prototype 2D nanomaterial) and nitrogenized graphene (C₃N) was investigated by molecular dynamics simulations. The influence of the IL film thickness on the amount of CO₂ adsorption was systematically analyzed. Our data clearly indicate that at the IL-gas interface the CO₂ accumulation is significantly enhanced. In contrast, at the IL-GRA and IL-C₃N interfaces, only slight enhancement was observed for CO₂ accumulation. Quantitative calculations of the adsorption-free energy for CO₂ inside the IL film further support the simulation results. Our present results also reveal that the sub-nanometer IL film possesses a considerably high CO₂ capture efficiency because of the formation of the reduced bulk IL region. Moreover, the nanomaterial substrate surfaces can effectively accelerate the diffusion of CO₂, which is beneficial for the CO₂ mass transfer. In general, our theoretical study provides a deep microscopic understanding of the CO₂ capture by nanomaterials and IL composites. These results could benefit the design and fabrication of a high-performance CO₂ capture and storage medium through the synthetic effects of ILs and nanomaterials.

---

1. Introduction

To slow the pace of global warming, considerable efforts have been made towards the development of novel strategies for mitigating CO₂ emission, particularly for the CO₂ capture and sequestration (CCS) in the past two decades.^1,2 Among numerous materials under study, ionic liquids (ILs) have received remarkable attention from both industry and academia since the first report of using ILs in CCS.³ The outstanding capability of ILs to attract CO₂ lies in their unique physical and chemical properties including the high affinity towards CO₂, the vanishingly low vapor pressure, the high thermal and chemical stability and the low energy consumption for regeneration.^4–6 Benefitting from these distinguished features, ILs are widely accepted as promising candidates for being economic and environmental-friendly capture-medium for CO₂. More importantly, the advantages of ILs are reflected in the abundant types of the composed cations and anions in ILs. Thus, the properties of ILs can be customized to match different CCS requirements.⁷

Despite the high performance of ILs for CO₂ capture, severe drawbacks still exist, which hinder their industrial usage. In particular, the high viscosity induced by the hydrogen-bond network inside ILs unavoidably slows down the kinetic speed for CO₂ uptake and separation.⁸ In addition, it is also reported that the absorption of CO₂ can in turn increase the viscosity of ILs.⁹ To overcome this drawback, one solution is building supported IL membranes (SILMs) by coating ILs on the surface of solid materials.^10,11 In previous studies, researchers found that the thickness of IL films is a key factor to attract CO₂. Earlier experiments have shown that when the thickness of an IL membrane decreases from the microscale to nanoscale, the CO₂ mass transfer rate has a hundred-fold increase.¹² This, according to the diffusion-reaction theory,¹³ is because the diffusion-controlled process in the thick layer is switched to the surface-reaction controlled process when IL thickness decreases. It is also reported that the thin IL film on the substrate can enhance the confinement effect, which accelerates its diffusion and the adsorption of CO₂.^14,15 Despite these previous studies, the molecular mechanism of how IL molecules interact with the CO₂ gas in the confined space is still poorly documented, which becomes the main hindrance for further experimental research.

Here, molecular dynamics (MD) simulations have been conducted to explore the molecular mechanism of SILM for CO₂ capture. The commonly used 1-ethyl-3-methyl-imidazole-tetrafluoroborate [EMIM][BF₄] was used as the representing IL model. However, for nanomaterial substrates, we choose graphene (GRA, a prototype 2D nanomaterial) and C₃N (a recently synthesized 2D structure,¹⁶ which was used for controllable CO₂ capture and release^17–19), as representative models. The influence of the IL film of various thicknesses in CO₂ capturing was systematically addressed.

Our results clearly indicate that the GRA/C₃N substrates introduced two interfacial faces in the IL film, the substrate-IL interface and IL-gas interface, with distinct IL ordering structures compared to that in the bulk IL. The CO₂ capturing efficiency of the supported IL increases with a decrease in the IL film thickness. Among two interfacial regions, the clear accumulation of CO₂ is found at the IL-gas interface region, which plays a key role in enhancing the CO₂ capturing capability than that with bulk IL. In contrast, the substrate-IL interfacial region plays a minor role in CO₂ capture. However, the substrate-IL interface is beneficial for the diffusion of IL and CO₂, particularly for graphene, which is beneficial for the CO₂ mass transfer. More importantly, our results also reveal that the sub-nanometer IL film possesses a significantly higher CO₂ capture efficiency because of the formation of a reduced bulk IL region. Our present findings provide a deep understanding of the interactions between the SILM components with CO₂, which are important to guide the design and optimization of SILMs for the CCS usage.

2. Simulation methods

2.1 Model setup

The GRA nanosheet has a dimension of 2.993 × 3.455 nm² and is composed of 392 carbon atoms. The C₃N nanosheet has a dimension of 3.104 × 3.584 nm² and is composed of 96 nitrogen and 288 carbon atoms. Previous studies have shown that compared to the monolayer substrate the multilayer substrates can affect the adsorption of molecules in a limited manner.²⁰ Thus, in our simulation, we only focused on the monolayer substrate. The nanosheet is placed parallel to the x–y plane of the simulation box and the cross-section of the simulation box is equal to the size of the nanosheet. A space of 10 nm is placed along the z direction. Five systems with different numbers of [EMIM][BF₄], which were 20, 50, 80, 110 and 140, coated on the surfaces of GRA and C₃N were prepared, respectively. The thickness of the IL film was estimated from the 2D layer to the IL-gas interface with the IL density being half of the bulk value. After equilibration of 10 ns, the IL film thicknesses were estimated to be 0.57 nm, 1.57 nm, 2.41 nm, 3.26 nm and 4.09 nm for the five GRA systems, and 0.57 nm, 1.43 nm, 2.21 nm, 3.00 nm and 3.76 nm for the five C₃N systems, as demonstrated in Fig. S1 in the ESI.† Then, 200 CO₂ molecules were randomly added to the vacuum region of each system, as depicted in Fig. 1c and d.

Fig. 1 Structures of cation [EMIM]⁺ (a) and anion [BF₄]⁻ (b), and schematic representations of the simulation box for IL coated on GRA (c) and C₃N (d).

2.2 Simulation parameters

All MD simulations were performed using the GROMACS 5.1.3 package.²¹ The force field parameters for [EMIM][BF₄],²² CO_2,¹³ and C₃N²³ were adopted from previous studies in which they were optimized by precise quantum chemistry calculations. During the simulation, a leap-frog algorithm was used to integrate the Newton equations of motion, and the time step was set to 1 fs. All bonds involving hydrogen atoms were constrained using the LINCS algorithm.²⁴ All the atoms in GRA/C₃N were position-restrained to maintain a planar surface during the whole simulation. The cutoff for both the electrostatic and van der Waals (vdW) interactions was 1.0 nm. Long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) method.²⁵ All the simulations were performed under the NVT ensemble, where the V-rescale algorithm²⁶ was applied to the system to maintain a constant temperature of 298 K. Totally 200 ns production trajectory for each system was collected for CO₂ molecules to dissolve into the IL film, and the last 30 ns trajectories were used for the data analysis.

2.3 Quantitative interpretation of the CO₂-binding affinity by umbrella sampling

To quantitatively evaluate the binding affinity of CO₂ by the IL film, the potential of mean force (PMF) of a single CO₂ molecule diffusing from the IL-gas interface into the substrate-IL interface was quantitatively calculated via the umbrella-sampling method.^27–29 The separation between CO₂ and the nanosheet is selected as a collective variable (CV). Totally 20 points were sampled along with the CV. The harmonic force with a constant of 1000 kJ mol⁻¹·nm⁻² was applied to the CO₂ molecule. Each point was sampled for 10 ns, and the PMF profile was generated using the Weighted Histogram Analysis Method.²⁸

3. Results and discussion

3.1 IL distributions on the substrate

To investigate the role of GRA and C₃N substrates on the IL film, we first characterized the structural features of the IL film on the two substrates. Fig. 2a–d depict the number density profiles of cation [EMIM]⁺ and ion [BF₄]⁻ along the normal direction of GRA and C₃N, respectively. Except for the system with an IL thickness of 0.57 nm, two clear peaks, which represented the first and second binding shells of cations (Fig. 2a and b) and anions (Fig. 2c and d) around the substrate, were detected. The high value indicates a strong attraction of the substrate towards IL molecules. Such phenomenon agrees well with the previous reports about the 2D material and IL interaction.^14,30 Beyond these two peaks, the density profiles of the ions reach a plateau with only slight fluctuation, until they drop to zero when being far away from the nanosheets. Based on the variation of the density profiles of ILs, the IL film is divided into three regions: (1) the substrate-IL interfacial region that covers the first two binding shell of IL, (2) the bulk IL region where IL is evenly distributed, and (3) the IL-gas interfacial region where the IL density starts to diminish. For the three regions, the thicknesses of the two interfacial regions are relatively constant, except that the bulk IL region varies in different systems. In the extreme case of the IL film with a thickness of only 0.57 nm (the black lines in Fig. 2a–d), the substrate-IL and IL-gas interfacial regions merge without the bulk IL region.

Fig. 2 Number density profiles of different compositions. (a and b) Represent the number density profiles of cations; (c and d) represent the number density profiles of anions; (e and f) represent the number density profiles of CO₂ molecules in the C₃N and GRA systems.

3.2 The CO₂ assembly in the IL film

The CO₂ molecules were initially placed outside the IL film. During the simulations, varying number of CO₂ molecules adsorbed and intruded into the IL film. Fig. 2e and f depict the number density profiles of CO₂ along the normal direction of GRA and C₃N. For each system, there are two detectable CO₂ peaks at the substrate-IL interface, which represent the direct attraction of GRA and C₃N towards CO₂. In addition, one dominant peak can be detected for the density profile, which was located at the edges of the IL film, particularly, at the IL-gas interface. When the thickness of the IL film decreases, this peak shifts according to the IL-gas interface. For the thinnest IL film with a thickness of 0.57 nm, the two accumulative regions of CO₂ merge into one. More importantly, CO₂ accumulation inside a 0.57 nm IL film is significantly enhanced than that in other systems, as indicated by the highest peak among the five systems. This reflects synergistic effects from both the substrate-IL and IL-gas interfaces in the thinnest IL film. Moreover, the CO₂ densities in the bulk IL regions are found to be constant, which are considerably lower than those found in the two interfacial regions, independent of the variation of the IL film thickness. Fig. 3a depicts the total number of captured CO₂ with respect to the thickness of IL film. In the studied thickness region, roughly a linear dependency is found for both the GRA and C₃N systems. This is because when the IL film thickness increases, more CO₂ can be stored inside. However, when the CO₂ storage is normalized by the number of IL molecules (storage efficiency of the IL), which is depicted in Fig. 3b, the CO₂ storage efficiency increases quickly with the decrease in IL thickness, which is well consistent with SILM experimental reports.³¹ These results are expected because a high-density layer of CO₂ is formed at the 0.7–1.1 nm from the substrates (Fig. 2e and f) where the IL density becomes considerably low (Fig. 2a–d). From these results, it is also deduced that the bulk IL region has lower CO₂ storage capability than two interfacial regions.

Fig. 3 (a) Total number of captured CO₂ with respect to IL thickness; (b) the CO₂ capture efficiency of the IL films with different thicknesses calculated by dividing the number of captured CO₂ by the total number of IL molecules.

From the above-mentioned results, the role of two substrates is highlighted in the introduction of the two interfacial regions, which are rather attractive to CO₂ than to bulk IL. To quantitatively distinguish the roles of interfacial regions in attracting CO₂, we classified the captured CO₂ according to their locations, in detail: the substrate-IL interfacial region, the IL bulk and the IL-gas interfacial regions. As summarized in Table S1† (GRA systems) and Table S2† (C₃N systems), CO₂ residing inside the two interfacial regions are relatively constant despite the IL film thickness. More importantly, the IL-gas interface plays a more dominating role in the CO₂ storage than does the solid-IL interface. For the bulk IL region, the CO₂ storage efficiency is also considerably lower than that of the IL-gas interface.

In addition, no significant differences were detected between the GRA and C₃N systems. This is because the regulation of the substrate to the IL structure is short-ranged, which is limited to the solid-IL interface, and the contribution of the solid-IL interface in CO₂ storage is rather small. This is consistent with the previous report on using hydrophilic TiO₂ as SILM¹⁴ where the CO₂ storage in the substrate-IL interfacial region was found to be small. Taken together, these results suggest that the capability of SILM is not sensitive to the specific type of supporting materials, which are mainly treated as substrates in the SILMs instead of amplifiers for the CO₂ capture.

3.3 Thin IL film facilitates quick CO₂ diffusion

Except for the gas capture process, the CO₂ mass transfer in the IL film is another determinative factor to evaluate the performance of the medium because a high rate of mass transfer is typically required. In addition, we studied the diffusion character of CO₂ molecules in the IL film. Technically, the mean square displacement (MSD) of CO₂ molecules along the transverse direction of the IL film was calculated. Then, the lateral diffusion coefficient (D) was fitted based on Einstein's relation , where t is the time, and r_i(t) and r_i(0) are the positions of CO₂i at time t and time 0, respectively.

Fig. 4a presents the specific values of D of the dissolved CO₂ in the GRA and C₃N systems with different IL film thicknesses. It is found that the diffusion of the dissolved CO₂ in both systems shows a decreasing trend with the increase in the IL thickness. According to the diffusion-reaction theory,¹³ the bulk IL region was deduced to act as major resistance for the CO₂ diffusion. A thinner IL film can enhance the CO₂ diffusion. Accordingly, it is notable that there is a clear increase in the CO₂ diffusion at the IL thickness of 0.57 nm where the bulk IL region totally vanishes. These results suggest that the introduction of the interfacial IL regions by substrates greatly facilities the CO₂ diffusion. This is consistent with previous experimental reports that CO₂ exhibits faster kinetics in IL with substrates than that in pure IL.^11,32,33

Fig. 4 Self-diffusion coefficients of the dissolved CO₂ molecules (a) and IL ions (b).

As in the IL medium, the diffusion of CO₂ is highly related to the cations and anions. We calculated the lateral diffusion coefficients of the cations and anions in both the GRA and C₃N systems (Fig. 4b) and found a similar vertiginous decrease tendency with the IL film thickness. Shi et al. reported that ions in IL diffuse faster with confinement from nano-materials at room temperature.¹⁵ In our simulations, the confinement rate of IL by nano-sheets decreases with the IL film thickness, thus the ion diffusion is the fastest in the systems with the smallest IL film thickness. Additionally, the diffusion coefficients of ions in the C₃N systems are found comprehensively lower than those in the GRA systems. This is attributed to the stronger surface-IL interactions, which slow down the self-diffusion of ions, as reported by Shi et al.¹⁵ We calculated the interaction energy of one ion binding onto GRA (−48.34 ± 11.24 kJ mol⁻¹ for [EMIM]⁺ and −13.53 ± 3.82 kJ mol⁻¹ for [BF₄]⁻) and C₃N (−49.02 ± 10.55 kJ mol⁻¹ for [EMIM]⁺ and −14.76 ± 3.32 kJ mol⁻¹ for [BF₄]⁻). For both the cation and anion, their interactions with C₃N are slightly higher than those with GRA. This explains the slower ion diffusion in the C₃N systems than that in the GRA systems.

3.4 Quantitative interpreting CO₂-binding affinity and structural mechanisms

To gain a deep understanding of the binding behavior of CO₂ inside the IL film, the binding affinity changes of a CO₂ molecule translocating from the gas region into the IL interior is studied by calculating the potential of mean force (PMF) by the umbrella-sampling method. As summarized in Fig. 5, the PMF profiles illustrate that the global minimum, around −6 kJ mol⁻¹, is located at the IL-gas interface in both GRA and C₃N systems, which is well consistent with the CO₂ distribution shown in Fig. 2c and f. In addition, there are two local minima near the two nanosheets, corresponding to the two density peaks of CO₂ at the solid-IL interface in Fig. 2c and f. These minima, particularly those at the IL-gas interface, reflect the strong attraction of interfacial regions to CO₂ molecules and the accumulation of CO₂. In the bulk IL, the PMF profiles are relatively flat and generally higher than those at the IL-gas interface, corresponding to the overall CO₂ storage capability of IL itself.

Fig. 5 Potential of mean force (PMF) for CO₂ diffusion from the IL-gas interface to the substrate-IL interface. The dashed lines represent the position of the two interfaces.

To probe the atomic mechanism of such patterned CO₂ distribution inside the IL film, we further accessed the microscopic structure of the molecule assembly in the IL film by calculating the orientations of the [EMIM]⁺ cations and CO₂. The [BF₄]⁻ anions were ignored because of their tetrahedron shapes. We defined two tilt angles for [EMIM]⁺ to describe the orientation of the cation: one is the angle between the ethyl chain and the normal direction of the substrate (Fig. 6a); the other is the dihedral angle between the imidazole ring plane and the substrate (Fig. 6e). The orientation of the CO₂ molecule is described by the angle between the linear CO₂ and the normal direction of the substrate (Fig. 6i).

Fig. 6 Definition and distribution of the tilt angles of the ethyl side chain of cations (a–d), imidazolium ring of cations (e–h) and CO₂ (i–l) in three IL films attached to GRA with the IL thickness of 0.57 nm, 2.41 nm and 4.09 nm, respectively.

We used the GRA system as a representative for discussions, and three IL films with thicknesses of 0.57, 2.41 and 4.09 nm were compared, as shown in Fig. 6. The corresponding results for the C₃N system are generally similar, as shown in Fig. S2 in ESI.† In a thin IL film, orientations of the cations are strongly influenced by the nanosurface. The orientations of the IL cations and CO₂ are highly biased at the substrate-IL and IL-gas regions. From the angle distributions, we observed two distinguishable spots for the ethyl chain tilt angle (Fig. 6b), which are labelled as s1 (around 50°) and s2 (90°), respectively. Two indistinguishable spots for the imidazole ring tilt angle (Fig. 6f), which are both labelled s3 (at equivalent 0° and 180°), can also be detected. From this information, it is seen that the cations stack to the substrate with the imidazole ring parallels to the substrate surface. However, for the ethyl chains, they are flexible enough to be alternatively attached to or detached from the surface, as illustrated by the two representative structures embedded in Fig. 6b. In addition, the CO₂ molecules at the substrate-IL interfacial region adopt a parallel binding pose to the nanosurface (as shown in Fig. 6j) because the tilt angle is located at around 90° (Fig. 6j, s4). Such a binding mode maximizes the intimate contacts between CO₂ and the substrates, representing the stronger attraction from the nano-surface than from IL.

It is notable that the three tilt angles are highly conserved in the substrate-IL interfacial region, regardless of the IL thickness. This is also true for the tilt angles of the ethyl chains, which are restrained between 30° and 60° (Fig. 6c and d, s5), indicating that the ethyl chains of [EMIM]⁺ point to the gas phase. This is in line with the finding that more ethyl groups are located at the IL-gas interfaces than at the imidazolium rings (Fig. S3†). Such structure is also consistent with previous experimental reports that, at the IL-gas interface, a longer alkyl chain of cations tends to point to the gas phase.³⁴ As for the CO₂ molecules at the IL-gas interface, a parallel orientation along the surface is also favored although the specific angle spreads in a range of 60–120° (Fig. 6j–l, s6). Such orientation also ensures intimate binding of CO₂ with IL.³⁵ In general, from these analyses, the assembled CO₂ in two interfacial regions are accompanied by the directionally-oriented CO₂ and cations, indicating structural compatibility among them.

From the above-mentioned analyses, the orientations of CO₂ in IL are highly regulated by the two interfacial regions. In addition, the dynamic characteristics of the CO₂ molecules in the two interfacial regions are also different, as revealed by the residence time of CO₂ (Fig. S4†). The CO₂ residence time at the IL-gas interface distributes mainly around 40–50 ps in both GRA systems and C₃N systems, which indicates the stable CO₂ adsorption at the IL-gas interfaces. In contrast, the residence time at the substrate-IL interface is more dispersed than that at the IL-gas interface (Fig. S4†). It is found that the binding at the substrate-IL interface can either last for more than 100 ps (mainly representing CO₂ molecules in the first adsorption layer of the substrate-IL interface) or less than 10 ps, which represent CO₂ molecules that cannot form a significant residence in the IL-gas interface region.

4. Conclusion

Using all atomic molecular dynamics simulations, we studied the structural and energetic characteristics of SILM with CO₂ capture. By coating different number of [EMIM][BF₄] pairs onto the nanosheet, the influences of the IL film thickness are systematically discussed. First, the CO₂ accumulation at the IL-gas interfacial region dominates the CO₂ capture by SILM, suggesting the minor role of specific substrate materials. In general, the CO₂ capturing efficiency of supported IL increases with a decrease in the IL film thickness. A sub-nanometer IL film possesses the highest CO₂ capture efficiency because of the formation of a reduced bulk IL region. The nanomaterial surfaces can effectively accelerate the diffusion of CO₂, which is beneficial for the CO₂ mass transfer. The current study provides deep understanding of the molecular interactions at the interfaces of SILMs, and paves the way for the design and fabrication of optimized supporting nanomaterials of SILM for the CCS usage.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (Grant No. 11874238), the Natural Science Foundation of Shandong Province (Grant No. ZR2018MA034, ZR2020MB074) and Fundamental Research Fund of Shandong University.

References

1. A. E. Creamer and B. Gao, Environ. Sci. Technol., 2016, 50, 7276–7289.
2. M.-O. Schach, R. d. Schneider, H. Schramm and J.-U. Repke, Ind. Eng. Chem. Res., 2010, 49, 2363–2370.
3. L. A. Blanchard, D. Hancu, E. J. Beckman and J. F. Brennecke, Nature, 1999, 399, 28–29.
4. M. Aghaie, N. Rezaei and S. Zendehboudi, Renewable Sustainable Energy Rev., 2018, 96, 502–525.
5. M. J. Earle, J. M. Esperanca, M. A. Gilea, J. N. Lopes, L. P. Rebelo, J. W. Magee, K. R. Seddon and J. A. Widegren, Nature, 2006, 439, 831–834.
6. M. Kosmulski, J. Gustafsson and J. B. Rosenholm, Thermochim. Acta, 2004, 412, 47–53.
7. M. Bui, C. S. Adjiman, A. Bardow, E. J. Anthony, A. Boston, S. Brown, P. S. Fennell, S. Fuss, A. Galindo, L. A. Hackett, J. P. Hallett, H. J. Herzog, G. Jackson, J. Kemper, S. Krevor, G. C. Maitland, M. Matuszewski, I. S. Metcalfe, C. Petit, G. Puxty, J. Reimer, D. M. Reiner, E. S. Rubin, S. A. Scott, N. Shah, B. Smit, J. P. M. Trusler, P. Webley, J. Wilcox and N. Mac Dowell, Energy Environ. Sci., 2018, 11, 1062–1176.
8. X. Liu, G. Zhou, S. Zhang and X. Yao, Fluid Phase Equilib., 2009, 284, 44–49.
9. B. F. Goodrich, J. C. de la Fuente, B. E. Gurkan, D. J. Zadigian, E. A. Price, Y. Huang and J. F. Brennecke, Ind. Eng. Chem. Res., 2011, 50, 111–118.
10. R. D. Noble and D. L. Gin, J. Membr. Sci., 2011, 369, 1–4.
11. X. Wang, N. G. Akhmedov, Y. Duan, D. Luebke and B. Li, J. Mater. Chem. A, 2013, 1, 2978–2982.
12. W. Xie, X. Ji, X. Feng and X. Lu, Ind. Eng. Chem. Res., 2015, 55, 366–372.
13. W. Xie, X. Ji, X. Feng and X. Lu, AIChE J., 2015, 61, 4437–4444.
14. Z. Tang, L. Lu, Z. Dai, W. Xie, L. Shi and X. Lu, Langmuir, 2017, 33, 11658–11669.
15. W. Shi and D. R. Luebke, Langmuir, 2013, 29, 5563–5572.
16. S. Yang, W. Li, C. Ye, G. Wang, H. Tian, C. Zhu, P. He, G. Ding, X. Xie, Y. Liu, Y. Lifshitz, S. T. Lee, Z. Kang and M. Jiang, Adv. Mater., 2017, 29, 1605625.
17. X. Li, T. Guo, L. Zhu, C. Ling, Q. Xue and W. Xing, Chem. Eng. J., 2018, 338, 92–98.
18. G. Qin, Q. Cui, W. Wang, P. Li, A. Du and Q. Sun, Chemphyschem, 2018, 19, 2788–2795.
19. X. Li, Y. Yin, X. Chang, Y. Xiong, L. Zhu, W. Xing and Q. Xue, Chem. Eng. J., 2020, 387, 123403.
20. J. Rafiee, X. Mi, H. Gullapalli, A. V. Thomas, F. Yavari, Y. Shi, P. M. Ajayan and N. A. Koratkar, Nat. Mater., 2012, 11, 217–222.
21. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
22. V. V. Chaban, I. V. Voroshylova and O. N. Kalugin, Phys. Chem. Chem. Phys., 2011, 13, 7910–7920.
23. Y. Deng, F. Wang, Y. Liu, Y. Yang, Y. Qu, M. Zhao, Y. Mu and W. Li, Nanoscale, 2020, 12, 5217–5226.
24. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
25. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
26. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
27. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
28. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1995, 16, 1339–1350.
29. B. Roux, Comput. Phys. Commun., 1995, 91, 275–282.
30. M. Atilhan and S. Aparicio, J. Phys. Chem. C, 2018, 122, 1645–1656.
31. N. Wu, X. Ji, W. Xie, C. Liu, X. Feng and X. Lu, Langmuir, 2017, 33, 11719–11726.
32. P. Fan, X. Qiu, F. U. Shah, Q. Ji and R. An, Phys. Chem. Chem. Phys., 2020, 22, 1097–1106.
33. X. Wang, N. G. Akhmedov, Y. Duan, D. Luebke, D. Hopkinson and B. Li, ACS Appl. Mater. Interfaces, 2013, 5, 8670–8677.
34. S. Rivera-Rubero and S. Baldelli, J. Phys. Chem. B, 2006, 110, 4756–4765.
35. M. E. Perez-Blanco and E. J. Maginn, J. Phys. Chem. B, 2010, 114, 11827–11837.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0na00875c

---