Data-driven molecular dynamics simulation of water isotope separation using a catalytically active ultrathin membrane

Abstract

Water isotope separation, specifically separating heavy from light water, is a technologically important problem due to the usage of heavy water in applications such as nuclear magnetic resonance, nuclear power, and spectroscopy. Separation of heavy water from light water is difficult due to very similar physical and chemical properties between the isotopes. We show that a catalytically active ultrathin membrane (e.g., a nanopore in MoS₂) can enable chemical exchange processes and physicochemical mechanisms that lead to efficient separation of deuterium from hydrogen. The separation process is inherently multiscale in nature with the shorter times representing chemical exchange processes and the longer timescales representing the transport phenomena. To bridge the timescales, we employ a deep learning methodology which uses short time scale ab initio molecular dynamics data for training and extends the timescales to the classical molecular dynamics regime to demonstrate isotope separation and reveal the underlying complex physicochemical processes.

---

Introduction

Heavy water, an isotope of H₂O, is of interest as many of its properties are quite similar to those of regular water. Heavy water is required for various applications, such as in nuclear power plants,¹ spectroscopy,² and nuclear magnetic resonance.³ Different production methods have been developed,^4,5 including chemical exchange,⁶ distillation,⁷ electrolysis,⁸ crystallization,⁹ and biological techniques;¹⁰ however, only a few of these methods are used for industrial purposes owing to their higher energy efficiency, capacity, and selectivity. The most popular method is the Girdler sulfide chemical exchange process. In the mixture of water and hydrogen sulfide, deuterium in the hydrogen sulfide is transferred to water and vice versa (H₂O + HDS ⇌ HDO + H₂S). The equilibrium constants depend on temperature (e.g., 2.33 and 1.82 at 303 K and 403 K, respectively), which is utilized for concentrating deuterium in water.⁴ Distillation exploits the boiling temperature of heavy water being higher than regular water. Water vapor is enriched with regular hydrogen, and deuterium is concentrated in liquid water. Electrolysis decomposes water into hydrogen (deuterium) and oxygen gases, where regular hydrogen is more likely to be decomposed than deuterium. However, these methods are expensive, and advances in separation technologies are needed to tackle isotope separations.

Owing to the recent advances in nanotechnology, promising alternative approaches have been developed. Mohammadi et al. experimentally studied the physicochemical properties of graphene oxide membranes that affect heavy-water filtration performance and discovered that a low oxidation level of the membrane gives rise to a high heavy-water rejection rate.¹¹ Furthermore, Ono et al. suggested activated carbon fiber (ACF) as a heavy water separation medium by generating adsorption and desorption differences controlled by external pressure.¹²

In this work, we investigate separation of water isotopes using ultrathin membranes with catalytic sites. While there are a number of ultrathin membranes, including two-dimensional materials and MXenes, that contain catalytic sites, as a model ultrathin membrane we study MoS₂. Owing to its selective transport characteristics, MoS₂ nanopore is widely studied for molecule separation research, such as water desalination,¹³ deoxyribonucleic acid (DNA) sequencing,¹⁴ and molecule separation.¹⁵ A single layer MoS₂ has good mechanical properties with an effective Young's modulus of 270 ± 100 GPa,¹⁶ exhibiting high structural stability. The presence of Mo atoms in the pore region can facilitate chemical exchange over atomic length scales, and we exploit this phenomenon for separation for water isotopes.

The physicochemical processes involved in water isotope separation include adsorption/desorption, molecule dissociation and formation of new molecules, chemical reactions, induced driving forces, and transport. The timescales associated with these processes can range from femtoseconds to nanoseconds and longer. Ab initio molecular dynamics (AIMD) calculations, combining density functional theory with classical molecular dynamics, can reasonably describe the physicochemical processes involved in water isotope separation. However, these calculations are limited to short timescales – typically, several tens of picoseconds. As a result, it is difficult to observe long timescale physicochemical processes such as complete translocation of molecules through the nanopore. In this work, to overcome this limitation, we employ a deep learning methodology which uses AIMD data for training. We demonstrate that the deep learning methodology can resolve the physicochemical processes involved in water isotope separation processes.

Methods

We conducted a quantum-level computational study on heavy water separation through a MoS₂ nanopore. A quantum-accurate method is required to explore the complete physicochemical transport^17,18 and the catalytic effect of the Mo atom on water.^19–21 AIMD simulation is one of the popular and quantum-accurate approaches in molecular simulation. The strength of AIMD simulation is its quantum-level precision that can help study reactive systems;^22–25 however, the method is impractical for large systems owing to significant computational costs. To reduce the computational cost of simulating reactive systems, reactive interatomic potentials based on quantum-accurate methods have been developed.^26,27 Neural network potentials^28–34 are especially popular owing to their universal reproducibility for a target system without ad hoc approximation and fitting functions. Here, we train neural network potentials for AIMD data to generate long-time reactive simulation trajectories, capturing physical and chemical details in the separation process.

AIMD simulation

AIMD simulations are used to solve Newtonian dynamics where the force is obtained from the density functional theory.³⁵ We use the Vienna ab initio simulation (VASP) package, in which the projector-augmented wave method with a plane-wave basis set is used.^36,37 The Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional from generalized gradient approximation is employed,³⁸ and the DFT-D3 method³⁹ is used for van der Waals interaction correction. Energy cutoff, Gaussian smearing width, and timestep size are taken to be 400 eV, 0.05 eV, and 0.5 fs, respectively. No external pressure is applied, and the temperature is maintained at 400 K. We perform AIMD simulations, having the same amount of H₂O and D₂O (81 for each) in the system. In this study, the MoS₂ nanopore is constructed using a repeating arrangement of molybdenum and sulfur atoms, referred to as the mixed-type MoS₂ nanopore.¹³ This nanopore is positioned at the system's center and is created by extracting 12 Mo atoms and 24 S atoms from a MoS₂ membrane. It has an approximate diameter of 13 Å, allowing room for molecules to pass through. The masses of Mo, S, O, H, and D atoms are 95.94, 32.066, 16, 1.013, and 2.013 Da, respectively.

Because the computational cost of running AIMD simulations for the target system is considerable, it is useful to generate data for the subsystems containing essential physics and use the data for pretraining. This is followed by fine-tuning where the pre-trained neural network is further trained with the whole system data. We consider three subsystems designed to highlight water–water interaction in bulk, water–MoS₂ wall interaction, and water–pore interaction (ESI,† S.2 for details).

DeePMD training

Using the AIMD data generated above, we train a DeePMD-SE model,³¹ where a neural network is trained to predict the potential energy for a descriptor containing information about atomic configuration in a local neighborhood. The descriptor is an embedded matrix that takes generalized local coordinates as input. The local coordinates of atom i, Rⁱ, comprise N_i row vectors, {x_ji, y_ji, z_ji}, which are then transformed to global coordinates R̃ⁱ the row vector of which can be written as {s(r_ij)x̂_ji, ŷ_ji, ẑ_ji}. Here, s(r_ij) is a differentiable scalar weighting function that is defined aswhere r_cs and r_c signify smooth cutoff and cutoff distances, which are 6.00 Å and 7.00 Å, respectively. Additionally, x̂_ji, ŷ_ji, and ẑ_ji are obtained by multiplying s(r_ij)/r_ij with x_ji, y_ji, and z_ji. The embedding network G(s(r_ij)) is used to map s(r_ij) to M₁ outputs, i.e., where we set M₁ as 64. This neural network maps a single value through several hidden layers, to the output layer on embedding feature dimension. Its difference to fitting layer is that embedding networks are atom-type specific and independent from each other and located at the head of the whole neural network. The fitting network, the tail part, takes the output of the embedding networks, and all types of atoms share the same neural network. Having kth number of s(r_ij) inputs, one can obtain . Finally, the feature matrix can be obtained by computing , which is fed into the fitting neural network for potential energy prediction. The fitting neural network or fitting net can be expressed aswhere E, N, and NN_s represent the total potential energy, the total number of atoms, and the feed-forward fitting net model, respectively. The fitting net comprises four fully connected layers containing 128 nodes followed by the tanh activation function and output layers . The loss function is defined aswhere symbolizes minibatch and stands for the minibatch size, which is two in our training procedure. E^(GT)_l, F^(GT)_l, and Ξ^(GT)_l refer to the potential energy, force, and virial of the system, which can be obtained using AIMD simulation. E^(pred)_l, F^(pred)_l, and Ξ^(pred)_l indicate the potential energy, force, and virial of the system, computed using the DeePMD model. p_f, and p_ξ are tunable prefactors that prioritize or underestimate each component in the loss function. Moreover, and p_f have 0.02 and 1000 for initial values, which gradually reach unity. We assign zero for p_ξ during the entire training process, meaning our neural network is mainly trained to prioritize potential energy and force reproducibility. Using the ADAM optimizer with an initial learning rate of 1 × 10⁻³ decaying until 3.51 × 10⁻⁸, neural network training is iterated for 1 000 000 steps. The neural network is initially trained with subsystem data (see ESI,† S.1. for details) and fine-tuned using the entire system data. The fine-tuned neural network has an energy error of 4.05 × 10⁻³ eV and a force error of 1.84 × 10⁻² eV Å⁻¹. For all of the above training procedures, we use DeePMD-kit³⁰ written in Tensorflow.⁴⁰ Further validation considering physical data reproducibility are discussed in the ESI,† S.3. Section.

DeePMD simulation

LAMMPS,^41,42 included in the DeePMD-kit,³⁰ is used to perform time integration and to generate trajectory in a DeePMD simulation. Force is computed by obtaining the negative potential energy gradient with respect to atomic coordinates. We assign the same atomic masses used in AIMD analysis for the DeePMD simulation. The Nosé–Hoover thermostat⁴³ with a time constant of 100 fs is applied, maintaining the system temperature at 400 K.

Results

Heavy-water separation is studied using a neural network potential-based molecular dynamics simulation for a single layer MoS₂ nanopore equilibrated with 81 H₂O and 81 D₂O molecules. We observe chemical reactions in confinement, reactive translocation of molecules, and water isotope separation phenomena in our system (Fig. 1A). A molecule in and near the pore undergoes various chemical reactions (Fig. 1B), such as water molecule adsorption to Mo atom (Fig. 1B(i)), hydrogen detachment from Mo-attached water molecule (Fig. 1B(ii)), proton donation to a water molecule to form a hydronium molecule (Fig. 1B(iii)), proton detachment to from a water molecule (Fig. 1B(iv)), and proton donation to a hydroxyl group (Fig. 1B(v)). To clarify the molecular forms of OH and OD, we use OH/OD attachment (to Mo atom), OH/OD bond (as a part of a water molecule), and OH⁻/OD⁻ ions (solvated in water). We note that these reactions can also occur for deuterium and deuterium-containing molecules. Continuous proton donation and detachment can change the isotopic composition of a water isotope. To study heavy water separation performance, we classify pore physics into inflow, backflow, and translocation (Fig. 1C), which denote the entry of a water isotope molecule into the pore, return of a molecule to the originating reservoir, and escape of the water molecule to the opposite reservoir.

Fig. 1 (A) The physical and chemical phenomena involved in water isotope separation using a MoS₂ nanopore. The deep learning approach is used to efficiently link the physicochemical phenomena ranging from chemical reactions to long-time dynamics. (B) Representative chemical reactions observed in the pore. While the schematics are shown using only the regular hydrogen, these reactions can occur for deuterium and a deuterium-containing molecule. (i) Adsorption of a water isotope to the Mo atom in the pore. (ii) Detachment of a proton (or a deuterium) from the water molecule attached to the Mo atom. An OH⁻ group is now attached to the Mo atom. (iii) The detached H⁺ or D⁺ ion (from step (ii)) combines with a water isotope molecule to form a hydronium ion isotope. (iv) A hydronium ion isotope has lower stability and one of the hydrogen isotopes detaches and is available as a free charge. The freely available H⁺ or D⁺ can combine with a nearby water isotope as in step (iii). (v) Freely available H⁺ or D⁺ ions can combine with the hydroxyl group on the Mo atom to form a water isotope molecule followed by step (ii). (C) Schematic illustrating a molecule undergoing reactive translocation. In this example, an H₂O molecule enters the pore, undergoes H⁺/D⁺ attachment/detachment processes and leaves the pore as a D₂O molecule. The schematic also shows backflow (from the pore to the reservoir) of molecules, which can be different from the entering molecules during the chemical exchange within the pore.

Fig. 2A shows the time history of each translocation event. We track the translocating molecules for 2 ns of simulation trajectories and observe that about 10 of the 23 translocated water isotopes eventually escaped the pore exhibiting a different isotopic state. In Fig. 2B, the number of translocations for each water isotope group is counted, where the translocations are highest for HDO, D₂O, and H₂O. Similarly, the number of hydrogen and deuterium atom translocations is visualized in Fig. 2C, where deuterium translocation is higher than hydrogen translocation by 16. This provides another indication of the separation efficiency, as the fractions of H₂O, HDO, and D₂O after the separation are determined from the number of translocated H and D atoms. The translocation of any hydrogen (or deuterium) atom is counted regardless of the molecular form. For example, the translocation of the H atom includes hydrogen translocation in H₂O, HDO, H⁺, H₃O⁺, H₂DO⁺, and HD₂O⁺. Separation performance is quantified with separation ratios of D₂O over H₂O and deuterium over hydrogen, which are 4.5 and 1.73, respectively. The deuterium over hydrogen separation ratio is smaller than that of D₂O over H₂O as the translocation of HDO, the molecule having the highest degeneracy and population, contributes to both hydrogen and deuterium translocations. Due to the catalytic site in the pore and generation of free H⁺/D⁺ ions, chemical reactions are active in the pore (Fig. 2D); thus, the translocating molecules undergo chemical reactions, changing their isotopic composition giving rise to reactive translocation. For example, in Fig. 2E, H₂O enters the pore, undergoes chemical reactions, and translocates to the opposite reservoir. Similarly, in Fig. 2F, a D⁺ ion detached from D₃O⁺ enters the pore and translocates with a different oxygen atom.

Fig. 2 (A) The history of translocation events. A translocating water molecule is denoted by two circles connected by a dashed line, where the color of each circle indicates the isotopic state at the pore entrance (beginning time) and exit (end time with molecule leaving the pore), respectively. (B) The number of translocations of each molecule visualized over time. (C) The number of translocations of the H and D atoms visualized over time. (D) The probability of reaction with respect to the z-direction. The chemical reactions are enhanced as a molecule enters the pore. Chemical exchanges in the pore lead to a translocating molecule change its isotopic composition. (E) An example of a reactive translocation trajectory with chemical exchange – H₂O enters the pore and combines with a D⁺ ion, forming H₂DO⁺. Since the product is unstable, one of the hydrogen isotopes detaches, which is regular hydrogen in this example, forming HDO. The molecule then escapes the pore as an HDO molecule. Notice that the oxygen atom is the same in this translocation. (F) An example of a hydrogen isotope translocation. Deuterium discharged from D₃O⁺ enters the pore and combines with H₂O, forming H₂DO⁺. The product is unstable and splits into H⁺ and HDO. The HDO escapes the pore. This is an example of a hydrogen isotope translocation with a different oxygen atom.

We further study the physicochemical mechanisms that contribute to the dominance of D₂O and deuterium atom during the translocation, i.e., Mo-attached molecules and the composition of net inflow that is defined as the number of each molecule in the inflow subtracted by that of backflow (see Fig. 1C for details). Considering three molecule types (water isotope group, hydronium isotope group, and H⁺/D⁺ ion group), we observe that a higher number of deuterium-dominant molecules are translocated than their counterparts in the water isotope group and the hydronium isotope group, while more H⁺ ions are translocated than D⁺ ions (Fig. 3A). Understanding the composition of net inflow is crucial, as it allows us to distinguish the contributions of two different translocation mechanisms: the non-reactive diffusion mechanism and the Grotthuss mechanism. The diffusion-based mechanism involves the simple movement of molecules across boundaries without any accompanying chemical reactions (as illustrated in Fig. 3B left). In this process, molecules like water isotopes and hydronium isotopes cross the boundary purely through translational motion. This mechanism is straightforward and does not involve any chemical changes within the molecules during translocation. On the other hand, the Grotthuss mechanism involves a chemical reaction where hydrogen or deuterium ions (H⁺/D⁺) are detached from one molecule and then transferred to a nearby molecule (as shown in Fig. 3B center and right). In this mechanism, a hydrogen or deuterium ion is first detached from a hydronium isotope and then attached to a nearby water isotope, resulting in the formation of a new hydronium isotope with a different oxygen atom (Fig. 3B center). Additionally, a free H⁺/D⁺ ion can directly combine with a hydroxyl group to form a water isotope, or it can be transferred from a nearby hydronium isotope to a hydroxyl group (Fig. 3B right). Our observations indicate that although there are more hydrogen ions (H⁺) than deuterium ions (D⁺) in the net flow, the overall inflow is dominated by deuterium-containing molecules, particularly deuterium-dominant water and hydronium isotopes. This dominance leads to a higher rate of deuterium atom translocation compared to hydrogen atoms. We also study the isotopic composition (Fig. 3C) and physical behavior of a hydroxyl group (Fig. 3D), since they are closely related to chemical reactions in the pore, affecting reactive translocation and isotope separation performance. OH attachment is preferred over OD attachment, which is quantified with the OH to OD attachment time ratio, where a value of 1.25 is obtained. To quantify the reliability of this value, with deuterium occupancy set at 2 and hydrogen at 1, the standard error is 0.00256, resulting in a 95% confidence interval width of 0.00502. Another interesting observation is that the oxygen atoms attached to Mo are almost stationary and do not escape the pore (Fig. 3D). This indicates that when a translocating hydrogen isotope combines with the O atom in a hydroxyl group, it will stay next to the oxygen, i.e., at the center of the pore, until another H⁺/D⁺ ion combines with the hydroxyl group. The OH/OD attachments act like a hindrance to hydrogen isotope translocation, and the H atom is more affected by this physicochemistry. As a result, H atoms are less translocated than that of D atoms. While it is counterintuitive that more deuterium based molecules are transported over hydrogen based molecules, similar phenomena have been observed in other works. Oh et al. demonstrated that D₂ could be separated using a nanopore with pyridine molecules as a size-adjustable aperture.⁴⁴ They reported that H₂ is larger than D₂, originating from the nuclear quantum effect. Niimura et al. conducted D₂ gas separation studies considering various nanoporous materials (e.g., ACF, carbon molecular sieve, zeolite, and carbon nanotubes) and observed more D₂ being transported.⁴⁵

Fig. 3 (A) The net inflow of molecules entering the pore (the number of molecules leaving the pore (backflow) is subtracted from the number entering the pore (inflow). The plots for inflow, backflow, and translocation are shown in ESI;† see Fig. 1(C) for details). D₂O, D-dominant hydronium isotopes, and H⁺ dominate the net inflow molecular group. (B) The mechanisms governing the net inflow. Translation-driven inflow, where a molecule crosses the pore boundary without a chemical reaction – water isotope and hydronium isotope molecules rely on this mechanism. H⁺ and D⁺ ions generally cross the pore boundary through the Grotthuss mechanism, where their acceptors are a water isotope molecule or a hydroxyl group attached to the Mo atom. (C) The visualization of the chemistry of the hydroxyl group attached to the Mo atom over time. OH attachment is preferred over OD attachment, quantified by the attachment time ratio (t_OH/t_OD) of 1.25. (D) The position of Mo-attached oxygen atoms, showing that they are almost stationary. That is, once H⁺ and D⁺ ions are attached to a hydroxyl group, they are trapped until another H⁺ or D⁺ ion combines with the hydroxyl group.

To understand the underlying physics, we further study the energy differences between H₂O and D₂O, OH/OD bond stability, and H⁺/D⁺ ion transfer time. Fig. 4A shows conceptual potential energy surfaces (PES) of OH and OD bonds, as well as the physical quantities related to the stability of OH/OD bond: average potential energies of OH/OD bonds, the fluctuations of OH/OD bond, the fluctuations of force acting on H/D atom. In general, OD bond is more stable than OH bond since the heavier mass of deuterium reduces the fluctuation of OD bond, which reduces the potential energy of OD bond and increases the dissociation energy of OD bond. We believe there is a contribution from quantum mechanical phenomena (e.g. kinetic isotope effect), in addition to classical mechanics (e.g. reaction rate difference). We analyze the potential energy differences between H₂O and D₂O, in bulk and confined environments to obtain physical insights. We perform ab initio molecular dynamics (AIMD) simulations for a bulk system comprising 32 H₂O molecules at 400 K and a separate bulk system comprising 32 D₂O molecules at 400 K (ESI,† S.1. and Fig. S1(A) top). Similarly, AIMD simulations were performed for a MoS₂ nanopore system with four H₂O molecules confined in the pore at 400 K and a separate system with four D₂O molecules confined in the pore at 400 K (ESI,† S.1. and Fig. S1(A) bottom). The total energies of bulk H₂O, bulk D₂O, confined H₂O, and confined D₂O are −14.454 eV, −14.4508 eV, −580.345 eV, and −580.706 eV, respectively. Note that they have the same averaged kinetic energy since the systems are equilibrated at the same temperature. To systematically compare the preference between the water isotopes depending on their environment, we compute the energy differences between isotopic systems and normalize by the number of substituted water isotopes: and obtaining −1.01 × 10⁻⁵ eV and 9.02 × 10⁻² eV, respectively. We present the magnitude of those values in Fig. 4B. In order to evaluate the significance of the values, we compare them with the magnitude of H₂O total energy (4.52 × 10⁻¹ eV, the blue horizontal line in Fig. 4B) and the standard deviation of H₂O kinetic energy (1.81 × 10⁻² eV, the red horizontal line in Fig. 4B). The magnitudes of Δε^bulk and Δε^conf are 0.0223% and 20.0% of the magnitude of H₂O energy, respectively, signifying that H₂O and D₂O are almost equally preferred in bulk, but D₂O is more preferred in nanopore, under the catalytic effect of MoS₂ nanopore edge. Also, the magnitudes of Δε^bulk and Δε^conf have different order from that of the standard deviation of H₂O kinetic energy: 0.0559% and 498%, which supports their difference is not caused by thermal energy fluctuation.

Fig. 4 (A) Conceptual illustration of OH (OD) bond PES with quantities of interest: force, position of zero force (r*), bond length fluctuation, average potential energy, and dissociation energy. The heavier mass of deuterium reduces the range of bond fluctuation, resulting in lower average potential energy and higher force fluctuation. The lower average energy of the OD bond requires more energy for dissociation, i.e., higher dissociation energy. (B) The stability difference between water isotopes is quantified by the absolute value of average total energy difference. In bulk, the total energy of H₂O and D₂O are nearly the same, indicating they have similar stability. In confinement, the total energy of H₂O is significantly higher than that of D₂O implying lower stability of H₂O. To assess the significance of the energy difference, they are compared with the magnitude of average H₂O energy (blue) and kinetic energy fluctuation (red). (C) The fluctuations of force acting on hydrogen and deuterium. Force fluctuations of H₂O and D₂O are higher in confinement than in bulk, indicating that both molecules are relatively unstable when they enter the pore. The increase in the force fluctuations for H₂O is significantly larger compared to that of D₂O, implying H₂O is relatively more unstable than D₂O. (D) The fluctuations in OH and OD bond lengths, which is another measure of stability. The bond fluctuations show a similar tendency to that of the force fluctuations: the fluctuations are smaller in bulk but larger in confinement, and H₂O has higher fluctuation than D₂O when they are in the same environment. (E) The time required for a H⁺ or D⁺ ion transfer to a nearby Mo-attached oxygen atom. H⁺ ion shows lower transfer time due to the higher mobility and the lower dissociation energy.

We also analyze the fluctuation of the force acting on H/D and the fluctuation of OH/OD bond length, which are different quantities related to the stability of OH/OD bond. As the mass of deuterium is heavier than that of hydrogen, it has a smaller vibration amplitude. This reduces deuterium deviating from the position of zero net force (r* in Fig. 4A), stabilizing the molecule and elevating the activation energy required for chemical reaction. The deviation of atom position from r* can be quantified by the fluctuation of force acting on hydrogen (deuterium) atom and the fluctuation of OH/OD bond length. Herein, we use the standard deviation of the force magnitude acting on hydrogen (or deuterium) atoms and the standard deviation of OH (or OD) bond length, respectively. We perform AIMD simulations for the two systems: 12 H₂O molecules that can freely enter and escape the pore confined by an unmovable water shell (ESI,† S.1. and Fig. S1(B) left) and the same system substituting freely movable H₂O with D₂O (ESI,† S.1. and Fig. S1(B) right). As shown in Fig. 4C, the computed force magnitude fluctuation of H₂O outside- and inside-pore are 0.889 eV Å⁻¹ and 1.56 eV Å⁻¹, respectively, and the counterparts of D₂O are 0.835 and 1.03 eV Å⁻¹. In bulk, the hydrogen (deuterium) force fluctuations of water isotopes are almost the same, but it is significantly increased for H₂O in confinement. These results reveal that H₂O in bulk is chemically as stable as D₂O in bulk, while H₂O inside the pore is chemically less stable than D₂O inside the pore, because of the catalytic effect of MoS₂ on water isotopes. Similar results are observed in the bond length fluctuation analysis. The calculated force magnitude fluctuation of H₂O outside- and inside-pore are 0.255 and 0.322 Å, respectively, and the corresponding fluctuations in force for D₂O are 0.217 and 0.289 Å (Fig. 4D). The bond length fluctuations are the highest for D₂O outside the pore, H₂O outside the pore, D₂O inside the pore, and H₂O inside the pore, respectively, consistent with the force magnitude analysis. The AIMD analyses on energy, hydrogen (deuterium) force fluctuation, and OH/OD bond length fluctuation indicate that an OH bond in the pore is relatively less stable compared to an OD bond in the pore. The molecules with rich OD bonds (D₂O, D₃O⁺, and HD₂O⁺) are less reactive, and H⁺ ion is more prevalent in the pore. A similar phenomenon is reported by Pasquini et al., where catalytic effect of Co-based amorphous-oxide catalysts on water is decreased by 82% for D₂O.⁴⁶ The more reactive OH bond increases the probability of H⁺ ions becoming trapped in Mo-active sites, and the molecules with rich OD bonds are favored, resulting in a lower total energy of the system.

To understand proton transfer kinetics within the nanopore, we consider an oxygen atom attached to the Mo-site and consider transfer of proton/deuterium from a H₃O⁺/D₃O⁺ ion inside the pore. AIMD simulations are performed for a H₃O⁺ ion at various distances (d_OO) with an oxygen atom attached to the Mo-site of the nanopore (ESI,† S.1 and Fig. S1(C) left) and for an identical system with hydrogen substituted by deuterium (ESI,† S.1 and Fig. S1(C) right). We calculate the transfer time, defined as the duration between the release of a proton isotope from a hydronium isotope and its subsequent attachment to the oxygen atom bonded to the molybdenum (Mo) atom; it is observed that the H⁺ transfer time is shorter than the D⁺ transfer time (Fig. 4E) because of the higher mobility of H⁺ ions and lower stability of the OH bond. Here, the mobility of each ion plays an important role as once the attached oxygen at Mo-site receives one of the H⁺ or D⁺ ions, it is less likely to accommodate more H⁺/D⁺ ions. This first-come-first-serve mechanism, driven by mobility difference, is reported by Weinrauch et al., where H₂ adsorption can dominate over D₂.⁴⁷ This corresponds to the higher occupancy time of OH attachment than that of OD attachment, shown in Fig. 3C. Higher occupation of OH attachment can be also explained by collision rate theory.⁴⁸ The rate equation in collision theory, r(t) = Zρ exp(−E_a/RT) explains that the reaction frequency is proportional to the collision frequency (Z), steric factor (ρ), and the Boltzmann factor (exp(−E_a/RT)). In this context, the collision frequencies for the hydrogen and deuterium systems are the only difference since hydrogen atoms are lighter than their counterparts, resulting in higher average velocities and increased collision probabilities. Similar phenomenon is observed by Phillips, where the chemical reaction of CH₃⁺ with H₂ is more active than and that of D₂.

Conclusions

We showed that a catalytically active MoS₂ nanopore can effectively separate D₂O and deuterium with a separation ratio of 4.5 and 1.73, respectively. To understand the underlying mechanism, we analyze the dominant molecules comprising the net inflow and their transport mechanisms. D₂O, D-dominant hydronium isotopes, and H⁺ ions constitute the net inflow compared to their isotopic counterparts. Translation-driven inflow mechanism governs the net inflow of water and hydronium isotopes that contain oxygen atoms, and the proton-jumping mechanism governs the net inflow of H⁺ and D⁺ ions. Furthermore, we observed that OH attachment to a Mo-active site is preferred to OD attachment and the attached oxygen atoms are almost stationary, temporarily trapping a hydrogen isotope.

To gain physical insights, we analyzed the fluctuation of the force acting on hydrogen isotopes and the fluctuations of OH (OD) bond length, finding that the OH bond is less stable than the OD bond in the pore, while they have nearly identical stability in bulk. A similar tendency is observed in energy analysis: the energy difference between H₂O and D₂O in bulk is negligible, and that in confinement is considerable. These analyses account for the high number of D₂O and D-dominant hydronium ions in the net inflow. The analysis also explains the dominance of H⁺ ions in the net inflow and the preference for OH attachment to Mo atoms. The higher fluctuation in OH bond length increases the chance of proton detachment, and the lower mass of a hydrogen atom causes higher mobility of H⁺ ions. These effects promote more frequent attachment and detachment of proton to an oxygen atom in a hydroxyl group. The separation mechanism and physical interpretation presented here can be exploited for other separation mechanisms.

Author contributions

Conceptualization: Jinu Jeong, Chenxing Liang, and Narayana Aluru. Data curation: Jinu Jeong and Chenxing Liang. Formal analysis: Jinu Jeong. Funding acquisition: Narayana Aluru. Investigation: Jinu Jeong and Chenxing Liang. Methodology: Jinu Jeong, Chenxing Liang, Narayana Aluru. Project administration: Jinu Jeong and Narayana Aluru. Resources: Narayana Aluru. Software: Jinu Jeong and Chenxing Liang. Supervision: Narayana Aluru. Validation: Jinu Jeong. Visualization: Jinu Jeong. Writing – original draft: Jinu Jeong, Chenxing Liang, and Narayana Aluru. Writing – review and editing: Jinu Jeong, Chenxing Liang, and Narayana Aluru.

Data availability

Data for this article, including [description of data types] are available at Open Science Framework at https://osf.io/b2p7y/?view_only=a2d7e5e5916b4ae79bfead328ef9eb8b.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The work on water was supported by the Center for Enhanced Nanofluidic Transport (CENT), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0019112. All other aspects of the work are supported by the National Science Foundation under grants 2140225 and 2137157. The authors also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing the computing resources on Stampede2, Frontera, and Lonestar6 under Allocation No. TG-CDA100010, DMR20002, and DMR22008, respectively. Additionally, the authors acknowledge the use of the Illinois Campus Cluster Program (ICCP) in conjunction with the National Center for Supercomputing Applications (NCSA).

References

1. M. Ho, E. Obbard, P. A. Burr and G. Yeoh, Energy Proc., 2019, 160, 459–466.
2. G. Zuber, S. J. Prestrelski and K. Benedek, Anal. Biochem., 1992, 207, 150–156.
3. J. J. Ackerman, C. S. Ewy, N. N. Becker and R. A. Shalwitz, Proc. Natl. Acad. Sci. U. S. A., 1987, 84, 4099.
4. H. K. Rae, Separation of Hydrogen Isotopes, American Chemical Society, 1978, vol. 68, pp. 1–26.
5. K. I. Sakodynskii, Sov. J. At. Energy, 1960, 6, 8–12.
6. B. M. Andreev, Sep. Sci. Technol., 2001, 36, 1949–1989.
7. G. M. Keyser, D. B. McConnell, N. Anyas-Weiss and P. Kirkby, Separation of Hydrogen Isotopes, American Chemical Society, 1978, vol. 68, pp. 126–133.
8. G. N. Lewis and R. T. Macdonald, J. Chem. Phys., 1933, 1, 341–344.
9. J. C. Posey and H. A. Smith, J. Am. Chem. Soc., 1957, 79, 555–557.
10. L. Carlbom, R. Skjöldebrand and N. N. Elsen, Nature, 1956, 177, 988.
11. A. Mohammadi, M. R. Daymond and A. Docoslis, ACS Appl. Mater. Interfaces, 2020, 12, 34736–34745.
12. Y. Ono, R. Futamura, Y. Hattori, T. Sakai and K. Kaneko, J. Colloid Interface Sci., 2017, 508, 14–17.
13. M. Heiranian, A. B. Farimani and N. R. Aluru, Nat. Commun., 2015, 6, 8616.
14. A. B. Farimani, K. Min and N. R. Aluru, ACS Nano, 2014, 8, 7914–7922.
15. L. Sun, H. Huang and X. Peng, Chem. Commun., 2013, 49, 10718–10720.
16. S. Bertolazzi, J. Brivio and A. Kis, ACS Nano, 2011, 5, 9703–9709.
17. D. Cohen-Tanugi and J. C. Grossman, Nano Lett., 2012, 12, 3602–3608.
18. K. Sint, B. Wang and P. Král, J. Am. Chem. Soc., 2008, 130, 16448–16449.
19. A. B. Laursen, S. Kegnæs, S. Dahl and I. Chorkendorff, Energy Environ. Sci., 2012, 5, 5577–5591.
20. D. Voiry, H. Yamaguchi, J. Li, R. Silva, D. C. B. Alves, T. Fujita, M. Chen, T. Asefa, V. B. Shenoy, G. Eda and M. Chhowalla, Nat. Mater., 2013, 12, 850–855.
21. D. Merki and X. Hu, Energy Environ. Sci., 2011, 4, 3878–3888.
22. A. Rayabharam and N. R. Aluru, Appl. Phys. Lett., 2022, 120, 211601.
23. Y. Fu, L. Bernasconi and P. Liu, J. Am. Chem. Soc., 2021, 143, 1577–1589.
24. X. Zhao and Y. Liu, J. Am. Chem. Soc., 2020, 142, 5773–5777.
25. C. Li and J. M. J. Swanson, J. Phys. Chem. B, 2020, 124, 5696–5708.
26. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409.
27. W. A. Goddard, in Computational Materials, Chemistry, and Biochemistry: From Bold Initiatives to the Last Mile: In Honor of William A. Goddard's Contributions to Science and Engineering, ed. S. Shankar, R. Muller, T. Dunning and G. H. Chen, Springer International Publishing, Cham, 2021, pp. 1079–1087.
28. T. W. Ko, J. A. Finkler, S. Goedecker and J. Behler, Nat. Commun., 2021, 12, 398.
29. L. Zhang, J. Han, H. Wang, R. Car and W. E, Phys. Rev. Lett., 2018, 120, 143001.
30. H. Wang, L. Zhang, J. Han and W. E, Comput. Phys. Commun., 2018, 228, 178–184.
31. L. Zhang, J. Han, H. Wang, W. Saidi, R. Car and W. E, Advances in Neural Information Processing Systems, Curran Associates, Inc., 2018, vol. 31.
32. C. Bajaj and M. Nguyen, Learning Optimal Control with Stochastic Models of Hamiltonian Dynamics, 2024, https://arxiv.org/abs/2111.08108v1, (accessed 5 June).
33. S. Schoenholz and E. D. Cubuk, Advances in Neural Information Processing Systems, Curran Associates, Inc., 2020, vol. 33, pp. 11428–11441.
34. A. Moradzadeh and N. R. Aluru, J. Phys. Chem. A, 2022, 126, 2031–2041.
35. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
36. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
37. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
38. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
39. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
40. M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, S. Moore, D. G. Murray, B. Steiner, P. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu and X. Zheng, arXiv, 2016, preprint, arXiv:1605.08695 DOI:10.48550/arXiv.1605.08695.
41. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
42. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
43. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
44. H. Oh, S. B. Kalidindi, Y. Um, S. Bureekaew, R. Schmid, R. A. Fischer and M. Hirscher, Angew. Chem., Int. Ed., 2013, 52, 13219–13222.
45. S. Niimura, T. Fujimori, D. Minami, Y. Hattori, L. Abrams, D. Corbin, K. Hata and K. Kaneko, J. Am. Chem. Soc., 2012, 134, 18483–18486.
46. C. Pasquini, I. Zaharieva, D. González-Flores, P. Chernev, M. R. Mohammadi, L. Guidoni, R. D. L. Smith and H. Dau, J. Am. Chem. Soc., 2019, 141, 2938–2948.
47. I. Weinrauch, I. Savchenko, D. Denysenko, S. M. Souliou, H.-H. Kim, M. Le Tacon, L. L. Daemen, Y. Cheng, A. Mavrandonakis, A. J. Ramirez-Cuesta, D. Volkmer, G. Schütz, M. Hirscher and T. Heine, Nat. Commun., 2017, 8, 14496.
48. L. F. Phillips, J. Phys. Chem., 1990, 94, 5265–5271.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp04020a

---