On assessing the carbon capture performance of graphynes with particle swarm optimization

Abstract

Tackling climate change is one of the greatest challenges of current times and therefore the development of efficient technologies to limit anthropogenic emissions is of utmost urgency. Recent research towards this goal has alluded to the use of carbon-based solid sorbents for carbon capture. Graphynes (GYs), an interesting class of porous carbon membranes, have recently proven their potential as excellent membranes for gas adsorption and separation. Herein, we explored the CO₂ and N₂ adsorption characteristics and CO₂/N₂ selectivities of a class of GYs, namely γ-GY-1, γ-GY-2 and γ-GY-4. We investigated the putative global minimum geometries of adsorbed unary (n = 2–10) and binary (n : m; n, m ∈ [1, 8]) clusters of CO₂ and N₂ by employing a stochastic global optimization method called particle swarm optimization in conjunction with empirical intermolecular force field formulations. The intervening interactions are modeled using various pairwise potentials, including Lennard-Jones potential, improved Lennard-Jones potential, Buckingham potential and Coulombic potential. The binding energies for both unary and binary clusters are highest for adsorption on γ-GY-1, followed by γ-GY-2. The putative global minimum geometries suggested that N₂ molecules preferred binding over the pore centres while CO₂ molecules showed higher clustering propensity than any binding site preference. The predicted interaction energies suggested higher selectivity for CO₂ over N₂ for all the three γ-GYs.

---

Introduction

With the rising concerns on meeting the target of the 2015 Paris Agreement,¹ wherein 2 °C above pre-industrial levels is defined as the upper limit for global warming, the need for concrete technologies to curb anthropogenic emissions is more apparent than ever. As well established through numerous research and agreed upon by various nations, going beyond this limit will result in catastrophic effects including heat waves, droughts, extreme precipitation and even extinction of several species.² One of the major sources of these emissions is the burning of fossil fuels and hence, efficient CO₂ capture from post-combustion mixtures is vital in realizing the goals set by the Intergovernmental Panel on Climate Change. When the combustion of fossil fuels occurs in ambient air, the flue gas, comprising primarily dinitrogen, carbon dioxide and water vapour, is emitted.³ The current method for separating CO₂ from flue gas in order to reduce the carbon content in atmosphere is via chemical absorption using aqueous amines.⁴ However, this method is costly and highly energy expensive, and researchers are on the lookout for solid sorbents as an alternative.

Through selective physisorption of CO₂ on the surface or within the porous networks, solid sorbents have proven to be highly energy efficient for carbon capture.⁵ A solid sorbent designed for this purpose should preferably be of high thermal and mechanical stability so as to withstand numerous desorption cycles and must possess good adsorption capacity even under humid conditions since flue gas contains vapor content. Presently, nanoporous materials such as zeolites and metal–organic frameworks are the front-runners in this area, albeit with certain drawbacks.^6,7 Despite being used industrially for carbon capture owing to its high CO₂ uptake, zeolite 13X has poor CO₂/N₂ selectivity.⁸ Carbon-based sorbents are an emerging class of promising materials for post-combustion CO₂ separation. The high stability, large surface area for adsorption, and scope of functionalization make graphene a substrate worth exploring. However, the small size of the hexagonal pores of graphene considerably limits its gas-uptake capacity. The introduction of large pores in the carbon frameworks can offer an interesting alternative. Pirani and co-workers demonstrated the excellent capacity of the graphtriyne multilayers for CO₂/N₂ separation using molecular dynamics simulations.⁹ Graphynes (GYs), with their uniform and well-architectured intrinsic pores are therefore an excellent choice for efficient carbon capture.

GYs exist in various structural patterns and are classified according to the number and arrangement of acetylenic linkages.¹⁰ Theoretical investigations have revealed GYs to be ideal substrates for gas storage and separation.¹¹ For instance, Bartolomei et al. have characterized the hydrogen storage capacity of GYs with the help of first-principles calculations.¹² Similarly, hydrogen purification from a mixture of methane and carbon monoxide was investigated by Jiao et al. using density functional theory (DFT).¹³ The inferences drawn from these theoretical studies regarding the physisorption strengths of gases on GYs further signify the promising nature of GYs for carbon capture. Interactions between CO₂ molecules and GYs have been investigated under two scenarios: (i) sensing of CO₂ through physisorption on GYs, and (ii) separation of CO₂ from multicomponent gaseous mixtures via permeation through GY membranes. The well-known GYs such as γ-GY-1¹⁴ and T_4,4,4-GY¹⁵ have been explored for the physisorption of CO₂, whereas γ-GY-2 and rhombic-GY have been investigated for separation and purification of CO₂ from gas mixtures based on permeability.¹⁶ Two recent studies employing first-principles and molecular dynamics simulations have reported the preferential physisorption of CO₂ molecules on γ-GY-3 over N₂ and H₂O.^9,17 Research on carbon capture using various metal-decorated GYs is also on the rise. CO₂ capture using alkali/alkaline earth metal-decorated γ-GY-1,¹⁸ transition metal-decorated γ-GY-1,^19–21 transition metal-decorated γ-GY-2²² and Li-decorated β1-GY²³ are some of the recent examples.

In view of the pressing need for reducing carbon emissions and taking into cognizance of the recent research on GYs, herein, we investigate a series of one-atom-thick carbon membranes, namely, γ-GY-1, γ-GY-2 and γ-GY-4 for CO₂ and N₂ adsorption. It is pertinent to note that all the three above-mentioned forms of GYs have been synthesised in recent times.^24–26 The different types of γ-GYs are considered in order to investigate the effects of variation in carbon density and pore size of the sheets on the energetics and adsorption pattern of the clusters. To explore the preferential binding sites and energetics of adsorption of unary and binary clusters of CO₂ and N₂ on the GYs, we employ a swarm intelligence technique called particle swarm optimization (PSO).^27–31 Swarm intelligence techniques can be considered as a sub-field of artificial intelligence that work on the principle of ‘emergence’, wherein the members of the swarm act with traits of intelligence.³² This property enables swarm intelligence techniques such as PSO to efficiently locate the global minimum structure of a high-dimensional chemical system. The population-based nature of the algorithm provides the researchers with an added benefit of parallel processing which can considerably reduce the computational complexity of the problem. Researchers have successfully applied PSO to predict minimum energy geometries of sizeable atomic and molecular clusters in their bare, adsorbed and intercalated forms.^33–42 Very recently, we looked into the binding of dinitrogen clusters on monolayer and bilayer graphene with the help of PSO.⁴² Majority of the theoretical studies related to adsorption of CO₂ on GYs reported till date are carried out using DFT, molecular dynamics or grand canonical Monte Carlo simulations.^{9,17–19,43–45} Our work offers a new strategy through which CO₂/N₂ adsorption features and selectivity of GYs can be ascertained.

Methodology

In the current study, we obtained the most favourable geometries of (i) bare CO₂ clusters, (ii) bare binary clusters of CO₂ and N₂ as well as (iii) CO₂ clusters, N₂ clusters, and their binary clusters adsorbed over GY model systems. For the bare and adsorbed unary clusters, we have considered 2–10 molecules in the cluster, while for binary clusters, we have studied two cases: (a) binary clusters with equal number of CO₂ and N₂ molecules ranging from 1 to 5, and (b) binary clusters in the ratio n : m where n + m = 10. We have carried out the adsorption studies on molecular models of three different γ-GYs, namely, γ-GY-1, γ-GY-2, and γ-GY-4 (Fig. 1). The geometries of the annulenic molecular models of the GYs are optimized at the ωB97X-D/6-311G(d,p) level of DFT using the Gaussian 16 software.⁴⁶ This level of theory was used previously in many studies for geometry optimization, especially in closely related systems including graphene.^47–49 The Cartesian coordinates of the molecular models are provided in the ESI.†

Fig. 1 Molecular model systems of GYs: (a) γ-GY-1 (C₁₉₈H₃₀), (b) γ-GY-2 (C₂₈₂H₃₀), and (c) γ-GY-4 (C₁₃₈H₁₈).

Modeling the clustering and adsorption of small molecules is a complex multidimensional problem and the use of ab initio techniques for this is computationally prohibitive. While, a metaheuristic technique such as PSO can track down the putative global minimum, i.e., the most stable geometry, comparatively faster when used with the right objective function. A study on predicting the geometry of carbon clusters by Chattaraj and co-workers has shown that PSO is less computationally expensive and has a faster convergence rate when compared to genetic algorithm. The same study pointed out the higher tendency of simulated annealing to get trapped in local minima when compared to PSO.³⁴ Recently, we have also reported the merits of PSO over various global optimization techniques including differential evolution, basin-hopping and dual annealing for evaluating the putative global minimum configurations of adsorbed and intercalated noble gas clusters on graphdiyne model systems.³⁸ In PSO, we begin by randomly initializing a set of particles, called swarm, in the desired search domain and the initial velocity of each particle is taken to be zero. These particles represent various points on the potential energy surface; the total intermolecular interaction energy constitutes the objective function. The position of each particle is an n-dimensional vector where n is the number of variables of the objective function. The velocity of a particle at a given iteration is also an n-dimensional vector that determines the direction in which the particle should move next. At each iteration t, each member of the swarm keeps a record of the best position that it has achieved (P_best,i) and the swarm has achieved (G_best) over the iterations. With this knowledge, the particles communicate with each other and update their positions (x_i) and velocities (v_i) using the equations,

x^t+1_i = x^t_i + v^t+1_i (1)

v^t+1_i = χ[v^t_i + c₁r^t₁(P^t_best,i − x^t_i) + c₂r^t₂(G^t_best − x^t_i)] (2)

where i is the index used to represent the particle. Here, c₁ and c₂ are the cognitive and the social acceleration constants, respectively, that determine the contributions of P_best and G_best to the velocity. The values of c₁ and c₂ are taken as 2.05, whereas r₁ and r₂ are random numbers generated between 0 and 1.^28,50 The constriction factor, χ that takes the value 0.729 is a modification to the basic PSO algorithm in order to improve the convergence.⁵⁰ The sharing of information among the members of the swarm is dependent on the neighbourhood topology. In the present study, we employ the star topology, where all the members of the swarm are connected to each other. i.e., the best solution of the entire swarm is known to all the members of the swarm.²⁹ As the iterations proceed, the particles converge to the putative global minimum, which is then considered as the starting point for a deterministic search using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method.⁵¹ All the PSO calculations are performed using an in-house code written in Python. The local optimizations using L-BFGS are performed using an in-built package in the SciPy module.⁵² Since we investigated both bare and adsorbed clusters of various sizes, we have considered a range of swarm sizes (500–4000) and iteration numbers (500–1000). For each system, we have carried out 20–25 independent PSO-L-BFGS runs (20 for CO₂ bare clusters and 25 for the remaining systems) and the lowest energy structure among them is reported. The rationale behind choosing 25 independent PSO-L-BFGS runs is provided in the Results and discussion section.

Since the study involves clusters of molecules, we use a rigid-body approximation to include both translational and rotational degrees of freedom of the molecules during the global optimization.⁵³ For a molecular cluster of size n, each particle of the swarm has 6n dimensions arising from the Cartesian coordinates of the centre of mass of the molecule and the three Euler angles (θ, φ and ϕ) describing the orientation of each molecule, i.e., the position of a swarm particle is a vector with these 6n variables as elements and at each iteration of PSO, all particles update their positions by updating these 6n variables. This step is then followed by the function evaluation, which in the current study is the total interaction energy. Once the centre of mass Cartesian coordinates (r_COM) and Euler angles of each molecule are defined, it is easy to obtain the Cartesian coordinates of each atom of the molecule which is essential for calculating the interaction energy. The Cartesian coordinate of the i^th atom of a molecule is given by

r_i = r_COM + R⁻¹(θ, φ, ϕ)r_rel,i (3)

where r_rel,i is the Cartesian coordinate vector of the i^th atom of the molecule relative to the centre of mass of the molecule (body-fixed coordinates) and R⁻¹(θ, φ, ϕ) is the inverse of the rotation matrix. Throughout the process, we bound the values of r_COM to the positive half-space of the GY model system in order to restrict the adsorption to a single side of the model system of GY. This restriction is prompted by many experimental studies on the applications of graphene, including adsorption studies, wherein graphene supported on a substrate is preferred over the free-standing graphene.^54–56 Furthermore, there are recent reports suggesting that the synthesis of GYs can be facilitated on supported substrates.^57,58 Here, the adsorption is restricted to a single side by clamping the Z-coordinate of centre of mass of the molecule between 0 Å and 5 Å. The position elements corresponding to the X and Y coordinates of centre of mass of the molecule were clamped between the maximum and minimum values of the X and Y coordinates of the molecular models of GYs. For the Euler angles, we clamped their values in the range of possible rotations which depends on the choice of the rotation matrix. For the current study, we considered the same rotation matrix and the bounds as in our previous work.⁴² Even though this approach works for adsorption over γ-GY-1 and γ-GY-2, the large pore size of γ-GY-4 can lead to putative global minima with the molecules being aligned vertically in the centre of the triangular pore with one of the atoms of the molecule in the negative half-space of the annulenic model of γ-GY-4. Hence, in the case of γ-GY-4, we additionally constrained the position of each atom to the positive half-space of the sheet. Apart from this, to improve the exploration of the algorithm, we have incorporated dynamic velocity clamping. Here, we initially define an upper bound for each velocity element which is the average of the maximum and minimum value for the corresponding position element. The same multiplied with (−1) serves as a lower bound for the velocity elements. Later, in each iteration, t, we multiply both the upper and the lower bound of each velocity component by a factor of (t − 1)/t_max, where t_max is the total number of iterations.

The objective function that needs to be minimized in our case is the total intermolecular interaction energy of the system. The details of the interaction energy terms of the various systems are given below:

Interactions in bare CO₂ clusters

For modelling CO₂ clusters, we used the TraPPE model,^59,60 wherein the intermolecular interactions in bare CO₂ clusters are represented by the sum of electrostatic and non-electrostatic interactions. The model was proposed to reproduce the vapour–liquid equilibria of CO₂ and it's mixture with alkanes and it has been validated recently against experimental data using molecular dynamics simulations.⁶¹ The total interaction energy is given byThe electrostatic interactions are modeled using the Coulombic pair potential between the partial charges residing on each atom of the two interacting CO₂ molecules. A representative model of the charge distribution in CO₂ is given in Fig. S1 of the ESI.† For a cluster of size n, the electrostatic interactions are calculated usingwhere C_el = 332 kcal Å mol⁻¹ e⁻², q_ki (or q_li) represents the partial charge on the k^th (l^th) atom of the i^th (j^th) molecule and r_kilj is the distance between the k^th and l^th atoms of the i^th and j^th molecules, respectively. The non-electrostatic contribution to the interaction energy is described using the Lennard-Jones (LJ) potential as follows:where k and l represent the atoms of the i^th and j^th molecules, respectively. ε_kl and σ_kl are the LJ parameters describing the well depth and the van der Waals distance of the k–l atom–atom interaction, respectively. The numerical values of the parameters corresponding to C–C interactions and O–O interactions are taken from the literature and are provided in Table 1⁵⁹ and the parameters for C–O interactions are obtained using the Lorentz–Berthelot mixing rules.⁶²

Table 1 Numerical values of the parameters used for calculating the intermolecular interaction energies

LJ potential
	ε (kcal mol^−1)	σ (Å)
C–C	0.0535	2.80
O–O	0.1570	3.05

---

Buckingham-type potential
	A (kcal mol^−1)	α (Å^−1)	C 6 (kcal Å^6 mol^−1)
N–N	29 995.9	3.4614	392.274

---

ILJ potential
	ε (kcal mol^−1)	σ (Å)	β
CCO2–N	0.0784	3.548	9.00
O–N	0.1038	3.699	9.00
CCO2–CGY	0.0820	3.564	6.75
O–CGY	0.1156	3.676	6.75
N–CGY	0.0916	3.828	7.50

---

Interactions in bare N₂ clusters

Recently, we predicted the energetics and geometries of bare N₂ clusters using PSO.⁴² For the sake of completeness, we provide the details here. The N₂–N₂ interactions are modeled as follows: the N₂–N₂ interaction energy (E_N2–N2) is described as the sum of electrostatic and non-electrostatic interaction terms, i.e.,The electrostatic interactions in E_N2–N2 arise from the quadrupole moment of the nitrogen molecule. It can be described using a four-point charge model (Fig. S1 of the ESI†). This model was previously used by Bertolus et al. for studying N₂ clusters using the Monte Carlo growth method.⁶³ Thus, the electrostatic contributions are given byHere, m represents the number of N₂ molecules in the cluster. The non-electrostatic part is described using a Buckingham-type potential, which is commonly employed for diatomic molecules.⁶⁴ Hence, the non-electrostatic contributions in N₂ clusters are expressed aswhere A, α and C₆ are the parameters of the Buckingham-type potential which are taken from literature and are provided in Table 1.⁶³ The Buckingham-type potential was previously used for studying the structural and thermodynamic properties of N₂ clusters.⁶⁵ A study by Calvo et al. compared the structural properties of N₂ clusters using both (i) Monte Carlo simulations with interactions described using Buckingham-type potential, and (ii) electron diffraction experiments, and observed a good agreement for diffraction patterns from lowest-energy structures.⁶⁶

Interactions in bare CO₂–N₂ binary clusters

In CO₂–N₂ binary clusters, there are three types of possible interactions: CO₂–CO₂, N₂–N₂ and CO₂–N₂. Therefore, the total interaction energy is given by the following equation:

E_binary = E_CO2–CO2 + E_N2–N2 + E_N2–CO2 (10)

The CO₂–CO₂ and N₂–N₂ interactions are modeled as described earlier. E_N2–CO2 corresponds to the interaction energy contributions from N₂–CO₂ interactions, which in itself contain electrostatic and non-electrostatic contributions. The N₂–CO₂ interaction energy is given aswhere the electrostatic contributions are described using the Coulomb's law as given byHere, m and n indicate the number of N₂ and CO₂ molecules in the cluster. The non-electrostatic contribution in the N₂–CO₂ interaction is modeled using an improved Lennard-Jones (ILJ) potential. The ILJ potential is a modified version of the LJ potential as it provides a better description of the long-range and short-range interactions.⁶⁷ The form of the ILJ potential and the total non-electrostatic N₂–CO₂ interaction energy are given by the following equations:where ε_kl, r_mkl and β_kl are the ILJ parameters for the k–l atom–atom interaction, whose values are also taken from literature and are also provided in Table 1.⁹

Interactions in adsorbed clusters

When a molecular cluster gets adsorbed on a molecular model of GY, an additional interaction term, arising from the interaction of the cluster with the GY, needs to be accounted for. Therefore, the total intermolecular interaction energy of an adsorbed unary/binary cluster is given asand

In the above, both E_CO2–GY and E_N2–GY are expressed using the ILJ potential⁹ as

andHere, p represents the number of atoms in the annulenic model of GY. The ILJ parameters of CO₂–GY and N₂–GY interactions (Table 1) are taken from a literature report that fits the potential against MP2C calculations for describing the interactions of CO₂ and N₂ with γ-GY-3.⁹ We assessed the transferability of these potentials to the GY models considered in our study, γ-GY-1, γ-GY-2, and γ-GY-4. For this, we used single-pore models of γ-GY-1, γ-GY-2, γ-GY-3, and γ-GY-4 (Fig. S2 of the ESI†), and performed single-point energy scans for the passage of both CO₂ and N₂ in three different orientations (along X, Y, and Z axes) towards the centre of the model systems. The calculations were performed at the B3LYP-D3/aug-cc-pVTZ (hereafter denoted as B3LYP-D3/aTZ) level of theory since the results obtained at this level of theory for the interactions of CO₂/N₂ with the model compound of γ-GY-3 were in close agreement with the MP2C results reported in the literature⁹ (Fig. S3–S5 of the ESI†). A comparison of the energy profiles for the CO₂–GY and N₂–GY complexes obtained using the ILJ potential with those of DFT is also shown in Fig. S3–S5 of the ESI.† The interaction energy profiles evaluated using the ILJ potential are in reasonable agreement with those obtained using DFT, especially at the equilibrium points, establishing the transferability of the ILJ parameters across the various GY systems. Additionally, we performed DFT calculations at the B3LYP-D3/aTZ level of theory to generate potential energy profiles for the CO₂–CO₂, N₂–N₂, and CO₂–N₂ dimers. When compared against the profiles evaluated using DFT, our corresponding empirical potential profiles showed good agreement as shown in Fig. S6 of the ESI,† validating the use of the above-mentioned empirical potentials for modeling the intermolecular interactions in molecular clusters.

Results and discussion

Global optimization of bare unary and binary clusters

To begin with, we investigated the putative global minimum geometries of bare CO₂ clusters using PSO. The obtained putative global minimum geometries are shown in Fig. 2. As we can see from Fig. 2, the CO₂ dimer has a slipped parallel structure with C_2h symmetry and the CO₂ trimer has a planar cyclic structure with C_3h symmetry, both of which are consistent with previous high-resolution spectroscopic studies.⁶⁸ The (CO₂)₆ cluster also showed a symmetric structure with S₆ point group. All clusters, except n = 2, 3 and 6, belongs to C₁ point group. A closer look at them, however, reveals that some of the geometries are close to high-symmetry structures when we consider the symmetry with respect to the centre of mass of the molecules. For example, the (CO₂)₄ and (CO₂)₅ clusters possess geometries very close to those of trigonal pyramidal and trigonal bipyramidal structures, respectively. The geometries are same as those predicted by Takeuchi in 2008, wherein he optimized CO₂ clusters of size 4–40 using a heuristic and unbiased method in which the CO₂ molecules were described using the Murthy–O’Shea–McDonald model.⁶⁹ The trends in energetics of these clusters are depicted in Fig. 3. Due to the larger number of possible interactions, the magnitude of interaction energy (Fig. 3a) as well as interaction energy per molecule (Fig. 3b) increase as we go from n = 2 to 10. In Fig. 3, for comparison, we have also included the interaction energies of bare N₂ clusters from our previous study.⁴² The interaction energies of CO₂ clusters are nearly four times those of the N₂ clusters, revealing the strong binding between CO₂ molecules. Fig. 3c represents the variation of the second difference in interaction energy as a function of cluster size. For a cluster of size n with an interaction energy E(n), the second difference of interaction energy is calculated as Δ₂E(n) = E(n + 1) + E(n − 1) − 2E(n). The second difference of interaction energy helps to identify clusters of high stability, known as magic number clusters. Magic number clusters are the clusters for which a maximum is observed in the plot. For bare CO₂ clusters, we obtained maxima at cluster sizes 4, 6 and 8 making them the magic number clusters, whereas for N₂, the magic number clusters obtained are clusters of size 4 and 7. As suggested by Bretón and co-workers in their study on rare gases adsorbed over graphene,⁷⁰ however, caution is warranted while we interpret these data. The putative global minimum geometries of N₂ clusters that we had obtained in our previous study are given in Fig. S7 of the ESI.† The obtained geometries were isomorphic to Ar_n clusters (n = 2–10), and exhibited symmetric structures when the centre of mass of each molecule is considered.

Fig. 2 Putative global minimum geometries of bare CO₂ clusters evaluated using PSO.

Fig. 3 A comparison of the (a) interaction energies, (b) interaction energies per molecule and (c) second difference of interaction energies of bare CO₂ and N₂⁴² clusters of size n = 2–10 evaluated using PSO.

We then looked into the putative global minimum geometries of bare binary clusters of CO₂ and N₂ molecules. There are no studies in literature reporting the global minimum geometries of bare binary clusters of CO₂ and N₂ and thus we explored the same using PSO. As mentioned in the methodology, we considered two cases: clusters with equal number of CO₂ and N₂ molecules and clusters with a total of 10 molecules with different ratios of CO₂ and N₂ molecules. As is evident from the putative global minimum geometries of these clusters given in Fig. 4, in most of the cases, CO₂ molecules form clusters while N₂ molecules prefer to surround the CO₂ molecules. For example, in a mixture containing 9 CO₂ molecules and 1 N₂ molecule, a cluster of CO₂ molecules with the N₂ molecule lying outside the cluster is observed, and for a mixture of 9 N₂ molecules and 1 CO₂ molecule, the single CO₂ molecule is surrounded by all the N₂ molecules. This observation can be explained on the basis of the interaction strengths of the various intermolecular components. The magnitude of the various interaction energy components follows the order: CO₂–CO₂ > CO₂–N₂ > N₂–N₂. The same accounts for the trends seen in Fig. 5 depicting the variation in the interaction energies and their various contributing energy terms of binary clusters. The various contributing energy terms include the interaction energies resulting from the interactions among the CO₂ molecules (CO₂–CO₂ energy), the interactions among the N₂ molecules (N₂–N₂ energy), and the interactions between the N₂ and the CO₂ molecules (N₂–CO₂ energy). For clusters with equal number of CO₂ and N₂ molecules, as the cluster size increases, the interaction energies showed a linear decrease, which is similar to the trend observed for unary clusters of CO₂ and N₂ molecules. For clusters of size 10 with varying CO₂ : N₂ ratios, with an increase in the number of CO₂ molecules, more negative interaction energies are obtained (see Fig. 5b). In order to validate the trends in energetics of the intervening interactions, we performed symmetry-adapted perturbation theory (SAPT)⁷¹ calculations. Using SAPT, the interaction energies and the various contributing terms including electrostatic, exchange, induction and dispersion of CO₂–CO₂, N₂–N₂ and CO₂–N₂ dimers obtained using PSO were evaluated. The calculations were performed at the SAPT0/AUG-cc-pVDZ level of theory using the Psi4 program package.⁷² The interaction energies and their components obtained are tabulated in Table S1 of the ESI.† With respect to all the contributions, the magnitudes of the derived interaction energies follow the same order as previously mentioned, i.e., CO₂–CO₂ > CO₂–N₂ > N₂–N₂. A major contribution from electrostatic (which is the leading contribution) and dispersion interactions makes the CO₂–CO₂ interaction most stable among all the three. For the N₂–N₂ interaction, the most significant contribution to the total interaction energy arises from the dispersion energy with the induction and electrostatic contributions being relatively small. Whilst for the CO₂–N₂ interaction, even though dispersion interaction forms a major contribution to the total interaction energy, electrostatic interactions also have a substantial role to play.

Fig. 4 Putative global minimum geometries of bare CO₂ and N₂ binary clusters of different N₂ : CO₂ (n : m) ratios evaluated using PSO.

Fig. 5 Variation of interaction energies and constituent contributions of bare CO₂ and N₂ binary clusters as a function of N₂ : CO₂ ratio (n : m) for (a) n = m = 1 to 5 and (b) n + m = 10, evaluated using PSO.

Global optimization of adsorbed unary clusters

Next, we investigated the structures and energetics of putative global minimum geometries of unary N₂ and CO₂ clusters of size 2 to 10 adsorbed over molecular models of three GYs namely γ-GY-1, γ- GY-2 and γ-GY-4. The adsorption features are strongly governed by the molecular size as well as the pore size of the GY substrate. Hence, we obtained the kinetic diameters of N₂ and CO₂ from literature,⁷³ and the pore diameters of the triangular pores of various γ-GYs using a Python library, namely pywindow.⁷⁴ The kinetic diameter of N₂ molecule is 3.64 Å while that of CO₂ is 3.30 Å. The triangular pore diameters of the GYs are 0.78 Å, 2.27 Å and 5.19 Å for γ-GY-1, γ-GY-2 and γ-GY-4, respectively. From this, we infer that the closest binding of the molecules will be observed in the case of γ-GY-4. Since the pore size of γ-GY-4 is larger than the diameters of N₂ and CO₂ molecules, the γ-GY-4 sheet allows easy passage of the molecules through its triangular pores. However, we restricted the current study to a single-sided adsorption by adding constraints in the algorithm, as mentioned in the Methodology section. This allowed us to probe the adsorption features, rather than the permeation features.

A few representative global minimum geometries of adsorbed N₂ and CO₂ clusters are shown in Fig. 6 and 7, respectively and the geometries of the rest of the adsorbed clusters are given in Fig. S8–S10 of the ESI.† Since we have carried out multiple PSO trials, we found other local minimum geometries which are very close in energy to the reported global minimum geometries. For example, in the ESI,† as a representative case, we have provided a few local minimum geometries of (N₂)₄ adsorbed on γ-GY-4 (Fig. S11 of the ESI†). The interaction energies corresponding to the various local minima geometries are very close to the interaction energy of the reported putative global minimum. The predicted geometries are also in good correspondence. The N₂ molecules, when adsorbed over γ-GY-1, showed an interesting pattern of adsorption by occupying the positions above both the hexagonal and triangular pore centres. The patterns followed by N₂ clusters of size n = 2–7, except n = 3, adsorbed over γ-GY-1 are similar to those of N₂ clusters adsorbed on circumcircumcoronene, which represents a model for graphene.⁴² Although the molecules are seen to adsorb as a monolayer on γ-GY-1 similar to that of adsorption on graphene, they adsorb much closer to the γ-GY-1 sheet (∼3.1 Å) when compared to graphene (∼3.4 Å).⁴² This can be attributed to the lower carbon atom density of γ-GY-1 when compared to graphene, which leads to reduced repulsive interactions in γ-GY-1 that impede the molecules from being proximal to the sheets. Upon moving to γ-GY-2, the pattern fades and the molecules even adsorb over the carbon skeleton along with the positions above the triangular pore centres. Unlike the case of adsorption on γ-GY-1 sheet where N₂ clusters nearly formed planar configurations, on γ-GY-2, the clusters form configurations in which some of the N₂ molecules are tilted towards the sheet. The adsorption configurations on γ-GY-4 are found to be different from those on both γ-GY-1 and γ-GY-2. This can be attributed to the larger pore size of γ-GY-4 which forces the molecules to cluster around the carbon skeleton of the sheet.

Fig. 6 Putative global minimum geometries (top and side views) of N₂ clusters of size n = 2, 4, 7 and 10 adsorbed on (a) γ-GY-1, (b) γ-GY-2 and (c) γ-GY-4 evaluated using PSO.

Fig. 7 Putative global minimum geometries (top and side views) of CO₂ clusters of size n = 2, 4, 7 and 10 adsorbed on (a) γ-GY-1, (b) γ-GY-2 and (c) γ-GY-4 evaluated using PSO.

The adsorption of CO₂ molecules over γ-GY-1 and γ-GY-2 sheets showed an entirely different behaviour. The CO₂ molecules showed more preference for clustering without any specific binding site preferences. This is a result of strong CO₂–CO₂ interactions present in adsorbed CO₂ clusters. Interesting growth patterns for CO₂ clusters over γ-GY-1 and γ-GY-2 can be observed in Fig. 7. For γ-GY-1, the molecules adsorbed at around 3 Å from the sheet forming a monolayer. As explained before for the case of N₂ adsorption over γ-GY-4, the CO₂ molecules also prefer to adsorb near the carbon skeleton of γ-GY-4. There are few reports in literature that have studied CO₂ adsorption on pristine γ-GYs. In 2017, Kwon et al. studied the adsorption of CO₂ and H₂ on γ-GY-1 using DFT-D3 (PBE/DNP) and they predicted that CO₂ adsorbs over the triangular pore site.¹⁴ Similar conclusion is drawn from our PSO calculations too. In 2018, Lu et al. also predicted a similar geometry in their DFT study.⁷⁵ Later, Reisi-Vanani and co-workers investigated three different sites on γ-GY-1 for CO₂ adsorption using DFT-D3 (PBE/DNP), namely, hexagonal pore centre, triangular pore centre and a site above the acetylenic linkage.²⁰ They considered both horizontal and vertical orientations of a single CO₂ with respect to GY for optimization and they found the most favourable geometry to be the one in which CO₂ is adsorbed horizontally above the triangular pore, which is concordant with our PSO calculations. In their successive studies on γ-GY-1 too, they have reported the same site with the horizontal orientation as the most favourable geometry for CO₂ adsorption.^19,44,76 A very recent study by the same group has investigated the adsorption of single CO₂ molecule on γ-GY-2 as well using a DFT-D2 method.²² The geometry obtained is consistent with the putative global minimum geometry we obtained wherein the CO₂ molecule occupies a position above the triangular pore. Similar results have been obtained by Fang et al. for γ-GY-2 in their DFT-D3 study.⁷⁷ In the same study, they have analysed single N₂ adsorption as well and have predicted that N₂ also prefers to adsorb above the triangular pore centre, which is again similar to our PSO result. Reisi-Vanani and co-workers, have also analysed the structure and energy of CO₂ trimer adsorbed on γ-GY-1²⁰ and have reported a geometry similar to the one we have obtained using PSO. All the above observations validate that PSO is a viable approach for investigating gas molecule adsorption on GYs.

We have also looked into the energetics of adsorbed N₂ and CO₂ unary clusters. The plots comparing the interaction energies for the adsorption over the model systems of three different GYs for both N₂ and CO₂ clusters are given in Fig. 8. Interaction energies for adsorbed CO₂ clusters are lower than those of N₂ clusters for all the three GYs, indicating that these GYs can be considered as potential candidates for the selective adsorption of CO₂ over N₂, enabling separation of CO₂ and N₂. Among the GYs, the magnitudes of interaction energies are higher for γ-GY-1 followed by γ-GY-2. A closer look at Fig. 8a and b reveals that the decrease in the magnitudes of interaction energies as we go from γ-GY-2 to γ-GY-4 is more pronounced for N₂ clusters than for CO₂ clusters. This is a consequence of the increase in pore size which reduces the possible interactions of N₂ molecules with the γ-GY-4 sheet. The trends obtained for the adsorption over GYs are related to the carbon densities of the sheet considered. The higher the carbon density of the sheet, the higher the magnitude of the interaction energy. In order to ascertain this, we further compared the energetics of adsorption of N₂ and CO₂ over GYs with that over graphene. For this, we modeled the adsorption of CO₂ unary clusters on circumcircumcoronene, where the interaction of CO₂ with circumcircumcoronene is described using the LJ potential.⁶⁰ The further details of the same are provided in the ESI.† The interaction energies associated with putative global minimum geometries obtained using PSO for CO₂ adsorption over circumcircumcoronene are also given in the ESI† (Fig. S12). The energetics and geometries of N₂ clusters adsorbed on circumcircumcoronene are already reported in our previous study.⁴² Since γ-GY-1 has the closest carbon density to graphene, we deemed it to be prudent to compare the interaction energies of the clusters adsorbed on graphene with those adsorbed on γ-GY-1. Fig. S12 of the ESI† clearly shows how close the interaction energies of N₂ clusters adsorbed on the graphene model system are with respect to γ-GY-1 for all cluster sizes. Graphene can indeed serve as a better adsorbent than γ-GY-1 for CO₂. Nevertheless, γ-GY-1 can be considered as a preferable alternative to graphene, since the selective separation of CO₂ from N₂ can be better achieved using γ-GY-1.

Fig. 8 A comparison of variation of interaction energies of (a) CO₂ clusters and (b) N₂ clusters adsorbed on γ-GY-1, γ-GY-2 and γ-GY-4 as a function of cluster size evaluated using PSO. Variation of second difference of interaction energies of (c) CO₂ clusters and (d) N₂ clusters adsorbed on γ-GY-1, γ-GY-2, and γ-GY-4 as a function of cluster size.

Then, in order to study the relative stabilities of the adsorbed N₂ and CO₂ clusters, we calculated the second difference of interaction energies. The plots of Δ₂E(n) for the adsorbed N₂ and CO₂ clusters are shown in Fig. 8(c and d). In our current study, the magic number clusters (maxima of the plot) in case of adsorbed N₂ are 4 and 7 for γ-GY-1, 5 and 7 for γ-GY-2 and 4 and 7 for γ-GY-4 whereas for adsorbed CO₂ clusters the magic numbers are 3, 5 and 8 for γ-GY-1, 3 and 7 for γ-GY-2 and 3, 5 and 8 for γ-GY-4.

We also analysed the components of the total interaction energy for a deeper understanding of the adsorption patterns. The total interaction energy comprises of contributions from the interactions within the gas molecules (cluster energy) and the interactions between the cluster and the molecular model of GY (cluster–GY energy). Fig. 9 depicts the contribution of each of these to the total interaction energies for the N₂ and CO₂ adsorption over GYs. For N₂ clusters, the contribution from cluster energy is so small that the total interaction energy is derived mainly from the cluster–GY interactions. The smaller slope in Fig. 8b for N₂ adsorption over γ-GY-4 when compared to γ-GY-1 and γ-GY-2 can also be attributed to the lower cluster-γ-GY-4 interaction energy contribution (Fig. 9b). While adsorbed N₂ clusters featured very less cluster energy, adsorbed CO₂ clusters feature a significant cluster energy contribution which is manifest in the form of clustering configurations of CO₂ observed in Fig. 7. In order to explicitly elucidate the same, we plotted the percentage relative contributions of cluster energy and cluster–GY energy to the total interaction energy in Fig. S13 of the ESI.† For adsorbed N₂ clusters, the cluster energy contributions are always less than 17% while in the case of adsorbed CO₂ clusters, a maximum contribution of up to 44% from the cluster energy was observed. We can also see that the increase in cluster energy contribution and the decrease in cluster–GY energy contribution with an increase in the cluster size is more pronounced in the case of CO₂ than N₂.

Fig. 9 Variation in the contributions of various terms to the total interaction energies for the adsorption of (a) CO₂ and (b) N₂ clusters on γ-GY-1, γ-GY-2, and γ-GY-4 as a function of cluster size.

All of the above analysis was performed for the putative global minimum geometries obtained using PSO. The reported energy minimum is the best value obtained from 25 independent PSO-L-BFGS runs. The choice of the number of independent runs was validated by comparing the putative global minimum energies obtained from 25 trial runs to those obtained after 100 trial runs for CO₂ and N₂ clusters of size 10 adsorbed on γ-GY-1, γ-GY-2, and γ-GY-4 (Fig. S14 of the ESI†). The interaction energies at minima predicted after 100 trials were in good agreement (deviations within 2 kcal mol⁻¹) with those obtained after 25 trials. As we can see, there are slight improvements in certain cases (for example, in the case of (CO₂)₁₀-γ-GY-1). However, achieving this slight improvement demands high computational resources. If we weigh the improvement in the results with respect to the computational resources required in running 100 vs. 25 trials, we note that the changes in geometries/energetics are not significant enough to warrant the use of these additional resources. Bearing these factors in mind, we concluded that 25 runs are adequate for reliably locating the structures of adsorbed CO₂ and N₂ clusters.

Next, we validated our choice of the annulenic model systems of γ-GY-1, γ-GY-2, and γ-GY-4 by comparing the interaction energies of molecular clusters of size 5 and 10 adsorbed on the currently employed model systems and a set of larger annulenic models of GYs, namely C₄₀₂H₄₂, C₅₈₂H₄₂, and C₄₅₀H₃₀. The geometries of the larger models of γ-GY-1, γ-GY-2, and γ-GY-4 are given in Fig. S15 of the ESI.† The interaction energies obtained for the adsorption of molecular clusters on the larger model systems showed good agreement with those obtained for our chosen model systems (Fig. S16 of the ESI†), implying the adequacy of our current model systems in representing GY sheets.

Global optimization of adsorbed binary clusters

Finally, we moved on to investigate the adsorption of binary mixtures of N₂ and CO₂ on γ-GYs. The putative global minimum geometries that we obtained for the binary mixtures of CO₂ and N₂ in various ratios are given in Fig. 10, 11 and Fig. S17 of the ESI.† In case of almost all the binary clusters considered, we observed that the N₂ molecules are scattered over different pore centres whereas CO₂ molecules prefer to cluster together as they did in their unary clusters but to a lesser extent. It is essentially the presence of the CO₂ molecules that scatter the N₂ molecules. The scattering of N₂ molecules is more apparent in the cases of adsorption on γ-GY-2 and γ-GY-4. The origin of these configurations can be explained with the help of the various contributing energy terms. For an adsorbed binary cluster, the total interaction energy includes energy contributions arising from: (i) the interactions between the CO₂ molecules and the model system of GY (CO₂–GY energy), (ii) the interactions between the N₂ molecules and the model system of GY (N₂–GY energy), (iii) the interactions among the CO₂ molecules (CO₂–CO₂ energy), (iv) the interactions among the N₂ molecules (N₂–N₂ energy), and (v) the interactions between the N₂ and the CO₂ molecules (N₂–CO₂ energy). Fig. 12 shows the total interaction energies and their components for all the adsorbed binary clusters we have considered in our study. An analysis of the interaction energy contributions for the mixtures with equal number of CO₂ and N₂ molecules revealed the strengths of the various interactions. Among the different terms, the CO₂–GY contribution is found to be the largest, which is followed by the N₂–GY term. The N₂–N₂ contribution is the lowest of them all. The CO₂–CO₂ interactions and CO₂–N₂ interactions compete with each other in order to arrive at the most favoured configuration. A similar analysis can be performed for the binary clusters of size 10 with unequal numbers of N₂ and CO₂ molecules. A crossover between the CO₂–GY and N₂–GY contributions can be observed in the plot at 6 : 4 (N₂ : CO₂) ratio (Fig. 12b). If we look at the leading terms of the total interaction energy for the 6 : 4 cluster, we can observe an interesting feature: even though there is a larger concentration of N₂ than CO₂, the contribution of CO₂–GY and N₂–GY interactions are almost the same. The contribution from N₂–N₂ interactions is always the lowest, even for mixtures with a larger concentration of N₂ molecules. For a deeper understanding, we have also plotted the variation of cluster energies and cluster–GY energies of CO₂ and N₂ in adsorbed binary mixtures (with equal number of N₂ and CO₂ molecules) with respect to adsorbed unary clusters in Fig. S18 of the ESI.† A relatively high reduction in cluster energies of N₂ molecules is observed upon forming binary mixtures, which is a direct consequence of the scattered arrangement of N₂ molecules seen in the putative global minimum geometries of the adsorbed binary clusters.

Fig. 10 Putative global minimum geometries (top and side views) of binary clusters of N₂ and CO₂ in the ratio n : m where n = m = 1–5 adsorbed on (a) γ-GY-1, (b) γ-GY-2 and (c) γ-GY-4 evaluated using PSO.

Fig. 11 Putative global minimum geometries (top and side views) of binary clusters of N₂ and CO₂ in the ratio 2 : 8, 4 : 6, 6 : 4, 8 : 2 adsorbed on (a) γ-GY-1, (b) γ-GY-2 and (c) γ-GY-4 evaluated using PSO.

Fig. 12 Variation in the contributions of various terms to the total interaction energies for the adsorption of binary clusters of CO₂ and N₂ adsorbed on γ-GY-1, γ-GY-2, and γ-GY-4 evaluated using PSO as a function of N₂ : CO₂ ratio, n : m where (a) n = m = 1–5 and (b) n + m = 10.

Conclusions

In conclusion, we have investigated three γ-GYs namely, γ-GY-1, γ-GY-2 and γ-GY-4, for the adsorption of unary and binary clusters of CO₂ and N₂ molecules up to a cluster size of 10. A swarm intelligence technique along with empirical force field formulations was implemented to obtain the putative global minimum geometries of the adsorbed clusters. The non-electrostatic interactions between CO₂–CO₂, N₂–N₂ and CO₂–N₂ were modeled using LJ, Buckingham and ILJ potentials, respectively, while the electrostatic interactions were described using the Coulombic potential. The interactions of the molecules with the molecular models of GYs were described using the ILJ potential. The putative global minimum geometries obtained for adsorbed unary clusters of CO₂ and N₂ indicated their higher affinity to γ-GY-1 followed by γ-GY-2. When adsorbed on γ-GY-1, CO₂ and N₂ molecules formed a monolayer at ∼3 Å from the sheet, while no such monolayer formation was observed in the case of γ-GY-2 and γ-GY-4. For γ-GY-1 and γ-GY-2, N₂ preferred adsorption over the pore centres, while CO₂ preferred clustering without any site preference. Clustering was prominently observed for CO₂ molecules due to the stronger CO₂–CO₂ interactions. The same was observed in case of adsorbed binary clusters, wherein, the clustering of CO₂ molecules scattered the N₂ molecules of the mixture. The predicted energetics demonstrate the high affinity of CO₂ in binding to all the three γ-GYs when compared to N₂ suggesting the higher adsorption selectivity of γ-GYs for CO₂ than for N₂. The ease of implementation of swarm intelligence technique along with the promising selectivity offered by GYs boosts our confidence towards designing efficient membranes for carbon capture in the near future.

Author contributions

The manuscript was written through contributions of all the authors. All authors have given approval to the final version of the manuscript.

Data availability

The data supporting this article have been included as part of the ESI.† The Python code executing global optimization of bare and adsorbed molecular cluster configurations using particle swarm optimization can be found at the GitHub repository https://github.com/SMMACG-IISER/PSO_N2_CO2_GY/tree/main.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge use of the Padmanabha cluster at the Centre for High-performance Computing at IISER TVM. R. S. S. acknowledges the Science and Engineering Research Board (SERB), Government of India for financial support of this work, through the SERB Core Research Grant (CRG/2022/006873). M. R. and C. J. thank IISER TVM for their fellowships.

References

1. UNFCCC, 1/Cp. 21, Adoption of the Paris Agreement, https://unfccc.int/resource/docs/2015/cop21/eng/l09r01.pdf.
2. Y. Gao, X. Gao and X. Zhang, Engineering, 2017, 3, 272–278.
3. J. G. Speight, Natural Gas (Second Edition), Gulf Professional Publishing, Boston, 2019, pp. 59–98.
4. M. R. M. Abu-Zahra, L. H. J. Schneiders, J. P. M. Niederer, P. H. M. Feron and G. F. Versteeg, Int. J. Greenhouse Gas Control, 2007, 1, 37–46.
5. P. J. E. Harlick and F. H. Tezel, Microporous Mesoporous Mater., 2004, 76, 71–79.
6. R. V. Siriwardane, M.-S. Shen, E. P. Fisher and J. Losch, Energy Fuels, 2005, 19, 1153–1159.
7. K. Sumida, D. L. Rogow, J. A. Mason, T. M. McDonald, E. D. Bloch, Z. R. Herm, T.-H. Bae and J. R. Long, Chem. Rev., 2012, 112, 724–781.
8. S. Cavenati, C. A. Grande and A. E. Rodrigues, J. Chem. Eng. Data, 2004, 49, 1095–1101.
9. Y. B. Apriliyanto, N. F. Lago, A. Lombardi, S. Evangelisti, M. Bartolomei, T. Leininger and F. Pirani, J. Phys. Chem. C, 2018, 122, 16195–16208.
10. R. H. Baughman, H. Eckhardt and M. Kertesz, J. Chem. Phys., 1987, 87, 6687–6699.
11. A. James, C. John, C. Owais, S. N. Myakala, S. Chandra Shekar, J. R. Choudhuri and R. S. Swathi, RSC Adv., 2018, 8, 22998–23018.
12. M. Bartolomei, E. Carmona-Novillo and G. Giorgi, Carbon, 2015, 95, 1076–1081.
13. Y. Jiao, A. Du, M. Hankel, Z. Zhu, V. Rudolph and S. C. Smith, Chem. Commun., 2011, 47, 11843–11845.
14. H. J. Kwon, Y. Kwon, T. Kim, Y. Jung, S. Lee, M. Cho and S. Kwon, AIP Adv., 2017, 7, 125013.
15. R. Majidi and U. Sarkar, Mol. Simul., 2020, 46, 1383–1389.
16. S. J. Mahdizadeh and E. K. Goharshadi, RSC Adv., 2020, 10, 24255–24264.
17. M. Bartolomei and G. Giorgi, ACS Appl. Mater. Interfaces, 2016, 8, 27996–28003.
18. S. Zhou, M. Wang, S. Wei, H. Xin, W. Zhai, S. Xu, S. Liu, S. Liu, Z. Wang, C.-M. Lawrence Wu and X. Lu, Mater. Today Phys., 2021, 21, 100539.
19. M. H. Darvishnejad and A. Reisi-Vanani, Int. J. Quantum Chem., 2020, 120, e26342.
20. M. H. Darvishnejad and A. Reisi-Vanani, Comput. Mater. Sci., 2020, 176, 109539.
21. C. He, R. Wang, D. Xiang, X. Li, L. Fu, Z. Jian, J. Huo and S. Li, Appl. Surf. Sci., 2020, 509, 145392.
22. Z. Tabandeh and A. Reisi-Vanani, Fuel, 2023, 333, 126295.
23. S. Xu, S. Wei, L. Wang, S. Liu, M. Wang, S. Liu, Z. Wang, T. Yang and X. Lu, Appl. Surf. Sci., 2022, 605, 154724.
24. G. Li, Y. Li, H. Liu, Y. Guo, Y. Li and D. Zhu, Chem. Commun., 2010, 46, 3256–3258.
25. Y. Hu, C. Wu, Q. Pan, Y. Jin, R. Lyu, V. Martinez, S. Huang, J. Wu, L. J. Wayment, N. A. Clark, M. B. Raschke, Y. Zhao and W. Zhang, Nat. Synth., 2022, 1, 449–454.
26. J. Gao, J. F. Li, Y. H. Chen, Z. C. Zuo, Y. J. Li, H. B. Liu and Y. L. Li, Nano Energy, 2018, 43, 192–199.
27. J. Kennedy and R. Eberhart, Particle Swarm Optimization, in Proceedings of ICNN'95 - International Conference on Neural Networks, 1995, vol. 4, pp. 1942–1948.
28. R. C. Eberhart and S. Yuhui, Particle Swarm Optimization: Developments, Applications and Resources, in Proceedings of the 2001 Congress on Evolutionary Computation (IEEE Cat. No. 01TH8546), 2001, vol. 1, pp. 81–86.
29. A. P. Engelbrecht, Computational Intelligence, John Wiley & Sons, Ltd, England, 2nd edn, 2007, pp. 289–358.
30. M. Clerc, Particle Swarm Optimization, ISTE Ltd, UK, 2006.
31. X.-S. Yang, Nature-Inspired Metaheuristic Algorithms, Luniver Press, Engelbrech, 2nd edn, 2008.
32. H. Iba, AI and Swarm: Evolutionary Approach to Emergent Intelligence, CRC Press, Boca Raton, 2019.
33. R. J. W. Hodgson, Partical Swarm Optimization Applied to the Atomic Cluster Optimization Problem, in Proceedings of the 4th Annual Conference on Genetic and Evolutionary Computation, 2002, pp. 68–73.
34. G. Jana, A. Mitra, S. Pan, S. Sural and P. K. Chattaraj, Front. Chem., 2019, 7, 485.
35. C. Owais, C. John and R. S. Swathi, Eur. Phys. J. D, 2021, 75, 16.
36. C. Owais, C. John and R. S. Swathi, Phys. Chem. Chem. Phys., 2020, 22, 20693–20703.
37. C. John, C. Owais, A. James and R. S. Swathi, J. Phys. Chem. C, 2021, 125, 2811–2823.
38. C. John, M. Rajeevan and R. S. Swathi, Chem. – Asian J., 2022, 17, e202200625.
39. S. T. Call, D. Y. Zubarev and A. I. Boldyrev, J. Comput. Chem., 2007, 28, 1177–1186.
40. L.-T. Shi, M. Tang, X.-R. Chen, C.-E. Hu and Y. Cheng, Comput. Theor. Chem., 2017, 1120, 102–111.
41. Y.-M. Wen, S.-K. Zhang, C.-E. Hu and Y. Cheng, Theor. Chem. Acc., 2019, 138, 83.
42. C. John and R. S. Swathi, J. Phys. Chem. A, 2023, 127, 4632–4642.
43. S. Zhou, M. Wang, S. Wei, S. Cao, Z. Wang, S. Liu, D. Sun and X. Lu, Mater. Today Phys., 2021, 16, 100301.
44. M. H. Darvishnejad and A. Reisi-Vanani, J. CO2 Util., 2021, 46, 101469.
45. L. Fu, R. Wang, C. Zhao, J. Huo, C. He, K.-H. Kim and W. Zhang, Chem. Eng. J., 2021, 414, 128857.
46. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg Williams, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, 2016.
47. J. Mawwa, S. U. D. Shamim, S. Khanom, M. K. Hossain and F. Ahmed, RSC Adv., 2021, 11, 32810–32823.
48. R. Amin, R. Ghiasi and A. Marjani, Russ. J. Phys. Chem. A, 2020, 94, 2760–2769.
49. P. K. Ramya and C. H. Suresh, J. Phys. Chem. A, 2023, 127, 2511–2522.
50. K.-L. Du and M. N. S. Swamy, Search and Optimization by Metaheuristics, Springer International Publishing, Switzerland, 2016.
51. D. C. Liu and J. Nocedal, Math. Program., 1989, 45, 503–528.
52. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa and P. van Mulbregt, Nat. Methods, 2020, 17, 261–272.
53. J. L. Llanio-Trujillo, J. M. C. Marques and F. B. Pereira, J. Phys. Chem. A, 2011, 115, 2130–2138.
54. H. R. Soni, J. Gebhardt and A. Görling, J. Phys. Chem. C, 2018, 122, 2761–2772.
55. M. Lafkioti, B. Krauss, T. Lohmann, U. Zschieschang, H. Klauk, K. v Klitzing and J. H. Smet, Nano Lett., 2010, 10, 1149–1153.
56. A. Chakradhar, N. Sivapragasam, M. T. Nayakasinghe and U. Burghaus, J. Vac. Sci. Technol., A, 2015, 34, 021402.
57. X. Gao, Y. Zhu, D. Yi, J. Zhou, S. Zhang, C. Yin, F. Ding, S. Zhang, X. Yi, J. Wang, L. Tong, Y. Han, Z. Liu and J. Zhang, Sci. Adv., 2018, 4, eaat6378.
58. J. Gao, J. Li, Y. Chen, Z. Zuo, Y. Li, H. Liu and Y. Li, Nano Energy, 2018, 43, 192–199.
59. J. J. Potoff and J. I. Siepmann, AIChE J., 2001, 47, 1676–1682.
60. G. M. Meconi and R. N. Zangi, Phys. Chem. Chem. Phys., 2020, 22, 21031–21041.
61. G. M. Meconi, R. Tomovska and R. Zangi, J. CO2 Util., 2019, 32, 92–105.
62. A. J. Stone, The Theory of Intermolecular Forces, Oxford University Press, Oxford, 2nd edn, 2013.
63. M. Bertolus, V. Brenner, P. Millie and J. B. Maillet, Z. Phys. D, 1997, 39, 239–250.
64. H.-J. Böhm and R. Ahlrichs, Mol. Phys., 1985, 55, 1159–1169.
65. J.-B. Maillet, A. Boutin, S. Buttefey, F. Calvo and A. H. Fuchs, J. Chem. Phys., 1998, 109, 329–337.
66. F. Calvo, G. Torchet and M.-F. de Feraudy, J. Chem. Phys., 1999, 111, 4650–4658.
67. F. Pirani, S. Brizi, L. F. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi, Phys. Chem. Chem. Phys., 2008, 10, 5489–5503.
68. M. J. Weida, J. M. Sperhac and D. J. Nesbitt, J. Chem. Phys., 1995, 103, 7685–7699.
69. H. Takeuchi, J. Phys. Chem. A, 2008, 112, 7492–7497.
70. E. García-Arroyo, M. I. Hernández, T. González-Lezana, J. H.-R. Campos-Martínez and J. Bretón, Theor. Chem. Acc., 2021, 140, 79.
71. B. Jeziorski, R. Moszynski and K. Szalewicz, Chem. Rev., 1994, 94, 1887–1930.
72. D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O'Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer, A. Y. Sokolov, K. Patkowski, A. E. DePrince, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford and C. D. Sherrill, J. Chem. Phys., 2020, 152, 184108.
73. N. Mehio, S. Dai and D. E. Jiang, J. Phys. Chem. A, 2014, 118, 1150–1154.
74. M. Miklitz and K. E. Jelfs, J. Chem. Inf. Model., 2018, 58, 2387–2391.
75. Z. S. Lu, P. Lv, D. W. Ma, X. W. Yang, S. Li and Z. X. Yang, J. Phys. D: Appl. Phys., 2018, 51, 065109.
76. M. H. Darvishnejad and A. Reisi-Vanani, Chem. Eng. J., 2021, 406, 126749.
77. L. Fang and Z. X. Cao, J. Phys. Chem. C, 2020, 124, 2712–2720.

---

Footnote

† Electronic supplementary information (ESI) available: Cartesian coordinates of the molecular models (text files), additional theoretical details, models, putative global minimum geometries and plots (PDF). See DOI: https://doi.org/10.1039/d4cp02843k

---