How does the shape and surface energy of pores affect the adsorption of nanoconfined fluids?

Funding information Engineering and Physical Sciences Research Council, Grant/Award Numbers: EP/E016340, EP/J014958, EP/P020194, EP/G036888/1 Abstract We report a systematic molecular simulation study of the behavior of Lennard–Jones fluids inside nanopores of diverse shapes, focusing on the effect that the pore geometry and the local energetic environment have on the adsorption isotherms. Infinitely long pores with polygon (triangle, square, pentagon, hexagon, octagon, decagon, and circle) cross sections are considered. Three different pore sizes commensurate with the molecular diameters along with three different values of fluid–solid energy interactions are chosen to perform Grand Canonical Monte Carlo simulations at a subcritical temperature. Overall, the effect of nanoconfinement on the adsorption of fluids is seen to be a delicate balance between the geometric packing restrictions imposed by the hard cores of the molecules and the surfaces, the excess adsorption induced by the presence (or absence) of energetically favored “hot spots” and the overall ratio of surface/bulk fluid volume present in the pore.


| INTRODUCTION
It is expected that the specific and unique details of the interactions between the molecules of a fluid have a big impact on their adsorption and diffusion under confinement. [1][2][3] However, in addition to the fluid-fluid interactions, the detailed characteristics of the adsorbate material play an equally important role on the ultimate behavior of the confined fluid. 4 Both natural and anthropogenic nanoporous materials have a plethora of pore geometries and sizes and it is precisely the subtleties in the local geometries that has a profound influence on the interfacial and transport properties of confined fluids. Fine-tuning the specifics of porous material geometry, something which is now in the realm of material sciences (cf. Figure 1), provides for avenues to tailor adsorbents to particular applications as has been exemplified in the rectification of ionic currents, 7,8 in the detection of nanoparticles, 9 for capillary condensation/evaporation in materials, 10 and in gas storage, 11 among others.
Theoretical and semiempirical models used to predict the adsorption isotherms of gases require assumptions to be made on the pore geometry and the topology of the existing pore network. The most common workhorse models are based on simplified scenarios, typically either cylindrical [12][13][14][15][16][17] or slit-like [18][19][20][21][22][23] pore geometries. In contrast with these idealized conformations, a detailed analysis of porous materials (e.g., montmorillonite, silica, carbons, clays, etc.) through techniques such as atomic force microscopy, scanning electron microscopy, and transmission electron microscopy (TEM) reveal the widespread presence of steps, angles and wedges. [24][25][26][27][28] In particular, the presence of angular pores (or wedges) in materials can affect the behavior of confined fluids as evidenced by changes in filling transitions, 29,30 the rise of a liquid and bubbles in a tube, 31 unexpected adsorption of liquids, 32 formation of solid phases, 33 and the nucleation kinetics due to orientational order effects. 28 In a similar way, the effect of the roughness of natural surfaces on adsorption has been considered in several molecular-based studies. [34][35][36][37] Prof. Keith E. Gubbins, to whom this manuscript is dedicated, has been a pioneer in the use of molecular simulations to study the behavior of confined fluids. In his early papers, he commented on how "The use of classical thermodynamics to investigate fluids in narrow pores is suspect, because of the highly inhomogeneous nature of the fluid. An understanding of such systems must come from the use of statistical mechanics." 38 Some of his seminal papers in the area (which in spite of being over 30 years old are still pertinent) describe the use of Grand Canonical Monte Carlo (GCMC) methods to study adsorption isotherms of Lennard-Jones fluids in cylindrical pores, 39 and the contrasting effects of pore geometry (between slit-pore and cylindrical models) on the adsorption isotherms. 11 Clearly, both then and now, molecular simulations provide for an unequivocal approach to studying the general traits of ultraconfinement in fluids. More recently, Sarkisov and Monson 40 quantified the effect of pore geometry in the hysteresis phenomena by comparing adsorption and desorption using molecular dynamics (MD) simulations. They considered a slit pore (with and without a closed end), a triangular pore (wedge), and an ink bottle geometry, concluding that different geometries lead to different adsorption isotherms types and different hysteresis loops. Coasne et al. 41 explored the adsorption of Argon on silica surfaces and nanopores using molecular simulations. In terms of pore geometry, their work compared the adsorption inside cylindrical and hexagonal pores, showing that both isotherms are similar in shape. However, the condensation and evaporation pressures for the cylindrical pore were smaller than the pressures for the hexagonal pore.
The work of Sedghi et al. 42 presented the study of wetting layer formation in shale rocks through MD simulations, where the pores structures are modeled as a circle with both smooth and rough wall, a square and a right triangle. From their results, it stems that the presence of angular pores can lead to metastable layer formations when water-flooding simulations were performed. Wallacher et al. 43 presented a study of the capillary condensation in linear mesopores of different shape. Surprisingly, they concluded that the individual local geometries such as the pore mouth, a constriction or adding a blind end do not affect the shape of the adsorption isotherm significantly. Sun et al. 44 employed elliptical cross sections, which gradually transition from circles to slit-like cross sections by manipulating the aspect ratio, finding that both geometry and surface energetics affect the transport of water through nanopores.
As seen above, the presence of irregular pores seems to play a role in the behavior of confined fluids, although we are not aware of any systematic study elucidating the structure-adsorption relationship. In spite of this, most researchers and modelers exclusively employ the cylindrical pore (or a slit pore) as an approximation of the more realistic randomshaped pores found in natural amorphous porous materials. There is an implicit presumption that the general shape and overall characteristics of the adsorption isotherm will not be affected by the shape of the pore, as long as the pore surface and volume are retained. We here question if this is a valid assumption. The core of this manuscript reports GCMC simulations of model fluids in well-defined pores in order to study the effect that the pore geometry has on the adsorption isotherms. The chosen pores have regular structures with different cross-sectional shapes but identical pore cross-sectional areas. The proposed pores have triangle, square, pentagon, hexagon, octagon, decagon, and circular crosssectional area of different sizes (probing different surface/bulk fluid ratios) and surface energies to quantify the effect that these energetical and morphological heterogeneities have on adsorbed fluid behavior.

| METHODOLOGY
The model for both the fluid-fluid and the fluid-solid interactions is based on the Mie potential, 45 which is a generalized form of the Lennard-Jones (LJ) potential: where σ is a length scale that corresponds loosely with an effective segment diameter, r is the distance between the centers of the particles, λ r and λ a are the repulsive and attractive exponents which the molecular diameter σ = 0.3 nm. 46 The values for the repulsive, λ r , and attractive, λ a , exponents are 12 and 6, respectively, corresponding to a single-sphere Lennard-Jones molecule. Lennard-Jones, and Mie fluids in general, show a conformal behavior. 47 In practice, this means that the thermophysical properties may be scaled (reduced) appropriately, for example, the temperature, T, may be represented as T* = k b T/ε, where k b is Boltzmann's constant, and the number density, ρ, may be scaled as ρ* = ρσ 3 . Other properties can be similarly scaled.
The behavior of a fluid, expressed in these reduced units, is accordingly universal (to within the validity of the corresponding states principle 48 ). In this sense, the results in this paper can be taken as representative of a whole family of model isotropic fluids.
The model pore geometries considered in this work are displayed in Figure 2, showing the corresponding triangle, square, pentagon, hexagon, octagon, decagon and circle cross sections. All the different pore geometries used have the same cross-sectional area in order to preserve equivalent pore volumes. The solid porous materials are built using spherical segments with σ = 0.2 nm, placed in a square lattice at a center-to-center distance of 0.1 nm in order to obtain a close-tosmooth surface in a multilayer formation. This guarantees that the density of the solid material is identical for all the pore structures.
The interaction between the fluid particles and the solid walls follows the Mie potential (Equation (1)), where the repulsive, λr, and attractive, λ a , exponents are 11 and 5, respectively. 35,49 For the fluid-solid interaction, three different values of potential depth, ε w , are employed, labeled weak, medium, and strong (Table 1). Furthermore, three different pore areas are employed for the simulations, as shown in Figure 3; 0.78, 3.14, and 12.56 nm 2 (from left to right). As a reference, the smaller circular pores have a diameter of 1 nm, while the larger ones, are of 2 and 4 nm in diameter. The length of each pore in the axial (z) direction was 4 nm; however, the system has periodicity in the z direction, essentially making the pores infinite in length. The number of beads used for the solid material depends on the pore geometry and the pore size. For the small, medium, and big pores, the numbers of beads are approximately 60,000; 85,000; and 150,000, respectively.
The number of particles in the fluid phase will depend on the chemical potential and the pore size, but the upper limit for the bigger pores at  Simulations are performed in the GCMC ensemble 50,51 to obtain the number of particles adsorbed in the system. In this ensemble, the number of particles is allowed to vary at a fixed temperature and chemical potential. Each type of move (translation, deletion, and insertion) is chosen with equal probability. Acceptance ratio for translations is targeted at 30%. Each system is allowed to equilibrate discarding the first 10 7 configurations, and the samples for the number of particles adsorbed inside the pore were taken over the following 10 7 configurations. For all the cases analyzed in this work, the temperature T is fixed to 100 K (corresponding to T* = 1, below the bulk critical point). The cut-off radius for both fluid-fluid and fluidsolid intermolecular potentials is taken to be 1.5 nm (ffi5σ).
It is customary to represent experimental data for adsorption isotherms as the adsorbed amount as a function of the pressure of the bulk gas in equilibrium. In GCMC simulations, it is the chemical potential that is specified instead of the equilibrium pressure. In order to relate the chemical potential to the external pressure P, the same chemical potential used for adsorption is used for simulations of a bulk phase at the same temperature and at an arbitrary volume. From the output of the latter, the bulk pressure is determined by calculating the virial. 51

| RESULTS
This section showcases the results obtained from GCMC simulations for the adsorption of a single-bead LJ fluid on pores of different geometries with constant cross-sectional area at the subcritical temperature of 100 K.

| Large pores
The large pores use as a reference a cylindrical pore with diameter   Table 1). For the weak wall interaction (Figure 5a) all isotherms have a Type III shape. 52 The adsorption into the triangular pore shows the highest uptake in the gas region; however, after the first-order transition, in the liquid phase, the adsorbed density is the lowest for all the geometries. The adsorption in the pore with square shape is lower than the triangular geometry in the gas region, but higher in the liquid region. The isotherms for the pentagon, hexagon, octagon, and decagon continue the trend gradually converging to curves that seem indistinguishable from those of the cylindrical pore. On the other hand, at higher wall energies (Figure 5b,c), a Type IV isotherm is observed, where the triangular-shaped pore sees a higher uptake across the whole pressure range. As the number of vertices increase, the curves tend to converge to a single result coinciding with the cylindrical pore. This behavior can be rationalized with help from the energy maps in Figure 4,

| Medium pores
The medium size pores are based on a cylindrical pore with 2.0 nm diameter as reference, leading to a cross-sectional area of 3.14 nm 2 .
The energy density distribution in these pores (Figure 7) for the strong fluid-solid interaction resembles that seen in the larger pores ( Figure 4): the triangular geometry has the highest energy anisotropy with the highest energy in the corners, while the distribution flattens out in the limiting case of the circular pore. While in the larger pores ( Figure 4), the center of the pore experiences almost no attractive energy from the wall (evidenced by the white regions where ε w /k b ffi 0), in the energy maps for the medium size pores, this central (white) area is seen to be significantly smaller (almost inexistent), which means that the effect of the wall potential extends over the whole volume inside the pore.  Figure 8c, we observe that the adsorption in the triangular pore is the highest, the uptake keeps descending as the number of sides is increased and the lowest limit is the cylindrical pore, following the same trend as the more energetic adsorption results in the big size pores (Figure 5c). We presume this is related to a favorable layering order induced by the underlying expected solid phase, which would have its highest packing in an FCC or BCC lattice.
These solid structures are incompatible with the increasingly curved surfaces, but are enhanced by the presence of flat surfaces, which are more prominent the less angles the pore has.  Figure 8c, where in the saturation region, the adsorption hierarchy for these systems was: triangle > hexagon > circle. Isolating the first adsorbed layer of fluid (adjacent to the solid surface), allows us to calculate the hexagonal in-plane order parameter, ψ, that quantifies the structure of the layer. This order parameter is defined as where θ j is the angle formed between the vector joining the centers of two adjacent beads and a randomly taken director in the plane and the sum is taken over all N beads. The order parameter can take values from 0 to 1, from the unordered to the fully ordered state, respectively. Figure 10 showcases  Figure 10, where we can observe that the order in the triangular pore is the highest, with around 71% of order, followed closely by the hexagonal pore with 69% of order. For the circular pore, the F I G U R E 7 Example energy maps, u(x,y), of the fluid-solid interaction for pore geometries used in this work. The pore size corresponds to the largest pore with area 3.14 nm 2 , and the potential depth is chosen as the strong fluid-solid interaction, ε w /k b = 12.0 K for each fluid particle interacting with a particle from the wall [Color figure can be viewed at wileyonlinelibrary.com] In the case of the weak interaction, Figure 8a, the adsorption in the triangular pore is the highest for the gas region (low pressure), while the adsorption in the cylindrical pore is the lowest.  Figure 8c). If the strength of the wall interaction is low, the wall will become increasingly nonwetting and the fluid molecules will prefer to be removed from there. At these conditions, the uptake inside the pore becomes a geometric problem, where the available free volume is lower for the triangular pore and higher for the decagonal cross section. For the medium interaction, Figure 8b shows that all the isotherms converge to a similar value in the liquid region, as for this condition both the filling imposed by the geometry and the wall order induced by the enhanced energy are equally important.

| Small pores
The small size pores use a cylindrical pore with 1 nm diameter as ref- erence, leading to a cross-sectional area of 0.78 nm 2 . For the small pores, only the weak fluid-solid interaction is presented, as for the medium and strong interactions, a solid phase is induced due to the ultra-confinement effect, which increases the effective strength of the fluid-solid interactions. Figure 11 displays the energy maps for the different geometries of small pores, and it is possible to observe the high energy inside the pore, a consequence of the superposition of the attractions between neighboring walls causing the whole internal space to be strongly affected by the attraction of the solid walls. Recently, Song et al. 53 presented simulation results relating to the effect of pore shape on the adsorption of methane in graphitic media.
Their results, which span only the very confined spaces up to 1.15 nm, are consistent with our findings for the smaller geometries.

| Patches
In order to illustrate the effect of high-energy vertices in isolation to the effects of having regions of flat surfaces, we employ a cylindrical pore and add three patches in a triangular configuration. This analysis F I G U R E 1 1 Example energy maps, u(x,y), of the fluid-solid interaction for pore geometries used in this work. The pore size corresponds to the largest pore with area 0.78 nm 2 , and the potential depth is chosen as the strong fluid-solid interaction, ε w /k b = 12.0 K for each fluid particle interacting with a particle from the wall [Color figure can be viewed at wileyonlinelibrary.com]  (Table 1), while the patches (blue particles in Figure 13b) have an even stronger interaction of ε w,patch /k b = 200 K for each particle in the patches. Figure 13a shows the energy map for the system used, where it is possible to observe that the energy is higher in the regions around the patches. The rest of the wall has a homogeneous energy, which corresponds to the same energy interaction that was used for the results presented in Figure 5c.
The adsorption isotherm obtained from this system is presented in Figure 13c, where three isotherms are compared, the cylinder with patches (blue crosses) with the adsorption in the triangular pore (black filled triangles), the regular cylindrical pore (red filled circles) (cf. Figure 5c). From this figure, we can observe that the isotherm obtained from the system with the patches lies in between the isotherms from the adsorption in the triangular pore and the one in the regular cylindrical pore. If we compare this with the results for the strong interaction for large pores, Figure 5c, the isotherm from the patchy system behaves similar to the adsorption into the square pore (green symbols), but the uptake is higher for low pressures, and as the pressure is increased, the isotherm converges to that of the square pore. Figure 13c clearly shows that the role of the energetical heterogeneities is only partially responsible for the change in adsorption behavior. The presence of edges, corners, and flat wall regions has an equivalently profound effect on the overall shape of the isotherm.

| CONCLUSIONS
This work shows the effect that different pore geometries, pore sizes, surface/bulk fluid volume present in the pore. This has an important implication in the modeling of adsorbents with very tight pores, for example, activated carbons and shale oils. In this case, and seen with this perspective, density functional kernels based on slit pores 54 and simulations based on collections of planar subunits 55 will be better posed to represent the acute angles and pore geometries seen in carbonaceous materials. 56 In the limit of ultraconfined fluids, the realm of pores of 1 nm of effective diameter (with only two or three molecular diameters in span), the wall energetics dominate the adsorption and a continuous smooth uptake is observed, leading to Type V isotherms. For these pores of 1 nm or less, the particular geometry is very relevant, as some of the pore volume becomes inaccessible by the acute angles formed by the walls. In these cases, it would be unrealistic to describe the adsorption complex geometries with the simple kernels discussed above.

ACKNOWLEDGMENTS
Computations were performed employing the resources of the Impe-