Leveraging metaheuristics to uncover water confinement in multilayer graphynes

Abstract

Global optimization is an effective approach to study the geometries and energetics of atomic or molecular confinement within nanostructures. The high computational cost associated with modeling such complex chemical systems calls for the adoption of stochastic global optimization techniques. Herein, we employ a swarm intelligence-based technique, namely, particle swarm optimization (PSO), to study the confinement of water clusters in monolayer and multilayer graphynes (GYs), including γ-GY-2, γ-GY-3, and γ-GY-4. The water molecules are described using the TIP4P model. The non-electrostatic part of GY–water and GY–GY interactions is modeled using the optimally fitted improved Lennard-Jones potential and the anisotropic Hod's interlayer potential, while the Coulombic potential is employed to account for the electrostatic interactions between GYs and water. Our PSO results reveal that the pore size of GYs is vital to the confinement of water clusters in multilayer γ-GYs. The γ-GY-2 multilayer tends to accommodate water as a monolayer between its two layers for large cluster sizes. A single-file confinement of water molecules is observed when water clusters were confined within the γ-GY-3 trilayer. In contrast, γ-GY-4, with the largest pore size, allowed clustering of water molecules within the triangular channels. Our findings established the importance of incorporating the twist features of GYs in the modeling formulation, as well as the accurate description of empirical formulations that can enable large-scale simulations. Our findings hold promise for extended research on water transport through twisted multilayer GYs.

---

Introduction

Properties of water under confinement have always been a hot topic in the scientific community. The structural features of confined water molecules differ from those of bulk water owing to the change in the hydrogen-bonding network upon confinement. This results in a change in their physicochemical properties from bulk water, spurring the incessant interest in confined water. For instance, the density and surface tension of water were found to be altered when confined in silica pores.¹ Similarly, the freezing point and melting point of water were found to be reduced on nano-scale confinement.² The influence of confinement on the dielectric constant of water has been discussed as well in the literature.³ Various experimental studies involving infrared spectroscopy have revealed the effect of water confinement on hydrogen bonding and dynamics of water.^4,5 A review by Bagchi and co-workers outlined the variation of dielectric relaxation and solvation dynamics of water molecules when confined within biological systems.⁶ Moreover, confinement brings in variation in the fluid mechanism when passing through nanoporous channels. In particular, water confined in carbon nanotubes exhibited enhanced flow compared to bulk water in macroscopic tubes.^7,8 However, many theoretical studies demonstrated that such a change in fluid behaviour, as well as the structural features of confined water molecules, is highly dependent on the characteristics of nanotubes.^9–12 For example, carbon nanotubes with a diameter of nearly 0.55 nm were experimentally demonstrated to have single-file water confinement (one-dimensional water chains).¹³ Water confinement in single-file arrangement has been of heightened interest as many of the characteristics of confined water are independent of the atomic structure of confining species, unlocking biomimetic and nanotechnological applications.¹⁴ The thermodynamic stability of single-file water confinement originates from the gain in rotational entropy of water molecules, which compensates for the increase in energy due to the reduction in the number of hydrogen bonds.¹⁵ For the case of single-file confinement in carbon nanotubes, water molecules possess strong anisotropic orientational relaxation, wherein slowing down of dipolar relaxation and speeding up of the HH vector (the vector joining the two hydrogen atoms in a molecule) relaxation were observed.¹⁶ Additionally, the water molecules exhibited angular jumps in which the two hydrogen atoms on a given water get interchanged preserving the H-bond with the same neighbour.¹⁷ Studies have shown that single-file water has high proton mobility along the water chain, and it happens via the Grotthuss mechanism.^18,19 Also, the length-dependent dielectric properties of the 1D water chains make them potential candidates for sensor applications.^20,21 As with the change in properties of water upon confinement, certain properties of the confining structures were also found to be influenced by the presence of water, especially when the confining medium is layered two-dimensional structures.²² In 2015, Wang et al. designed a seamless, thermally transparent, and electrically insulating interface with the intercalation of water between graphene and metal substrates.²³ Similarly, the presence of water has been shown to repress charge transfer in graphene/mica heterostructures.²⁴

The interaction of water molecules with polycyclic aromatic hydrocarbons (PAHs) and other model compounds of graphene is extensively studied owing to the hydrophobicity of the carbon materials. Many ab initio studies have been reported discussing the interaction of a single water molecule with PAHs or graphene.^25–30 Most of them theoretically analysed the energetics with respect to various orientations of the adsorbed water molecule and the dependence of the size of PAHs on the stability of the adsorbed water. Moreover, absolutely localized molecular orbital energy decomposition analysis (ALMO-EDA) predicted that the interactions significantly contributing to the stability of PAH–water complexes are dispersion and electrostatics.²⁷ Apart from these, various theoretical studies, including molecular dynamics, global optimization, and electronic structure calculations, have been employed to study the lowest energy geometries of water clusters adsorbed on PAHs and graphitic surfaces.^31–39 These investigations probed various features such as the energetics, thermodynamic properties, dynamics, charge transfer, and spectroscopy of the adsorbed water clusters,^34,36–39 and the influence of water molecules on the hydrophobicity and the electronic structure of graphene.^32,33,39 The electronic state of graphene was found to have a very minimal effect on water adsorption.^33,39 Similarly, researchers have also explored the structures and energetics of water clusters intercalated between PAHs and graphene layers.^40–43 Beyond the intercalation, the transport of water molecules through the slit pores of graphene bilayer has also been assessed via molecular dynamics.⁴⁴ The study has shown that the formation of hydrogen-bonded networks by water molecules between the graphene layers is hindered due to the high van der Waals interaction between the graphene layers. This issue no longer persisted when graphene oxide layers were considered.⁴⁴ Besides the slit pores, the permeation of water molecules through graphene layers is impeded due to the small pore size of hexagonal pores in graphene. However, incorporation of a nanopore in graphene overcomes this challenge. For a single layer of graphene, such a nanopore brings in the possibility of water desalination, which largely depends on the pore size and chemical substitution of the nanopore.^45,46 Similarly, nanoporous graphene multilayers exhibit water permeation; however, the water transport and desalination efficiency additionally depend on the interlayer separation and the pore alignment.^47,48

Building on the well-established results of water confinement within graphene/nanoporous graphene layers, graphynes (GYs), another class of hydrophobic carbon membranes, offer promising yet underexplored potential for water confinement. GYs are planar carbon sheets composed of sp² and sp carbon atoms.⁴⁹ Various classifications of GYs are made based on the arrangement of acetylenic linkages in the carbon network. The classification includes α-GY, β-GY, γ-GY, rhombic-GY, 6,6,12-GY, etc.⁵⁰ Each of these GYs can be again classified based on the number of acetylenic linkages (N) present between two sp² carbon atoms (denoted as GY-N). This structural variation in GYs enables the possibility of finetuning the porosity of GYs. The remarkable structural, electronic, thermal, and mechanical properties of GYs have resulted in several intriguing applications that have received a great deal of attention.^51–54 One among them is water desalination.^55,56 The majority of investigations on water desalination using GYs employ γ-GYs.^57–62 Within γ-GYs, the pore size of γ-GY-3 makes it most appropriate for salt rejection; meanwhile, considering both salt rejection and permeability, γ-GY-4 is established as a better purification membrane.^61,62 For a better understanding of the mechanism of water desalination, various studies have been performed analysing the confinement and transport of water through the layers of γ-GYs.^63–66 In a recent molecular dynamics study by Li et al., the authors reported the order of diffusion rates of water across γ-GYs as γ-GY-4 > γ-GY-5 > γ-GY-3, a trend which is non-linear with respect to the pore size.⁶⁵ In the same study, they mentioned the possibility of the single-file mode of transport through γ-GY-3 multilayers at a lower time scale and the Fickian mode of transport at a larger time scale. A similar observation was also reported for the case of fast transport of water through a unit pore of single layer γ-GY-3.⁶³ However, the Fickian mode of transport observed at a larger time scale might also be a consequence of the absence of twist in γ-GY-3 multilayers. It is, therefore, essential to assess the expected role of twisting of γ-GYs while studying water transport in multilayer γ-GYs. Moreover, molecular dynamics studies investigating water–GY interactions employed empirical formulations that are not parametrized for the systems under consideration. Thus, it is prudent to revisit the confinement of water in multilayer γ-GYs and to ascertain if single-file water arrangements can still be observed in multilayer γ-GY-3 after taking into account the twisting features.

In the present study, we employ a swarm intelligence-based technique, namely particle swarm optimization (PSO), to assess the single-file confinement of water in multilayer γ-GYs. We evaluate the structures and energetics of water molecules confined in monolayers, bilayers, and trilayers of γ-GY-2, γ-GY-3, and γ-GY-4. To model the water confinement in bilayers and trilayers more accurately, our modeling formalism included the twisting features of GY layers by way of an empirical anisotropic interlayer potential. The water–water interactions and GY–water interactions are also described using empirical formulations. The empirical formulations describing GY–water interactions and GY–GY interactions are parametrized against reference electronic structure calculations. Subsequently, the putative global minimum geometries of water molecules up to clusters of size 20 confined on (within) monolayer (multilayer) γ-GYs are predicted using PSO in conjunction with empirical formulations.

Methodology

Herein, we explore the water confinement in monolayers, bilayers, and trilayers of three different γ-GYs, namely γ-GY-2, γ-GY-3, and γ-GY-4. The annulenic model systems of all the three GYs are optimized using density functional theory (DFT) at the B3LYP/6-311G(d,p) level of theory incorporating D3 correction with Becke–Johnson damping (B3LYP-D3(BJ)) (Fig. 1).^67–69 The choice of the B3LYP-D3 functional for DFT calculations was based on the literature reports suggesting the reliability of this functional in modelling non-covalent interactions in closely related systems.^70–72 The geometry optimizations were performed using the Gaussian 16 suite of programs.⁷³ The acquired geometries were confirmed as actual minima by verifying the absence of imaginary frequencies and obtaining convergence for maximum force, RMS force, maximum displacement, and RMS displacement, as specified in the Gaussian suite of programs. We considered the confinement of 1–10, 1–12, and 10–20 water molecules, respectively, in monolayer, bilayer, and trilayer GYs. The number of molecules considered for confinement in the bilayer and trilayer GYs are chosen based on the estimated availability of pores for accommodating them. The putative global minimum geometries of various water–GY complexes are obtained by performing global optimization of the total intermolecular interaction energy using PSO. Note that, in the literature thus far, free water clusters in the size range of 1–37 are investigated using the state-of-the-art global optimization techniques.^74–78 In the current study, we consider up to a cluster size of 20 for intercalated clusters, given the additional computational costs associated with the incorporation of the carbon membranes for adsorption and intercalation into the modeling formalism. The total intermolecular interaction energy includes contributions from water–water interactions, GY–water interactions, and GY–GY interactions. The details regarding the employed empirical formulations, their parametrization, and the PSO algorithm are provided in the following subsections:

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

Water–water interactions

The interactions among water molecules are described using the well-known four-site transferable intermolecular potential (TIP4P) model.⁷⁹ The model, as mentioned, represents a water molecule with four sites, three atomic sites of water molecules, and a fictitious M-site. The M-site is a charged fictitious site placed along the bisecting plane of the HOH angle. In the TIP4P model, the intermolecular interaction between two water molecules is a sum of the Lennard-Jones (LJ) interaction between two O atoms and the Coulombic interaction between the rest of the three sites (two H atoms and M-site) of the interacting water molecules. For an n-molecule cluster of water, the TIP4P interaction energy is given by

indicates the LJ potential between the oxygen atoms of k^th and l^th molecules. V_el(r_ik,jl) refers to the Coulombic interaction between the i^th and j^th sites of the k^th and l^th molecules, respectively. The LJ potential⁸⁰ and the Coulombic potential take the forms

where r_ij is the distance between the two interacting sites i and j. ε and σ are the LJ parameters representing the well depth and the distance at which the pair potential is zero. q_i indicates the charge of the i^th site, and C_el is a constant that takes the value of 331.934 kcal Å mol⁻¹ e⁻². The structural and potential parameters of the TIP4P model are adopted from ref. 79 and tabulated in Table S1 of the ESI.†

GY–water interactions

The interaction between a water molecule and a GY model compound is described in terms of a sum of electrostatic and non-electrostatic components as given by

The non-electrostatic interactions are described using an isotropic potential, namely improved Lennard-Jones (ILJ) potential.⁸¹ As the name suggests, the ILJ potential is an improved version of the LJ potential with an additional parameter accounting for the better description of long-range and short-range interactions. For two atoms, i and j, separated by a distance of r_ij, the ILJ potential is expressed as

The first and the second terms of the potential describe the short-range repulsion and the long-range attraction, respectively. The potential description involves three parameters, ε, r_m, and β, which represent the well depth, equilibrium distance, and a factor governing the hardness of the interaction between the two atoms i and j. In the present study, we consider water as a pseudoatom positioned at the oxygen atom that interacts with the atoms of the GY model systems, making the total non-electrostatic interaction energy as

where p and m indicate the number of carbon and hydrogen atoms in the model system of GY, and n denotes the number of water molecules. The same description was previously employed by Bartolomei et al. for describing the non-electrostatic interactions of water molecules with γ-GY-2.⁶⁶

For the electrostatic part of the interaction energy, we employ the Coulombic potential wherein the partial charges residing on atoms of GYs are evaluated using the Merz–Singh–Kollman (MK) scheme.⁸² The partial charges on water are considered the same as in the TIP4P model. Thus, the total electrostatic component of the interaction energy is defined as

here, w indicates the total number of atoms in a model system of GY (including carbon atoms and hydrogen atoms), n refers to the number of water molecules, and l denotes the two H atoms and the M-site of a water molecule. q_w and q_lk represent the partial charges of the w^th atom of GY and the l^th site of the k^th water molecule, respectively.

GY–GY interactions

The GY–GY interactions in the present study are modeled using an anisotropic interlayer potential proposed by Hod and co-workers.⁸³ Herein, we denote the potential as H–ILP. Though the H–ILP was first proposed for layered hexagonal boron nitride, the parameters for graphitic systems were subsequently developed.^84–86 Utilising the H–ILP parameters of the graphitic systems as the initial parameters, we have recently reparametrized the H–ILP for twisted bilayers of various GY-1 systems against dispersion-corrected DFT calculations.⁸⁷ Owing to its anisotropic nature, the potential described the interlayer features with reasonable accuracy, superseding isotropic potentials like the ILJ potential. The H–ILP, in its original form, has three terms: (i) long-range attraction, (ii) short-range repulsion, and (iii) electrostatic interaction. As the carbon atoms in GYs do not possess sizable effective charges, we excluded the electrostatic term from our H–ILP description, making the current H–ILP formulation as given below:

V_H–ILP(r_ij,ρ_ij) = V_rep(r_ij,ρ_ij) + V_att(r_ij) (9)

The repulsive and attractive terms are given by

here, r_ij is the distance between atoms i and j from two different layers of GYs. For the same pair of atoms, ρ_ij denotes the lateral distance between them, which can be evaluated using the equations ρ_ij² = r_ij² − (n_i·r_ij)² and ρ_ji² = r_ij² − (n_j·r_ij)². n_i and n_j represent the normal vectors to the plane where the respective atoms and their nearest neighbours reside. The ρ_ij brings in the anisotropy in the repulsive component of H–ILP. The repulsive and attractive contributions include a long-range taper cut-off function given bywhich limits the interactions to within an interatomic distance of R_cut. In the present study, R_cut was assigned as 16.0 Å.⁸⁷ In the repulsive part of the potential, the parameters ε_ij and C_ij scale the isotropic and the anisotropic parts of the repulsion, respectively, α_ij sets the slope of the isotropic repulsion wall, and the parameters β_ij and γ_ij determine the range of the respective interactions. The attractive term incorporates a damped LJ-type attraction adopted from the Tkatchenko–Scheffler dispersion correction scheme.⁸⁸d_ij and S_R,ij are the unitless parameters governing the slope and onset of the damping function with numerical values 15.0 and 1.0, respectively.⁸⁷C_6,ij and r^eff_ij are the dispersion coefficient and the effective atomic radius of the two interacting atoms.

In the current study, we consider only the interactions between the carbon atoms of the two layers owing to the minimal influence of the hydrogen atoms on the total interaction energy between two GY layers in a bilayer, as observed in our previous study.⁸⁷ For interaction between two GY layers with p carbons each, the total interaction energy is calculated using the equation

Parametrization of the empirical potentials

For an accurate description of the energetics of water confinement, we parametrize the ILJ potential for water–GY interactions and the H–ILP for GY–GY interactions. Since the water–GY interaction includes electrostatic interaction as well, the same was also included along with ILJ potential during parametrization. The parametrization of a potential is performed by minimizing the root mean square deviation (RMSD) of the interaction energy profile obtained using the empirical potential (E_int) against the reference data (E_ref). The geometry optimizations and other reference electronic structure calculations reported here are carried out using the Gaussian 16 suite of programs.⁷³ The RMSD is expressed aswhere M represents the number of data points in the interaction energy profile. The minimization is performed in Python using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm.⁸⁹ During the process of fitting, whenever necessary, we have added a weightage to the minima by multiplying e^−αE_ref,i with the squared error. α is a constant whose value is varied between 0.01 and 0.5 as per the requirement. We initially parametrize the ILJ and the H–ILP potentials for the case of γ-GY-3 and assess the transferability for γ-GY-2 and γ-GY-4. In cases where the parameters are not transferable, they are further finetuned. The details of the reference data that have been employed for the ILJ and the H–ILP parametrization are given below:

ILJ potential. For parametrizing the ILJ potential, we considered the MP2C data reported by Bartolomei et al. as the reference data.⁶⁶ The data correspond to the vertical motion of the water molecule towards the center of the triangular pore of the single-pore models of γ-GYs. However, the study reported the MP2C data for the interaction of water molecules with γ-GY-1, γ-GY-2, and γ-GY-3. Hence, in order to check the transferability for the water-γ-GY-4 system, we obtained similar data at the B3LYP-D3(BJ)/6-311G(d,p) level of theory with the incorporation of counterpoise correction.⁹⁰ The single-pore models of γ-GYs utilized for empirical calculations are also optimized at the same level of theory and are provided in Fig. S1 of the ESI.†

H–ILP. We employed the twisting profile of the γ-GY-3 bilayer kept at a stacking distance of 3.6 Å as the reference data for parametrizing H–ILP. The twisting profiles were generated using dispersion-corrected DFT at the B3LYP-D3(BJ)/6-311G(d,p) level of theory (incorporating counterpoise correction) by varying the twist angle from 0° to 30°. The twisting profiles of the γ-GY-2 and γ-GY-4 bilayers, as well as the stacking profiles of all the three γ-GY bilayers (by varying the vertical distance between the two GYs for AA stacking) were assessed for the transferability of the obtained potential parameters.

PSO: algorithm and implementation

PSO is a swarm intelligence-based metaheuristic algorithm developed by Kennedy and Eberhart in 1995.⁹¹ The simplicity and effectiveness of the algorithm made PSO the most cited (as per Google Scholar) metaheuristic algorithm for the period 2000–2022.⁹² The algorithm is initiated by a random distribution of a group of particles (swarm) with zero velocity across the desired search space. Basically, at a given iteration, a particle is defined by two factors, namely position and velocity. Position represents a vector with the variables that need to be optimized as the elements, and thus, the number of variables gives the dimension of the position vector. Hence, each particle represents a point in the parameter space. The velocity of the particle is also a vector with the same dimension as the position and represents the direction of motion of the particle. The velocity of the particle, i, at an iteration, t + 1, is determined based on three pieces of information – the previous velocity (v^t_i) and the best solution that the particle (P^t_best,i) and the entire swarm (G^t_best) have achieved until then, expressed as

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)] (15)

c₁ and c₂ are the acceleration constants governing the contributions of individual and group memory to the velocity. r₁ and r₂ are random numbers generated within the range [0,1], thereby incorporating stochasticity in the algorithm. The inclusion of the constriction factor, χ, is an improvisation to the basic PSO algorithm to enhance the performance.⁹³ The numerical value of χ (taken as 0.729) is dependent on the numerical values of c₁ (2.05) and c₂ (2.05).⁹⁴ Knowing the velocity, the position of the particle at (t + 1)^th iteration is calculated as

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

Thus, as the iterations pass, the position and the velocity of each particle get updated, and the entire swarm converges to the global minimum solution. As mentioned, the global best solution at each iteration is evaluated considering the entire swarm, thereby employing a star topology wherein all particles are connected to each other.⁹⁵ Additionally, we have incorporated dynamic velocity clamping to improve the exploration as performed in our previous study.⁹⁶

In the current study, the objective function to be minimized using PSO is the total intermolecular interaction energy of the systems considered, including water–water, GY–water, and GY–GY interactions. For the adsorption of water on monolayer GYs and, intercalation within bilayer and trilayer GYs, the following forms of intermolecular interaction energies are employed:

E_{monolayer–H2O} = E_GY–H2O + E_H2O–H2O (17)

E_{bilayer–H2O} = E_GY–GY + E_GY–H2O(layer 1) + E_GY–H2O(layer 2) + E_H2O–H2O (18)

E_{trilayer–H2O} = E_GY–GY(layer 1, layer 2) + E_GY–GY(layer 1, layer 3) + E_GY–GY(layer 2, layer 3) + E_GY–H2O(layer 1) + E_GY–H2O(layer 2) + E_GY–H2O(layer 3) + E_H2O–H2O (19)

We employ the rigid body approximation for water.⁹⁷ Within the approximation, we consider the water molecules to be rigid bodies with only rotational and translational degrees of freedom. Hence, a water molecule is defined by its centre of mass coordinates (r_COM) and Euler angles (θ, φ, and ϕ) describing its orientation. If r_rel,i is the Cartesian coordinate of each atom in a water molecule with respect to its centre of mass, the Cartesian coordinate of the same atom in the search space is evaluated as

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

here, R⁻¹(θ, φ, ϕ) is the inverse of the rotation matrix. Thus, for the adsorption of n water molecules on monolayer GY, the dimension of the problem is 6n as the position coordinates of the atoms of GYs are not varied. For the case of intercalation of water within bilayer and trilayer GYs, the twist angles between the layers and the interlayer distances are varied as well, making the dimensions of the problems (6n + 2) and (6n + 4), respectively. Note that, for intercalation within trilayer GYs, variations of interlayer distances and twist angles of second and third layers with respect to the first layer contribute to the four dimensions in addition to the degrees of freedom of water.

In our PSO implementation, the number of total iterations performed is 1000, except for adsorption of water molecules on the monolayer GYs wherein we performed 2000 iterations as well. For the adsorption of water molecules on the monolayer GYs, we performed PSO calculations with 500 particles, while for the bilayer and trilayer GYs, we considered 200–300 particles. To ensure the optimality of the solution, we performed a local optimization using the L-BFGS algorithm at the end of each PSO run. Owing to the stochastic nature of the algorithm, we performed 25 trial runs of PSO for every system. For better convergence and efficiency, it is necessary to define the search space effectively. Hence, we provided bounds for each of the variables in the position vector. The rotation matrix and the bounds for Euler angles used are the same as in our previous study.⁹⁶ Owing to the symmetry of GYs, the interlayer twist angles are varied in the range 0° to 30°. We employed various limits for the interlayer distances, with a lower bound of 2.5 or 3 Å and an upper bound of 3.5 or 4 or 7 Å, depending on the GY sheet under consideration. For the trilayer GY, the lower and upper bounds for the interlayer distance between the first and the third layers were chosen as 5 or 6 or 8 Å and 7 or 8 or 12 Å, respectively. Whenever multiple limits have been used to analyse a single case, the best result obtained is reported. The bounds for the centre of mass coordinates of water molecules for each case are intuitively chosen so as to correspond to the adsorption or intercalation within the GY framework. This requires the introduction of a constraint in our PSO implementation by way of limiting the Z-coordinate of the centre of mass of water molecules in each iteration. However, in the local optimization step that is implemented after each PSO run, we perform the optimization with and without this constraint for the Z-coordinate of the centre of mass of water (hereafter mentioned as unconstrained local optimization), leading to two sets of results for each local optimization. Therefore, for 25 trials of PSO, we obtain 50 solutions upon local optimization. Out of the 50 solutions, the undesired configurations in which the water molecules move out of the GY framework are eliminated, and the best solutions among the rest are reported. All the PSO calculations are performed using in-house developed Python codes.

Results and discussion

The present study aims to obtain the minimum energy configurations of water-confined monolayers, bilayers, and trilayers of various γ-GYs, including γ-GY-2, γ-GY-3, and γ-GY-4. The objective is achieved in two steps: (i) parametrizing the empirical formulations describing the intervening intermolecular interactions against reference electronic structure data, and (ii) employing the optimized empirical potentials in PSO to obtain the putative global minimum geometries of multilayer GYs confining water molecules. The findings from both the steps are provided in the following subsections:

Parametrization of empirical potentials

ILJ potential. As mentioned in the Methodology section, we initially parametrized the ILJ potential for the interaction of a water molecule with γ-GY-3 using the MP2C data of Bartolomei et al.⁶⁶ The newly parametrized ILJ potential showed good correspondence to the MP2C data with an RMSD of 0.0483 kcal mol⁻¹. The obtained potential could describe the water-γ-GY-4 interactions accurately, establishing the transferability of the parameters for higher members of the γ-GY family. However, considerable deviation was observed for the water-γ-GY-2 interaction. Hence, we finetuned the parameters to arrive at a new set of parameters that could describe the water-γ-GY-2 interaction effectively. The obtained potential energy profiles are provided in Fig. 2 and the corresponding ILJ parameters are tabulated in Table 1. Additionally, the figure reveals that γ-GY-3 is better than the other GYs in stabilising water as it has the largest binding strength. The presence of an energy barrier in the case of γ-GY-2 indicates an impedance to water permeation through the triangular pore of γ-GY-2. Interestingly, a study by Bartolomei et al. has shown a reduction in the permeation barrier for a water molecule when another water molecule is placed at the other end of the γ-GY-2 pore, making the permeation more feasible.⁶⁶ The reduction in the binding strength for water at the center of the pore towards γ-GYs beyond γ-GY-3 can be attributed to their large pore size.

Fig. 2 Interaction energy profiles evaluated using the reference electronic structure method and fitted intermolecular analytic potential with electrostatic and non-electrostatic components (ILJ potential) for the interaction of a single water molecule with various GYs.

Table 1 Numerical values of optimal parameters of the ILJ potential and H–ILP

Interaction	GY	ε (kcal mol^−1)	r m (Å)	β
ILJ potential
C–H2O	γ-GY-2	0.233	3.508	8.423
	γ-GY-N (N = 3, 4)	0.200	3.632	8.500
H–H2O	γ-GY-N (N = 2, 3, 4)	0.150	3.500	8.500

---

Interaction	GY	α	β (Å)	γ (Å)	ε (kcal mol^−1)	C (kcal mol^−1)	r ^eff (Å)	C 6 (kcal Å^6 mol^−1)
H–ILP
C–C	γ-GY-N (N = 2, 3, 4)	8.863	2.717	0.921	0.150	3.000	3.519	559.395

---

Moreover, we assessed the transferability of the parametrized ILJ potential in describing the lateral motion of a water molecule over the GY models (Fig. 1). The corresponding DFT (B3LYP-D3(BJ)/6-311G(d,p)) data are generated by fixing the position of a water molecule at a distance of 3.5 Å above the GY plane and then varying its position from the center of the hexagonal ring along a direction bisecting the triangular pore. The resulting potential energy profiles demonstrated an overall reasonable agreement and a good agreement around the pore region with the corresponding DFT profiles. However, we reparametrized the ILJ potential (denoted as ILJ-I) by additionally including the lateral potential energy profiles in the reference data set. The obtained potential energy profiles are shown in Fig. S2 (ESI†), and the reparametrized ILJ parameters are tabulated in Table S2 (ESI†). Though the ILJ-I potential described the lateral motion better than the ILJ potential, we observed a slight decrease in the accuracy of the ILJ-I potential in describing vertical motion. As the primary objective of our study is to model water confinement in multilayer GYs involving the water molecules in the pore regions, we consider the first set of parameters (ILJ) for further analysis.

H–ILP. In order to capture the twisting features of bilayer GYs, it is essential to incorporate anisotropy in the intermolecular potential. Hence, we resort to the H–ILP for describing GY–GY interactions. The parametrization of H–ILP was initially performed to capture the twist features of the bilayer γ-GY-3, which yielded a reasonable fit with an RMSD of 1.239 kcal mol⁻¹. The obtained parameters (see Table 1) were found to be transferable for describing interlayer interactions pertaining to both twisting and stacking features of bilayers of γ-GY-2 and γ-GY-4, as well as the stacking feature of bilayer γ-GY-3. The potential energy profiles corresponding to the twisting and stacking motion of bilayers of γ-GY-N (N = 2–4) obtained using DFT and H–ILP are shown in Fig. 3. As expected, the bilayer binding strengths follow the order: γ-GY-2 < γ-GY-3 < γ-GY-4. The parametrized H–ILP predicted the nature of twisting and stacking profiles of all the γ-GYs with reasonable accuracy. However, it is necessary to ensure that the deviation in the energetics of twisting profiles does not arise from the lack of electrostatics in the intermolecular description. Thus, we parametrized H–ILP along with electrostatics (using point charge description) against DFT for the case of the twisting profile of the bilayer of the γ-GY-3, and subsequently, assessed the transferability of the parameters for the bilayers of γ-GY-2 and γ-GY-4, as well as for the stacking feature of bilayer γ-GY-3. The electrostatic interactions are accounted for by the Coulombic potential, where the partial charges residing on various atoms are evaluated using the MK scheme⁸² at the B3LYP-D3(BJ)/6-311G(d,p) level of theory. The obtained parameters and the corresponding interaction energy profiles are provided in Table S3 and Fig. S3 of the ESI.† The nature of twisting profiles remained the same even after the inclusion of electrostatics, and there is no significant improvement in the fits of stacking and twisting profiles upon the incorporation of electrostatics. Hence, throughout the remainder of the current study, we used only H–ILP (without electrostatics) to describe interlayer interactions in bilayer and trilayer GYs.

Fig. 3 The (a) twisting and (b) stacking interaction energy profiles evaluated using the reference electronic structure method and fitted H–ILP for bilayer γ-GYs.

Water confinement in GY layers

Now that we have optimized the empirical formulations for modeling intermolecular interactions of interest, we employed PSO to obtain the putative global minimum configurations of water molecules adsorbed on monolayer GYs, and intercalated into bilayer and trilayer GYs. The putative global minimum geometries for the adsorption of 1–10 water molecules are given in Fig. S4–S6 of the ESI.† A few representative geometries are provided in Fig. 4. A single water molecule occupied the position above the triangular pore center of γ-GY-2 with the hydrogen atoms pointing towards the pore. While, in the case of the adsorption of a water molecule on γ-GY-3 and γ-GY-4, the molecule gets adsorbed over the hexagonal pore center with the oxygen atom facing the pore, unlike the case of the water–graphene complex in which both the hydrogen atoms of water face the hexagonal pore upon adsorption.²⁹ This difference in the occupied positions mainly arises from the difference in pore sizes of the GYs. The pore size of γ-GY-2 being small, the water molecule tends to occupy positions above triangular pores to maximize the interactions, whilst for γ-GY-3 and γ-GY-4, the maximum binding energies are attained above hexagonal pore centers. Meanwhile, the orientation of the water molecule occupied above triangular or hexagonal pores is governed by the effective partial charges of the atoms constituting the pores. The electrostatic potential energy surfaces (see Fig. S7, ESI†) indicate that the hexagonal pore center exhibits an effective positive charge, while the triangular pore center displays an effective negative charge. Consequently, when a water molecule is positioned above these pore centers, it orients in such a way that the positively charged hydrogen atoms face the triangular pore center, while the negatively charged M-site faces the hexagonal pore center. The water molecules of cluster sizes 2–7 exhibited similar geometries upon adsorption on all the three GY monolayers. Clusters of size 3–5 featured ring-shaped/regular polygonal geometries with H-bonds constituting the sides of the polygon. The predicted geometries are similar to those reported for water confinement on various PAHs, including coronene.^35,36 For adsorbed (H₂O)₆, we obtained a fused ring structure with two 4-membered rings (flat book structure), similar to the case of a stable isomer for (H₂O)₆ adsorbed on coronene arrived at using a density functional tight binding (DFTB) approach.³⁶ The (H₂O)₇ cluster displayed a fused ring structure with 5 and 4 membered rings when adsorbed on all the three GYs. Global optimization studies using basin hopping have reported the same geometries of water clusters for (H₂O)₆ and (H₂O)₇ adsorbed on circumcoronene.⁹⁸ For water clusters of sizes 8–10, adsorption on γ-GY-2, γ-GY-3, and γ-GY-4 yielded geometries possessing a 3D network of water molecules rather than nearly planar configurations. Such 3D adsorbed configurations for larger water clusters have been previously reported for water binding on nanographene models.^36,38

Fig. 4 Putative global minimum geometries (top and side views) obtained using PSO for the adsorption of 3, 5, 7, and 9 water molecules on monolayers of (a) γ-GY-2, (b) γ-GY-3, and (c) γ-GY-4.

As noted, the adsorption features of water molecules on the three γ-GYs do not differ significantly as the geometries of adsorbed water molecules are similar. However, the intercalation features of water within bilayers and trilayers of γ-GYs differ significantly. The variation in pore sizes of γ-GY bilayers resulted in interesting geometries of water under confinement. The putative global minimum geometries for the intercalation of 1–12 water molecules within bilayer γ-GYs are given in Fig. S8–S10 of the ESI.† A few representative geometries are provided in Fig. 5. For intercalation of water into bilayer γ-GY-2, until a cluster size of 8, the bilayer accommodated water molecules as a single water molecule or as a dimer within the triangular channels (including interlayer regions and the triangular pores) of two sheets maintaining the optimal interlayer distance of 3.2–3.4 Å. For higher-order clusters ((H₂O)_n, n = 9–12) the molecules preferred to stay together by forming a layer within the interlayer region of bilayer γ-GY-2, and hence the interlayer distance between the two GY layers increased to nearly 6 Å. For the bilayer γ-GY-3, the relatively large pore size allowed the easy confinement of water dimers within the triangular channels without any considerable alteration in the optimal interlayer distance. Such a behaviour suggested the possibility of single-file confinement of water molecules within multilayer γ-GY-3. On the other hand, a further increase in the pore size in the case of γ-GY-4 resulted in the entrapment of water as small clusters with 3–4 water molecules within the triangular channels of the bilayer γ-GY-4. The interlayer distances and the twist angles of GY bilayers under water confinement are tabulated in Table 2. γ-GY-3 and γ-GY-4 bilayers with confined water molecules exhibited an interlayer separation of 3.1–3.2 Å, which is lower when compared to that of γ-GY-2 (3.2–5.9 Å). In all the cases, the GY layers undergo a twist, establishing the importance of incorporating the twist features of multilayers in intercalation studies. Moreover, the obtained twist angles for the putative global minima are closer to the minima observed in the bilayer twisting energy profiles (Fig. 3), demonstrating the significance of accurate parametrization of the H–ILP in describing the intercalation features. This results in similar structural features for some of the water-intercalated bilayer GYs. The bilayer GYs with similar structural features are indicated by similar cell colors for the interlayer distances and twist angles in Table 2. The twist angles and interlayer distances of GY bilayers with deviations within 0.50° and 0.11 Å are color-coded the same.

Fig. 5 Putative global minimum geometries (top and side views) obtained using PSO for the intercalation of 3, 6, 9, and 12 water molecules within bilayers of (a) γ-GY-2, (b) γ-GY-3, and (c) γ-GY-4.

Table 2 The interlayer distances and twist angles of water-intercalated GY bilayers in the putative global minimum configurations obtained using PSO. Similarity in cell colors for the interlayer distances and twist angles indicates similarity in the structural features of bilayer GYs

---

Subsequently, we analysed the putative global minimum geometries of water intercalated into trilayer γ-GYs to have a comprehensive understanding of water confinement in multilayer γ-GYs. Additionally, we increased the number of water molecules under study to a range of 10–20 molecules. The putative global minimum geometries for the intercalation of 10–20 water molecules within trilayer γ-GYs are given in Fig. S11–S13 of the ESI.† A few representative geometries are provided in Fig. 6. Similar to our observation for the intercalation within bilayer γ-GY-2, in all the cases except for clusters of size 10–12, the water molecules formed a layer between two γ-GY-2 layers in the trilayer γ-GY-2. For cluster sizes 10–12, we observed dimers and single-file trimers of water molecules within the triangular channels in the trilayer γ-GY-2. However, as the number of water molecules increases, the single-file dimers and trimers no longer exist. The smaller pore size of γ-GY-2 restricts any further addition of water molecules in such single-file configurations. This could be a consequence of emerging repulsive interactions within the system as the cluster size increases. To overcome this, water molecules form a layer between two γ-GY-2 layers, thereby maximizing water–water interactions. On the other hand, the geometries obtained for the intercalation within trilayer γ-GY-3 exhibited single-file water confinement within all the six triangular channels. The single-file arrangement of water observed during the intercalation of 16 water molecules within trilayer γ-GY-3 is shown in Fig. 7. The geometries for the intercalation of 10–20 water molecules within trilayer γ-GY-3 are provided in Fig. S14 of the ESI.† In a majority of these cases, the water chains followed either a nearly linear or zigzag orientation to accommodate the interlayer twist observed in the trilayer γ-GY-3. For clusters of size 18–20, we observed the clustering of water molecules in one or two triangular channels, which is a consequence of constraining the molecules between the first and third layers. This is evident from the observation of a more stable single-file configuration of these clusters when unconstrained local optimization is implemented. The single-file configurations here are facilitated by the water molecules moving beyond the GY framework. Such putative global minimum geometries and corresponding interaction energies of 19 and 20 water clusters confined in the trilayer γ-GY-3 are provided in Fig. S15 of the ESI.† This suggests the possibility of finetuning the length of water chain by including more layers of γ-GY-3. Additionally, unlike the case of single-file motion of water in carbon nanotubes, the confinement in multilayer GYs has an added advantage of accommodating multiple single-file water chains in their triangular porous frameworks. Meanwhile, the γ-GY-4 trilayer accommodates water molecules as small clusters, the largest being (H₂O)₆, within the triangular channels. The interlayer distances and twist angles estimated for trilayers of γ-GYs with water molecules entrapped are provided in Table 3. The obtained distances and angles provide insights into the stacking of GY layers when water molecules are confined. Interestingly, we observed that some of the trilayer GYs adopt similar structures for the confinement of clusters of different sizes (see the color scheme in Table 3). The twist angles and interlayer distances of GY trilayers with deviations within 0.50° and 0.11 Å are color-coded the same. A possible explanation to this lies in the notable contribution of interlayer interactions resulting in the favourable stacking of GY layers. These stacking configurations do not directly relate to the number of molecules under confinement. For instance, the γ-GY-2 trilayer shows similar stacking configurations for the intercalation of water clusters of sizes 10 and 11, 13 and 17, 15 and 19, as well as 18 and 20. Moreover, for the case of γ-GY-3 trilayers, the obtained configurations result in different kinds of water chains. The γ-GY-3 trilayer exhibited three unique stacking configurations for the intercalation of (10, 12, 16), (11, 14, 17), as well as (13, 15, 18) water clusters. Here, the first and the third configurations exhibited nearly linear geometries for three-membered water chains confined inside them, while the second configuration resulted in a zigzag water chain. This is attributed to the fact that, in the second configuration, the twist with respect to the first layer displayed by the third layer is less than that of the second layer. For single-file water confinement in carbon nanotubes, it was reported that the OH bonds involved in the hydrogen bonds are nearly parallel with the nanotube axis.⁹⁹ Whilst in our case, most of the OH bonds are aligned with the vector joining two pore centers of the triangular channel within which the water molecule resides. Since the first and the third stacking configurations mentioned above have these pore-to-pore vectors connected almost linearly, the water chains confined within bear resemblance to the ones confined within the carbon nanotubes. However, because the confinement occurs in a different environment with a distinct pore design, the dynamics of the single-file water chain are expected to differ from that of carbon nanotubes, a topic worth exploring in the future. Additionally, we compared the water occupancies in carbon nanotubes and multilayer γ-GY-3. A 13.4 Å long nanotube of diameter 8.1 Å can accommodate five water molecules on an average, making it an available length of 2.68 Å for a single water molecule.⁹⁹ In parallel, the γ-GY-3 trilayers with an average separation of 6.4 Å between the first and the third layers can host three water molecules resulting in an available length of 2.13 Å per water, establishing that γ-GY-3 trilayers can accommodate water molecules more efficiently than carbon nanotubes. The interlayer twist features in multilayer GYs enable more facile passage pathways for water.

Fig. 6 Putative global minimum geometries (top and side views) obtained using PSO for the intercalation of 10, 13, 17, and 20 water molecules within trilayers of (a) γ-GY-2, (b) γ-GY-3, and (c) γ-GY-4.

Fig. 7 The geometries of 16 water molecules (top and side views) confined within trilayer γ-GY-3 (after removing the GY framework) as seen in the putative global minimum geometries obtained using PSO.

Table 3 The interlayer distances and twist angles of water-intercalated trilayer GYs in the putative global minimum configurations obtained using PSO. The interlayer distances provided refer to the distances between the first and the second layers, and between the second and the third layers. The twist angle provided refer to the twist of the second and the third layers with respect to the first layer. Similarity in cell colors for the interlayer distances and twist angles indicates similarity in the structural features of trilayer GYs

---

In order to ensure the universality of the presence of nearly linear and zigzag water chains intercalated in trilayer γ-GY-3, we performed PSO calculations for the intercalation of 10–17 water molecules within trilayer γ-GY-3 by employing the extended simple point charge (SPC/E) model for water.¹⁰⁰ Prior to PSO calculations, we verified the transferability of the parametrized ILJ potential in describing GY–water non-electrostatic interactions with water described using the SPC/E model (see Fig. S16, ESI†). The obtained global minimum geometries also suggested the single-file conformations of intercalated water molecules and the zigzag or nearly linear orientations of the resulting water chains (Fig. S17, ESI†). However, with the SPC/E model, we obtained zigzag water chains for most of the cases, while with the TIP4P model, it was mostly near-linear configurations. Apart from this, we also performed PSO calculations for the intercalation of water clusters of sizes 12 and 16 within trilayers of γ-GY-2, γ-GY-3, and γ-GY-4, with the water–GY non-electrostatic interactions being described by the ILJ-I potential. We obtained similar putative minimum geometries (Fig. S18, ESI†) as those obtained using the ILJ potential, suggesting the reliability of our findings. While the ILJ potential captures the key features of water–GY interactions across the three GYs studied, further improvement could be achieved by incorporating an anisotropic interfacial potential such as the one recently employed by Ouyang and co-workers in describing the interaction of water with graphene and hexagonal boron nitride surfaces.^101,102 Such a potential may also provide an advantage of having a single parameter set to represent water–GY interactions across all the three GY systems investigated in our study. This may serve as a potential avenue for future investigations.

Next, we analysed the energetics of adsorption of water on monolayer γ-GYs and intercalation within bilayer and trilayer γ-GYs. The variation in total intermolecular interaction energy for adsorption and intercalation as a function of cluster size of water is shown in Fig. 8. The interaction energies for the adsorption of water on all three monolayer GYs are nearly the same, with adsorption on γ-GY-2 giving rise to slightly more favourable intermolecular interactions. This is a consequence of the higher carbon density of the γ-GY-2 membrane. On the other hand, for intercalation within bilayers and trilayers, the total interaction energies show a distinct trend: γ-GY-4 > γ-GY-3 > γ-GY-2. The interaction energies substantially increase from adsorption on monolayers to intercalation within trilayers owing to the significant increase in the number of intervening interactions. For a better description, we analysed various intermolecular interaction energy components that contribute to the total intermolecular interaction energy. These components include: (i) cluster energy, the interaction energy among water molecules modeled using the TIP4P model, (ii) GY-cluster energy, the interaction energy between water molecules and GY layers described by the ILJ potential and the Coulombic potential, and (iii) GY–GY energy, which sums the total interaction between the GY layers described by the H–ILP. The variation in total intermolecular interaction energies and the components for adsorption and intercalation as a function of cluster size of water are shown Fig. 9. For adsorption of water molecules over monolayer γ-GYs, the leading contribution is from the water–water interactions. The contribution of water–GY interactions is significantly low for all the γ-GYs. This dominance of cluster energy over GY-cluster energy results in the 3D network structures of larger adsorbed water clusters. Similar to the total intermolecular interaction energies, both cluster energies and cluster-GY energies are nearly the same for adsorption of water on all three monolayer GYs. The GY–GY interactions make a dominant contribution for intercalation involving bilayers and trilayers of γ-GY-3 and γ-GY-4. Moreover, this contribution is nearly the same irrespective of the number of water molecules. For intercalation within bilayer and trilayer γ-GY-3, except for a few large clusters intercalated within trilayers, the next leading contribution is from the GY-cluster interactions, which is marginally higher than cluster energies. The deviation for the large clusters is correlated to the putative global minimum geometries described earlier wherein large clusters resulted in aggregation of water in the triangular channels instead of a single-file arrangement of water seen as in the case of mid-sized clusters. On the other hand, for bilayers and trilayers of γ-GY-4, the cluster energies were found to be higher than the cluster-GY energies, as is also evident from their confined geometries as described earlier. A different scenario arises for the confinement of water within trilayer and bilayer γ-GY-2. After a cluster size of 8 in the γ-GY-2 bilayer and a cluster size of 12 in the γ-GY-2 trilayer, a crossover in the dominance of the intermolecular components is observed between the cluster energies and GY–GY energies, which is a consequence of the formation of a water layer within the interlayer region of bilayer and trilayer γ-GY-2.

Fig. 8 A comparison of variation in total intermolecular interaction energies evaluated using PSO for the adsorption (intercalation) of water on (within) (a) monolayers, (b) bilayers, and (c) trilayers of γ-GY-2, γ-GY-3, and γ-GY-4 as a function of cluster size.

Fig. 9 Variation in the contributions of various terms to the total intermolecular interaction energies for the adsorption (intercalation) of water on (within) (a) monolayers, (b) bilayers, and (c) trilayers of γ-GY-2, γ-GY-3, and γ-GY-4 as a function of cluster size.

Conclusions

The current manuscript reports the structures and energetics of adsorbed and intercalated water clusters in monolayers, bilayers, and trilayers of a series of GYs, namely γ-GY-2, γ-GY-3, and γ-GY-4, with an explicit incorporation of the twist features. The putative global minimum geometries of confined water clusters were obtained by minimizing the total intermolecular interaction energies using a metaheuristic global optimization technique, namely PSO. The water molecules were modeled using the TIP4P model and the interactions between GY layers were described using an anisotropic interlayer potential that accounts for the twist features in multilayer GYs. The electrostatic and non-electrostatic interactions between water and GYs were expressed using the Coulombic and ILJ potentials, respectively. A better description of interactions was ensured by proposing potential parameters for the specified interactions by parametrizing the potentials against electronic structure data. The monolayer confinement of water molecules resulted in similar adsorbed geometries of clusters for up to cluster size 7, and the associated binding energies exhibited no significant variation across various GYs. However, the water molecules confined in multilayer γ-GYs displayed distinct confined configurations depending on the pore size of GYs. At higher water occupancies, bilayer and trilayer γ-GY-2 tend to squeeze in water as a layer between the GY layers. Meanwhile, the multilayer γ-GY-3 enabled single-file confinement of water molecules, with the length tuneable by increasing the number of layers. The large pore size of the multilayer γ-GY-4 resulted in the clustering of water molecules within the triangular channels. The single-file confinement of water obtained for transport through multilayer γ-GY-3 opens up the possibility of further fundamental exploration of its dynamics and potential applications. We believe that the methods employed are state-of-the-art for the problem in hand, and the findings are important for the general chemistry community as many experimental groups are exploring the studied carbon membranes for water permeation and desalination applications.

Author contributions

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

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the ESI.†

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. thanks IISER TVM for her fellowship.

References

1. T. Takei, K. Mukasa, M. Kofuji, M. Fuji, T. Watanabe, M. Chikazawa and T. Kanazawa, Colloid Polym. Sci., 2000, 278, 475–480.
2. A. W. Knight, N. G. Kalugin, E. Coker and A. G. Ilgen, Sci. Rep., 2019, 9, 8246.
3. S. Senapati and A. Chandra, J. Phys. Chem. B, 2001, 105, 5106–5109.
4. J. B. Brubach, A. Mermet, A. Filabozzi, A. Gerschel, D. Lairez, M. P. Krafft and P. Roy, J. Phys. Chem. B, 2001, 105, 430–435.
5. S. Le Caër, S. Pin, S. Esnouf, Q. Raffy, J. P. Renault, J. B. Brubach, G. Creff and P. Roy, Phys. Chem. Chem. Phys., 2011, 13, 17658–17666.
6. N. Nandi, K. Bhattacharyya and B. Bagchi, Chem. Rev., 2000, 100, 2013–2046.
7. M. Majumder, N. Chopra, R. Andrews and B. J. Hinds, Nature, 2005, 438, 44.
8. M. Whitby and N. Quirke, Nat. Nanotechnol., 2007, 2, 87–94.
9. M. Shaat, U. Javed and S. Faroughi, Sci. Rep., 2020, 10, 17167.
10. A. Alexiadis and S. Kassinos, Chem. Rev., 2008, 108, 5014–5034.
11. J. Hernández-Rojas, F. Calvo, J. Bretón and J. M. Gomez Llorente, J. Phys. Chem. C, 2012, 116, 17019–17028.
12. M. K. Tripathy, D. K. Mahawar and K. R. S. Chandrakumar, J. Chem. Sci., 2019, 132, 7.
13. S. Cambré, B. Schoeters, S. Luyckx, E. Goovaerts and W. Wenseleers, Phys. Rev. Lett., 2010, 104, 207401.
14. J. Köfinger, G. Hummer and C. Dellago, Phys. Chem. Chem. Phys., 2011, 13, 15403–15417.
15. H. Kumar, B. Mukherjee, S.-T. Lin, C. Dasgupta, A. K. Sood and P. K. Maiti, J. Chem. Phys., 2011, 134, 124105.
16. B. Mukherjee, P. K. Maiti, C. Dasgupta and A. K. Sood, ACS Nano, 2008, 2, 1189–1196.
17. B. Mukherjee, P. K. Maiti, C. Dasgupta and A. K. Sood, J. Phys. Chem. B, 2009, 113, 10322–10330.
18. C. Dellago, M. M. Naor and G. Hummer, Phys. Rev. Lett., 2003, 90, 105902.
19. R. Pomès and B. Roux, Biophys. J., 2002, 82, 2304–2316.
20. J. Köfinger and C. Dellago, Phys. Rev. Lett., 2009, 103, 080601.
21. J. Köfinger and C. Dellago, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 205416.
22. P. Bampoulis, K. Sotthewes, E. Dollekamp and B. Poelsema, Surf. Sci. Rep., 2018, 73, 233–264.
23. Y. Wang and Z. Xu, ACS Appl. Mater. Interfaces, 2016, 8, 1970–1976.
24. J. Shim, C. H. Lui, T. Y. Ko, Y.-J. Yu, P. Kim, T. F. Heinz and S. Ryu, Nano Lett., 2012, 12, 648–654.
25. E. Voloshina, D. Usvyat, M. Schütz, Y. Dedkov and B. Paulus, Phys. Chem. Chem. Phys., 2011, 13, 12041–12047.
26. S. Abe, Y. Nagoya, F. Watari and H. Tachikawa, Jpn. J. Appl. Phys., 2010, 49, 01AH07.
27. A. O. Ajala, V. Voora, N. Mardirossian, F. Furche and F. Paesani, J. Chem. Theory Comput., 2019, 15, 2359–2374.
28. D. Ž. Veljković, J. Mol. Graphics Modell., 2018, 80, 121–125.
29. E. M. Cabaleiro-Lago, J. A. Carrazana-García and J. Rodríguez-Otero, J. Chem. Phys., 2009, 130, 234307.
30. J. G. Brandenburg, A. Zen, M. Fitzner, B. Ramberger, G. Kresse, T. Tsatsoulis, A. Grüneis, A. Michaelides and D. Alfè, J. Phys. Chem. Lett., 2019, 10, 358–368.
31. C. Melios, C. E. Giusca, V. Panchal and O. Kazakova, 2D Mater., 2018, 5, 022001.
32. A. Akaishi, T. Yonemaru and J. Nakamura, ACS Omega, 2017, 2, 2184–2190.
33. O. Leenaerts, B. Partoens and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 235440.
34. L. F. L. Oliveira, J. Cuny, M. Morinière, L. Dontot, A. Simon, F. Spiegelman and M. Rapacioli, Phys. Chem. Chem. Phys., 2015, 17, 17079–17089.
35. E. Michoulier, N. Ben Amor, M. Rapacioli, J. A. Noble, J. Mascetti, C. Toubin and A. Simon, Phys. Chem. Chem. Phys., 2018, 20, 11941–11953.
36. A. Simon and F. Spiegelman, J. Chem. Phys., 2013, 138, 194309.
37. A. Simon, M. Rapacioli, J. Mascetti and F. Spiegelman, Phys. Chem. Chem. Phys., 2012, 14, 6771–6786.
38. B. S. González, J. Hernández-Rojas, J. Bretón and J. M. Gomez Llorente, J. Phys. Chem. C, 2007, 111, 14862–14869.
39. R. R. Q. Freitas, R. Rivelino, F. de Brito Mota and C. M. C. de Castilho, J. Phys. Chem. A, 2011, 115, 12348–12356.
40. K. Rohini, D. M. R. Sylvinson and R. S. Swathi, J. Phys. Chem. A, 2015, 119, 10935–10945.
41. S. McKenzie and H. C. Kang, Phys. Chem. Chem. Phys., 2014, 16, 26004–26015.
42. A. Zamir, E. R. Molina, M. Ahmed and T. Stein, Phys. Chem. Chem. Phys., 2022, 24, 28788–28793.
43. S. K. Kim, W. Chen, S. Pourasad and K. S. Kim, J. Phys. Chem. C, 2016, 120, 19212–19224.
44. N. Raghav, S. Chakraborty and P. K. Maiti, Phys. Chem. Chem. Phys., 2015, 17, 20557–20562.
45. D. Cohen-Tanugi and J. C. Grossman, Nano Lett., 2012, 12, 3602–3608.
46. D. Cohen-Tanugi and J. C. Grossman, Desalination, 2015, 366, 59–70.
47. D. Cohen-Tanugi, L.-C. Lin and J. C. Grossman, Nano Lett., 2016, 16, 1027–1033.
48. J. Muscatello, F. Jaeger, O. K. Matar and E. A. Müller, ACS Appl. Mater. Interfaces, 2016, 8, 12330–12336.
49. R. H. Baughman, H. Eckhardt and M. Kertesz, J. Chem. Phys., 1987, 87, 6687–6699.
50. B. G. Kim and H. J. Choi, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 115435.
51. 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.
52. G. Narang, D. Bansal, S. Joarder, P. Singh, L. Kumar, V. Mishra, S. Singh, K. Tumba and K. Kumari, FlatChem, 2023, 40, 100517.
53. J. Kang, Z. Wei and J. Li, ACS Appl. Mater. Interfaces, 2019, 11, 2692–2706.
54. T. He, Y. Kong, A. R. P. Santiago, M. A. Ahsan, H. Pan and A. Du, Mater. Chem. Front., 2021, 5, 6392–6412.
55. K. H. Lasisi, O. K. Abass, K. Zhang, T. F. Ajibade, F. O. Ajibade, J. O. Ojediran, E. S. Okonofua, J. R. Adewumi and P. D. Ibikunle, Front. Chem., 2023, 11 DOI:10.3389/fchem.2023.1125625.
56. H. Qiu, M. Xue, C. Shen, Z. Zhang and W. Guo, Adv. Mater., 2019, 31, 1803772.
57. S. Lin and M. J. Buehler, Nanoscale, 2013, 5, 11801–11807.
58. M. Akhavan, J. Schofield and S. Jalili, Phys. Chem. Chem. Phys., 2018, 20, 13607–13615.
59. M. Mehrdad and A. Moosavi, Polymer, 2019, 175, 310–319.
60. M. Xue, H. Qiu and W. Guo, Nanotechnology, 2013, 24, 505720.
61. J. Kou, X. Zhou, H. Lu, F. Wu and J. Fan, Nanoscale, 2014, 6, 1865–1870.
62. C. Zhu, H. Li, X. C. Zeng, E. G. Wang and S. Meng, Sci. Rep., 2013, 3, 3163.
63. J. Kou, X. Zhou, Y. Chen, H. Lu, F. Wu and J. Fan, J. Chem. Phys., 2013, 139, 064705.
64. S. Bagheri, A. Shameli, M. Darvishi and G. Fakhrpour, Physica E, 2017, 90, 123–130.
65. L. Li, F. Fang, J. Li, G. Zhou and Z. Yang, Phys. Chem. Chem. Phys., 2022, 24, 2534–2542.
66. M. Bartolomei, E. Carmona-Novillo, M. I. Hernández, J. Campos-Martínez, F. Pirani, G. Giorgi and K. Yamashita, J. Phys. Chem. Lett., 2014, 5, 751–755.
67. J. Tirado-Rives and W. L. Jorgensen, J. Chem. Theory Comput., 2008, 4, 297–306.
68. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
69. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
70. D. G. A. Smith and K. Patkowski, J. Chem. Theory Comput., 2013, 9, 370–389.
71. K. Ohno, H. Satoh and T. Iwamoto, Chem. Phys. Lett., 2019, 716, 147–154.
72. B. Brauer, M. K. Kesharwani, S. Kozuch and J. M. L. Martin, Phys. Chem. Chem. Phys., 2016, 18, 20905–20925.
73. 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, D. Williams-Young, 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.
74. D. J. Wales and M. P. Hodges, Chem. Phys. Lett., 1998, 286, 65–72.
75. H. Kabrede and R. Hentschke, J. Phys. Chem. B, 2003, 107, 3914–3920.
76. H. Takeuchi, J. Chem. Inf. Model., 2008, 48, 2226–2233.
77. S. Kazachenko and A. J. Thakkar, Chem. Phys. Lett., 2009, 476, 120–124.
78. B. Hartke, Phys. Chem. Chem. Phys., 2003, 5, 275–284.
79. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
80. J. E. Jones and S. Chapman, Proc. R. Soc. London, Ser. A, 1997, 106, 463–477.
81. F. Pirani, S. Brizi, L. F. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi, Phys. Chem. Chem. Phys., 2008, 10, 5489–5503.
82. U. C. Singh and P. A. Kollman, J. Comput. Chem., 1984, 5, 129–145.
83. I. Leven, I. Azuri, L. Kronik and O. Hod, J. Chem. Phys., 2014, 140, 104106.
84. T. Maaravi, I. Leven, I. Azuri, L. Kronik and O. Hod, J. Phys. Chem. C, 2017, 121, 22826–22835.
85. I. Leven, T. Maaravi, I. Azuri, L. Kronik and O. Hod, J. Chem. Theory Comput., 2016, 12, 2896–2905.
86. W. Ouyang, I. Azuri, D. Mandelli, A. Tkatchenko, L. Kronik, M. Urbakh and O. Hod, J. Chem. Theory Comput., 2020, 16, 666–676.
87. A. Melekamburath, A. James, M. Rajeevan, C. John and R. S. Swathi, Phys. Chem. Chem. Phys., 2021, 23, 27031–27041.
88. A. Tkatchenko, R. A. DiStasio, R. Car and M. Scheffler, Phys. Rev. Lett., 2012, 108, 236402.
89. R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, SIAM J. Sci. Comput., 1995, 16, 1190–1208.
90. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
91. J. Kennedy and R. Eberhart, Particle Swarm Optimization, Proceedings of ICNN'95 – International Conference on Neural Networks, 1995, vol. 4, pp. 1942–1948.
92. K. Rajwar, K. Deep and S. Das, Artif. Intell. Rev., 2023, 56, 13187–13257.
93. K.-L. Du and M. N. S. Swamy, Search and Optimization by Metaheuristics, Springer International Publishing, Switzerland, 2016.
94. R. C. Eberhart and S. Yuhui, Particle Swarm Optimization: Developments, Applications and Resources, Proceedings of the 2001 Congress on Evolutionary Computation (IEEE Cat. No. 01TH8546), 2001, vol. 1, pp. 81–86.
95. A. P. Engelbrecht, in Computational Intelligence, John Wiley & Sons, Ltd, England, 2nd edn, 2007, pp. 289–358.
96. M. Rajeevan, C. John and R. S. Swathi, Phys. Chem. Chem. Phys., 2024, 26, 23152–23167.
97. J. L. Llanio-Trujillo, J. M. C. Marques and F. B. Pereira, J. Phys. Chem. A, 2011, 115, 2130–2138.
98. B. S. González, J. Hernández-Rojas, J. Bretón and J. M. Gomez Llorente, J. Phys. Chem. C, 2008, 112, 16497–16504.
99. G. Hummer, J. C. Rasaiah and J. P. Noworyta, Nature, 2001, 414, 188–190.
100. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
101. Z. Feng, Y. Yao, J. Liu, B. Wu, Z. Liu and W. Ouyang, J. Phys. Chem. C, 2023, 127, 8704–8713.
102. Z. Feng, Z. Lei, Y. Yao, J. Liu, B. Wu and W. Ouyang, Langmuir, 2023, 39, 18198–18207.

---

Footnote

† Electronic supplementary information (ESI) available: Potential parameters, models, putative global minimum geometries and plots (PDF). See DOI: https://doi.org/10.1039/d5cp02507a

---