Alternative graph-4-yne stacking fashion: toward selective CO₂ capture

Abstract

In the present study, we investigate the stacking configurations of bilayer and trilayer graphtetrayne, a promising carbon-rich two-dimensional material derived from graphene, by integrating literature-reported structures with a synergistic computational approach that combines classical molecular dynamics and ab initio calculations. The classical molecular dynamics simulations are performed adopting the improved Lennard-Jones (ILJ) formulation of the potential, a crucial choice to easily explore the configurational space identifying low-energy arrangements followed by a more precise investigation performed through ab initio calculations for accurate energetic evaluation and refinement of structures. Through this combined methodology, we identified a stacking mode present in both isolated and stacked bilayer structures, which, although previously reported in aqueous environments, is here characterized in the gas phase and found to be significantly more stable than the configurations proposed in the literature. Furthermore, molecular dynamics simulations of an equimolar CO₂/N₂ mixture revealed a marked preference for carbon dioxide uptake by the bilayer structures. These results indicate that the bilayer structure preferentially accommodates CO₂ due to favorable quadrupole–π interactions, linear geometry, and accessible pore morphology, providing initial quantitative evidence of its potential for selective CO₂ capture. The present findings not only contribute to the fundamental understanding of graphtetrayne bilayers, but also highlight the potential of carbon-based nanomaterials for the design of efficient membranes for environmental and energy-related applications.

---

1 Introduction

The increasing level of atmospheric greenhouse gases (GHGs) is considered one of the major causes of global warming and climate change. Addressing this issue requires innovative technologies that can efficiently capture, separate, and store CO₂ – one of the most ubiquitous and abundant GHGs – from industrial emissions in the atmosphere. A broad array of materials technologies have been implemented for this purpose, including zeolites, metal–organic frameworks (MOFs) and covalent organic polymers (COFs/COPs).^1–4 Among them, membrane-based gas separation has emerged as a promising technique,^5–8 due to its modularity, low energy consumption, and potential for continuous operation. In the framework of such so called covalent organic nanosheets (CONs),⁹ structures based on carbon atoms have proven to be potential candidates as separation and storage materials. These systems are thin two-dimensional lattices with size-tunable pores and structural and electronic properties that are suitable for selective gas transport and adsorption. A good starting point for designing such materials is graphene. Although monolayer graphene membranes appear to be impermeable to most gases,¹⁰ the introduction of pores and/or defects on the surface creates pathways for selective gas permeation and adsorption.^11–13 For example, nitrogen-functionalization of graphene layers significantly enhances the CO₂/N₂ capture selectivity, due to the increased interaction strength between CO₂ molecules and the nitrogen dopant.¹⁴ In addition, nitrogen-doped graphene has attracted considerable interest in electrocatalysis and biosensing due to its unique electronic structure.^15–17 Nitrogen atoms can be placed in a graphitic, pyridinic, or pyrrolic configuration, each contributing differently to the material's chemical reactivity and conductivity.¹⁸ These recent results indicate a leading role for carbon material design and functionalization in the development of new technologies for gas storage and separation. Among graphene derivatives, carbon nanotubes (CNTs) have received significant attention for their potential to adsorb small gas molecules such as H₂, H₂O, CH₄, and CO₂, due to their exceptional surface area, tunable electronic properties, and structural versatility.^19–21

A new and promising class of post-graphene CONs is represented by graph-N-ynes,²² 2D carbon layers consisting of hexagonal rings linked by chains of N acetylenic units, resulting in natural porosity and flexibility, since the carbon chains define triangle-like pores of sub-nanometer dimensions in the layers, which are uniformly distributed and can be adjusted in size, by changing the number of acetylenic bonds in the structure.

Graphynes consist of benzene rings interconnected by a single (N = 1) acetylenic linkage, forming a planar network of sp- and sp²-hybridized carbon atoms. Its various predicted allotropes (including α-, β-, and γ-) exhibit tunable electronic features, including direction-dependent charge transport and the presence of Dirac cones in specific configurations.²³ These unique anisotropic properties, along with mechanical stability and π-conjugation pathways, make graphynes suitable platforms for applications requiring tailored electronic behavior, such as field-effect transistors, sensors, and nanoelectronic components.⁹ In addition to its CO₂ adsorption capabilities, γ-graphyne has been shown to intercalate sodium ions within its multilayer structure,²⁴ forming stable compounds with a high theoretical capacity. This highlights its potential not only for gas separation but also as a versatile material for energy storage applications. Graphdiyne (GDY, N = 2), characterized by diacetylene bridges between phenyl units, has been successfully synthesized and experimentally validated as a highly conjugated two-dimensional semiconductor.²⁵ Its hybrid sp/sp² carbon lattice introduces natural porosity and a moderate, direct bandgap, making it suitable for use in optoelectronics, photocatalysis, and membrane-based separation technologies.²⁵ Moreover, GDY exhibits excellent charge mobility and chemical tunability, positioning it as a strong candidate for organic electronic devices and energy storage systems.^26–28 The material's structural versatility and compatibility with bottom-up synthesis approaches further support its integration into scalable, sustainable device architectures.⁹

Among all variants identified, graphtriyne (N = 3), in particular, has triangular pores of sub-nanometer dimensions whose diameter efficiently allows for selective diffusion of gas molecules.^29,30 The possibility of forming multilayer structures represents a key aspect in the design of membrane-based technologies. In detail, molecular dynamics (MD) simulations show how membranes consisting of graphtriyne layers³¹ exhibit high permeability and selectivity for CO₂. The stacking actually enables synergistic effects in gas capture, depending on the number of layers, a parameter which is likely to be crucial for practical applications.^21,32 In addition, first-principle studies³³ have shown that graphtriyne exhibits highly favorable physisorption properties for hydrogen due to its large pore size and low surface density, further emphasizing the versatility of graphynes as 2D materials tunable for gas adsorption.

Graphtetrayne (N = 4), a less studied variant, is obtained by introducing an additional acetylenic bond, leading to larger pores and a more flexible framework. Because of the large pores, the possible applications of graphtetrayne for gas separation and storage depend on the stacking configurations of multilayer systems, which arguably determine the actual pore size. In this respect, Pan et al.³⁴ reported the first successful direct synthesis of crystalline graphtetrayne, using a bottom-up synthetic strategy that leads to the formation of a highly ordered crystalline material. Detailed structural analysis revealed that the graphtetrayne layers exhibit a well-defined stacking arrangement. Specifically, the material adopts an AB-stacking configuration—referred to as 9-fold stacking in the nomenclature of Pan et al.³⁴—which differs from the more commonly observed AA stacking (also reported as 1-fold stacking) in related graphyne structures. The study not only introduces a new allotrope to the family of graphyne materials but also highlights the importance of stacking order in modulating their physical behavior, with potential applications in nanoelectronics, membrane technology, and energy storage.

Experimental and computational investigations on graphtetrayne membranes are currently limited, but their theoretical structure suggests that they could exhibit even higher CO₂ permeability, provided mechanical stability is retained. Future studies will be crucial in determining its viability as a next-generation membrane material. Various experimental approaches have been proposed for the design of carbon-based nanostructures,^35,36 but fully exploring, characterizing, and assessing the performance of all potential configurations remains a significant experimental challenge. Despite the advances in the synthesis of two-dimensional carbon materials, the experimental realization of graphyne-based systems remains largely confined to monolayer or few-layer structures. For graphene as well, after all, the controlled fabrication of multilayer architectures with a well-defined stacking sequence along the out-of-plane direction is highly challenging, while for more complex graphyne derivatives, such control—especially for thicknesses exceeding a few layers—is currently arguably beyond experimental reach.^4,37–39 Computational techniques, particularly MD simulations combined with targeted quantum chemical calculations, can offer a powerful alternative by providing insights into the structure and the reactive properties of the diverse systems. An appropriate intermolecular force field—which describes the interplay of the most fundamental interaction components in governing the stacking of layers and the molecular adsorption—besides ensuring both the accuracy and the reliability of the simulations, is essential as a guide in the search for stationary points in the configuration space (for both stacking and adsorption), avoiding often long and unfeasible ab initio calculations.⁴⁰ This approach can effectively support the experimental work, providing predictive insights that may guide future experimental efforts as suitable fabrication strategies become available.

In the present work, in line with the objectives outlined above, we carried out a theoretical investigation of the bilayer and trilayer stacking configurations of graphtetrayne, using a semi-empirical ILJ formulation of the potential to perform MD simulations.⁴¹ The results revealed a new stacking structure, not found in previous studies, which appears to be more stable than the AB mode proposed by Pan et al.,³⁴ as also confirmed by accurate electronic structure calculations. The same force-field-based approach (FF-based) has also been used to characterize the trilayer structure. A preliminary analysis performed on the two structures revealed comparable performances in terms of selective adsorption of CO₂.

2 Theoretical methods

2.1 Molecular dynamic simulations

To verify the stability of the bilayer structure, molecular dynamics simulations were performed using the DL_POLY software⁴² in the canonical (NVT) ensemble with periodic boundary conditions (PBC) set in all directions. A simulation box of 92 × 84 × 70 ų was chosen. A Nose–Hoover thermostat with a relaxation constant of 0.5 ps was employed to keep the system temperature T at 10 K, while the pressure was set at 1 atm. The cutoff distance for the non-electrostatic and electrostatic components was set to 15 Å, while the Ewald method⁴³ was applied for the calculation of electrostatic contributions.

In detail, the total potential energy includes both van der Waals and Coulombic contributions. For the pristine graphtetrayne bilayer, carbon atoms were assigned zero partial charges, and thus electrostatic interactions do not contribute. When gas molecules, such as CO₂ and N₂, are included (see Section 3.2 for details), partial charges are assigned according to ref. 31, and Coulombic interactions are evaluated via the Ewald method. During the simulations, the graphtetrayne membranes, composed of two layers of 480 atoms each (960 atoms in total), have been considered as rigid structures. A schematic representation of the initial configuration of the system is shown in Fig. 1.

Fig. 1 Schematic representation of the initial configuration of the bilayer system (which corresponds to the AB-stacking configuration, also referred to as 9-fold stacking in the nomenclature of Pan et al.³⁴) [For the sake of visualization the two graphtetrayne layers are in different colors, green and brown, respectively].

Each simulation was performed for 5 ns after a properly long equilibration period (equal to 1 ns) with a fixed time step of 1 fs.

To model the intermolecular interaction potential, a force field based on a formulation developed in our laboratory within a phenomenological framework was employed, namely the Improved Lennard-Jones (ILJ) potential.^41,44–46 This formulation of the interaction has been demonstrated to be highly versatile and applicable to systems involving carbon nanostructures as well as atmospheric components such as CO₂, N₂, H₂O, CO, CH₄, etc… (see ref. 47–54 and references therein). Within this approach, the total interaction energy is expressed as a sum of pairwise additive contributions. Each term explicitly accounts for the interplay between attractive interactions arising from induced dipole–induced dipole forces and short-range repulsive effects associated with the finite size of the interacting species. Accordingly, the intermolecular interaction potential can be written as:

where ε and r^m describe the depth of the well and its position, respectively, associated with each interaction pair, while m assumes the value of 6 in the case of neutral–neutral systems. The values are estimated on the basis of the polarizability, evaluated considering each interaction center at a certain distance r. The term n(r) shows a dependence on the distance r, defined as:where the parameter β represents the hardness of the two interacting fragments and can be properly modulated. The following ILJ parameters have been used to describe C–C interactions: ε = 2.84 meV, r_m = 4.07 Å and β = 7. A similar computational strategy has been used in different cases to understand the role of non covalent interactions and how long-range complexes can drive chemical reactivity,^55,56 in contrast to approaches limited to electronic structure calculations.^57,58

2.2 Ab initio calculations

Ab initio calculations have been performed in the framework of density functional theory (DFT) as implemented in the pseudoatomic orbital (PAO)-based SIESTA code.^59,60 The generalized gradient approximation (GGA) for the exchange and correlation functional parametrized by Perdew, Burke, and Ernzerhof (PBE)⁶¹ was exploited along with norm-conserving pseudopotentials of the Troullier–Martins-type⁶² for the description of the core electrons. In detail, a non-relativistic pseudopotential with 4 electrons in the valence (2s, 2p) has been used to describe C atoms associated with a triple-ζ plus polarization (TZP) basis set. The structures were fully relaxed (cell and ions) until the residual forces were less than 0.02 eV Å⁻¹. In the optimization of the structures, an 8 × 8 × 12 Γ-centered k-point sampling of the Brillouin Zone (BZ) was used. The real-space grid mesh cut-off parameter used in the calculations is 200 Ry. In order to explicitly assess the role or dispersion in stabilizing the layered systems, additional calculations with inclusion of such a contribution (D3) were performed, see Section3.1.

3 Results

3.1 Stacking configurations

We adopted the 9-fold stacking structure reported by Pan et al.³⁴ as the initial configuration for our molecular dynamics simulations, which corresponds to the AB-stacking arrangement of the bilayer system. Both layers were kept rigid during the calculations to properly analyze the effect of the non-covalent interactions in the coupling of the two layers. Upon structural optimization, the bilayer system relaxed into a more energetically favorable configuration, in which the final stacking geometry deviates from the initial 9-fold (AB) one and converges toward an intermediate arrangement between AB (9-fold) and AA (1-fold) configurations. More specifically, the two layers exhibit a slight lateral shift with respect to each other, resulting in a misalignment that is neither fully eclipsed (AA) nor fully staggered (AB). MD calculations performed at 10 K revealed that this intermediate stacking arrangement (hereafter also int) corresponds to a local energy minimum, suggesting enhanced stability compared to the original configuration. Interestingly, similar stacking behavior has also been reported in the liquid phase. In particular, Ying et al.⁶³ investigated the influence of stacking on proton penetration in two-dimensional graphtetrayne bilayers under aqueous conditions. Although their study focused on a solvated bilayer system, the current work reveals an analogous stacking motif under gas-phase isolated conditions, suggesting a potentially intrinsic structural preference for graphtetrayne layers, regardless of the surrounding environment. In Fig. 2, the configurational (potential) energy is reported as a function of the simulation time.

Fig. 2 Configurational (potential) energy as a function of the simulation time obtained through MD calculations for the bilayer system.

As can be seen, the decrease in the energy corresponds to the identification of a stable structure. The simulation parameters are enough to keep the fluctuations of the total energy below 10⁻⁴ kcal mol⁻¹. The lateral shift between the two layers observed in the optimized structure not only alters the stacking geometry but also significantly affects the pore dimensions. In the AA stacking configuration, the pores exhibit a diameter of approximately 7.9 Å, which is relatively large and can reduce interactions with small molecules, such as CO₂ due to reduced confinement effects. After optimization, the shift between layers leads to a contraction of the pore size to approximately 6.5 Å. This reduction enhances the potential for van der Waals interactions and confinement, making the modified structure more suitable for the adsorption of small molecules. Notably, the optimized pore size is larger than the kinetic diameters of CO₂ (3.3 Å) and N₂ (3.6 Å), suggesting that both species can diffuse through the pores, but with confinement effects more pronounced for CO₂. This size relationship highlights the potential for selective adsorption arising from differences in molecular interactions within the confined space. In Fig. 3 the structure of the bilayer is reported, together with the diameter of the pores and the interlayer distance.

Fig. 3 Main structure of the new stacking mode (int) for the bilayer configuration (the yellow circle represents the dimension of the effective pore in the new configuration).

As for the separation of the layers, the MD results show that the system tends to organize itself in a structure characterized by an interlayer distance of 3.7 Å, in accordance with the ILJ formulation of the interaction, where the r_m appears to be equal to 4.07 Å. As already demonstrated in the case of graphtriyne,^30,31 the interlayer region plays a central role in the adsorption process, giving the structure important properties in terms of permeation and selectivity.

To assess the reliability of the ILJ and molecular dynamics predictions, the results were subsequently examined within an ab initio framework through a systematic computational campaign based on density functional theory (DFT) calculations (see Section 2.2). Periodic boundary conditions were applied along all three spatial directions, with the bilayer stacked along the z axis. Within this framework, three distinct structural models were considered and compared (see Fig. 4), corresponding to different stacking arrangements, namely the AA configuration, the stack-9 structure³⁴ and the int configuration predicted by the initial MD analysis.

Fig. 4 Top and lateral view of each of the three optimized structures of the configurations discussed in the text. (a) AA, (b) stack-9, and (c) int stacking. [For the sake of visualization the two graphtetrayne layers are in different colors, green and brown, respectively].

DFT results support the previously observed trend, with the latter fashion which is thermodynamically ∼3 kcal mol⁻¹ (∼0.13 eV) more stable than the stack-9 one. The least stable configuration is, as expected, the eclipsed AA one. Additional DFT calculation with inclusion of dispersion (D3) shows an enhanced energy difference up to 6.9 kcal mol⁻¹ (∼0.3 eV). In Table 1, we report the main energetics and geometric parameters for the three configurations.

Table 1 Main thermodynamic (ΔE in eV) and structural (interlayer distance IL d, in Å) parameters for the three configurations, with energies computed both including and excluding dispersion corrections

	AA	stack-9	int
ΔE (DFT) (eV)	+0.15	+0.13	—
ΔE (DFT-D3) (eV)	+0.29	+0.30	—
IL d (Å)	3.839	3.658	3.724

---

Interestingly enough, one may notice the perfect agreement between the interlayer distance value derived from MD and that obtained at the DFT level of theory (∼3.7 Å).

The energetic stability of the newly identified bilayer stacking configuration was further assessed by comparing classical MD simulations and density functional theory (DFT) results, with the aim of validating the MD predictions against a more accurate electronic-structure description. Within the PBE approximation, DFT calculations predict an energy difference of 0.13 eV between the AB stacking and the new configuration, corresponding to an energy density of approximately 0.69 meV Å⁻² when normalized by the unit-cell area. Upon inclusion of dispersion corrections (D3), the energy difference increases to 0.30 eV, i.e., about 1.59 meV Å⁻², highlighting the crucial role of van der Waals interactions in determining the interlayer energetics. Classical MD simulations reproduce the same energetic ordering, with the new stacking being more stable than the AB configuration, although a significantly smaller energy density of about 0.19 meV Å⁻² is obtained. The larger discrepancy observed with respect to dispersion-corrected DFT reflects the different levels of theory and the intrinsic differences between the two approaches. In particular, the DFT-D3 method provides a more complete description of long-range dispersion interactions within an electronic-structure framework, while the MD simulations are based on an ILJ potential, which represents an effective and computationally efficient description of interlayer interactions. Within this context, the ILJ potential is primarily designed to reproduce structural and energetic trends across extended systems rather than absolute binding energies. Despite the quantitative differences, both approaches consistently predict the same energetic ordering, with the new stacking being more stable than the AB configuration. This agreement confirms the robustness of the ILJ-based MD results in capturing the relative stability of the different stacking configurations, while the DFT-D3 results provide a more accurate reference for the magnitude of the stacking energy. In addition, a quantitative comparison between the MD-derived bilayer structure and the corresponding ab initio-optimized geometry reveals excellent agreement in terms of both interlayer separation and lateral displacement. This consistency indicates that the adopted force field reliably reproduces the essential features of the underlying interlayer interaction potential of graphtriyne with sufficient accuracy to correctly describe the stacking configurations, despite the presence of quantitative differences in the absolute stacking energies. Overall, this agreement provides a practical validation of the classical approach and, in particular, demonstrates its suitability for describing stacking modes in extended graphyne frameworks. Crucially, this validation enables the use of MD-based sampling to investigate more complex assemblies, such as the trilayer system, for which full ab initio optimization becomes computationally prohibitive.

Guided by the bilayer benchmark, we subsequently determined the equilibrium configuration of the trilayer system, which provides a realistic structural model for assessing the adsorption and confinement properties of small penetrant species. Notably, multilayer graphyne-like membranes have been shown to retain high performance in CO₂ and light-gas adsorption and transport,³¹ thereby further supporting the relevance of such extended architectures (Fig. 5).

Fig. 5 Main structure of the new stacking mode for the trilayer configuration (the yellow circle represents the dimension of the effective pore in the new configuration).

The optimized trilayer structure displays a stacking motif that is fully consistent with the one identified for the bilayer. In particular, the three sheets spontaneously adopt a registry pattern that minimizes the interlayer configurational energy, leading to a stable sequence of relative shifts that preserves the characteristic triangular pore arrangement. The resulting structure is remarkably ordered: the lateral displacement between adjacent layers converges toward the same intermediate stacking mode established for the bilayer, indicating that the interlayer potential-energy landscape governing the two-layer system extends naturally to multilayer assemblies. As a direct consequence of this registry, the effective aperture of the pores contracts, with the largest inscribed sphere decreasing from the bilayer value to approximately 5.5 Å, while the average interlayer separation remains close to 3.7 Å. Such a reduction in pore dimension is expected to influence the confinement and adsorption behavior of small penetrants, making trilayer architectures potentially relevant for selective gas uptake and transport, in line with previous observations for multilayer graphyne-like membranes.

3.2 Gas adsorption properties

In order to preliminarily assess the capability of the newly designed bilayer structure for selective atmospheric CO₂ adsorption, MD simulations have been performed. In detail, classical MD simulations have been performed with the same DL-POLY software, considering the canonical ensemble (NVT) with periodic boundary conditions (PBC) in all directions. A Nose–Hoover thermostat with a relaxation constant of 0.5 ps was used to keep the temperature of the system at 300 K. During the simulations, the graphtetrayne membrane structure was kept frozen, and the gas molecules were treated as rigid bodies. Each simulation was performed for 5 ns after a sufficiently long equilibration period with a fixed time step of 2 fs. A mixture gas of 50 CO₂ and 50 N₂ molecules has been considered and the value of the gas pressure for the NVT simulations at 300 K was obtained from the Peng–Robinson equation of state, corresponding to 1.26 atm. The temperatures in the simulations had relative fluctuations of around 5%. The force-field description adopted for gas–gas and gas–membrane interactions follows the same methodological framework reported in ref. 31 and 63. In particular, electrostatic interactions were modeled using a distributed charge scheme, employing three interaction centers for N₂ and five interaction sites for CO₂, for a better reproduction of the molecular quadrupole moment. Non-electrostatic interactions were described using the Improved Lennard-Jones (ILJ) potential, with the corresponding interaction parameters reported in Table 2.

Table 2 List of the interaction parameters describing the ILJ interaction used in the present study (adapted from ref. 31 and 63). G4 stands for graphtetrayne

System	Interacting pairs	ε (meV)	r^m (Å)	β
N2–G4	NN2–C	3.97	3.83	7.50
CO2–G4	CCO2–C	3.55	3.56	6.75
	OCO2–C	5.01	3.67	6.75
N2–N2	NN2–NN2	3.52	3.77	9.00
CO2–N2	CCO2–NN2	3.40	3.55	9.0
	OCO2–NN2	4.50	3.70	9.0
CO2–CO2	CCO2–CCO2	3.42	3.32	9.0
	OCO2–CCO2	5.52	3.36	8.5
	OCO2–OCO2	5.33	3.65	9.0

---

By systematic analysis of the simulation results, the structural and adsorption properties of the CO₂/N₂ mixture were thoroughly characterized in the vicinity of the graphtetrayne membrane. In particular, radial distribution functions were computed to characterize the local organization of the gas molecules and their interactions with the membrane, while number density profiles along the direction perpendicular to the membrane were evaluated to quantify adsorption and interfacial structuring effects. Finally, the relative affinity of the membrane for CO₂ and N₂ was assessed by calculating the adsorption selectivity.

Fig. 6 shows the radial distribution functions (RDFs), g(r), between the gas molecules (CO₂ and N₂) and the carbon atoms of the graphtetrayne-based frameworks. These calculations were performed to evaluate the selective adsorption behavior of the new porous structures: as can be seen, the RDF profiles clearly indicate a preferential interaction of CO₂ molecules with the carbon atoms of the frameworks compared to N₂. This is evidenced by a significantly higher peak intensity for pairs CO₂–C_{graphtetrayne}, particularly in the short-range region (around 4 Å), which corresponds to the interlayer spacing of the material.

Fig. 6 Radial distribution functions, g(r), between the gas molecules (CO₂ and N₂) and the carbon atoms of the graphtetrayne bilayer in the new stacking configuration.

In Fig. 7, the mean number density profiles along the direction perpendicular to the Graph-4-yne plane, ρ(z), are reported for CO₂ and N₂. The positions of the layers are indicated by vertical orange lines.

Fig. 7 Z-density profile for the CO₂/N₂ mixture with the bilayer membrane, calculated at 1.26 atm.

As can be seen, far from the membrane (|z| ≥ 15 Å), both molecules exhibit nearly constant density values, representative of bulk-like behavior. Upon approaching the surface of the Graph-4-yne, the number density shows well-defined peaks, indicating that molecules preferentially reside in the vicinity of the material as an effect of gas–surface interactions. Notably, CO₂ displays higher and sharper density maxima compared to N₂, suggesting stronger interactions with Graph-4-yne and a higher affinity for adsorption. The presence of symmetric density peaks on both sides of the plane reflects the two-dimensional nature of the material and the equivalence of the two interfaces. Overall, the density profiles highlight the preferential accumulation of CO₂ relative to N₂ near Graph-4-yne, pointing to a potential selectivity of the membrane toward CO₂. The interlayer peaks follow the same trend, indicating the preferential accumulation of CO₂ in the inner space.

To gain microscopic insight into the origin of the selective adsorption observed in molecular dynamics simulations, interaction energy profiles were computed for CO₂ and N₂ approaching the γ-graphyne–type pore along the direction perpendicular to the membrane plane. The interaction energy was evaluated as a function of the distance from the pore center using the improved Lennard-Jones potential (Fig. 8).

Fig. 8 Interaction potential (ILJ) profiles for CO₂ and N₂ approaching the center of a γ-graphyne (G4) pore along the direction perpendicular to the membrane plane. Top: single-layer G4. Bottom: G4 bilayer in the new int stacking.

Two structural models were considered. In the first case (Fig. 9, top), the calculations were performed on a single γ-graphyne layer, modeled as an isolated triangular pore unit. Due to the relatively large aperture of the pore, the interaction with both molecules is weak and characterized by a shallow minimum at the center of the pore; nevertheless, CO₂ systematically exhibits a deeper potential well than N₂, consistent with its larger quadrupole moment and stronger dispersive interactions with the carbon framework. In the second case (Fig. 9, lower panel), the same analysis was carried out for a bilayer adopting the new intermediate stacking, previously identified as the most stable configuration based on combined classical and ab initio calculations. In this stacked geometry, the effective pore size is significantly reduced by the relative alignment of the two layers, leading to a pronounced enhancement of confinement effects. As a consequence, the interaction potentials become markedly deeper, with a substantially stronger stabilization of CO₂ at the pore center compared to N₂. This increased energetic discrimination in the bilayer provides a direct molecular-level explanation for the enhanced CO₂ selectivity observed in the mixture simulations and confirms the crucial role of stacking-induced pore modulation in determining the molecular sieving properties of the membrane.

Fig. 9 Interaction potential profile evaluated using the ILJ formulation (circle points) and at the DFT level (B2PLYP-D3, solid lines) for CO₂ and N₂ approaching the center of a single layer γ-graphyne (G4) pore along the direction perpendicular to the membrane plane.

To validate the reliability of the ILJ-based interaction model, additional benchmark calculations were performed on reduced monolayer + gas systems (CO₂ and N₂). In these configurations, interaction energies were computed both using the ILJ potential and at the DFT level using the B2PLYP-D3 functional, which contains a MP2 correction and explicitly accounts for dispersion interactions. These results confirm the accuracy of the ILJ description in reproducing the relative interaction strengths between CO₂, N₂, and the carbon surface, while enabling its application to the full bilayer systems at a tractable computational cost.

Fig. 10 shows a representative snapshot taken from the final stage of the molecular dynamics trajectory, illustrating the spatial distribution of CO₂ and N₂ molecules in the vicinity of the graphtetrayne bilayer. The membrane preserves its structural integrity throughout the simulation and provides two equivalent adsorption interfaces, one on each side of the bilayer. Both gas species are observed to interact with the membrane surface; however, CO₂ molecules preferentially accumulate closer to the graphtetrayne layers, whereas N₂ molecules remain more homogeneously distributed in the surrounding gas phase. Several CO₂ molecules are found in proximity to the membrane and within the interlayer region, suggesting favorable interactions with the carbon framework. This qualitative picture is consistent with the enhanced CO₂ affinity inferred from the structural and density analyses discussed below.

Fig. 10 Snapshots of the final configuration obtained for the CO₂/N₂ mixture with the bilayer membrane, calculated at 1.26 atm. The orientation of the structure is explicitly indicated in the figure for clarity. The system is modeled using fully periodic boundary conditions in all directions, ensuring the absence of edge effects.

The composition of the adsorbed phase was quantified in terms of the adsorbed mole fraction (x^ads_i) and the adsorption selectivity (S_CO2/N2). The adsorbed mole fractions and can be defined as the ratio between the numbers of the single species adsorbed (N_CO2 and N_N2) and the total number of adsorbed molecules N, both estimated at equilibrium, as follows:

The adsorption selectivity S_CO2/N2 was used to assess the preferential uptake of CO₂ over N₂ and was calculated as the ratio between the CO₂/N₂ mole ratio in the adsorbed phase and that in the bulk gas phase,

where y denotes the mole fraction of the species in the gas phase. For an equimolar feed CO₂/N₂ (y_CO2 = y_N2 = 0.5), the selectivity reduces to S_CO2/N2 = N_CO2/N_N2.

In line with the emerging picture, the adsorption results confirm a preferential uptake of CO₂ over N₂ by the bilayer. Of a total of 25 adsorbed molecules, 17 were CO₂ and 8 were N₂, corresponding to adsorbed mole fractions of and . These values highlight a significant enrichment of CO₂ in the adsorbed phase in relation to the equimolar composition of the bulk gas. The adsorption selectivity, S_CO2/N2 = 2.125, quantifies this preference, indicating that the bilayer adsorbs more than twice as many CO₂ molecules as N₂ under the same conditions. In general, the high adsorbed mole fraction and selectivity values demonstrate the strong affinity of the bilayer for CO₂, consistent with its potential for selective gas capture applications.

Even though the selectivity value obtained here is lower than those reported for other porous carbon-based and hybrid membranes (see Table 3, where data and references for comparison are reported for a series of materials), it should be emphasized that the present results are intended as a proof-of-concept assessment of the separation capabilities of G4. A direct quantitative comparison—most notably with G3 systems—must indeed be made with caution, as different definitions of selectivity are involved. In particular, in the study of G3 membranes,³¹ selectivity was inferred from permeation rates or from free-energy barriers associated with molecular transport across the membrane, rather than from adsorption statistics directly extracted from MD simulations. In contrast, in the present work selectivity is defined in terms of the relative populations of CO₂ and N₂ molecules transiently confined within the pores, providing a simpler and more conservative descriptor of preferential affinity. On this basis, the non-negligible preference toward CO₂ already observed for G4 is likely to represent a lower bound, and further optimization—such as extended simulation times, pressure-dependent investigations, and free-energy-based permeation analyses analogous to those employed for G3—may reveal significantly enhanced separation performance in future studies.

Table 3 Adsorbed mole fraction x and selectivity S values obtained from the MD calculations for the CO₂/N₂ mixture with bilayer graphtetrayne membranes, calculated at 1.26 atm. For comparison, selectivities of different materials are included, such as graphtriyne, single- and trilayer (G3, G3tri), fluorine-modified porous graphene (FMod-Gr), covalent organic polymers (COP) and metal–organic frameworks (MEF)

Materials	xCO2	xN2	S
G4	0.68	0.32	2.1 (this work)
G3	—	—	5.4^31
FMod-Gr	—	—	4.0–11.0^64
G3-tri	—	—	11.3–12.5^31
COP	—	—	8.4–13.7^31
MEF	—	—	5-0–40.0^4,65

---

These findings suggest that the newly designed bilayer structures possess promising CO₂-selective adsorption characteristics, making them suitable materials for potential applications in membrane-based gas separation processes aimed at capturing CO₂ from highly polluted or industrial exhaust streams. The pronounced affinity of CO₂ is likely mediated by quadrupole–π interactions between the linear CO₂ molecule and the electron-rich pore walls, in combination with the accessible pore morphology that facilitates a favorable positioning of CO₂ near the interlayer carbon sites. Importantly, the pore diameter is only moderately larger than the kinetic diameters of small gas molecules such as CO₂ and N₂, providing a balance between adequate diffusion and confinement effects. This balance appears to preferentially enhance the adsorption of CO₂ over N₂, as evidenced by the observed adsorbed mole fractions and the adsorption selectivity (S_CO2/N2 = 2.125).

The higher population of CO₂ within the pores suggests that the bilayer architecture not only provides stronger confinement for CO₂ but also offers a greater number of energetically favorable adsorption sites compared to N₂. The combination of molecular shape, linearity, and quadrupole moment of CO₂ enables more effective interactions with the electron-rich carbon sites, while the slightly smaller polarizability and spherical geometry of N₂ limit its ability to achieve similarly favorable contacts. Moreover, the interlayer spacing and local pore environment appear to promote optimal packing of CO₂ molecules, minimizing steric hindrance and maximizing van der Waals contacts. This selective confinement effect contributes directly to the observed enrichment of CO₂ in the adsorbed phase, demonstrating that the bilayer structure inherently favors CO₂ uptake even in the presence of competing gases. Overall, these results highlight the potential of the bilayer as a selective adsorbent and provide a detailed microscopic understanding of the factors governing CO₂/N₂ separation, which can inform the rational design of membrane materials with enhanced selectivity and capacity.

4 Conclusions

By combining classical molecular dynamics with quantum-mechanical calculations, we have investigated the structural and adsorption properties of an element of an emerging class of carbonaceous materials, namely graphtetrayne. Building on recent literature data, we show that the landscape of stable stacking configurations is broader than previously suggested by purely theoretical predictions. In particular, we identify an intermediate stacking arrangement between the AA configuration and the previously reported (nominally nine-fold) AB stacking that emerges as the most thermodynamically stable among all investigated structures. This conclusion is consistently supported for both isolated bilayers and extended stacked systems.

The excellent agreement between molecular dynamics simulations and ab initio calculations confirms the robustness of the ILJ-based force field in describing interlayer interactions, thereby justifying its use for more complex architectures. On this basis, we confidently extended the analysis to the trilayer configuration, obtaining structurally consistent and physically stable models.

From an application-oriented perspective, we compared the adsorption behavior of the newly identified intermediate bilayer with that of the AB-stacked system. Both configurations exhibit an enhanced preference for CO₂ over N₂, indicating that alternative stacking arrangements of graphtetrayne can be regarded as promising candidates for the design of carbon-based membranes for gas separation. Molecular dynamics simulations of an equimolar CO₂/N₂ mixture further corroborate this conclusion: starting from an initially equimolar gas phase, both bilayer systems display a clear enrichment of CO₂ within the adsorbed phase at equilibrium, resulting in adsorbed mole fractions that markedly deviate from the initial composition and translate into a significant adsorption selectivity in favor of carbon dioxide.

The preferential uptake of CO₂ can be attributed to a combination of favorable quadrupole–π interactions with the electron-rich carbon framework, the linear geometry of the CO₂ molecule, and the accessible pore morphology, which together promote efficient confinement and optimal packing within the interlayer region. These results provide a clear assessment of the intrinsic adsorption selectivity of graphtetrayne bilayers in a mixed-gas environment.

Overall, this work represents a proof-of-concept study aimed at elucidating the structure–property relationships governing gas adsorption in graphtetrayne-based systems. Future investigations will systematically explore the effects of temperature, pressure, gas composition, and bilayer flexibility, as well as adsorption from pure gases, in order to establish pressure-dependent selectivity trends and to further clarify the role of structural dynamics. Such studies will be essential for informing the rational design of next-generation carbon-based membranes for CO₂ capture and separation applications.

Moreover, the present study establishes the validity of the proposed methodology through its application to representative systems that capture the essential physical features of the problem. Although extended comparison to a broader range of experimentally characterized structures could further generalize the analysis, such an effort would require a dedicated and systematic investigation beyond the scope of this work. This direction will be explored in future studies to further expand the applicability of the approach.

Author contributions

L. Mancini: investigation, visualization, writing – original draft. G. Giorgi: conceptualization, supervision, formal analysis, writing – review. Y. B. Apriliyanto: investigation, writing – review draft. A. Lombardi: conceptualization, supervision, formal analysis, funding acquisition, writing – review and editing. N. Faginas-Lago: conceptualization, supervision, formal analysis, funding acquisition, writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are available within the article and its supplementary information (SI). This includes optimized structural models of the bilayer and trilayer graphtetrayne systems, interaction parameters used in the molecular dynamics simulations, and representative analysis data underlying the figures and tables. No external datasets were generated or deposited in a public repository as part of this work, and no custom software or source code beyond standard, publicly available simulation packages was developed. Supplementary information is available. See DOI: https://doi.org/10.1039/d6cp00427j.

Acknowledgements

This work has been funded by the European Union – NextGenerationEU under the Italian Ministry of University and Research (MUR) National Innovation Ecosystem grant ECS00000041 – VITALITY – CUP J97G22000170005, GG similarly acknowledges the HYdrogen TRAceability (HYTRA) Project (CUP F57G25000150006). GG is similarly grateful to the ISCRA B and C initiatives for awarding access to computing resources on Leonardo at CINECA SuperComputer Center, Italy.

References

1. E. Braun, A. F. Zurhelle, W. Thijssen, S. K. Schnell, L.-C. Lin, J. Kim, J. A. Thompson and B. Smit, Mol. Syst. Des. Eng., 2016, 1, 175–188.
2. Z. Xiang, R. Mercado, J. M. Huck, H. Wang, Z. Guo, W. Wang, D. Cao, M. Haranczyk and B. Smit, J. Am. Chem. Soc., 2015, 137, 13301–13307.
3. A.-H. Lu and G.-P. Hao, Annu. Rep. Prog. Chem., Sect. A:Inorg. Chem., 2013, 109, 484–503.
4. X. Zhou, W. Huang, J. Miao, Q. Xia, Z. Zhang, H. Wang and Z. Li, Chem. Eng. J., 2015, 266, 339–344.
5. H. Du, J. Li, J. Zhang, G. Su, X. Li and Y. Zhao, J. Phys. Chem. C, 2011, 115, 23261–23266.
6. Y. Tao, Q. Xue, Z. Liu, M. Shan, C. Ling, T. Wu and X. Li, ACS Appl. Mater. Interfaces, 2014, 6, 8048–8058.
7. M. Alaghemandi, Chem. Phys. Lett., 2015, 629, 65–69.
8. Z. Meng, X. Zhang, Y. Zhang, H. Gao, Y. Wang, Q. Shi, D. Rao, Y. Liu, K. Deng and R. Lu, ACS Appl. Mater. Interfaces, 2016, 8, 28166–28170.
9. M. Palummo, K. Yamashita and G. Giorgi, Sustainable Strategies in Organic Electronics, Elsevier, 2022, pp. 391–422.
10. J. S. Bunch, S. S. Verbridge, J. S. Alden, A. M. Van Der Zande, J. M. Parpia, H. G. Craighead and P. L. McEuen, Nano Lett., 2008, 8, 2458–2462.
11. D.-e Jiang, V. R. Cooper and S. Dai, Nano Lett., 2009, 9, 4019–4024.
12. S. P. Koenig, L. Wang, J. Pellegrino and J. S. Bunch, Nat. Nanotechnol., 2012, 7, 728–732.
13. S. P. Surwade, S. N. Smirnov, I. V. Vlassiouk, R. R. Unocic, G. M. Veith, S. Dai and S. M. Mahurin, Nat. Nanotechnol., 2015, 10, 459–464.
14. K.-J. Hsu, S. Li, M. Micari, H.-Y. Chi, L. F. Villalobos, S. Huang, L. Zhong, S. Song, X. Duan and A. Züttel, et al., Nat. Energy, 2024, 9, 964–974.
15. H. Wang, T. Maiyalagan and X. Wang, ACS Catal., 2012, 2, 781–794.
16. H. Liang, F. Ming and H. N. Alshareef, Adv. Energy Mater., 2018, 8, 1801804.
17. Y. Wang, Y. Shao, D. W. Matson, J. Li and Y. Lin, ACS Nano, 2010, 4, 1790–1798.
18. L. Mancini, G. Vanuzzo, P. Recio, A. Caracciolo, N. Faginas-Lago, M. Rosi, P. Casavecchia and N. Balucani, J. Phys. Chem. A, 2024, 128, 7177–7194.
19. E. Singh, M. Meyyappan and H. S. Nalwa, ACS Appl. Mater. Interfaces, 2017, 9, 34544–34586.
20. N. Faginas-Lago, Y. B. Apriliyanto and A. Lombardi, Eur. Phys. J. D, 2021, 75, 161.
21. N. Faginas-Lago, Y. B. Apriliyanto and A. Lombardi, International Conference on Computational Science and Its Applications, 2019, pp. 374-387.
22. R. Baughman, H. Eckhardt and M. Kertesz, J. Chem. Phys., 1987, 87, 6687–6699.
23. D. Malko, C. Neiss, F. Vines and A. Görling, Phys. Rev. Lett., 2012, 108, 086804.
24. M. Bartolomei and G. Giorgi, ACS Mater. Lett., 2024, 6, 4682–4689.
25. G. Li, Y. Li, H. Liu, Y. Guo, Y. Li and D. Zhu, Chem. Commun., 2010, 46, 3256–3258.
26. M. Long, L. Tang, D. Wang, Y. Li and Z. Shuai, ACS Nano, 2011, 5, 2593–2600.
27. M. Bartolomei, E. Carmona-Novillo, M. I. Hernández, J. Campos-Martnez, F. Pirani and G. Giorgi, J. Phys. Chem. C, 2014, 118, 29966–29972.
28. M. Bartolomei, E. Carmona-Novillo, M. I. Hernández, J. Campos-Martnez, F. Pirani, G. Giorgi and K. Yamashita, J. Phys. Chem. Lett., 2014, 5, 751–755.
29. M. Bartolomei and G. Giorgi, ACS Appl. Mater. Interfaces, 2016, 8, 27996–28003.
30. Y. B. Apriliyanto, N. Faginas-Lago, S. Evangelisti, M. Bartolomei, T. Leininger, F. Pirani, L. Pacifici and A. Lombardi, Molecules, 2022, 27, 5958.
31. Y. B. Apriliyanto, N. Faginas Lago, A. Lombardi, S. Evangelisti, M. Bartolomei, T. Leininger and F. Pirani, J. Phys. Chem. C, 2018, 122, 16195–16208.
32. N. Faginas-Lago, Y. B. Apriliyanto and A. Lombardi, International Conference on Computational Science and Its Applications, 2020, pp. 489-501.
33. M. Bartolomei, E. Carmona-Novillo and G. Giorgi, Carbon, 2015, 95, 1076–1081.
34. Q. Pan, S. Chen, C. Wu, F. Shao, J. Sun, L. Sun, Z. Zhang, Y. Man, Z. Li and L. He, et al., CCS Chem., 2021, 3, 1368–1375.
35. Q. Li, Y. Li, Y. Chen, L. Wu, C. Yang and X. Cui, Carbon, 2018, 136, 248–254.
36. J. Gao, J. Li, Y. Chen, Z. Zuo, Y. Li, H. Liu and Y. Li, Nano Energy, 2018, 43, 192–199.
37. V. G. Desyatkin, W. B. Martin, A. E. Aliev, N. E. Chapman, A. F. Fonseca, D. S. Galvão, E. R. Miller, K. H. Stone, Z. Wang and D. Zakhidov, et al., J. Am. Chem. Soc., 2022, 144, 17999–18008.
38. X. Gao, Y. Zhu, D. Yi, J. Zhou, S. Zhang, C. Yin, F. Ding, S. Zhang, X. Yi and J. Wang, et al., Sci. Adv., 2018, 4, eaat6378.
39. X.-H. Wang, Z.-C. Zhang, J.-J. Wang, X.-D. Chen, B.-W. Yao, Y.-X. Hou, M.-X. Yu, Y. Li and T.-B. Lu, ACS Appl. Mater. Interfaces, 2020, 12, 33069–33075.
40. J. M. Marques, F. V. Prudente and F. Pirani, Intermolecular forces: From atoms and molecules to nanostructures, 2022.
41. F. Pirani, S. Brizi, L. F. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi, Phys. Chem. Chem. Phys., 2008, 10, 5489–5503.
42. W. Smith, C. Yong and P. Rodger, Mol. Simul., 2002, 28, 385–471.
43. P. P. Ewald, Annalen der Physik, 1921, 369, 253–287.
44. F. Pirani, M. Albert, A. Castro, M. M. Teixidor and D. Cappelletti, Chem. Phys. Lett., 2004, 394, 37–44.
45. F. Pirani, G. S. Maciel, D. Cappelletti and V. Aquilanti, Int. Rev. Phys. Chem., 2006, 25, 165–199.
46. G. Liuti and F. Pirani, Chem. Phys. Lett., 1985, 122, 245–250.
47. M. Bartolomei, F. Pirani, A. Lagana and A. Lombardi, J. Comput. Chem., 2012, 33, 1806–1819.
48. A. Lombardi, A. Laganà, F. Pirani, F. Palazzetti and N. F. Lago, International Conference on Computational Science and Its Applications, 2013, pp. 17-31.
49. A. Lombardi, N. Faginas-Lago, L. Pacifici and G. Grossi, J. Chem. Phys., 2015, 143, 034307.
50. A. Lombardi, F. Pirani, A. Lagana and M. Bartolomei, J. Comput. Chem., 2016, 37, 1463–1475.
51. A. Lombardi, N. Faginas-Lago, G. Gaia, P. Federico and V. Aquilanti, International Conference on Computational Science and Its Applications, 2016, pp. 246-257.
52. A. Lombardi, F. Pirani, M. Bartolomei, C. Coletti and A. Laganà, Front. Chem., 2019, 7, 309.
53. L. D. Pietanza, O. Guaitella, V. Aquilanti, I. Armenise, A. Bogaerts, M. Capitelli, G. Colonna, V. Guerra, R. Engeln and E. Kustova, et al., Eur. Phys. J. D, 2021, 75, 237.
54. Q. Hong, M. Bartolomei, C. Coletti, A. Lombardi, Q. Sun and F. Pirani, Molecules, 2021, 26, 7152.
55. L. Mancini, E. V. F. de Aragão, F. Pirani, M. Rosi, N. Faginas-Lago, V. Richardson, L. M. Martini, L. Podio, M. Lippi and C. Codella, et al., Astron. Astrophys., 2024, 691, A83.
56. M. Michielan, L. Mancini, D. Ascenzi, M. Rosi, P. Tosi, F. Pirani, N. Balucani, D. Skouteris and C. Ceccarelli, Astron. Astrophys., 2025, 698, A205.
57. M. Rosi, D. Skouteris, N. Balucani, L. Mancini, N. F. Lago, L. Podio, C. Codella, B. Lefloch and C. Ceccarelli, International Conference on Computational Science and Its Applications, 2019, pp. 306-315.
58. L. Mancini, M. Rosi, D. Skouteris, G. Vanuzzo, G. Pannacci, P. Casavecchia and N. Balucani, Comput. Theor. Chem., 2023, 1229, 114341.
59. J. M. Soler, E. Artacho, J. D. Gale, A. Garca, J. Junquera, P. Ordejón and D. Sánchez-Portal, J. Phys.: Condens. Matter, 2002, 14, 2745.
60. E. Artacho, E. Anglada, O. Diéguez, J. D. Gale, A. Garca, J. Junquera, R. M. Martin, P. Ordejón, J. M. Pruneda and D. Sánchez-Portal, et al., J. Phys.: Condens. Matter, 2008, 20, 064208.
61. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
62. N. Troullier and J. L. Martins, Phys. Rev. B:Condens. Matter Mater. Phys., 1991, 43, 1993.
63. L. Mancini, F. Pirani, Y. B. Apriliyanto, A. Lombardi and N. Faginas-Lago, A Snapshot of Molecular Electronic Structure Theory and its Applications, J. Mol. Struct., 2025, 92, 179.
64. T. Wu, Q. Xue, C. Ling, M. Shan, Z. Liu, Y. Tao and X. Li, J. Phys. Chem. C, 2014, 118, 7369–7376.
65. S. Ghosh, M. Sevilla, A. B. Fuertes, E. Andreoli, J. Ho and A. R. Barron, J. Mater. Chem. A, 2016, 4, 14739–14751.

---