On the D₂/H₂ adsorption selectivity of sodium decorated polycyclic aromatic hydrocarbons. I. Accurate modeling approaches and low pressure selectivity

Abstract

In the search for an economical approach to separate deuterium and protium isotopes, we have theoretically investigated the equilibrium adsorption of these diatoms onto Na-decorated polycyclic aromatic hydrocarbons: protonated coronene and, as a graphene prototype, circumcircumcoronene. We have developed accurate interaction potentials between H₂ and these substrates (within pseudo-atom or rigid-body approximations) and have solved the Schrödinger equation for the adsorbed H₂ and D₂ in order to compute the corresponding partition functions. In this way, adsorption isotherms and D₂/H₂ selectivity have been estimated providing a complete characterization of the low pressure behavior. Our results show that the pseudo-atom model is insufficient to correctly predict the mentioned properties and that neglecting the contribution of thermally excited states provides a lower bound of the exact selectivity. As for quantitative aspects, both materials have been found to afford a deuterium selectivity well above 6, the required value for industrial applications. Also, separation of spin isomers of H₂ is shown to be effective. In this respect, these substrates appear to be suitable candidates as chromatographic stationary phases, warranting additional studies on their separation capabilities at higher pressure. The theoretical framework developed for such a task and the results thus obtained are presented in a second article in this series.

---

1 Introduction

Heavy hydrogen isotopes represent a valuable commodity for today applications in science and technology as well as potential fuels for controlled fusion power plants.^1–3 As such, it is a pitiable situation that both deuterium and tritium only represent minuscule fractions of the naturally occurring Z = 1 species. One thus often faces the task of separating D and T atoms from H ones starting from extremely diluted mixtures. Besides, many physico-chemical properties of the isotopologues are quite similar, a fact that makes the separation processes quite cumbersome. In the end, one must rely on a limited difference in the quantum-mechanical nature of their motion, a fact that, in turn, imposes the usage of low temperatures to maximize the impact of the latter.

One possible way of exploiting the limited differences in properties between isotopomers is chemical affinity quantum sieving,^4–11 where we can distinguish two limiting cases, chemisorption or physisorption based on differences in the interaction strengths. Obviously, there may be advantages and disadvantages in both cases. For instance, higher temperatures may be used when chemisorption takes place, as the desorption requires a higher amount of energy (e.g. see the cases of environmentally harmless and easily recycled nanoscopic alkali-earth oxides¹²). If so, less work is required to reduce the system temperature in order to make the separation processes effective. Conversely, physisorption leaves intact the chemical nature of the adsorbed species, so that surface catalyzed isotopic scrambling could be avoided. In this respect, a goal worth pursuing would be to increase the interaction strength between active sites and physisorbing species, as this should also induce larger differences in their adsorption energies. Attempts in this direction have already been made^13–17 by decorating carbonaceous materials with heteroatoms to improve hydrogen storage, exploiting substrates such as graphene^18–20 and its derivatives,^21,22 or carbon nanotubes (CNTs).^6,23–27 This approach is even more attractive by the possibility of using vastly available metals such as alkali (Na and K) or alkali-earths (Mg and Ca), for which theoretical studies have suggested potentially useful performances.^16,20,28 In this arena, we have previously compared the hydrogen adsorption capability of coronene and sodium decorated coronene and the corresponding stabilities of the clusters formed inside superfluid helium droplets, evidencing a substantial increase in the uptake for the decorated system.¹⁷

Albeit the adsorption energy on decorated carbonaceous materials may still not be sufficiently large to allow H₂ storage for automotive applications, we must stress that it may suffice to differentiate between D₂ and H₂ adsorption. This capability may, de facto, emerge as a consequence of relatively small energy differences provided the operating temperature is sufficiently low; thus, even a small difference (e.g. ΔH_ads/k_B ⋍ 30 K over amorphous ice²⁹) in adsorption enthalpies may be sufficient to generate a measurable selectivity for the heavier isotopologue at 20 K, as it is substantially defined by the exponential function exp[ΔH_ads/(k_BT)].³⁰

In this work, we explore the effectiveness in D₂/H₂ separation that decorated carbonaceous supports may afford. We do so by building sufficiently accurate force fields to describe the pairwise interaction energy between X₂ (X representing H or D) and the carbonaceous materials, describing X₂ either as a pseudo-atom (PsAt) or as a rigid-body (RB). In addition, the change in chemical potential is estimated as a function of the temperature at low pressure via numerical eigenstate determination for the adsorbed species, comparing exact results with low temperature approximations for the same quantities.³⁰ The protocol is applied to progressively larger polycyclic aromatic hydrocarbons (PAHs) decorated with Na atoms: protonated coronene ([H–Cor–Na]⁺) and circumcircumcoronene ([CCCor–Na]). On the one hand, [H–Cor–Na]⁺ is interesting since it is one PAH for which there are accurate mass spectrometry measurements of the attachment of hydrogen in comparison with theory¹⁷ (the substrate is modeled with an extra proton due to its abundance in the experiments). On the other hand, computations of H₂ adsorption on neutral [CCCor–Na] have been recently accomplished²⁸ to provide a more accurate model of an extended support such as graphene, so here we aim to find out the dependence of the D₂/H₂ selectivities with the size and charge of the substrate.

To summarize the main contents of this paper (part I of a two-parts work), we mention that

(1) we have built a force field for the X₂/[H–Cor–Na]⁺ interactions employing the RB approach;

(2) using this force field and previously developed ones for [H–Cor–Na]⁺ (PsAt-type) and [CCCor–Na] (PsAt and RB), we have solved the Schrödinger equation for the vibrational motion of a single adsorbed species on the decorated materials;

(3) we have predicted the adsorption isotherm^30,31 and the low pressure selectivity toward D₂ molecules by estimating the partition functions for adsorbed H₂ and D₂, also assessing the extent of some approximations such as low-temperature partition functions and the PsAt approach;

(4) the predicted RB selectivities are well above the industrially required standard (ref. 32) and outperform a few of the most promising material candidates for an effective and economical separation;^4,33

(5) data produced in this investigation are also used to explore the possibility of separating the different spin isomers of either H₂ or D₂.³⁴

All results concerning operando conditions are, instead, described in a subsequent part of this investigation (part II), where the assumption that a single molecule is adsorbed per metal site (or that a few non-interacting molecules adsorb on the site) is removed. In this way, effects on the selectivity due to an increase in molecular crowding around the metal will be analyzed. To that aim, a diffusion Monte Carlo (DMC) framework will be applied to relate the energetics of adsorbed aggregates composed of D₂ molecules and H₂ ones with the equilibrium gas phase mixture. Some of the conclusions reached in the present work will be useful to assess the impact of some limitations of the DMC method, specifically, that it only provides the ground state of a system.

The outline of the manuscript is as follows. In Section 2, we derive the calculation of adsorption isotherms and adsorption selectivity (2.1) and describe the interaction potential model (2.2) as well as the numerical approach to compute vibrational eigenstates (2.3). Results on the adsorption of D₂ and H₂ as well as on the D₂/H₂ separation are presented in Section 3, first for the [H–Cor–Na]⁺ substrate (3.1) and following with the [CCCor–Na] one (3.2). Finally, Section 4 reports our conclusions on the suitability of the investigated materials for the chosen task, together with additional discussions relevant for future developments.

2 Models and methods

2.1 Theoretical approach

Assuming independent behavior for surface adsorption sites binding at most a single molecule, then the statistical mechanics approach commonly employed to derive Langmuir-type isotherms can be followed³¹ provided the appropriate corrections are introduced to account for the description of gas phase molecular rotation when necessary.³⁰ Discarding the chance that strongly cooperative behavior between adsorbed species is present, a similar analysis could be used also in case of multi-molecule adsorbing sites when the pressure is low, even though several species may be simultaneously present.

We begin recalling that estimating adsorption isotherms for a diatomic species X₂ requires two pieces of information: first, the definition of their gas phase chemical potential (μ_X2), which would depend on the composition of the gaseous phase, x_X2, and its total pressure, p_tot, so that the partial pressure is p_X2 = x_X2p_tot; second, an estimate of the chemical potential for the adsorbed species , which would depend on the fractional amount of occupied sites . In the case of ideal gases composed of linear species, the analytical expression ³⁰ can be used. Here, is the standard chemical potential at temperature T for a point-like particle with a mass equal to X₂, while q_r(X₂) is the rotational partition function of X₂. Equating the last expression with allows one to estimate the fractional occupation of the sites for each pure component when the total pressure is low. Thus, a straightforward modification of the demonstration provided by Hill³¹ to include the description of rotational motion leads to the Langmuir-like expression:

Here,

with , and q_ads(T,X₂) being the partition function for the adsorbed linear rotor. This is, as usual, written as “sum over states” for the adsorbed species (β ≡ (k_BT)⁻¹):

Here, E_l are the energies of the quantum librational states of the adsorbed molecules.

As for the adsorption selectivity between two species, this can be defined in terms of the relative amount of adsorbed species with respect to the ones present at equilibrium in the gas phase. Being the chemical potential a function of intensive composition variables, it may be more convenient to use the molar fraction of each adsorbed species (i.e. ), or the fraction of occupied adsorbent sites to define the selectivity which, in the case involving D₂ and H₂, can be defined as:^7,35

Substituting eqn (1) into eqn (4) when p_X2χ_r(T,X₂) ≪ 1, one obtains

where ΔE₀ = E₀(D₂) − E₀(H₂) is the difference in ground state (l = 0) energy for the two adsorbed species, and ΔE_(l,0)(X₂) = E_l(X₂) − E₀(X₂), the energy gap between the ground state and the l-th excited state. Note that the condition p_X2χ_r(T,X₂) ≪ 1 sets an upper limit to the pressures, i.e., pressures should be lower than χ_r(T,X₂)⁻¹, the half coverage pressure. For example, for para D₂, these limiting pressures are about 10⁻²⁶ and 10⁻¹² bar at 20 and 40 K, respectively.

Notice that, differently from the case of CNT selectivity,^25,36,37 all translational degrees of freedom of molecules are quantized when a molecule is adsorbed, as the interaction potential with the site restrains the center of mass motion along all of them. From this fact, it descends that the factor, related to the motion of the free molecules, remains unmodified. Moreover, all ΔE_(l,0)(X₂) values are positive defined, so that the contribution due to terms with large l decreases upon increasing l. In other words, it may be possible to neglect terms in the sums when ΔE_(l,0) ≫ k_BT. Alternatively, one may altogether neglect the summations ratio if the numerator and denominator grow at a similar pace upon increasing l, so that the overall value does not deviate substantially from unity. In fact, one would expect the exact estimate for the ratio to be just above one, as the density of states for the adsorbed D₂ molecules ought to be higher due to the potentially smaller energy gaps. If that is the case, neglecting the ratio would produce a lower bound for the selectivity, the tightness of which may somewhat depend on the specific system under investigation.

The approach just briefly summarized may be extended to the case of a simultaneous co-adsorption of more than two species, if they can be considered independent in both phases. This is, indeed, important if one wishes to estimate selectivity toward “normal” D₂ and H₂ gases, for which conversion between the nuclear spin states is hindered by the lack of an interaction source coupling them. As pointed out in ref. 38, the appropriate approach to deal with non-convertible ortho (o) and para (p) X₂ is to consider them as independent and with constant relative populations defined by the spin degeneracy factors. Hence, one has to consider the contribution for both o and p states of D₂ and H₂ to the total amount of adsorbed species to derive an estimate for the normal-D₂ selectivity. In doing so, however, a few details emerged as important, given the fact that the four species afford different adsorption energies. In order to clarify such aspects, we begin by deriving a more general formula for estimating the amount of adsorbed molecules within the low pressure regime when in presence of multi-species adsorption. In fact, the derivation is also correct if independence between sites and the adsorption of a single molecule per site are valid hypotheses, as in the derivation of the Langmuir isotherm. From the latter result, one easily derives adsorption selectivity for normal gases adding the contribution from both spin isomers of D₂ or H₂.

Assuming molecular species as non-interacting, the partition function for N_s different chemical species that may be adsorbed onto M independent sites, each bearing at most a single molecule of the n_s totally adsorbed for species s (1 ≤ s ≤ N_s), may be written as:

where N = (n₁, n₂, …, n_Ns) collectively indicates the number of adsorbed molecules for each species, q_ads(s) is the partition function for an adsorbed s molecule, and the combinatorial factor gives the total number of possible ways the molecules are placed on the M sites once their population (i.e. n_s) is defined. Using the Stirling formula n! = (n/e)ⁿ and computing the logarithm of Q one obtains:

The chemical potential for the s-th adsorbed species is easily obtained using the statistical mechanics relationship μ_s = −k_BT(∂Q/∂n_s) and reads:

where θ_s ≡ n_s/M is the fraction of the M total sites occupied by species s due to adsorption. At equilibrium, the chemical potential for the adsorbed s species is equal to its gas phase one (, with p_s being the partial pressure for species s), so that one obtains:where we have used the definition for χ_r(s) given previously. From the latter results, one realizes that θ_s/θ_t = p_sχ_r(s)/(p_tχ_r(t)) or θ_t = θ_sp_tχ_r(t)/(p_sχ_r(s)). Substituting in the denominator of eqn (9) the sum over all species for the latter result, and solving for θ_s, one gets:

From eqn (10), it becomes clear that the fraction of sites occupied by two species (e.g., o and p D₂) is easily estimated as ; a similar results is obtained also for the hydrogen diatoms, so that one has:

with x_s being the gas phase molar fraction of species s. Defining, for instance, x_oH2 = g_oH2x_nH2, with gx_oH2 being the statistical weight for oH₂ in normal H₂ (nH₂) and x_nH2, the overall molar fraction of H₂-like species in the gas mixture, one gets:from which the low pressure selectivity for normal D₂ and H₂ is found to bewhich allows to estimate once χ_r(s) are computed.

To carry out this task for the species allowed to populate only odd J states, one should recall that their gas phase rotational ground state energy differs from zero as E_J = BJ(J + 1), with B being the rotational constant. In this case, the exact q_ads(s)/q_r(s) ratio can be written as:

which, if T is sufficiently low to allow one neglecting the thermal excitation to states with J ≥ 3, becomeswhere E_des(s) = E₀(s) − 2B. In other words, the low-T approximation requires to use the difference between adsorbed and gas phase ground state energy rather than the simple ground state energy for the adsorbed species.

The derivation presented above makes clearly apparent that information on the quantum states of adsorbed diatoms, E_l(X₂), are needed in order to evaluate χ_r(s) and, subsequently, estimate both adsorption isotherms and selectivity. We have approached this requirement numerically tackling the solution of the Schrödinger's equation as described in Section 2.3, representing the molecules either as spherical pseudo-atoms (PsAt) or as rigid rotors (RR). The force fields used for those calculations are presented in Section 2.2.

We conclude this section by noting that our approach deviates other rigorous methods such as quantum Grand Canonical Monte Carlo, where is directly estimated.³⁹ This, however, is suitable when adsorption takes place inside porous materials (e.g. CNT^25,26,39 or MOF), as the molecular insertion moves employed to equilibrate chemical potential in the two phases are well defined in terms of the insertion location, with a limited volume for the attempt. Our cases are more complicated due to the lack of natural confining walls allowing the straightforward labeling of adsorbed molecules. In principle, one may place an arbitrarily confining wall parallel to the adsorption surface; this choice would, however, impose an undesired quantization on the translational motion of molecules, potentially biasing their chemical potential. Besides, investigating the low temperature regime that may be needed for the separation could become computationally demanding if quantum effects are important, as many replicas ought to be included in the “necklace” representing the delocalized distribution of each molecule.

2.2 Interaction potentials

In this study we have considered the Na-decorated systems as rigid structures and the surrounding hydrogen molecule either as a rigid body (RB) or as a point particle (PsAt model). Both PsAt and RB analytic interaction potentials have been already published for [CCCor–Na],²⁸ while for [H–Cor–Na]⁺, only the PsAt potential has been reported.¹⁷ Thus, we here present the X₂–[H–Cor–Na]⁺ interaction model within the RB approach, following a procedure analogous to that employed for X₂–[CCCor–Na]²⁸ as well as other systems,^16,40 with satisfactory comparisons with high level electronic structure calculations. Note that within the standard Born–Oppenheimer approximation, applied here, the electronic Hamiltonian does not depend on the masses or other properties of the nuclei, hence the interaction potentials employed are the same for H₂ and D₂.

The RB global potential energy V for the interaction between X₂ (X = H or D) and [H–Cor–Na]⁺ is expressed as a sum of pairwise interactions between X₂ and every atom j of the substrate, each term, in turn, having an electrostatic plus a non-covalent contribution:

The electrostatic contribution V^j_elec is a Coulomb interaction between point charges of X₂ and the partial charge of atom j in the substrate. For X₂, with an interatomic distance r_e = 0.7666 Å, three point charges were considered to reproduce its quadrupole moment:^40,41 q₁ = q₂ = 0.45955e on the X nuclei and q₃ = −2q₁, at the midbond site. Coordinates and partial charges associated to [H–Cor–Na]⁺ are given in Table S1 of the SI, corresponding to an optimized geometry obtained as described in ref. 17 with atomic charges estimated by means of a CM5 approach.⁴² The non-covalent component involves both induction and van der Waals interactions and it is expressed as a function of (R_j,γ_j), where R_j is the distance from the center of mass of X₂ to the atomic site j and γ_j, the angle formed between the X₂ molecular axis and the vector joining the site and the X₂ center of mass. Removing the subscript j to alleviate the notation, such a contribution is represented by the Improved Lennard-Jones (ILJ) formula,^43,44

wherewhich depend on the parameters β, R^⊥_e, R^‖_e, ε^⊥ and ε^‖. Finally, m is set to 4 for the X₂–Na interacting pair, for which the long range induction component is dominant, and to 6 for the remaining pairs(X₂–C and X₂–H).

The above mentioned parameters were fine-tuned (from an initial guess based on the polarizabilities of the partners) to reproduce electronic structure calculations carried out as detailed in ref. 17. Their values are reported in Table 1. Fig. 1 shows a comparison between the analytical representation and the electronic structure calculations of the X₂–[H–Cor–Na]⁺ interaction. They correspond to three different directions of approach of X₂ toward the decorated PAH and, for each of them, two different orientations of the X₂ molecular axis with respect to the PAH plane. Globally, a good agreement can be appreciated, with a very attractive interaction noticed for the X₂ monomer perpendicularly approaching the Na atom (red curves in upper panel of Fig. 1), mainly due to a maximization of the electrostatic interaction between the X₂ quadrupole and the partial positive charge associated to the alkali atom.

Table 1 Optimized parameters for the non-covalent contribution to the X₂–[H–Cor–Na]⁺ RB analytical potential (X = H and D). ε^⊥,‖ are in meV, R^⊥,‖_e are in Å and β and m are dimensionless

Pair	X2–[H–Cor–Na]^+
	m	β	R^⊥e	R^‖e	ε^⊥	ε^‖
X2–Na	4	4.70	2.700	2.820	59.00	50.00
X2–C	6	6.25	3.510	3.565	4.519	4.650
X2–H	6	6.25	3.153	3.228	2.778	2.756

---

Fig. 1 Interaction potentials for the X₂–[H–Cor–Na]⁺ dimer (X = H and D). The substrate is depicted in the insets, with Na, C and H atoms shown in blue, gray and red, respectively. Lines with circles correspond to reference electronic structure calculations¹⁷ while dashed lines, to the analytical curves. Red and black colors refer to parallel and perpendicular orientations of the H₂ molecular axis with respect to the PAH plane, respectively. Upper panel: Approach of X₂ on top of the Na atom as a function of the X₂–Na distance. Intermediate panel: X₂ approaching the opposite side of the PAH as a function of the distance from X₂ to the center of the coronene plane. Lower panel: X₂ approach towards an outer C–C bond, R being the distance between X₂ and the midpoint of the bond.

2.3 Bound states of H₂/D₂ attached to Na-decorated carbonaceous materials

To correctly estimate H₂/D₂ adsorption isotherms and selectivities, we have computed the bound states of an X₂ molecule attached to [H–Cor–Na]⁺ and [CCCor–Na]. Within the RB formulation, the Hamiltonian of the system is written aswhere (x,y,z) are the coordinates of the center of mass of the diatomic molecule, m_X2 = 2m_X (m_X = 1.0078 and 2.0141 u for H and D), , the angular momentum operator associated to the rotation of X₂ and B = ħ²/(2μ_X2r_e²)), the rotational constant, with values 81.89 and 40.98 K for H₂ and D₂, respectively. Finally, V is the interaction potential of eqn (16), depending on (x,y,z) and (θ,ϕ), the angles determining the orientation of the diatomic molecule with respect to the reference system. The bound states were obtained by building and diagonalizing a Hamiltonian matrix using, as a basis set, a direct product of fixed-node discrete variable representation (DVR) functions⁴⁵ for each Cartesian coordinate, and, for the angular degrees of freedom, spherical harmonics Y_JM(θ,ϕ). It is worth noting that this Hamiltonian does not couple rotational states with different parities, so in the end, the adsorbed molecular states will consist in linear combinations of spherical harmonics with either even or odd J values. In other words, within the present model, ortho and para X₂ molecules act as distinct species, as anticipated in a previous paragraph.

These calculations have been also carried out within the PsAt approach, where the procedure becomes simplified since the rotational term disappears in the Hamiltonian and the interaction potential, corresponding to the average over the orientations of the diatomic molecule, only depends on the position of the center of mass of X₂. Diagonalization of the Hamiltonian is carried out similarly, limiting the employed basis set to the DVR functions for the (x,y,z) coordinates.

3 Results and discussion

Following the steps indicated in Section 2, in this section we show the results concerning two different PAHs: [H–Cor–Na]⁺ and [CCCor–Na]. For simplicity, we use the parity of the rotational states (even-J or odd-J) to name the spin isomers of H₂ and D₂, bearing in mind that odd-J and even-J correspond to the ortho and para states of H₂, respectively, and the reverse for D₂.

3.1 The case of [H–Cor–Na]⁺

The ground state energies of the H₂/D₂–[H–Cor–Na]⁺ complex are presented in Table 2 within the PsAt and RB models. As expected from the larger mass of D₂, its energies are lower than the H₂ energies. This implies that the difference appearing in the last column of Table 2 is positive and ultimately, will lead to a selectivity favoring the D₂ species (see the exponential factor in eqn (5)). In addition, it can be seen that the PsAt ground state energies are less negative than the RB ones. This is due to the anisotropy of the H₂–[H–Cor–Na]⁺ interaction (i.e. see upper panel of Fig. 1), which makes the molecules orient along the most attractive regions of the RB potential, whereas the PsAt potential, being an average over all the orientations, is less attractive. Also, it has to be noted that, although the odd-J energies are smaller in absolute value than even-J ones when they are referred to the same asymptote (E₀), the odd species are more stable because their desorption energies, E_des = E₀ − 2B, are larger in absolute value than the ones of the even species, which coincide with E₀. This feature will favor the separation of the spin isomers of a given isotopologue, as discussed below.

Table 2 Ground state energies (in K), of H₂/D₂ attached to [H–Cor–Na]⁺ within the PsAt approach as well as for the even- and odd-J RB model. Energies E₀ are referred to a rotationless (J = 0) molecule infinitely separated from the substrate and coincide with the desorption energies for the PsAt and even-J RB cases (E_des = E₀). For the odd-J RB block, the ground state of the free molecule is J = 1 so E_des = E₀ − 2B. −ΔE(D₂ − H₂) is the difference between the D₂ and H₂ energies of the same row

		H2	D2	−ΔE(D2 − H2)
PsAt	E0	−871.1	−940.3	69.2
RB, even-J	E0	−1116.3	−1250.7	134.4
RB, odd-J	E0	−1016.8	−1208.1	191.3
RB, odd-J	Edes	−1180.5	−1290.0	109.5

---

Sets of bound state energies for H₂ and D₂ are displayed in Fig. 2 for the PsAt (a), RB even-J (b) and RB odd-J (c) approaches. Corresponding numerical values are provided in Tables S2 (PsAt), S4 (RB even-J) and S5 (RB odd-J) of the SI. They are given as energy gaps, ΔE_(l,0) = E_l − E₀, as these are the quantities used to compute the last factor in eqn (5). Recall that with a view toward a simplified modeling for multiply adsorbed diatoms, in a subsequent work, we aim to assess the impact of neglecting this factor in the determination of the selectivity ratios. To understand the structure of these states, we have first studied the probability distributions of the bound states within the PsAt approximation. Note that all the states are localized onto the Na-decorated face of the adsorbent. We have found that they can be approximately assigned to excitations along nearly separable modes in cylindrical coordinates: z, the axis perpendicular to the carbon support, ρ = (x² + y²)^1/2, the radial coordinate nearly centered on the sodium decoration, and the angle φ = tan⁻¹(y/x). Most of the low-lying states of Fig. 2a correspond to excitations along the angular coordinate. Probability distributions D_l(x,y) for some selected states are displayed in Fig. 3. Notice that the wavefunctions of the lowest levels (l = 0–3) are somewhat localized on a limited range of φ values: this is due to the extra proton in [H–Cor–Na]⁺ which leads to a more attractive interaction for some specific values of φ (see Fig. S1 of the SI). However, the distributions of higher states (l > 6) become quite similar to those of the motion of a free rotor around the z axis (solutions of the “particle on a ring” problem with cylindrical symmetry). For instance, panels e and f of Fig. 3 show two states (l = 9 and 10) with an equal number of nodes (only differing in their location); besides, these states are nearly degenerated (ΔE_(l,0) = 55.6 and 55.7 K). Other pairs of degenerated states can be identified in Fig. 2a, with the difference in energy between consecutive pairs growing nearly linearly, as expected when a free rotational motion is present. For the highest levels, some of the states are rather assigned to excitations along the z or ρ coordinates, therefore leading to a more complex spectrum.

Fig. 2 Energy levels of H₂ (yellow) and D₂ (green) adsorbed on [H–Cor–Na]⁺, as differences with respect to the ground state energy, ΔE_(l,0) = E_l − E₀ (in K), for the PsAt (a), RB even-J (b) and RB odd-J (c) approaches. See text for discussion.

Fig. 3 Probability densities projected onto the coronene plane, , for the l = 0, 1, 2, 3, 9 and 10 states (panels a–f) of the H₂–[H–Cor–Na]⁺ complex within the PsAt approximation.

For the RB approach and within the J = even block, the bound states are nearly a product of a node-less function of the angular degrees of freedom (θ,ϕ) and a function of the (x,y,z) coordinates, the latter with probability distributions resembling those of the PsAt approach; indeed, notice the similarity of the energy gaps of panels a and b of Fig. 2. Regarding the J = odd block, all the studied levels can be approximately assigned to J = 1, however, the structure of the bound states becomes more complex due to the contribution of states corresponding to the three projections of J = 1. Examples of these states are given in Fig. S2 of the SI; there it can be seen that the two lowest states have equivalent configurations on (x,y) but different distributions along (θ,ϕ), suggesting two distinct orientations of the H₂ axis. Similar behaviors are found for the D₂–[H–Cor–Na]⁺ bound states.

We have examined the sum involving the energy gaps , as it is one of the terms contributing to the D₂/H₂ selectivity. Fig. S3 of the SI shows these sums as functions of the total number of levels considered, for H₂ and D₂ and various temperatures. A sufficiently large number of states have been included to converge these sums with errors smaller than one part per hundreds. Notice that the sum is larger for D₂ than for H₂, therefore, the ratio appearing in the last factor of eqn (5) is larger than one at all temperatures, indicating that neglecting this factor leads to an underestimated value of the selectivity ratio .

Turning to the adsorption process of H₂ and D₂, Fig. 4 presents the Langmuir-like adsorption isotherms (eqn (1)) for the PsAt and RB approaches and including all relevant excited states (isotherms computed neglecting this contribution are provided in Fig. S4 of the SI). As one may notice, increasing the temperature exponentially shifts the adsorption isotherms toward higher partial pressures by several orders of magnitude, an interesting detail in the context of the possible usage of metal-decorated carbonaceous materials as hydrogen diatom adsorbents. Also, a clear separation between isotherms of the two isotopologues at each temperature is noticed, albeit the effect is reduced upon increasing the temperature. For instance, for the RB even-J block at T = 20 K, the half coverage pressure is about 300 times lower for D₂ than for H₂, whereas the disparity reduces to a factor of about ten at 40 K. In addition, comparing the PsAt and RB approaches, it is clear that the RB model leads to a considerable displacement of the isotherm curves towards lower pressures.

Fig. 4 Adsorption isotherms on [H–Cor–Na]⁺ for both H₂ and D₂ at 20, 30, and 40 K. Top panel (a): within the PsAt approach. Medium and lower panels (b and c): even-J and odd-J states, respectively, within a RB description for the adsorbed molecules.

Selectivities for adsorbing D₂ vs. H₂ on [H–Cor–Na]⁺ in the 20–40 K temperature range are reported in Fig. 5, both within the PsAt (upper panel) and RB (middle and lower panels) approaches. The PsAt selectivities are relatively small: they only become larger than six – the acceptable standard for an industrial application³² – for the lowest temperature. We also compare the exact PsAt selectivity (eqn (5)) with the low-temperature approximation where only the ground state energies are considered, the fractional or relative error involved being shown in the inset of the panel. It is found that the approximate selectivities are always lower than the exact ones, as anticipated above, with relative errors that increase with the temperature. In other words, we numerically demonstrate that estimating the selectivity by neglecting the thermal excitations in q_ads(T) produces a reasonably tight lower bound of this quantity.

Fig. 5 Exact and low-temperature approximate selectivities for D₂ versus H₂ on [H–Cor–Na]⁺. The insets show the fractional error of the approximate estimate. Upper panel (a): PsAt approximation. Middle panel (b): RB for a normal mixture of even and odd states of H₂ and D₂. Lower panel (c): RB approach, even-J states of H₂ and D₂.

Regarding the D₂ adsorption selectivities within the accurate RB approach and assuming that there is not interconversion between even and odd-J molecular states, we have computed selectivities following eqn (13) and considering a normal mixture of the rotational states (relative populations of 1 : 3 and 2 : 1 for the even:odd J states of H₂ and D₂, respectively). Results are presented in the middle panel of Fig. 5, also showing approximate selectivities where thermal excitations are neglected. Compared with the adsorption of PsAt-like species, the RB-normal description of the molecular gases significantly increases , by a factor of 2–4 depending on the temperature. Therefore, the PsAt approach is not accurate for predicting the D₂/H₂ selectivities on [H–Cor–Na]⁺, due to the anisotropy in this system. Also, notice that the RB-normal selectivity remains either well above, or very close to the industrially acceptable standard for the temperature range investigated.

Albeit the results just presented are clearly supportive for the capability of [H–Cor–Na]⁺ to chromatographically separate the two isotopomers, even better results could be obtained if each gas was allowed to reach the thermodynamic equilibrium for the interconversion between odd and even J rotational states. In order to have an estimate of the selectivity in these conditions, we have assumed that the H₂ and D₂ species are just composed of even-J states and have computed the corresponding selectivities, which are presented in the lower panel of Fig. 5. It can be seen that the selectivities become about 7 and 3 times larger than those for the normal mixture, for T = 20 and 40 K, respectively. While it may seem counter-intuitive to find a higher selectivity for the spin-purified species compared to a normal mixture, this is justified by realizing that, in the normal gas case (eqn (13)), the two species with the highest relative populations are even-J for D₂ and odd-J for H₂, whose ground state desorption energies differ the least (69 K), whereas the difference between the desorption energies of even-J D₂ and H₂ is much larger (134 K).

Within the assumption that there is not interconversion between the spin isomers of H₂ or D₂, states, one may further wonder whether [H–Cor–Na]⁺ would be capable of separating, e.g., odd-H₂ from even-H₂. Fig. 6 displays the selectivities for the separation of odd-H₂ and odd-D₂, which are the species most tightly adsorbed in comparison with their even-J counterparts. Separation of the spin isomers is clearly more efficient for H₂ than for D₂; this is mainly due to the larger difference between the odd and even ground energies in H₂ (64 K) as compared with D₂ (39 K). We also show the selectivities resulting of neglecting thermal excited states as well as the error of this approximation. As in the study of the D₂/H₂ separation, these approximate selectivities are lower bounds of the accurate ones; however, the associated errors are larger. This is explained by the small energy gaps of the excited states of the odd-J species (related to the threefold degeneracy of J = 1) which lead to relatively large sums-over-states contributing to the exact selectivity.

Fig. 6 Exact (solid lines) and approximate (dashed lines) selectivities for separation between odd-J and even-J isomers of H₂ (in brown) and D₂ (in pink), on [H–Cor–Na]⁺; the inset shows the fractional errors of the approximate estimates.

3.2 The case of [CCCor–Na]

This section presents the study of the adsorption of H₂ and D₂ on a different PAH, sodium-decorated circumcircumcoronene, [CCCor–Na]. The goal has been to investigate whether the conclusions reached for Na-decorated protonated coronene hold for [CCCor–Na], a prototype of an extended neutral system, i.e., Na-decorated graphene. Comparing with the smaller PAH, the X₂–[CCCor–Na] interaction is somewhat less attractive²⁸ mainly due to a weakening of the H₂–Na interaction, which is compensated by a stronger interaction with the more extended carbon plane. Thus, it is worth to examine how these differences may affect the detailed features of the H₂/D₂ adsorption and selectivities.

Table 3 presents the ground state energies of H₂ and D₂ adsorbed on [CCCor–Na] for the different models considered. Comparing with the [H–Cor–Na]⁺ support (Table 2), the energies are somewhat smaller in absolute value because of the mentioned weakening of the X₂–[CCCor–Na] interaction. Excited states spectra are presented in Fig. S5 as energy gaps, ΔE_(l,0). The general structure of these states is similar to that of decorated coronene, however, the gaps are relatively smaller, hence we can expect a larger contribution of the excited states to the selectivities. In addition, there are some differences in the structure and probability distributions of the excitations, mainly due to the C₆ symmetry of [CCCor–Na], a symmetry that is lacking in [H–Cor–Na]⁺ due to the extra proton in the periphery of the PAH. We present some of these probability distributions in Fig. 7, for the J = odd block of the RB approach. For the ground state (l = 0) (Fig. 7a and b), it can be seen that the center of mass of H₂ is completely delocalized on a ring around the Na atom and, regarding the orientations distribution , the state is assigned to J = 1 with the molecule nearly perpendicular to the PAH plane (i.e., a p_z-like orbital shape). The first excited state (Fig. 7c and d) has a node along the y axis, while the distribution is the same as that of the l = 0 state (on the contrary, the l = 0, 1 states of H₂–[H–Cor–Na]⁺, in Fig. S2, have different molecular orientation). Counterparts of the l = 0 and 1 states are the l = 7 and 9 states (Fig. 7e–h) which have the same spatial distributions but differ in the molecular orientation. Then, the splittings of the three-fold J = 1 states are larger than those of the smaller PAH. Finally, note that the l = 2 state (not shown) is nearly degenerate with l = 1 (gaps of 1.9 and 1.6 K, respectively), the only difference being that the direction of the node in is along the x axis. This scheme is repeated for most of the excited states, which appear as nearly degenerated pairs. In this way, the H₂–[CCCor–Na] bound states closely resemble those of a free rotor around the axis perpendicular to the carbon support.

Table 3 As in Table 2 for the ground state energies (in K), of H₂/D₂ attached to [CCCor–Na]

		H2	D2	−ΔE(D2 − H2)
PsAt	E0	−788.5	−849.0	60.5
RB, even-J	E0	−1097.9	−1239.2	141.3
RB, odd-J	E0	−1011.5	−1206.1	194.6
RB, odd-J	Edes	−1175.3	−1288.0	112.7

---

Fig. 7 Probability densities of selected states of H₂ adsorbed on [CCCor–Na], within the RB model, J = odd block (l = 0, 1, 7, 9). Densities projected onto the (x,y) plane are presented in the left panels (a, c, e, g) whereas, in the right panels (b, d, f, h), they have been projected onto the orientation of the molecule (θ,ϕ) and are displayed as polar plots (i.e., the density is the radius in a representation using spherical coordinates).

With regard to the properties gauging the adsorption of [CCCor–Na], isotherms for both isotopologues are depicted in Fig. S6 of the SI. It can be noticed that, except for some differences within the PsAt approach, [CCCor–Na] isotherms are very similar to the ones for the [H–Cor–Na]⁺ support (Fig. 4).

Turning to the adsorption selectivity, Fig. 8 reports D₂/H₂ selectivities within the same approaches and approximations presented for [H–Cor–Na]⁺. General features of these selectivities are similar to those encountered in the smaller prototype, in particular, the PsAt model also fails to provide reasonable values for the selectivities and the low-temperature approximation also leads to lower bounds of the selectivities. Nevertheless, it can be noticed that the RB selectivities (normal mixture and even-J block) are larger for [CCCor–Na]. This is due in part to the larger differences between the ground state energies of D₂ and H₂ (see last columns of Tables 2 and 3). In addition, the role of thermally excited states in enhancing the selectivities with respect to the low temperature approximation is more important in the present substrate, because of the narrower energy gaps of its spectrum of adsorbed states. Therefore, [CCCor–Na] is somehow more effective than [H–Cor–Na]⁺ for our separation purpose, a result that could seem surprising given the somewhat weaker binding in the case of the larger prototype.

Fig. 8 As in Fig. 5 for the exact and approximate selectivities for D₂ versus H₂ on [CCCor–Na].

Finally, Fig. 9 reports the selectivities for the separation of odd-J H₂ from even-J H₂, and equivalently for D₂. Again, a similar behavior is found as compared with the smaller substrate (Fig. 6) but with rather larger selectivities involved: varies from 59 at 20 K to 10 at 40 K, while spans the range 14-5 over the same temperature interval. Thus, with the assumption that the kinetics of converting odd-J species into even-J ones is very slow, it seems that [CCCor–Na] may be able to effectively separate ortho-H₂ from para-H₂.

Fig. 9 As in Fig. 6 for the exact and approximate selectivity for odd-J versus even-J states adsorbed on [CCCor–Na].

4 Conclusions

In this work, we have investigated whether it is possible to exploit sodium decorated polycyclic aromatic hydrocarbons (PAHs) for the adsorption and separation of hydrogen isotopomers, specifically D₂ and H₂. In particular, we have considered Na-decorated protonated coronene and Na-decorated circumcircumcoronene. We have approached this task using electronic structure calculations to calibrate analytical model potentials, describing the diatoms either as pseudo-atoms (PsAt) or as rigid bodies (RB) species. The Statistical Mechanics description of these adsorption processes requires the determination of the eigenenergies of the adsorbed species, which we have obtained by numerically solving the involved Schrödinger equation. Compared to previous discussions in the literature limited to the usage of the harmonic oscillator approximation, our approach allowed us to highlight some details playing an important role in the studied systems, among which the facile excitation of X₂ orbital motion around the metallic decoration, resembling the ones obtained for the “particle on a ring” model system. The latter aspect would have been entirely missing if the quantum-mechanical description was limited to the level of localized harmonic oscillator states. Moreover, we have explicitly evaluated the impact of having a mixture of odd and even-J species determining their respective adsorption properties, a feature that is not commonly addressed; in fact, only recently Path Integral Monte Carlo approaches have been extended to effectively sample odd-J states in H₂.⁴⁶

Within the assumption of a single molecule adsorbed per metal site, we estimated the low pressure selectivity toward D₂ of the decorated PAHs at low temperatures (20–40 K). For a normal mixture of ortho and para states of H₂ and D₂ the PAHs studied afford a deuterium selectivity well above six – the required value for industrial applications – for most of the temperature range, reaching a value of 48 at 20 K. Actually, the sodium-decorated circumcircumcoronene displays higher selectivity toward normal D₂ despite the lower evaporation energy, suggesting that it may not be necessary for an adsorbent to bear a net charge to be an effective stationary phase for the chromatographic separation of the isotopologues. These selectivity values can be compared with experimental D₂/H₂ adsorption selectivities reported by Scholl⁴ and Chen et al.⁷ for various porous materials and temperatures (20–100 K). While some high selectivity values (about 20–35) are found for some zeolites or MOFs, selectivities for the reported carbon materials are particularly low (below 7). In this way, the Na-decorated substrates here studied constitute promising carbon-based materials for effective separation of hydrogen isotopologues.

Even higher selectivities may be obtained if the normal gas mixtures were allowed to reach the appropriate equilibrium, i.e. if the transformation odd-JX₂ → even-JX₂ could take place. In principle, the equilibration would require a spin source capable of flipping the nuclear spin states, which may be provided by either the sodium atom nucleus or by its unpaired valence electron. Whether this is feasible in practice would fundamentally depend on the interaction strength between the spin sources, the exploration of this possibility being left for future work. In addition, we have studied the separation of the ortho and para species of a given isotopologue, finding that this is indeed a feasible process for hydrogen.

The validity of some approximated treatments for the calculation of adsorption isotherms and low pressure selectivities has been also examined in this work. First, it is found that the PsAt model is far from being adequate to estimate any of the two quantities, due to the considerable anisotropy of the interactions. In addition, neglecting the thermal excitations in the calculation of the partition functions provides a lower bound to the accurate selectivities, while the isotherms are shifted to higher pressure regions. These results may be useful in undertaking more complex studies, where either the RB approach or the calculations of excited bound states become computationally challenging. This is the case of the study of the adsorption of multiple H₂ and D₂ molecules on the metal decorated PAHs, where it will be interesting to investigate the conditions for the selectivity keeping the values here predicted or becoming increased, due to cooperative effects, or reduced, due to solvation of the decorating atom (part II in this series).

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Tables S1–S7. See DOI: https://doi.org/10.1039/d6cp00873a.

Acknowledgements

This work was supported by the Spanish Agencia Estatal de Investigación, project reference PID2023-149406NB-I00/AEI/10.13039/501100011033 and Fondo di Ateneo per la Ricerca (FAR) dell'Universitá degli Studi dell'Insubria. COST Action CA21101 (COSY) and allocation of computing time by CESGA (Spain) are also acknowledged.

References

1. S. Kopf, F. Bourriquen, W. Li, H. Neumann, K. Junge and M. Beller, Recent Developments for the Deuterium and Tritium Labeling of Organic Molecules, Chem. Rev., 2022, 122(6), 6634–6718,  DOI:10.1021/acs.chemrev.1c00795.
2. M. Rethinasabapathy, S. M. Ghoreishian, S.-K. Hwang, Y.-K. Han, C. Roh and Y. S. Huh, Recent Progress in Functional Nanomaterials towards the Storage, Separation, and Removal of Tritium, Adv. Mater., 2023, 35(48), 2301589,  DOI:10.1002/adma.202301589.
3. L. Shere, A. K. Hill, T. J. Mays, R. Lawless, R. Brown and S. P. Perera, The next generation of low tritium hydrogen isotope separation technologies for future fusion power plants, Int. J. Hydrogen Energy, 2024, 55, 319–338,  DOI:10.1016/j.ijhydene.2023.10.282.
4. D. S. Sholl, Quantum Sieving for Isotopic Separations of Gases Using Porous Materials–30 Years of Progress, J. Phys. Chem. C, 2025, 129(24), 10793–10801,  DOI:10.1021/acs.jpcc.5c02981.
5. J. Beenakker, V. Borman and S. Y. Krylov, Molecular transport in subnanometer pores: zero-point energy, reduced dimensionality and quantum sieving, Chem. Phys. Lett., 1995, 232(4), 379–382.
6. Q. Wang, S. Challa, D. Sholl and J. Johnson, Quantum sieving in carbon nanotubes and zeolites, Phys. Rev. Lett., 1999, 82(5), 956–959,  DOI:10.1103/PhysRevLett.82.956.
7. Y. Chen, X. Bai, D. Liu, X. Fu and Q. Yang, High-Throughput Computational Exploration of MOFs with Open Cu Sites for Adsorptive Separation of Hydrogen Isotopes, ACS Appl. Mater Interfaces, 2022, 14(21), 24980–24991,  DOI:10.1021/acsami.2c06966.
8. J. Y. Kim, R. Balderas-Xicohténcatl, L. Zhang, S. G. Kang, M. Hirscher, H. Oh and H. R. Moon, Exploiting diffusion barrier and chemical affinity of metal–organic frameworks for efficient hydrogen isotope separation, J. Am. Chem. Soc., 2017, 139(42), 15135–15141.
9. R. Muhammad, S. Jee, M. Jung, J. Park, S. G. Kang, K. M. Choi and H. Oh, Exploiting the specific isotope-selective adsorption of metal–organic framework for hydrogen isotope separation, J. Am. Chem. Soc., 2021, 143(22), 8232–8236.
10. S. Kumar, J. Park, H. Kim, S. Jung, T. Park, J. Park, H. J. Cho, C. W. Park and H. Oh, Quantum sieving in robust microporous zeolite SSZ–13 for efficient hydrogen isotope separation, Int. J. Hydrogen Energy, 2025, 100, 184–190,  DOI:10.1016/j.ijhydene.2024.12.298 , cited by: 4.
11. S. H. So and H. Oh, A mini-review of the current progress and future challenges of zeolites for hydrogen isotopes separation through a quantum effect, Int. J. Hydrogen Energy, 2024, 50, 539–560,  DOI:10.1016/j.ijhydene.2023.08.241.
12. S. Puricelli, G. Bruno, C. Gatti, A. Ponti and M. Mella, Viability of hydrogen isotopes separation via heterolytic dissociation-driven Chemical Affinity Quantum Sieving on inexpensive alkali-earth oxides, Appl. Surf. Sci., 2024, 657, 159596,  DOI:10.1016/j.apsusc.2024.159596.
13. R. C. Lochan and M. Head-Gordon, Computational studies of molecular hydrogen binding affinities: The role of dispersion forces, electrostatics, and orbital interactions, Phys. Chem. Chem. Phys., 2006, 8, 1357–1370.
14. E. Anikina, S. R. Naqvi, H. Bae, H. Lee, W. Luo, R. Ahuja and T. Hussain, High-capacity reversible hydrogen storage properties of metal-decorated nitrogenated holey graphenes, Int. J. Hydrogen Energy, 2022, 47(19), 10654–10664.
15. H. Li, J. H. Lim, Y. Lv, N. Li, B. Kang and J. Y. Lee, Graphynes and Graphdiynes for Energy Storage and Catalytic Utilization: Theoretical Insights into Recent Advances, Chem. Rev., 2023, 123(8), 4795–4854,  DOI:10.1021/acs.chemrev.2c00729.
16. A. M. Reider, S. Kollotzek, P. Scheier, F. Calvo, E. Yurtsever, F. Pirani, M. Bartolomei, M. I. Hernández, T. González-Lezana and J. Campos-Martínez, Experimental and theoretical assessment of the enhanced hydrogen adsorption on polycyclic aromatic hydrocarbons upon decoration with alkali metals, Int. J. Hydrogen Energy, 2024, 58, 525–535,  DOI:10.1016/j.ijhydene.2024.01.244.
17. E. García-Arroyo, A. M. Reider, S. Kollotzek, F. Foitzik, J. Campos-Martínez, M. Bartolomei, F. Pirani, M. I. Hernández, M. Mella and P. Scheier, The role of Na decoration on the hydrogen adsorption on coronene: A combined experimental and computational study, Int. J. Hydrogen Energy, 2024, 83, 387–395,  DOI:10.1016/j.ijhydene.2024.07.425.
18. I. Cabria, M. J. López and J. A. Alonso, Enhancement of hydrogen physisorption on graphene and carbon nanotubes by Li doping, J. Chem. Phys., 2005, 123(20), 204721.
19. P. Pei, M. B. Whitwick, W. L. Sun, G. Quan, M. Cannon and E. Kjeang, Enhanced hydrogen adsorption on graphene by manganese and manganese vanadium alloy decoration, Nanoscale, 2017, 9, 4143–4153,  10.1039/C6NR09545C.
20. S. Ghotia, T. Rimza, S. Singh, N. Dwivedi, A. K. Srivastava and P. Kumar, Hetero-atom doped graphene for marvellous hydrogen storage: unveiling recent advances and future pathways, J. Mater. Chem. A, 2024, 12, 12325–12357,  10.1039/D4TA00717D.
21. R. K. Sahoo, B. Chakraborty and S. Sahu, Reversible hydrogen storage on alkali metal (Li and Na) decorated C₂₀ fullerene: a density functional study, Int. J. Hydrogen Energy, 2021, 46(80), 40251–40261.
22. V. Mahamiya, A. Shukla and B. Chakraborty, Potential reversible hydrogen storage in Li-decorated carbon allotrope PAI-Graphene: A first-principles study, Int. J. Hydrogen Energy, 2023, 48(96), 37898–37907.
23. M. P. de Lara-Castells and A. O. Mitrushchenkov, Spectroscopy of a rotating hydrogen molecule in carbon nanotubes, Phys. Chem. Chem. Phys., 2019, 21, 3423–3430,  10.1039/C8CP04109A.
24. M. P. de Lara-Castells, A. W. Hauser, A. O. Mitrushchenkov and R. Fernández-Perea, Quantum confinement of molecular deuterium clusters in carbon nanotubes: ab initio evidence for hexagonal close packing, Phys. Chem. Chem. Phys., 2017, 19, 28621–28629,  10.1039/C7CP05869A.
25. G. Garberoglio and J. Johnson, Hydrogen isotope separation in carbon nanotubes: Calculation of coupled rotational and translational states at high densities, ACS Nano, 2010, 4(3), 1703–1715.
26. S. R. Challa, D. S. Sholl and J. K. Johnson, Light isotope separation in carbon nanotubes through quantum molecular sieving, Phys. Rev. B:Condens. Matter Mater. Phys., 2001, 63, 245419,  DOI:10.1103/PhysRevB.63.245419.
27. G. Garberoglio, M. M. DeKlavon and J. K. Johnson, Quantum sieving in single-walled carbon nanotubes: effect of interaction potential and rotational-translational coupling, J. Phys. Chem. B, 2006, 110(4), 1733–1741,  DOI:10.1021/jp054511p.
28. E. García-Arroyo, M. Bartolomei, M. I. Hernández, T. González-Lezana, M. Mella and J. Campos-Martínez, Adsorption of Molecular Hydrogen on Extended Graphene Prototypes Decorated with Sodium: A Quantum-Mechanical Study, Phys. Scr., 2024, 100(1), 015411.
29. F. Dulieu, L. Amiaud, S. Baouche, A. Momeni, J.-H. Fillion and J. Lemaire, Isotopic segregation of molecular hydrogen on water ice surface at low temperature, Chem. Phys. Lett., 2005, 404, 187–191,  DOI:10.1016/j.cplett.2005.01.044.
30. M. Mella and E. Curotto, Quest for inexpensive hydrogen isotopic fractionation: Do we need 2D quantum confining in porous materials or are rough surfaces enough? The case of ammonia nanoclusters, J. Phys. Chem. A, 2016, 120(41), 8148–8159,  DOI:10.1021/acs.jpca.6b08005.
31. T. L. Hill, An Introduction to Statistical Thermodynamics, Dover, New York, NY, 1986.
32. Z. Zhu, Permeance should be used to characterize the productivity of a polymeric gas separation membrane, J. Membr. Sci., 2006, 281(1), 754–756,  DOI:10.1016/j.memsci.2006.04.040.
33. D. Cao, J. Ren, Y. Gong, H. Huang, X. Fu, M. Chang, X. Chen, C. Xiao, D. Liu, Q. Yang, C. Zhong, S. Peng and Z. Zhang, Quantum sieving of H₂/D₂ in MOFs: A study on the correlation between the separation performance, pore size and temperature, J. Mater. Chem. A, 2020, 8(13), 6319–6327,  10.1039/c9ta14254a.
34. T. A. Grinev, A. A. Buchachenko and R. V. Krems, Separation of ortho-and para-hydrogen in van der Waals complex formation, ChemPhysChem, 2007, 8(6), 815–818.
35. A. L. Myers and J. M. Prausnitz, Thermodynamics of mixed-gas adsorption, AIChE J., 1965, 11(1), 121–127,  DOI:10.1002/aic.690110125.
36. S. R. Challa, D. S. Sholl and J. K. Johnson, Adsorption and separation of hydrogen isotopes in carbon nanotubes: Multicomponent grand canonical Monte Carlo simulations, J. Chem. Phys., 2002, 116(2), 814–824,  DOI:10.1063/1.1423665.
37. B. C. Hathorn, B. G. Sumpter and D. W. Noid, Contribution of restricted rotors to quantum sieving of hydrogen isotopes, Phys. Rev. A:At., Mol., Opt. Phys., 2001, 64, 022903,  DOI:10.1103/PhysRevA.64.022903.
38. G. Colonna, A. D'Angola and M. Capitelli, Statistical thermodynamic description of H₂ molecules in normal ortho/para mixture, Int. J. Hydrogen Energy, 2012, 37(12), 9656–9668,  DOI:10.1016/j.ijhydene.2012.03.103.
39. G. Garberoglio, M. DeKlavon and J. Johnson, Quantum sieving in single-walled carbon nanotubes: Effect of interaction potential and rotational-translational coupling, J. Phys. Chem. B, 2006, 110(4), 1733–1741,  DOI:10.1021/jp054511p.
40. J. Ortiz de Zárate, M. Bartolomei, T. González-Lezana, J. Campos-Martínez, M. I. Hernández, R. Pérez de Tudela, J. Hernández-Rojas, J. Bretón, F. Pirani, L. Kranabetter, P. Martini, M. Kuhn, F. Laimer and P. Scheier, Snowball formation for Cs⁺ solvation in molecular hydrogen and deuterium, Phys. Chem. Chem. Phys., 2019, 21, 15662–15668.
41. M. Bartolomei, T. González-Lezana, J. Campos-Martínez, M. I. Hernández and F. Pirani, Complexes of Alkali Metal Cations and Molecular Hydrogen: Potential Energy Surfaces and Bound States, J. Phys. Chem. A, 2019, 123(39), 8397–8405,  DOI:10.1021/acs.jpca.9b05937.
42. A. V. Marenich, S. V. Jerome, C. J. Cramer and D. G. Truhlar, Charge Model 5: An extension of Hirshfeld population analysis for the accurate description of molecular interactions in gaseous and condensed phases, J. Chem. Theory Comput., 2012, 8(2), 527–541.
43. F. Pirani, S. Brizi, L. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi, Beyond the Lennard–Jones Model: A Simple and Accurate Potential Function Probed by High Resolution Scattering Data Useful for Molecular Dynamics Simulations, Phys. Chem. Chem. Phys., 2008, 10, 5489–5503.
44. F. Pirani, M. Albertí, A. Castro, M. M. Teixidor and D. Cappelletti, Atom-Bond Pairwise Additive Representation for Intermolecular Potential Energy Surfaces, Chem. Phys. Lett., 2004, 394, 37–44.
45. J. T. Muckerman, Some Useful Discrete Variable Representations for Problems in Time-Dependent and Time-Independent Quantum Mechanics, Chem. Phys. Lett., 1990, 173, 200–205.
46. G. Garberoglio and A. H. Harvey, First-Principles Calculation of ortho–para Effects in the Second Virial Coefficients of H₂ and D₂ at Low Temperatures, J. Phys. Chem. A, 2025, 129(36), 8453–8463,  DOI:10.1021/acs.jpca.5c04224.

---