Enhanced Selective Hydrogen Permeation through Graphdiyne Membrane: A Theoretical Study

Graphdiyne (GDY), with uniform pores and atomic thickness, is attracting widespread attention for application in H2 separation in recent years. However, the challenge lies in the rational design of GDYs for fast and selective H2 permeation. By MD and DFT calculations, several flexible GDYs were constructed to investigate the permeation properties of four pure gas (H2, N2, CO2, and CH4) and three equimolar binary mixtures (H2/N2, H2/CO2, and H2/CH4) in this study. When the pore size is smaller than 2.1 Å, the GDYs acted as an exceptional filter for H2 with an approximately infinite H2 selectivity. Beyond the size-sieving effect, in the separation process of binary mixtures, the blocking effect arising from the strong gas–membrane interaction was proven to greatly impede H2 permeation. After understanding the mechanism, the H2 permeance of the mixtures of H2/CO2 was further increased to 2.84 × 105 GPU by reducing the blocking effect with the addition of a tiny amount of surface charges, without sacrificing the selectivity. This theoretical study provides an additional atomic understanding of H2 permeation crossing GDYs, indicating that the GDY membrane could be a potential candidate for H2 purification.


Introduction
As an attractive alternative fuel source, hydrogen (H 2 ) could eliminate the use of polluting fossil fuels in industry and transport in the future [1]. So, the H 2 energy is critical to reduce global greenhouse gases and promote sustainable development because of its natural abundance and high efficiency of combustion [2,3]. Today, in many industrial and drilling streams, H 2 is mainly produced from natural gas, hydro-electrolysis, and the combustion of hydrides [4,5]. However, these processes release several million tonnes of by-products (such as carbon dioxide, nitrogen, and methane) per year [6]. How to separate the target product of H 2 from the undesirable species is crucial to improve production efficiency and reduce the cost. Traditional separation techniques, such as pressure swing adsorption and cryogenic separation, would consume a considerable amount of energy to collect H 2 [7]. With the improved performances and lower operating conditions, the advanced membrane-based separation is seen as an alternative to significantly improve energy efficiency [8,9]. At room temperature and low transmembrane pressure, it is more instructive for experiments to develop membranes with both good H 2 permeance and selectivity. [10].
Recently, two-dimensional (2D) carbon-based membranes have sparked global attention due to their atomic thickness [11]. Among various membranes, graphene-based materials are, assuredly, one of the powerful candidates for membrane separation according to our previous experimental [12,13] and simulation studies [14,15]. According to the separation mechanism of the size sieving effect, controllable pores in the 2D-material membrane are imperative for gas separation. However, it is Based on the optimized membranes, MD calculations were performed to simulate the gas separation. Figure 2 shows the simulation system, where two chambers were isolated by the GDY membrane. The left chamber was a gas reservoir, comprising 2000 molecules of pure gases (or binary mixtures with the mixing ratio of 1000:1000, in volume ratio). The right one treated as a vacuum is the permeate side, which is the most common setting in MD simulations for collecting the permeated molecules. [20,23,30] At both ends of each chamber, rigid graphene was placed to prevent molecules from roaming between the periodic boxes. Considering that carbon-based membranes usually have good flexibility in experiments [32], all pores in GDYs were treated as flexible in our present work so that the atoms around pores could perform small displacements. Thus, the aperture would be enlarged to allow the passage of H2 even though the molecular size is larger than the pore size. On the contrary, the atoms on the edge of membranes were imposed with a position restriction to ensure the uniform lattice of GDYs and decrease the impact on gas permeation [33]. Furthermore, no collapse appeared in our system, suggesting good stability of the flexible structures. The framework of GDYs was described by the all-atom optimized potential. The four gases were represented by Lennard-Jones (LJ) and electrostatic potentials originating from our previous work [14]. Between dissimilar atoms, the interactions were assembled by the Lorentz-Berthelot combination rule [34]. Based on the optimized membranes, MD calculations were performed to simulate the gas separation. Figure 2 shows the simulation system, where two chambers were isolated by the GDY membrane. This left chamber was a gas reservoir, comprising 2000 molecules of pure gases (or binary mixtures with the mixing ratio of 1000:1000, in volume ratio). This right one treated as a vacuum is the permeate side, which is the most common setting in MD simulations for collecting the permeated molecules. [20,23,30] At both ends of each chamber, rigid graphene was placed to prevent molecules from roaming between the periodic boxes. Considering that carbon-based membranes usually have good flexibility in experiments [32], all pores in GDYs were treated as flexible in our present work so that the atoms around pores could perform small displacements. Thus, the aperture would be enlarged to allow the passage of H 2 even though the molecular size is larger than the pore size. On the contrary, the atoms on the edge of membranes were imposed with a position restriction to ensure the uniform lattice of GDYs and decrease the impact on gas permeation [33]. Furthermore, no collapse appeared in our system, suggesting good stability of the flexible structures. This framework of GDYs was described by the all-atom optimized potential. This four gases were represented by Lennard-Jones (LJ) and electrostatic potentials originating from our previous work [14]. Between dissimilar atoms, the interactions were assembled by the Lorentz-Berthelot combination rule [34]. In this study, MD simulations were all carried out using the GROMACS package (version 4.5.5) [35]. A static minimization with the steepest descent method was first carried out to remove the unreasonable contacts among each atom. Then, the separation system was pre-equilibrated in the isobaric-isothermal ensemble (constant temperature and constant pressure ensemble, named as NPT) for 2 ns only in the z-direction. Following this, another 20 ns MD simulations were performed in the canonical ensemble (constant temperature and constant volume ensemble, named as NVT) for data collection and further analysis, and the trajectories were stored every 1 ps. During the simulations, the classical equations of motion were integrated with a time step of 1 fs by using the leapfrog algorithm. The temperature was maintained at 300 K by coupling with the velocity rescaling [36] thermostats. The long-range electrostatic interactions were computed by using the method of particle mesh Ewald [37]. While the short-range van der Waals interactions were truncated at a cut-off distance of 1.2 nm, the periodic boundary conditions were implemented in all three directions. With the partial pressure gradient as the driving force [38], the gases would permeate through the GDYs. For a better understanding of the separation process, an animation is provided in the supplementary materials (Video S1).

Gas Permeation Behavior
The relationship between gas flux ( , mol·m −2 ·s −1 ) and permeance ( , mol·m −2 ·s −1 ·Pa −1 ) was described as Equation (1): where N refers to the number of permeated molecules, is the Avogadro constant, and ∆ is the transmembrane gas partial pressure estimated by Equation (2): in which refers to the number of gases adsorbed on the membrane surface, and are the volumes of the left and right chambers, respectively, and k represents the Boltzmann constant. Therefore, the time evolution of the number of permeated molecules is integrated as Equation (3), where R is the gas constant, In this study, MD simulations were all carried out using the GROMACS package (version 4.5.5) [35]. A static minimization with the steepest descent method was first carried out to remove the unreasonable contacts among each atom. Then, the separation system was pre-equilibrated in the isobaric-isothermal ensemble (constant temperature and constant pressure ensemble, named as NPT) for 2 ns only in the z-direction. Following this, another 20 ns MD simulations were performed in the canonical ensemble (constant temperature and constant volume ensemble, named as NVT) for data collection and further analysis, and the trajectories were stored every 1 ps. During the simulations, the classical equations of motion were integrated with a time step of 1 fs by using the leapfrog algorithm. This temperature was maintained at 300 K by coupling with the velocity rescaling [36] thermostats. This long-range electrostatic interactions were computed by using the method of particle mesh Ewald [37]. While the short-range van der Waals interactions were truncated at a cut-off distance of 1.2 nm, the periodic boundary conditions were implemented in all three directions. With the partial pressure gradient as the driving force [38], the gases would permeate through the GDYs. For a better understanding of the separation process, an animation is provided in the Supplementary Materials (Video S1).

Gas Permeation Behavior
The relationship between gas flux (J, mol·m −2 ·s −1 ) and permeance (S, mol·m −2 ·s −1 ·Pa −1 ) was described as Equation (1): where N refers to the number of permeated molecules, N A is the Avogadro constant, and ∆P is the transmembrane gas partial pressure estimated by Equation (2): in which N ad refers to the number of gases adsorbed on the membrane surface, V l and V r are the volumes of the left and right chambers, respectively, and k represents the Boltzmann constant. Therefore, the time evolution of the number of permeated molecules is integrated as Equation (3), where R is the gas constant, Membranes 2020, 10, 286

of 11
As presented in Figure 3a, the number of permeated H 2 is increased exponentially with the simulation time, which agrees well with the above mathematical analysis. This permeance of pure H 2 is shown in Figure 3b. Obviously, it is remarkably enhanced with the increase of porosity and pore size. With the smallest pore of 1.5 Å, the porosity of 0.03 and 0.07 have the H 2 permeance of 1.34 × 10 5 and 2.55 × 10 5 GPU (1 GPU = 3.35 × 10 −10 mol·m −2 ·s −1 ·Pa −1 ), respectively. This permeation of other gases (i.e., CO 2 , N 2 and CH 4 ) through the GDYs was also simulated. It was shown that only the biggest pore with 2.5 Å allowed the passage of N 2 , CO 2 , and CH 4 , as presented in Figure 3c. This incompatibility of the kinetic diameter of gases (i.e., H 2 , 2.89 Å; N 2 , 3.64 Å; CO 2 , 3.30 Å; CH 4 , 3.8 Å) [39], and the pore size is ascribed to the flexible structures, implying that the membrane of GDY_2.5Å is not suitable for H 2 purification. Moreover, the comparable pore with 2.1 Å diameter can not only completely block the other three gases but also process the good H 2 permeance of 5.94 × 10 5 GPU, which implies that the selectivities of H 2 over CO 2 , N 2 , and CH 4 can be extremely high in the membrane of GDY_2.1Å. We noted that the molecules of N 2 , CO 2 , and CH 4 can not be detected in the permeation side as long as the pore diameter was smaller than 2.1 Å, indicating the infinitely low permeance of N 2 , CO 2 , and CH 4 .
Membranes 2020, 10, x FOR PEER REVIEW 5 of 11 As presented in Figure 3a, the number of permeated H2 is increased exponentially with the simulation time, which agrees well with the above mathematical analysis. The permeance of pure H2 is shown in Figure 3b. Obviously, it is remarkably enhanced with the increase of porosity and pore size. With the smallest pore of 1.5 Å, the porosity of 0.03 and 0.07 have the H2 permeance of 1.34 × 10 5 and 2.55 × 10 5 GPU (1 GPU = 3.35 × 10 −10 mol·m −2 ·s −1 ·Pa −1 ), respectively. The permeation of other gases (i.e., CO2, N2 and CH4) through the GDYs was also simulated. It was shown that only the biggest pore with 2.5 Å allowed the passage of N2, CO2, and CH4, as presented in Figure 3c. The incompatibility of the kinetic diameter of gases (i.e., H2, 2.89 Å; N2, 3.64 Å; CO2, 3.30 Å; CH4, 3.8 Å) [39], and the pore size is ascribed to the flexible structures, implying that the membrane of GDY_2.5Å is not suitable for H2 purification. Moreover, the comparable pore with 2.1 Å diameter can not only completely block the other three gases but also process the good H2 permeance of 5.94 × 10 5 GPU, which implies that the selectivities of H2 over CO2, N2, and CH4 can be extremely high in the membrane of GDY_2.1Å. We noted that the molecules of N2, CO2, and CH4 can not be detected in the permeation side as long as the pore diameter was smaller than 2.1 Å, indicating the infinitely low permeance of N2, CO2, and CH4. To evaluate the ideal selectivities of H2 over other three gases in GDY_2.1 Å, the DFT calculations were performed to calculate permeation barriers as per our previous study [15]. Figure  4a illustrates the minimum energy pathway (MEP) of four gases crossing the membrane. The inset configurations are the energetically stable states of CO2 permeation at different locations. By searching the saddle point in MEPs, the energy barriers of permeation can be calculated. Evidently, the permeation barrier of H2 crossing the flexible GDY_2.1 Å is drastically reduced to 1.55 kJ/mol (Figure 4a, inset), which results in extraordinary H2 permeance. On the contrary, the permeation barriers of CO2, N2, and CH4 are all extremely high. In Figure 4b, the temperature-dependent H2 selectivity has an inverse correlation with the temperature according to the Arrhenius Equation (4), whereas Pi is the permeation rate, A is the permeation prefactor that can be assumed as 10 11 s −1 [38], and T is the temperature. For H2/CO2 and H2/N2, the ideal selectivities of H2 can be up to 10 17 at 300 K, which is the same order with the modified GDYs [25]. It remains very high (>10 8 ) even at 600 K, further suggesting the extraordinary H2 separation performance through GDY_2.1 Å. To evaluate the ideal selectivities of H 2 over other three gases in GDY_2.1 Å, the DFT calculations were performed to calculate permeation barriers as per our previous study [15]. Figure 4a illustrates the minimum energy pathway (MEP) of four gases crossing the membrane. This inset configurations are the energetically stable states of CO 2 permeation at different locations. By searching the saddle point in MEPs, the energy barriers of permeation can be calculated. Evidently, the permeation barrier of H 2 crossing the flexible GDY_2.1 Å is drastically reduced to 1.55 kJ/mol (Figure 4a, inset), which results in extraordinary H 2 permeance. On the contrary, the permeation barriers of CO 2 , N 2 , and CH 4 are all extremely high. In Figure 4b, the temperature-dependent H 2 selectivity has an inverse correlation with the temperature according to the Arrhenius Equation (4), whereas P i is the permeation rate, A is the permeation prefactor that can be assumed as 10 11 s −1 [38], and T is the temperature. For H 2 /CO 2 and H 2 /N 2 , the ideal selectivities of H 2 can be up to 10 17 at 300 K, which is the same order with the modified GDYs [25]. It remains very high (>10 8 ) even at 600 K, further suggesting the extraordinary H 2 separation performance through GDY_2.1 Å.

Transport Mechanism: Blocking Effect
For binary gas mixtures, the MD calculations were also performed to simulate the separation process. To visualize H 2 purification properties through the membrane of GDY_2.1Å, the equilibrium configurations of final frames are presented in Figure 5a-c. As seen, the GDY_2.1Å membrane acts as an effective filter for H 2 separation while completely blocking the passage of CO 2 , N 2 , and CH 4 , so that the H 2 is largely gathering in the vacuum chamber. This corresponding H 2 permeance of the three binary mixtures is presented in Figure 5d. Interestingly, the mixture of H 2 /N 2 exhibits the highest H 2 permeance of 4.71 × 10 5 GPU, which is a little higher than that of H 2 /CH 4 and almost twice as many as the mixture of H 2 /CO 2 . This primary reason is that on the membrane surface, the strongly interacting gas (CO 2 ) preferentially adsorbs, blocking the transport pores of H 2 molecules, thus resulting in a relatively low H 2 permeance.

Transport Mechanism: Blocking Effect
For binary gas mixtures, the MD calculations were also performed to simulate the separation process. To visualize H2 purification properties through the membrane of GDY_2.1Å, the equilibrium configurations of final frames are presented in Figure 5a-c. As seen, the GDY_2.1Å membrane acts as an effective filter for H2 separation while completely blocking the passage of CO2, N2, and CH4, so that the H2 is largely gathering in the vacuum chamber. The corresponding H2 permeance of the three binary mixtures is presented in Figure 5d. Interestingly, the mixture of H2/N2 exhibits the highest H2 permeance of 4.71 × 10 5 GPU, which is a little higher than that of H2/CH4 and almost twice as many as the mixture of H2/CO2. The primary reason is that on the membrane surface, the strongly interacting gas (CO2) preferentially adsorbs, blocking the transport pores of H2 molecules, thus resulting in a relatively low H2 permeance.  Further analyses are carried out to understand this blocking effect. This radial distribution function (RDF) was calculated to present the affinity of GDY_2.1Å to four gases by using Equation (5): where N ij (r,r + ∆r) is the number of species j around i within a shell from r to r + ∆r, r is the distance between two species, and N i and N j refer to the numbers of atom types i and j, respectively. As presented in Figure 6a, there is an increasing trend of the gas-membrane interaction as H 2 < N 2 < CH 4 < CO 2 . This weakest interaction together with the smallest molecular size endows H 2 with exceptional permeance. Meanwhile, the strong interaction promotes the gases, particularly CO 2 to preferentially adsorb on the membrane surface. To offer more intuitionistic information, the density contours of the distribution of N 2 , CH 4 , and CO 2 were plotted on the GDY_2.1Å surface. This general distribution behavior of the three gases is similar as shown in Figure 6b-d, where the pores that are largely clogged by these gases are larger than the pore diameters. Nevertheless, the intensity of gas accumulation is quite different. As seen in Figure 6d, the pores are the most clogged, with the number of adsorbed CO 2 molecules exceeding 6.5 N w /uc in every pore. Similar to the affinity analysis in Figure 6a, the extent of the blocking effect follows the increasing trend of N 2 < CH 4 < CO 2 , decreasing the placeholders of H 2 on the membrane surface in Figure 7a.
Membranes 2020, 10, x FOR PEER REVIEW 7 of 11 Further analyses are carried out to understand this blocking effect. The radial distribution function (RDF) was calculated to present the affinity of GDY_2.1Å to four gases by using Equation (5): where Nij (r,r + Δr) is the number of species j around i within a shell from r to r + Δr, r is the distance between two species, and Ni and Nj refer to the numbers of atom types i and j, respectively. As presented in Figure 6a, there is an increasing trend of the gas-membrane interaction as H2 N2 CH4 CO2. The weakest interaction together with the smallest molecular size endows H2 with exceptional permeance. Meanwhile, the strong interaction promotes the gases, particularly CO2 to preferentially adsorb on the membrane surface. To offer more intuitionistic information, the density contours of the distribution of N2, CH4, and CO2 were plotted on the GDY_2.1Å surface. The general distribution behavior of the three gases is similar as shown in Figure 6b-d, where the pores that are largely clogged by these gases are larger than the pore diameters. Nevertheless, the intensity of gas accumulation is quite different. As seen in Figure 6d, the pores are the most clogged, with the number of adsorbed CO2 molecules exceeding 6.5 Nw/uc in every pore. Similar to the affinity analysis in Figure 6a, the extent of the blocking effect follows the increasing trend of N2 CH4 CO2, decreasing the placeholders of H2 on the membrane surface in Figure 7a.  placeholders of H2, ascribing to the lower blocking effect of N2, promote the faster movement of H2, as presented in Figure 7b. H2 in three binary mixtures of H2/N2, H2/CH4, and H2/CO2 exhibits the diffusion coefficients of 0.57 × 10 −2 , 0.44 × 10 −2 , and 0.33 × 10 −2 cm 2 /s, respectively. That is, the H2 permeation is indeed impeded by the preferentially adsorbed CO2 molecules. In other words, the GDY_2.1Å can separate H2 faster and more selectively from such a kind of binary mixture that the other species has a weak interaction toward GDY surfaces, such as H2/N2.

Surface Charge Effect
Two important mechanisms dominated the gas separation. The first one is the size sieving effect, which favors the permeation of small-sized molecules of H2. While beyond the size sieving effect, in the separation process of binary mixtures, the blocking effect greatly impedes H2 permeation. The above understanding of the separation mechanism is beneficial to further guide the improvement of H2 permeance particularly for the mixture of H2/CO2 without sacrificing its selectivity. An effective way proposed here is to reduce the most severe blocking effect of CO2 by surface charge modification. As depicted in the inset figure in Figure 8a, the positive and negative charges are uniformly imposed on the network of GDY_2.1Å from ±0.00001 to ±0.35 e/atom, whereas the net charge of the whole membrane is zero. The CO2 molecules still cannot cross the membrane regardless of surface charges. Compared to the pristine GDY membrane, the increased H2 permeance of the binary mixture of H2/CO2 was observed when surface charges are lower than ±0.1 e/atom in Figure 8b. It can be up to 2.84 × 10 5 GPU by imposing a tiny amount of surface charge of ±0.035-0.050 e/atom. This exceptionally high permeance is several orders of magnitude greater than the existed experiments [40][41][42], which is ascribed to the ultimate pore size and atomic thickness of GDYs. The following decreasing trend in Figure 8b was ascribed to the surface overcharge, where According to the first peak in Figure 6a, the mean square displacement (MSD) of H 2 molecules was analyzed within 0.5 nm of the membrane surface (Figure 7b, inset) from where r i (t) 2 − r i (t 0 ) 2 refers to the distance traveled by the atom (i) over the time interval of t − t 0 . This diffusion coefficient (D) is determined by the linear slope of the MSD curve with Equation (7). This diffusion behavior in this confined region cannot sustain a normal movement with a large time interval, and it is sufficient to calculate the diffusion coefficient within 20 ps. This greater placeholders of H 2 , ascribing to the lower blocking effect of N 2 , promote the faster movement of H 2 , as presented in Figure 7b. H 2 in three binary mixtures of H 2 /N 2 , H 2 /CH 4 , and H 2 /CO 2 exhibits the diffusion coefficients of 0.57 × 10 −2 , 0.44 × 10 −2 , and 0.33 × 10 −2 cm 2 /s, respectively. That is, the H 2 permeation is indeed impeded by the preferentially adsorbed CO 2 molecules. In other words, the GDY_2.1Å can separate H 2 faster and more selectively from such a kind of binary mixture that the other species has a weak interaction toward GDY surfaces, such as H 2 /N 2 .

Surface Charge Effect
Two important mechanisms dominated the gas separation. This first one is the size sieving effect, which favors the permeation of small-sized molecules of H 2 . While beyond the size sieving effect, in the separation process of binary mixtures, the blocking effect greatly impedes H 2 permeation. This above understanding of the separation mechanism is beneficial to further guide the improvement of H 2 permeance particularly for the mixture of H 2 /CO 2 without sacrificing its selectivity. An effective way proposed here is to reduce the most severe blocking effect of CO 2 by surface charge modification. As depicted in the inset figure in Figure 8a, the positive and negative charges are uniformly imposed on the network of GDY_2.1Å from ±0.00001 to ±0.35 e/atom, whereas the net charge of the whole membrane is zero. This CO 2 molecules still cannot cross the membrane regardless of surface charges. Compared to the pristine GDY membrane, the increased H 2 permeance of the binary mixture of H 2 /CO 2 was observed when surface charges are lower than ±0.1 e/atom in Figure 8b. It can be up to 2.84 × 10 5 GPU by imposing a tiny amount of surface charge of ±0.035-0.050 e/atom. This exceptionally high permeance is several orders of magnitude greater than the existed experiments [40][41][42], which is ascribed to the ultimate pore size and atomic thickness of GDYs. This following decreasing trend in Figure 8b was ascribed to the surface overcharge, where the generated strong electrostatic repels not only CO 2 but also H 2 from approaching. Thus, achieving the perfect balance between these two mechanisms, a tiny amount of surface charges is qualified to maximize the H 2 permeance by reducing the blocking effect of strong interlaced molecules of CO 2 , meanwhile ensuring the placeholders of small-sized molecules of H 2 as well.
Membranes 2020, 10, x FOR PEER REVIEW 9 of 11 the generated strong electrostatic repels not only CO2 but also H2 from approaching. Thus, achieving the perfect balance between these two mechanisms, a tiny amount of surface charges is qualified to maximize the H2 permeance by reducing the blocking effect of strong interlaced molecules of CO2, meanwhile ensuring the placeholders of small-sized molecules of H2 as well.

Conclusions
In summary, a multiscale study combining MD and DFT calculations was performed to investigate the gas permeation through carefully designed flexible GDY membranes with different pore structures and surface charges. Four single gases and three equimolar binary mixtures (H2/other gas) were simulated to study the selective gas permeation through GDYs. Approximately infinite selectivities of H2 over N2, CO2, and CH4 in GDY_2.1Å membranes were demonstrated by DFT calculations. The underlying mechanism indicated that the blocking effect impeded H2 permeation and the GDY_2.1Å is prone to separate H2 from such binary mixtures in which the other species has a weak gas-membrane interaction. Moreover, by imposing a tiny amount of surface charges, the H2 permeance of the binary H2/CO2 mixture was further enhanced up to 2.84 × 10 5 GPU without sacrificing the selectivity. These excellent transport properties make the GDYs a promising candidate for efficient H2 purification. Although the present work is a theoretical study, it is believable that the GDYs can be elaborately designed for realistic separations with the improvement of a well-controlled synthesis strategy.