Isotopic separation of helium through graphyne membranes: a ring polymer molecular dynamics study

Abstract

Microscopic-level understanding of the separation mechanism for two-dimensional (2D) membranes is an active area of research due to potential implications of this class of membranes for various technological processes. Helium (He) purification from the natural resources is of particular interest due to the shortfall in its production. In this work, we applied the ring polymer molecular dynamics (RPMD) method to graphdiyne (Gr2) and graphtriyne (Gr3) 2D membranes having variable pore sizes for the separation of He isotopes, and compare for the first time with rigorous quantum calculations. We found that the transmission rate through Gr3 is many orders of magnitude greater than Gr2. The selectivity of either isotope at low temperatures is a consequence of a delicate balance between the zero-point energy effect and tunneling of ⁴He and ³He. In particular, a remarkable tunneling effect is reported on the Gr2 membrane at 10 K, leading to a much larger permeation of the lighter species as compared to the heavier isotope. RPMD provides an efficient approach for studying the separation of He isotopes, taking into account quantum effects of light nuclei motions at low temperatures, which classical methods fail to capture.

---

1 Introduction

The use of two-dimensional (2D) materials is an open field for many technological applications as well as for basic science.^1,2 Graphene, one of the pioneering materials that started this very active area of research, acts as an almost impermeable sheet for most atomic and molecular species due to its high electron density surrounding the aromatic rings.^3,4 Some years ago, experiments by Lozada-Hidalgo et al.^5,6 and recently by Creager et al.⁷ have shown that small charged particles such as the proton and its isotopes can penetrate through a monolayer of pristine graphene. There is no complete explanation despite many theoretical calculations,^8–10 although some kind of chemical interaction seems to be responsible for this process. Even more recently, it has also be found that hydrogen can permeate pristine graphene.¹¹ However, the penetration energy barrier for other neutral atoms is significantly higher and therefore supports the widely recognized notion of graphene impermeability.¹² Introducing defects and moderate annealing into the membrane, and nanopores of different sizes, may enable graphene and other 2D layers to act as atomic and molecular sieves.^12–20 The use of 2D materials as a filter at the molecular level is one of the most interesting applications.¹⁴ Besides using 2D membranes with fabricated nanopores, some compounds contain nanopores in their structures at different positions and sizes; therefore, they are more naturally suited for filtering at the molecular level.

An interesting class of 2D nanoporous material is graphyne,^21,22 which is composed of sp–sp² hybridized C atoms and can be considered as a graphene derivative. In graphyne, one-third of C–C bonds have been replaced by mono- and poly-acetylenic units (–C≡C–)_N. Baughman et al.,²³ theoretically proposed the first stable structures of graphyne, which was two decades later synthesized on a copper substrate.^24,25N defines the number of acetylenic linkages and, consequently, the size of the uniformly distributed and repeating sub-nanometer triangular pores. They are termed as graph-N-yne membranes, such as graphdiyne (Gr2), graphtriyne (Gr3), etc., for N = 2 and 3, respectively.²⁶ The chemical and mechanical properties of these graphyne membranes have many useful features, such as they are chemically inert and stable at ambient temperatures,^27–29 and flexible enough to withstand deformations induced by high pressures.^30,31 The above-mentioned properties of graphynes, coupled with their unique geometrical structure, make them an excellent candidate to be utilized in gas separation and water filtration technologies. Indeed, one can find many reports on the purification of gases such as H₂,³² N₂,³³ and O₂³⁴ from a mixture of gases and desalination and filtration of water^26,31,35 in the literature.

There is a growing worldwide demand for helium purification from its natural resources due to the shortfall in its production.³⁶ It has numerous industrial and scientific applications such as superconducting magnets, space rockets, arc welding, etc. In particular, the lighter isotope, ³He, also plays a pivotal role in fundamental research, for example, in neutron-scattering centers, ultracold physics and chemistry, etc.³⁷ The relative abundance of ³He is low (≈1.34 × 10⁻⁶%)³⁸ in comparison to its heavier isotope, and its extraction from natural gas is usually done by means of expensive cryogenic distillation and pressure-swing adsorption methods.³⁹ An alternate and more energy-efficient process is to use the 2D porous membranes in isotopic gas separation since they usually do not involve costly liquefaction of the gases.⁴⁰ Many theoretical works have been reported in the last decade on the separation of He isotopes using graphene derivatives, such as: polyphenylene (2D-PP),^37,41–44 functionalized graphene pores,^16,45–47 nanoporous multilayers,^46,48 Gr2,^37,49,50 Gr3,⁵¹ holey graphene⁵⁰, and graphenylene membranes,⁵²etc.

The light element He within the vicinity of subnanometer pores is an obvious environment for observing the important role of the quantum mechanical effects such as zero-point energy (ZPE) and tunneling effects. A wise approach for an effective isotopic separation maybe is to exploit these quantum properties that could run in the opposite directions. The heavier isotope with smaller ZPE will diffuse faster, while quantum mechanical (QM) tunneling favours the lighter species. Previously, some of us, using quantum three-dimensional wave packet calculations have shown that the ⁴He/³He selectivity increases with decreasing temperature for Gr2 and holey graphene membranes.^49,50 It has also been reported that the effect of ZPE is more dominant than tunneling at low temperatures (20–40 K).⁴⁹ On the contrary, within a low but acceptable gas flux and at low temperatures (10–30 K), separation on various functionalized membranes indicates increased selectivity for the lighter isotope due to QM tunneling.^16,45,46,52

The ring polymer molecular dynamics (RPMD) method,^53–56 based on the imaginary-time path integral formalism, is an efficient approach that can accurately and reliably describe the ZPE⁵⁷ and deep quantum tunneling effects.^53,58,59 RPMD method is essentially a classical molecular dynamics method in an extended ring polymer phase space. It can provide reliable estimates of thermal rate coefficients since the RPMD partition function rigorously converges to the QM partition function.⁶⁰ The long-time limit of the ring polymer flux-side correlation functions is independent of the choice of the dividing surface that separates reactants from products, while its short-time limit is related to various quantum transition state theories.⁵⁶ RPMD was introduced in an ad hoc manner by Craig and Manolopoulos to study the dynamics of the condensed phase processes,^53,54 owing to its simplicity and efficiency (scales favorably with the dimensionality of the system). Examples of its successful application include diffusion in and inelastic neutron scattering from liquid para hydrogen,^61,62 the translational and orientational diffusion in liquid water,⁶³ proton transfer in water,⁶⁴ diffusion of H and μ atoms in liquid water, hexagonal ice,⁶⁵ and on Ni surface,⁶⁶ electron transfer⁶⁷ and proton-coupled electron transfer,⁶⁸ enzyme catalysis,⁶⁹etc. However, the RPMD method is not restricted to the condensed phases and also found wide application in calculating rate coefficients for the gas-phase bimolecular reactions^70–72 as explored by Suleimanov and co-workers, see for instance, a review by one of us⁵⁶ and a recent paper⁷³ and references mentioned therein. The RPMD method has recently been used to describe phenomena such as atom scattering or sticking on surfaces^74,75 and a related problem on permeation of H⁺ and D⁺ on pristine graphene.⁸ Nonetheless, in all cases, the authors point out, again, the importance of considering both ZPE and tunneling for isotopic separation.

In this work, we test the accuracy of the RPMD method compared to pure quantum results, with a similar model to make a fair and righteous comparison, providing evidence that RPMD can be used down to very low temperature and in conditions where quantum phenomena dominate, which has not been reported before. To this end, coupled with the major scientific and industrial appeal, in this paper, we examine the feasibility of separation of He isotopes using 2D Gr2 and Gr3 membranes at low temperatures (10–250 K) using the RPMD method. It is worth noting that the quantum rate of permeation of He isotopes through the pores of Gr2 at cold temperature regime (down to 10 K) has not been reported previously, probably because wave packet calculations are often difficult to converge at very low temperature. A significant part of the success of the RPMD method is attributed to the fact that it gives the exact quantum-mechanical rate coefficient for the transmission through a parabolic barrier,⁵³ which is advantageous for the present investigation. The primary objective of this study is to provide a reliable estimate of the selectivity and to show that RPMD is a necessary alternative to study these processes where quantum effects are expected to be very important.

The paper is organized as follows: in the next Section (2), we provide the details of the RPMD approach along with the potential energy surface (PES) used in the present study. The results of RPMD rate coefficients and selectivity have been compared with earlier studies^49–51 in Section 3. Concluding remarks are provided in Section 4.

2 Computational details

2.1 The simulation setup and potential energies

We investigate the rate of transmission of both ³He and ⁴He through 2D graphdiyne (Gr2) and graphtriyne (Gr3) membranes using the RPMD method. The unit cell of Gr2 and Gr3 has the dimensions (in x- and y-coordinates) of (16.37 Å, 9.45 Å) and (20.82 Å, 12.02 Å), respectively (see Fig. 1(a) and (b)). A comprehensive study on the geometry of the membranes can be found in ref. 37. All molecular dynamics simulations have been performed on graphyne membranes using same periodic unit as in ref. 50.

Fig. 1 Schematic representation of the unit cell of two-dimensional (a) graphdiyne and (b) graphtriyne membranes employed in this study. The x- and y-dimensions are in the unit of Å. The solid lines represent bonds between two neighboring C atoms. The filled gray circle denotes He atom penetrating through one of the four energetically most favorable and equivalent site(s) of the membrane (x = 0; y = 0). The four equivalent sites can be identified in the inset plot of Fig. 2b of ref. 50.

The interaction potential between He–graphyne membranes is obtained as an additive improved Lennard-Jones (ILJ) He–C pair potentials, as described in detail in an earlier work.⁷⁶ The optimized values of the parameters of the ILJ potentials have been determined from “coupled” supermolecular second-order Møller–Plesset perturbation theory (MP2C) theory⁷⁷ using aug-cc-pVTZ and aug-cc-pV5Z basis set for C- and He-atoms respectively. These IJL potentials correspond to the same potentials as employed in previous studies of He transmission through graphyne membranes,^37,49–51 and a detailed description can be found elsewhere.^37,49 The accuracy of the MP2C method as well as the analytical potentials can be ascertained from rigorous comparison with previous force fields²⁶ and ab initio calculations,^37,49 respectively.

The potential energy curves of He along the transmission pathway (z-coordinates) to both membranes are illustrated in Fig. 2, while the contour plots in Fig. 3(a) and (b) highlight the in-pore displacements of He along x- and y-coordinates. It is quite evident from both these plots that the minimum energy path (MEP) for an effective He transmission may correspond to a straight line perpendicular to the center of the pore. For the transmission through the Gr2 membrane, the He atom lying at the center of a pore is the saddle point. The maximum potential barrier height (i.e. energy difference between the transition state, TS, (z = 0.0 Å) and at the asymptote (for example, z = 6.0 Å)) for the MEP is 37.85 meV, which is similar to the previous studies.^37,49,50 In contrast, MEP for the He–Gr3 PES is devoid of any potential barrier; instead, a large well (of depth 17.74 meV) appears within the vicinity of the membrane.⁵¹ These results suggest that He penetration of larger Gr3 pores should be much easier than that of smaller Gr2 pores, analogous to the report on water permeation by Bartolomei et al.²⁶ Furthermore, along this MEP and He atom–membrane perpendicular distance of around 6.0 Å, the interaction potential reaches a plateau and therefore can be characterized as an asymptotic reactant site. It is interesting to note that the interaction potential steeply rises to a high value for any movement along the in-pore x – and y – degrees of freedom [see Fig. 3(a) and (b)].

Fig. 2 Variation of the potential energies (in meV) for (a) He atom–graphdiyne and (b) He atom–graphtriyne interaction along the minimum energy transmission pathway (MEP) (in Å). The interaction energies are calculated by using improved Lennard-Jones (ILJ) potentials. The MEP corresponds to He atom perpendicularly approaching the geometric center of graphdiyne and graphtriyne pores with reaction coordinates (0,0, z).

Fig. 3 Improved Lennard-Jones (ILJ) interaction potentials (in meV) for the displacement of He atom along x and y directions from the center of the pore (0,0,0) on (a) graphdiyne and (b) graphtriyne membranes. The z coordinate of He atom was held at the origin. The energy contour step value is 50 meV and 10 meV for (a) and (b), respectively.

2.2 Ring polymer molecular dynamics (RPMD) method

The classical Hamiltonian for the system composed of ⁴He or ³He atom under the influence of external potential arising from M carbon atoms of fixed Gr2 or Gr3 membrane can be written in the atomic unit as:where, p̂ and r̂_He are the momentum and position vector of the He atom of mass m, and r̂_C(i)'s are the position vectors of the C atoms of Gr2 or Gr3. In the RPMD method, the He atom is treated as n classical replicas of the original particle, each connected with its nearest neighbor by harmonic springs. The modified ring polymer Hamiltonian has the form:where,n is the number of classical beads representing quantum He atom which are connected by a harmonic potential with force constant ω_n (= β_nħ). β_n ≡ β/n is the reciprocal temperature of the system, β = 1/k_BT. k_B is the Boltzmann constant, T is the system temperature, and ħ is the Dirac constant (ħ = 1 in atomic units). p^(j) and r^(j) are the momentum and position vectors of the jth bead in the ring polymer necklace of He atom, respectively. The ring polymer trajectory now evolves in f = 3n total degrees of freedom (in atomic Cartesian coordinates and including translation and rotational degrees of the freedom of the entire system which is a convenient method to propagate RPMD trajectories⁷¹).

We perform the RPMD simulations at temperatures T below 250 K because, previously, it was observed that the maximum selectivity was obtained at low temperatures.^49,50 We choose eight different T = 10, 20, 30, 50, 100, 150, 200, and 250 K, to study the variation of selectivities with T and to find an optimal temperature range for isotopic separation. Furthermore, the number of beads is a measure of the resolution of the path integral calculations, i.e., RPMD calculations scale linearly with n.⁵⁶ However, since RPMD calculations are approximately n times slower than the purely classical calculations,⁵⁶ after several trials, we have meticulously chosen the number of beads to be 128. This value of n gives a fair compromise between computational cost and accuracy (quantum effects of ZPE and tunneling) at low temperatures. Since the crossover temperature T_c(= ħω_i/2πk_B, where ω_i is the imaginary frequency at the TS) on the Gr2–He PES is found to be 17 K, in order to obtain the reliability of the RPMD simulations at very low temperature, calculations at 10 K have been performed with a greater number of beads (n = 256). The value of n = 1, will correspond to a classical calculation.

Since RPMD is simply classical molecular dynamics in an extended (n-bead imaginary time path integral) phase space, the ring polymer rate coefficient can be expressed as:^53,54,71

where, Q⁽ⁿ⁾_r(T) is the n-bead path integral approximation to the quantum mechanical partition function of the reactants per unit volume, and c̃⁽ⁿ⁾_fs(t) is a ring polymer flux-side correlation function⁷⁸

Here, subscript 0 and t indicates time, δ(r₀) is a delta function centered at r₀, v(r₀,p₀) is the velocity, and h(r_t) is the Heaviside step function. The RPMD rate coefficient in eqn (4) is not straightforward to solve numerically. Therefore, we introduce the Bennett–Chandler factorization scheme^79,80 to simplify eqn (4) that can be solved numerically without compromising its generality.⁸⁰ The method has been extensively discussed previously^54,70,71,81 and will not be repeated here. Briefly, in this approach, a reaction coordinate s(r) is defined, which monitors the progress of a reaction from the reactant (s > 0) to the product (s < 0) site. For the He transmission through 2D membranes, a reasonable reaction coordinate is s(r) = z̄, where z̄ is the z-component of the He atom centroid, which follows the minimum energy path of Fig. 2. Within the Bennett–Chandler approach, the RPMD rate coefficient for a process in which reactants and products separated by a dividing surface at s (in this work, at the center of the pore, s = 0.0 Å) can be expressed as a product of two terms:

k_RPMD(T) = k_QTST(T)κ(t_p). (6)

Here, the first factor, k_QTST(T), is the centroid-density quantum transition-state theory (QTST) rate coefficient.⁵⁴k_QTST(T) can be calculated from the centroid potential of mean force (PMF),^70–72W(s) along the reaction coordinate. If s_∞ is the asymptotic distance in which the He – membrane interaction potential is negligible and then introduce a dividing surface carefully placed at the TS region, s^‡, then k_QTST(T) can be calculated as:^66,81

W(s)'s can be computed by employing the umbrella integration procedure of Kästner and Thiel.^82–84 To calculate the PMF profiles, the reaction coordinate s of the He atom has been divided into 130 equally spaced windows (of width 0.05 Å) within the range from 6.00 Å (reactant site, s_∞) to −1.00 Å (product site). The PMFs in eqn (7) is then calculated as:with

Here, N_windows is the number of biasing windows placed along the reaction coordinate. The strength of the force constant (k_i) of the harmonic biasing potential was chosen to be 2.72 × 10⁻³T (K) eV a₀⁻². N_i is the total number of steps sampled for window i, s̄_i and σ_i² are the mean value and variance for the trajectory calculated for the ith window. In each umbrella sampling windows, 100 trajectories with different initial configurations were propagated for 100 ps following an initial equilibration period of 20 ps in the presence of an Andersen thermostat.⁸⁵ The ring polymer equation of motion were integrated using velocity Verlet integrator with a step size of 0.1 fs that involves alternating momentum updates and free ring polymer evolutions.⁶⁰

The second term in eqn (6), κ(t_p), is the long-time limit of a time-dependent ring polymer transmission coefficient or the ring polymer recrossing factor^65,71,72,81 and is a dynamic correction to k_QTST. Typically this factor is calculated at the top of the free energy barrier on the PMF profile so as to minimize the time required to reach plateau value.⁸⁰ This factor ensures that the final k_RPMD value is independent of the choice of the dividing surface.^54,56 Therefore, finding an optimum location for transmission coefficient calculation is not essential within the RPMD formalism; a close approximation of the TS (for example, s = 0.0 Å in Fig. 4 and 5) will suffice for a smooth and efficient convergence of κ(t_p). The mathematical expression for κ(t_p) can be written as:^66,72

where, δ[s(r₀)] constrains the initial configurations to the dividing surface, ṡ(r₀) is the velocity factor that accumulates the flux through the dividing surface, and h[s(r_t)] is a Heaviside function that gathers trajectories that have crossed over to the product side of the dividing surface. h[ṡ(r₀)] is basically a normalization factor that ensures κ(t → 0₊) = 1. To calculate the recrossing factor, a long “parent” trajectory for He atom of length 2 ns has been carried out after an initial thermalization period of 20 ps in the presence of Andersen thermostat with its centroid pinned at the dividing surface using RATTLE algorithms.⁸⁶ After each 2 ps propagation period of the parent trajectory, 100 trajectories have been generated that have initial position of the parent trajectory, but their momenta is randomly generated from a Boltzmann distribution. These “child” trajectories are then propagated for 1 ps in the absence of thermostat and dividing surface constraints.

Fig. 4 Variation of the RPMD potential of mean force, W(s), (in meV) for (i) ³He (ii) ⁴He atom along the reaction coordinate s (in Å) perpendicular to the graphdiyne membrane within the temperature range (A) 10–250 K, (B) 20–30 K, and (C) 200 K. The legends correspond to (A)–(C).

Fig. 5 Variation of the RPMD potential of mean force, W(s), (in meV) for (i) ³He (ii) ⁴He atom along the reaction coordinate s (in Å) perpendicular to the graphtriyne membrane within the temperature range (A) 10–250 K, (B) 20–30 K, and (C) 250 K. The legends correspond to (A)–(C).

3 Results and discussion

3.1 Centroid potentials of mean force and ring polymer recrossing factor

The variation of the RPMD potential of mean force W(s) for ³He and ⁴He along the reaction coordinate s at T = 10–250 K are plotted in Fig. 4 (Gr2) and Fig. 5 (Gr3), and the corresponding k_QTST values are reported in the Table S1 of ESI.† These PMF profiles include both potential energy and temperature-dependent entropic contributions. The barrier height at the TS, identified around the membrane plane, increases with increasing temperature. The free energy profiles enter a shallow well at around s = 2.0–3.4 Å from the reactant site before it steadily increases up to the TS. For the PMFs obtained for the transmission through Gr2, the calculated TS barrier height falls within 37–127 meV. It is evident that the PMF profiles at 10 K are markedly distinct from those obtained at other temperatures. At 10 K, the free energy profiles for both isotopes are almost flat within the vicinity of the pore (s = ±0.6 Å). The corresponding barrier heights at the TS are also significantly low (37–40 meV) compared to those obtained at 20 K (64–66 meV). The PMF profiles on Gr3, on the other hand, are more interesting from thermodynamic point of view. We recall that the He–Gr3 interaction potential does not have any barrier;⁵¹ rather, a well is built around the transmission zone. This behavior is reflected in the free energy profiles, particularly within the temperature range 10–30 K, in which the thermodynamic barrier in the TS region (at s = 0.0 Å) has a negative (or marginally positive) value compared to the reactant site. This energy barrier gradually takes the shape of a parabola with increasing temperature and reaches up to a height of 44 meV at 250 K.

Since the RPMD method averages points in configurational space on both sides of the potential energy barrier, modeling the effects of the tunneling and ZPE and leaving a mark on the free energy profiles,⁵⁹ we can roughly attempt to identify these properties. Comparing the PMF profiles of ⁴He at a particular temperature with the corresponding ³He's on Gr2 indicate the existence of a prominent ZPE effect; heavier isotope has a lower free energy at the TS than the lighter ones except at 10 and 20 K (see Fig. 4(B)), in which quantum tunneling may play a crucial role. In particular, at 10 K, the barrier height for ³He transmission is around 3 meV lower than that for ⁴He. Interestingly, at that temperature, the free energy barrier (37 meV) for ³He is marginally smaller than the pure potential energy at 0 K (38 meV), suggesting that the transmission of the lighter isotope may be driven by quantum mechanical tunneling. However, the difference between them at TS decreases with increasing temperature. For example, the difference at 30 K is 2.3 meV (see Fig. 4(B)), which decreases to 0.1 meV at 200 K (see Fig. 4(C)). Similarly, the ³He and ⁴He free energy comparison for Gr3 follows the same trend observed for Gr2, i.e., they decrease with increasing temperature. Within 10–30 K, ⁴He has a smaller TS free energy than ³He (≈1.0 meV, Fig. 5(B)), and with the increase in temperature in the range 50–150 K, the heavier isotope progressively requires more free energy to reach TS, i.e., PMF profiles for both isotope are almost identical. Within 200–250 K, the free energies of ⁴He are consistently larger than those obtained for ³He (see Fig. 5(C)). These observations can be notionally interpreted as follows: within 10–30 K, the free energy barrier is too broadened for effective tunneling, and consequently, the ZPE effect takes precedence.⁵⁷ There is a competition between them at moderate temperatures (50–150 K). At still more elevated temperature (200–250 K), the greater tunneling probability of ³He may facilitate a smaller TS free energy compared to ⁴He. Note that this analysis is conceptual because these two quantum effects cannot be rigorously separated within the RPMD formalism.

There are considerable differences observed between the PMF profiles obtained by the RPMD method with the classical ones at low temperatures (see Fig. S1–S4 of ESI†). This, therefore, reinforces our earlier argument of the existence of quantum effects in these transmission processes, which the classical calculation fails to capture. For the He transmission through Gr2 membrane, the free energy barrier height in the classical calculations is always lower than those obtained by the RPMD method, except at 10 K. Moreover, the flat region in the free energy profile in the RPMD simulations found around the TS at 10 K is absent in the classical calculations. The maximum difference between classical and RPMD barriers is ≈10 meV for ³He (at 100 K) and ≈12 meV for ⁴He (obtained at 20 K). Note that the free energy values exponentially contribute to the rate coefficients. As expected, these PMF plots tend to merge with increasing temperature as the quantum effects fade. Comparison between the classical and RPMD method on Gr3 shows that, for both isotopes, the maximum free energy difference at the TS is found at the lowest temperature (≈4 meV), and this difference smoothly decreases with increasing temperature. It is worth noting that the classical PMFs within the 10–20 K temperature range do not have any barrier around the center of the pore. For ⁴He within 200–250 K, the difference between the RPMD and classical barrier heights becomes 1 meV.

The time-dependent transmission coefficients, κ(t), for all temperatures considered in this study are calculated at s^‡ = 0.0 Å. The plots of κ(t) against time (t) are moved to the ESI† (Fig. S5 and S6). The corresponding plateau value of the transmission coefficients, κ(t → ∞), are provided in Table S1 (ESI†). The plateau value κ(t) on both membranes is always close to unity (0.98–1.00), except at 10 K on Gr2, where κ(t) falls within 0.63 (³He)–0.78 (⁴He). The greater recrossings of the lighter isotope at 10 K may be attributed to the fact that the slowly moving ³He can more readily tunnel back and forth through the potential barrier. On the other hand, at 20 K, κ(t) for ³He is marginally smaller than that obtained for ⁴He (0.98 and 0.99 respectively). Additionally, a longer propagation time is required to converge the transmission coefficients at 10 K (2 ps) than at any other temperatures (around 500 fs). Clearly, for temperatures above 10 K, recrossing dynamics do not play any significant role for the He transmission on both membranes, particularly from 50 K and higher temperature, where κ(t) value is always 1.00. This implies that most of the He trajectories that reach the center of the pore of either Gr2 or Gr3 would overcome the free energy barrier and transport to the other side of the dividing surface. On Gr3, if the recrossing trajectories are propagated for a more extended period of time (for example, 10 ps), then after around 1.1 ps (³He) or 1.2 ps (⁴He), κ(t) decays rapidly and eventually settle to a value close to zero (see Fig. S7 of ESI†). This decay is more prominent at low temperatures (10–20 K) than at higher temperatures (for example, 250 K). This phenomenon is apparently due to a “secondary event”, such as the He atom recrosses the dividing surface while performing “another” transmission.⁶⁶ This exponential decay in the κ(t) values is missing for the transmission through Gr2; therefore, they may also be linked to the shape of the PMFs on Gr3 at low temperatures, which is devoid of a substantial free energy barrier and consequently easier to perform such motions. The signature of the “secondary event” on the κ(t) profile can be effectively eliminated by placing the dividing surface, s^‡ around the asymptotic reactant site. However, the convergence of κ(t) at such separation would require a very large number of “child” trajectories to be propagated for an extended period of time (hundreds of ps for the present cases), thus incurring a substantial computational burden. As discussed previously, since the RPMD rate coefficient is independent of the position of this dividing surface; therefore, determination of κ(t) at the asymptotic reactant site would be disadvantageous in the present case. Note that the classical recrossing factor is always 1, even at low temperatures.

3.2 Thermal rate coefficient

Since the plateau of κ(t) is always close to unity, there is practically no significant difference between the values of k_QTST and k_RPMD on both membranes, except at 10 K on Gr2. The variation of k_RPMD with temperature is plotted in Fig. 6 and Fig. S8 of ESI,† while the numerical values of k_RPMD and k_QTST are reported in Table S1 of ESI.† It is evident that k_RPMD's on Gr3 is many orders of magnitude greater than that obtained on Gr2. The maximum difference was observed at the lowest temperature (by a factor of 10²³–10²⁰ at 10 K), and with the increase in temperature, this difference decreases rapidly (to 10² at 250 K) as the k_RPMDvs. T curves start to converge. This is due to the fact that the value of k_RPMD on Gr3 does not change drastically with temperature (confined within 3.63 × 10¹⁰–1.06 × 10¹¹ s⁻¹). It starts to decrease slightly with temperature up to 50 K and then increases successively with increasing T. On Gr2, however, for both isotopes, k_RPMD increases manifold with temperature, in agreement with the previous calculations.^37,49,50 This increment in the rate coefficient is more pronounced at the modest rise in T in the low temperature regime than at the high temperature range. The increase in the value of k_RPMD is in the order of 10⁴–10⁶ for the temperature rise from 10 K to 20 K, while for other cases (for example, 20 K to 30 K or 30 K to 50 K), the same raises by five orders of magnitude. However, k_RPMD increases much less than an order for the consecutive 50 K temperature jumps, starting at 150 K.

Fig. 6 Variation of the ring polymer molecular dynamics, RPMD (k_RPMD) and classical (k_cl) rate coefficients (in s⁻¹) for the transmission of ³He [red circle (RPMD) and green diamond (classical)] and ⁴He [blue square (RPMD) and orange pentagon (classical)] through (a) graphdiyne and (b) graphtriyne membranes with temperature T (in K). The legends correspond to both (a) and (b).

The above observations can be explained by inspecting the corresponding PMF profiles. At low temperatures (10–50 K), there is practically no free energy barrier for the transmission of either isotope on Gr3. However, on Gr2, the free energy barrier is already about 40 meV at 10 K. The slowly moving He atoms do not have enough kinetic energy to overcome this barrier to reach the other side of the membrane. Therefore, it is not surprising that k_RPMD's on Gr3 is many orders of magnitude greater than those obtained on Gr2, and the actual He flux through the pores of Gr2 membranes will be extremely slow.^49,50 Similar arguments can be presented at higher temperatures, although the increased kinetic energy of the incoming He atom will contribute to smaller differences in the k_RPMD values. We also point out the analogous quantum wave packet observation reported on the Gr2 and holey graphene (P7) sheet having a more diffused pore.⁵⁰ Finally, when comparing the classical and RPMD rate coefficients, it is obvious that the classical method overestimates the rate coefficient within the whole temperature regime studied in this work, except at 10 K for Gr2, where they underestimate by factor of 10⁴–10⁶. This discrepancy is particularly more apparent at low temperatures, when the quantum effects dominate, and on the Gr2 membrane. For example, at 20 K and for Gr2 membrane, the classical rate coefficient is more than three orders of magnitude greater than the corresponding k_RPMD. However, they do closely follow a similar temperature dependence as obtained by the RPMD calculations.

3.3 ⁴He/³He selectivity

The ⁴He/³He selectivity, defined as the ratio between the RPMD rate coefficient of the heavier He isotope to the lighter one, k_RPMD(⁴He)/k_RPMD(³He), is plotted as a function of temperature in Fig. 7, and the corresponding values are reported in Table S1 of ESI.† For Gr2 membrane, ⁴He/³He selectivity increases for the temperature rise from 10 K to 30 K. At 10 K, the calculated selectivity corresponds to the lowest value found in this study (0.013). Conversely, it signifies that the rate of permeation of the lighter isotope is about 77 times faster than the heavier isotope due to the predominant tunneling of ³He atoms, as elaborated in Section 3.1. Therefore, it represents an optimum temperature for isotope separation on Gr2. Note that such high selectivity on either Gr2 or Gr3 membrane has not been reported in the literature. There is no significant preference for either isotope found at 20 K. While, at 30 K, the maximum selectivity for the heavier isotope is obtained (≈2), indicating pronounced quantum ZPE effect as discussed previously. However, these increased selectivities at low temperatures are usually accompanied by a decreased permeability.⁸⁷ With the increase in temperature, this selectivity ratio progressively becomes smaller and starts to flatten out, commencing from 100 K. Although at higher temperatures (>100 K), the transmission of the lighter isotope is favoured to some extent. We note the analogous selectivity profile was obtained on the P7 sheet.⁵⁰ Similarly, for the Gr3 membrane, the maximum ⁴He/³He selectivity obtained at 10 K (1.25), which gradually decreases till 50 K, and then remains almost constant (0.87 to 0.83) marginally preferring ³He for higher temperature regime. The variation of the selectivity is fully consistent with the nature of their PMF profiles that inherits quantum effects. It is interesting to note that the selectivity plot for both Gr2 and Gr3 approach each other with increasing T, one from the top and another from the bottom (see Fig. S9 of ESI†), and at 250 K, they almost become equal (≈0.8 and partially in favour of ³He). This is reasonable since as the quantum effects diminish with temperature and with fewer constraints, the lighter isotope transmission will be slightly favoured, governed by the factor (≈ 0.87).

Fig. 7 Variation of the ⁴He/³He rate coefficient ratio, k(⁴He)/k(³He), calculated using (i) classical (red filled circle) (ii) ring polymer molecular dynamics, RPMD (blue filled square) and (iii) three-dimensional wave packet propagation method, WP3D (green solid line) for the He transmission through (a) graphdiyne and (b) graphtriyne membranes with temperature T (in K). The legends correspond to both (a) and (b).

The selectivity ratios obtained by the RPMD method show good agreement with those obtained by the quantum three-dimensional wave packet propagation (WP3D) results obtained by Gijón et al.⁵⁰ and Hernández et al.⁵¹ for the whole temperature range studied in this work. In general, the RPMD selectivity profile closely resembles the WP3D results; however, they do overestimate and underestimate to some degree. For example, on Gr2, the RPMD selectivity ratio overestimates the quantum calculations for the whole temperature range. The maximum difference between these two methods is found at 30 K (0.2). Similarly, on Gr3, the RPMD method seems to overestimate the selectivity ratio within 20–250 K, but have an excellent agreement at 10 K (1.25 and 1.30, respectively). Overall, the difference between RPMD and WP3D results for the selectivity ratio lies within 9–23% on Gr2 and 4–36% on Gr3 in comparison to the quantum calculation. The discrepancy between RPMD and WP3D selectivity ratios on Gr3, particularly at low temperatures, may be attributed to the formation of long-lived resonances within the potential well, which the RPMD method is not completely able to capture.⁵⁶ On the other hand, the classical selectivity shows no significant change with temperature and remains almost constant (0.80–0.91), slightly favoring ³He transmission. This is valid for both membranes. As expected, with the rise in temperature, the classical selectivity follows the RPMD and WP3D ones.

A natural question arises regarding the effect of the motion of the carbon atoms of Gr2/Gr3 on the ³He/⁴He permeation rates and the selectivity ratio. One way to foresee the importance of the membrane motion is to study the modifications that the penetrating atom produces on the membrane degrees of freedom at the atom–membrane transition state. It is known that the He atoms typically induce minute changes in the geometries of the pores, producing a small decrease in the barrier height upon relaxation of the membrane (of the order of few cm⁻¹) as found, for instance, in the electronic energy calculations of the interaction between He with nanoporous graphene¹⁶ or for this very same system (see Fig. S1 of ref. 37). For this reason, it is not expected that the inclusion of the motion of the C atoms will significantly alter the results reported in this paper. However, it would be interesting to evaluate the effect of the layer motion on the RPMD rate coefficient provided that an adequate C–C force field for these membranes is available. Other processes at the gas-membrane interface do require an explicit treatment of the carbon atoms dynamics since they involve large deformations of the carbon framework, including rehybridization and formation of new bonds.^8,88

4 Conclusions

In this work, we have calculated the thermal rate coefficient for the transmission of He isotopes through the pores of one atom thick graphdiyne (Gr2) and graphtriyne (Gr3) membranes using the ring polymer molecular dynamics (RPMD) method within the temperature range 10–250 K. Transmission through Gr2 has a substantial free energy barrier even at 10 K, and this barrier height increases with increasing temperature. Although, the free energy barrier at 10 K is smaller than the potential energy, signifying the presence of quantum mechanical tunneling. On the other hand, transmission through Gr3 can either have marginally negative (up to 30 K) or small positive (T ≥ 50 K) free energy barrier. The extent of the barrier height directly impacts the calculated rates, as evident from the fact that the rate coefficient on Gr3 is at least 10²⁰ order of magnitude greater than on Gr2 at 10 K. The rate coefficient on Gr2 increases rapidly with temperature till 150 K (∼10⁸ s⁻¹); thereafter, it increases moderately. However, the rate coefficients on Gr3 do not vary appreciably with temperature and remain almost constant for the whole temperature range (3.63 × 10¹⁰–1.06 × 10¹¹ s⁻¹). In general, the rate coefficient calculated for the transmission through Gr3 is always greater than the corresponding one on Gr2 over the whole temperature regime considered in this work. Moreover, we found that the recrossing dynamics have little or no effect on the final value of the rate coefficients for temperatures above 10 K.

The selectivity ratio, which indicates the preference of either isotope for its permeation through the membranes, has been calculated as a function of temperature. From the values of the selectivity ratio, the quantum effects driving the He transmission, particularly at low temperatures, i.e., the zero-point energy (ZPE) favoring the heavier isotope and the tunneling of the more mobile lighter isotope, has been coarsely deduced. The maximum selectivity ratio found in this study was obtained on Gr2 membrane at 10 K in which the rate of transmission of ³He is almost 77 times the rate of ⁴He, seemingly due to a more dominant quantum mechanical tunneling. On the other hand, presumably as a result of the enhanced ZPE effect, the transmission of the heavier isotope is favored by a factor of ≈2 at 30 K. A similar conclusion can be derived on Gr3, where the ⁴He/³He selectivity is around 1.2 within the temperature range 10–20 K. However, on the Gr3 membrane, the selectivity ratio does not vary considerably with the temperature (0.8–1.2) as compared to what was observed for the transmission through the Gr2 membrane (0.01–2.0). With the increase in temperature, the permeation of the lighter isotope is marginally favored. The RPMD selectivity ratio is consistent with the quantum calculations for the entire temperature range studied in this work.^50,51 On the other hand, the classical method failing to capture either of the quantum effects demonstrated a considerable discrepancy with the RPMD and quantum results, particularly at low temperatures. Therefore, it is of paramount importance to use robust and accurate methods, such as RPMD, that can correctly describe quantum effects such as ZPE and tunneling when studying physical processes for which strong quantum nature is expected.

In conclusion, in the present study, we have corroborated the efficient and rigorous nature of the RPMD method in a quantum sieving process within low to ambient temperature regime by determining thermal rate coefficients of physical processes of broad industrial and scientific significance. Previously unreported, an optimum temperature (10 K) has been found in this study that can be exploited for an efficient isotopic separation. We also hope that this work will stimulate future experimental measurements of the rate coefficients for the separation of He isotopes using graphene derivatives, taking advantage of quantum effects at low temperatures. As an extension to this work, we would like to investigate the influence of explicit motions of the membrane atoms as well as of surrounding He atoms on the RPMD rate coefficient. In such cases implementing QM methods is practically futile. It is to be noted that since the RPMD method does not take into account of the nuclear spin statistics, the bosonic and fermionic nature of the He isotopes will not be considered for a system consisting of many atoms. Despite this drawback, the superfluidity (formation of Bose–Einstein condensate) of the bosonic and the fermionic He initiates at far lower critical temperature or λ point (2.17 K⁸⁹ and below 2.7 mK,⁹⁰ respectively) than the proposed temperature regime (10–250 K). Recently, Parrinello and co-workers have shown that one can overcome the distinguishable nature and the explicit disregard of bosonic and fermionic exchange effects in path integral methods by avoiding the permutations of all identical particles.^91,92 We also envisage to apply the RPMD method for isotopic He separation on other graphene derivatives such as polyphenylene, functionalized graphene pores, holey graphene, graphenylene membranes, etc. that were previously reported to serve as excellent atomic sieves.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Y. V. S. and S. B. acknowledge the support of the European Regional Development Fund and the Republic of Cyprus through the Research Promotion Foundation (Projects: INFRASTRUCTURE/1216/0070 and Cy-Tera NEA YPODOMH/STPATH/0308/31). Y. V. S. was also supported by the Russian Foundation for Basic Research grant number 20-03-00833. J. C.-M. and M. I. H. were supported by Spanish MICINN with Grant FIS2017-84391-C2-2-P and PID2020-114957GB-I00. We are also indebted to “Centro de Supercomputación de Galicia” (CESGA) and specially to the people in the help desk.

References

1. K. Khan, A. K. Tareen, M. Aslam, R. Wang, Y. Zhang, A. Mahmood, Z. Ouyang, H. Zhang and Z. Guo, J. Mater. Chem. C, 2020, 8, 387–440.
2. C. Huang, Y. Li, N. Wang, W. Xue, Z. Zuo, H. Liu and Y. Li, Chem. Rev., 2018, 118, 7744–7803.
3. J. S. Bunch, S. S. Verbridge, J. S. Alden, A. M. Van Der Zande, J. M. Parpia, H. G. Craighead and P. L. McEuen, Nano Lett., 2008, 8, 2458–2462.
4. V. Berry, Carbon N. Y., 2013, 62, 1–10.
5. S. Hu, M. Lozada-Hidalgo, F. C. Wang, A. Mishchenko, F. Schedin, R. R. Nair, E. W. Hill, D. W. Boukhvalov, M. I. Katsnelson, R. A. W. Dryfe, I. V. Grigorieva, H. A. Wu and A. K. Geim, Nature, 2014, 516, 227–230.
6. M. Lozada-Hidalgo, S. Hu, O. Marshall, A. Mishchenko, A. N. Grigorenko, R. A. Dryfe, B. Radha, I. V. Grigorieva and A. K. Geim, Science, 2016, 351, 68–70.
7. S. Bukola and S. E. Creager, Electrochim. Acta, 2019, 296, 1–7.
8. J. W. Mazzuca and N. K. Haut, J. Chem. Phys., 2018, 148, 224301.
9. M. Miao, M. B. Nardelli, Q. Wang and Y. Liu, Phys. Chem. Chem. Phys., 2013, 15, 16132–16137.
10. M. Bartolomei, M. I. Hernández, J. Campos-Martínez and R. Hernández-Lamoneda, Carbon N. Y., 2019, 144, 724–730.
11. P. Z. Sun, Q. Yang, W. J. Kuang, Y. V. Stebunov, W. Q. Xiong, J. Yu, R. R. Nair, M. I. Katsnelson, S. J. Yuan, I. V. Grigorieva, M. Lozada-Hidalgo, F. C. Wang and A. K. Geim, Nature, 2019, 579, 229–232.
12. L. Tsetseris and S. T. Pantelides, Carbon N. Y., 2014, 67, 58–63.
13. X. Zhao, S. Villar-Rodil, A. J. Fletcher and K. M. Thomas, J. Phys. Chem. B, 2006, 110, 9947–9955.
14. S. P. Koenig, L. Wang, J. Pellegrino and J. S. Bunch, Nat. Nanotechnol., 2012, 7, 728–732.
15. D. E. Jiang, V. R. Cooper and S. Dai, Nano Lett., 2009, 9, 4019–4024.
16. A. W. Hauser and P. Schwerdtfeger, J. Phys. Chem. Lett., 2012, 3, 209–213.
17. R. K. Joshi, P. Carbone, F. C. Wang, V. G. Kravets, Y. Su, I. V. Grigorieva, H. A. Wu, A. K. Geim and R. R. Nair, Science, 2014, 343, 752–754.
18. C. Owais, A. James, C. John, R. Dhali and R. S. Swathi, J. Phys. Chem. B, 2018, 122, 5127–5146.
19. E. Griffin, L. Mogg, G.-P. Hao, G. Kalon, C. Bacaksiz, G. Lopez-Polin, T. Zhou, V. Guarochico, J. Cai, C. Neumann, A. Winter, M. Mohn, J. H. Lee, J. Lin, U. Kaiser, I. V. Grigorieva, K. Suenaga, B. Özyilmaz, H.-M. Cheng, W. Ren, A. Turchanin, F. M. Peeters, A. K. Geim and M. Lozada-Hidalgo, ACS Nano, 2020, 14, 7280–7286.
20. O. Leenaerts, B. Partoens and F. M. Peeters, Appl. Phys. Lett., 2008, 93, 193107.
21. J. Yeo, G. S. Jung, F. J. Martín-Martínez, J. Beem, Z. Qin and M. J. Buehler, Adv. Mater., 2019, 1805665, 1–24.
22. J. Kang, Z. Wei and J. Li, ACS Appl. Mater. Interfaces, 2019, 11, 2692–2706.
23. R. H. Baughman, H. Eckhardt and M. Kertesz, J. Chem. Phys., 1987, 87, 6687–6699.
24. G. Li, Y. Li, H. Liu, Y. Guo, Y. Li and D. Zhu, Chem. Commun., 2010, 46, 3256–3258.
25. J. Zhou, X. Gao, R. Liu, Z. Xie, J. Yang, S. Zhang, G. Zhang, H. Liu, Y. Li, J. Zhang and Z. Liu, J. Am. Chem. Soc., 2015, 137, 7596–7599.
26. 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.
27. S. W. Cranford and M. J. Buehler, Carbon N. Y., 2011, 9, 4111–4121.
28. A. N. Enyashin and A. L. Ivanovskii, Phys. Status Solidi Basic Res., 2011, 248, 1879–1883.
29. S. W. Cranford, D. B. Brommer and M. J. Buehler, Nanoscale, 2012, 4, 7797–7809.
30. Y. Yang and X. Xu, Comput. Mater. Sci., 2012, 61, 83–88.
31. S. Lin and M. J. Buehler, Nanoscale, 2013, 5, 11801–11807.
32. Y. Jiao, A. Du, M. Hankel, Z. Zhu, V. Rudolph and S. C. Smith, Chem. Commun., 2011, 47, 11843–11845.
33. L. Zhao, P. Sang, S. Guo, X. Liu, J. Li, H. Zhu and W. Guo, Appl. Surf. Sci., 2017, 405, 455–464.
34. Z. Meng, X. Zhang, Y. Zhang, H. Gao, Y. Wang, Q. Shi, D. Rao, Y. Liu, K. Deng and R. Lu, ACS Appl. Mater. Interfaces, 2016, 8, 28166–28170.
35. C. Zhu, H. Li, X. C. Zeng, E. G. Wang and S. Meng, Sci. Rep., 2013, 3, 3163.
36. W. J. Nuttall, R. H. Clarke and B. A. Glowacki, Nature, 2012, 485, 573–575.
37. M. Bartolomei, E. Carmona-Novillo, M. I. Hernández, J. Campos-Martínez, F. Pirani and G. Giorgi, J. Phys. Chem. C, 2014, 118, 29966–29972.
38. Atomic Weights and Isotopic Compositions for Helium, https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Heisotype=all.
39. N. K. Das, H. Chaudhuri, R. K. Bhandari, D. Ghose, P. Sen and B. Sinha, Curr. Sci., 2008, 95, 1684–1687.
40. P. Bernardo, E. Drioli and G. Golemme, Ind. Eng. Chem. Res., 2009, 48, 4638–4663.
41. S. Blankenburg, M. Bieri, R. Fasel, K. Müllen, C. A. Pignedoli and D. Passerone, Small, 2010, 6, 2266–2271.
42. J. Schrier, J. Phys. Chem. Lett., 2010, 1, 2284–2287.
43. J. Schrier and J. McClain, Chem. Phys. Lett., 2012, 521, 118–124.
44. A. M. Brockway and J. Schrier, J. Phys. Chem. C, 2013, 117, 393–402.
45. A. W. Hauser, J. Schrier and P. Schwerdtfeger, J. Phys. Chem. C, 2012, 116, 10819–10827.
46. M. Lalitha, S. Lakshmipathi and S. K. Bhatia, J. Phys. Chem. C, 2015, 119, 20940–20948.
47. X. Liu, X. Chang, L. Zhu and X. Li, Comput. Mater. Sci., 2019, 157, 1–5.
48. S. Mandrà, J. Schrier and M. Ceotto, J. Phys. Chem. A, 2014, 118, 6457–6465.
49. M. I. Hernández, M. Bartolomei and J. Campos-Martínez, J. Phys. Chem. A, 2015, 119, 10743–10749.
50. A. Gijón, J. Campos-Martínez and M. I. Hernández, J. Phys. Chem. C, 2017, 121, 19751–19757.
51. M. I. Hernández, M. Bartolomei and J. Campos-Martínez, Nanomaterials, 2021, 11, 73.
52. Y. Qu, F. Li and M. Zhao, Phys. Chem. Chem. Phys., 2017, 19, 21522–21526.
53. I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2005, 122, 084106.
54. I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2005, 123, 034102.
55. S. Habershon, D. E. Manolopoulos, T. E. Markland and T. F. Miller, Annu. Rev. Phys. Chem., 2013, 64, 387–413.
56. Y. V. Suleimanov, F. Javier Aoiz and H. Guo, J. Phys. Chem. A, 2016, 120, 8488–8502.
57. R. Pérez De Tudela, F. J. Aoiz, Y. V. Suleimanov and D. E. Manolopoulos, J. Phys. Chem. Lett., 2012, 3, 493–497.
58. J. O. Richardson and S. C. Althorpe, J. Chem. Phys., 2009, 131, 214106.
59. R. Pérez De Tudela, Y. V. Suleimanov, J. O. Richardson, V. Sáez Rábanos, W. H. Green and F. J. Aoiz, J. Phys. Chem. Lett., 2014, 5, 4219–4224.
60. I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2004, 121, 3368–3373.
61. T. F. Miller and D. E. Manolopoulos, J. Chem. Phys., 2005, 122, 184503.
62. I. R. Craig and D. E. Manolopoulos, Chem. Phys., 2006, 322, 236–246.
63. T. F. Miller and D. E. Manolopoulos, J. Chem. Phys., 2005, 123, 154504.
64. J. W. Mazzuca and C. P. Schultz, J. Phys. Chem. A, 2017, 121, 819–826.
65. T. E. Markland, S. Habershon and D. E. Manolopoulos, J. Chem. Phys., 2008, 128, 194506.
66. Y. V. Suleimanov, J. Phys. Chem. C, 2012, 116, 11141–11153.
67. A. R. Menzeleev, N. Ananth and T. F. Miller, J. Chem. Phys., 2011, 135, 074106.
68. J. S. Kretchmer and T. F. Miller, J. Chem. Phys., 2013, 138, 134109.
69. N. Boekelheide, R. Salomón-Ferrer and T. F. Miller, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 16159–16163.
70. R. Collepardo-Guevara, Y. V. Suleimanov and D. E. Manolopoulos, J. Chem. Phys., 2009, 130, 174713.
71. Y. V. Suleimanov, R. Collepardo-Guevara and D. E. Manolopoulos, J. Chem. Phys., 2011, 134, 044131.
72. Y. V. Suleimanov, J. W. Allen and W. H. Green, Comput. Phys. Commun., 2013, 184, 833–840.
73. S. Bhowmick, D. Bossion, Y. Scribano and Y. V. Suleimanov, Phys. Chem. Chem. Phys., 2018, 20, 26752–26763.
74. Q. Liu, L. Zhang, Y. Li and B. Jiang, J. Phys. Chem. Lett., 2019, 10, 7475–7481.
75. H. Jiang, M. Kammler, F. Ding, Y. Dorenkamp, F. R. Manby, A. M. Wodtke, T. F. Miller, A. Kandratsenka and O. Bünermann, Science, 2019, 364, 379–382.
76. M. Bartolomei, E. Carmona-Novillo, M. I. Hernández, J. Campos-Martínez and F. Pirani, J. Phys. Chem. C, 2013, 117, 10512–10522.
77. M. Pitoňák and A. Heßelmann, J. Chem. Theory Comput., 2010, 6, 168–178.
78. W. H. Miller, S. D. Schwartz and J. W. Tromp, J. Chem. Phys., 1983, 79, 4889–4898.
79. C. H. Bennett, Molecular Dynamics and Transition State Theory: The Simulation of Infrequent Events, 1977, vol. 46, pp. 63–97.
80. D. Chandler, J. Chem. Phys., 1978, 68, 2959–2970.
81. R. Collepardo-Guevara, I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2008, 128, 144502.
82. J. Kästner and W. Thiel, J. Chem. Phys., 2005, 123, 144104.
83. J. Kästner and W. Thiel, J. Chem. Phys., 2006, 124, 234106.
84. J. Kästner, J. Chem. Phys., 2009, 131, 034109.
85. H. C. Andersen, J. Chem. Phys., 1980, 72, 2384–2393.
86. H. C. Andersen, J. Comput. Phys., 1983, 52, 24–34.
87. L. M. Robeson, J. Membr. Sci., 2008, 320, 390–400.
88. H. Jiang, X. Tao, M. Kammler, F. Ding, A. M. Wodtke, A. Kandratsenka, T. F. Miller and O. Bünermann, J. Phys. Chem. Lett., 2021, 12, 1991–1996.
89. P. Kapitza, Nature, 1938, 141, 74.
90. D. D. Osheroff, R. C. Richardson and D. M. Lee, Phys. Rev. Lett., 1972, 28, 885–888.
91. B. Hirshberg, V. Rizzi and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 21445–21449.
92. B. Hirshberg, M. Invernizzi and M. Parrinello, J. Chem. Phys., 2020, 152, 171102.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1cp02121d

---