Polarizable Molecular Simulations Reveal How Silicon-containing Functional Groups Govern the Desalination Mechanism in Nanoporous Graphene

We report a molecular dynamics (MD) simulation study of reverse osmosis desalination using nanoporous monolayer graphene passivated by SiH$_2$ and Si(OH)$_2$ functional groups. A highly accurate and detailed polarizable molecular mechanics force field model was developed for simulating graphene nanopores of various sizes and geometries. The simulated water fluxes and ion rejection percentages are explained using detailed atomistic mechanisms derived from analysis of the simulation trajectories. Our main findings are: (1) The Si(OH)$_2$ pores possess superior ion rejection rates due to selective electrostatic repulsion of Cl$^-$ ions, but Na$^+$ ions are attracted to the pore and block water transfer. (2) By contrast, the SiH$_2$ pores operate via a steric mechanism that excludes ions based on the size and flexibility of their hydration layers. (3) In absence of ions, water flux is directly proportional to the solvent accessible area within the pore; however, simulated fluxes are lower than those inferred from recent experimental work. We also provide some hypotheses that could resolve the differences between simulation and experiment.


Introduction
Seawater desalination plays a crucial role in the research for overcoming water scarcity, one of the most serious challenges humans are facing. 1,2 Multiple approaches have been investigated, including multi-stage flash, 3 membrane distillation, 3 mechanical vapor compression, 4 and reverse osmosis (RO). 5 Among them, reverse osmosis technology is by far the most promising method, as it requires the least amount of energy input to purify the same amount of water.
Conventional RO membranes made of composite materials, like polymers, has improved over decades to reach an acceptable energy efficiency; 6 this has led the growth of RO plants all over the world, but improvements in membrane permeability have been limited since the 1990s. 7 Therefore, people seek new forms of materials, such as carbon nanotubes, 8 and nanoporous graphene 9,10 that hold promise for desalination applications.
Single-layer nanoporous graphene is a promising desalination membrane due to its unmatched thinness, good mechanical strength, and chemical stability. Experimental studies on nanoporous graphene involve various treatments to create nanopores in graphene, including electrical pulse, 11 ion bombardment, 12 and O 2 plasma etching. 10 Ion passage through the prepared nanoporous graphene was observed to be highly selective in several studies, [11][12][13] with cation/anion selectivity ratio of over 100; ion selectivity is positively correlated with salt rejection. 11 Recently, one experimental study has shown nanoporous graphene can act as a desalination membrane with surprisingly high water flux and nearly perfect salt rejection. 10 The reported 10 6 g m −2 s −1 water flux (3000 water molecules per pore per ps) under only 0.17 bar driving pressure is several orders of magnitude higher than conventional membranes characterized by fluxes of about 12 g m −2 s −1 under 83 bar. 6 Scanning transmission electron microscopy (STEM) images suggested that the pores produced by O 2 plasma were around 1 nm in diameter. Lee et al. 14 proposed that graphene nanopores could be stabilized by silicon based on STEM imaging and density functional calculations; these results suggest that the edges of nanopores produced by oxygen plasma may also contain silicon. These studies have contributed essential insights into the nanopore structure; however, the molecular mechanism of water transport and ion rejection responsible for the high performance remain unknown, creating an opportunity for theory and simulation to provide detailed insights into the desalination process in nanoporous graphene.
Several published theoretical studies have focused on the possibility of nanoporous graphene as a desalination membrane. In 2012, Cohen-Tanugi and Grossman investigated the use of hydrogen-or hydroxyl-passivated graphene nanopores to separate NaCl ions from water. 9 Their molecular dynamics (MD) study suggested a maximum diameter of 5.5Å for maintaining a high percentage of salt rejection. They also followed up the studies with reduced graphene oxide 15 and multilayer porous graphene. 16 Strong and Eaves investigated the hydrodynamics of the water transferring through porous graphene 17 and suggested that the non-equilibrium dynamics of the water transfer is controlled by two competing mechanismstranslocation mechanism and evaporation-condensation mechanism. Ebrahimi 18 modeled the nanoporous graphene passivated by Si atoms and found that the membrane curvature reduces the water flux. The relation between pore ion selectivity and surface charge was studied by Zhao et al., 19 who predicted that negative charges at pore edges could impede the passage of Cl − . Very recently, Li et al. 20 found that the ion selectivity of multi-layer graphene nanopores may decrease if the surface charge becomes too strong due to counterions entering the pore and causing a charge inversion. These studies all used fixed-charge models that do not account for the environmental dependence of molecular dipole moments, raising the possibility that a more detailed physical model is needed to explain the water transfer phenomenon at the interface of the nanopore.
The inclusion of polarizability greatly improves on the physical descriptions of molecular mechanics (MM) force field models by explicitly modeling induced polarization from external electric fields or intermolecular interactions. 21 In particular, the AMOEBA model includes polarizable point dipoles and also incorporates atomic multipole moments up through the quadrupole, 22-25 providing a more detailed description of molecular electrostatics compared to other widely implemented polarizable force fields such as the Drude oscillator 26,27 and the fluctuating-charge models. 28,29 As a result, AMOEBA models can reach quantitative agreement against gas-phase ab initio interaction energies at various distances and orientations, and produce reliable predictions when the models are used to simulate condensed-phase phenomena. [30][31][32][33][34] Because molecular polarization depends on changes in the electrostatic environment, we expect a polarizable model to yield a more accurate description of water transfer and ion rejection as water molecules and ions move from the bulk and through the nanopore.
In this paper we describe an accurate "AMOEBA-type" polarizable force field for siliconpassivated nanoporous graphene, and its application in MD simulations on model nanopores that reveal their performance in desalination. We find that hydrophobic nanopores, passivated by SiH 2 groups (also denoted as Si-H), has superior water transfer rates compared to hydrophilic Si(OH) 2 nanopores (also denoted as Si-OH) ; on the other hand, hydrophilic Si-OH nanopores selectively block the negatively charged Cl − ions resulting in an improved salt rejection. We also provide new atomistic insights into the physical mechanisms that govern water transfer rates and salt rejection percentages on the atomic scale; one of our novel observations is that ions attracted by hydrophilic functional groups may block the pores, thereby significantly reducing water transfer rates. However, under the same assumed pore density as in the experimental work, our calculated transfer rates are still two orders of magnitude lower than experiment, 10 suggesting there still exist differences between the experimental result and our molecular picture.
The rest of the paper is organized as follows. First, we describe the development of our polarizable nanopore model and force field, the setup of desalination simulations, and the simulation analysis methodology. Next, we describe the optimized parameters and simulation results, including our simulated water transfer and ion rejection rates and discussions of the molecular mechanisms. Finally, we provide some hypotheses that may explain the remaining differences between simulation and experiment.

Polarizable Model
The general functional form of the AMOEBA polarizable force field can be divided into four terms: The first term describe valence interactions which include bond stretching, angle bending, and torsional rotations; they may also include out-of-plane bending terms and higher-order couplings. The last three terms describe the non-bonded interactions, including the van der Waals (vdW) force and the electrostatic contributions from permanent and induced multipoles. Two types of local molecular frames used for defining the permanent dipole and quadrupole moments are shown in the Supporting Information. ( Figure S1 and Table S6).
We constructed an atomistic model of nanoporous graphene passivated by Si atoms in the form of five-membered rings according to previous STEM imaging and density functional studies. 14 Figure 1 shows the molecular models used in QM calculations to develop of force field parameters for the Si-H pore. We constructed two medium-sized molecules, C 26 H 16 Si 2 ( Figure 1a) and C 24 H 14 Si 2 (Figure 1b), to characterize both the Si-C-C-C-Si and the Si-C-C-Si structures in the complete pores. We also created a small model incorporating only one Si atom (SI Figure S2), but it was not included in the final data set because it lacks the interaction between Si passivating groups. For carbon atoms, five different atom types, namely Ca1, Ca2, Cb, Cc, Cs, were created based on the topology and their distance to the pore center, as shown in different colors. Figure 1c illustrates how a complete Si 12 H 24 pore can be built by assembling the two types of medium-sized models. This large model contains 120 carbon, 12 silicon, and 60 hydrogen atoms.  Table 1 summarizes the reference data set for the Si-H pore model. In order to generate data describing the intramolecular forces, we carried out ab initio MD simulations at the B3LYP/6-31G* level of theory 35 using the TeraChem software package. 36,37 Under the simulated temperature of 1500 K, 300 snapshots were generated with an 0.5 fs time interval, from which single-point energies and gradients were computed at the DF-MP2/cc-pVTZ level of theory using the Psi4 software package. 38 We then characterized the non-bonded interactions by computing interaction energies at various model geometries including some The water-"large pore" model contains 195 atoms and 4567 basis functions with def2-TZVP basis set, which limited the number of calculations that we could afford to only four. We included them mainly as a sanity check to ensure the results did not deviate significantly from the medium models. The entire procedure was repeated for the Si-OH pore models, in which all of the SiH 2 terminating groups were replaced by Si(OH) 2 groups. Total 2247 A. DF-MP2/cc-pVTZ on geometries from B3LYP/6-31G* AIMD simulations. B. DLPNO-CCSD(T)/def2-TZVP on geometries optimized at B3LYP/6-311G* level of theory. The graphene carbon atoms more than four bonds away from Si atoms were described using one extra atom type, Ca0, which was not included in these models. The interaction parameters of Ca0 were fitted to coronene-water binding energies (Supporting Information Section S3). Although the modeling of graphene with finite polycyclic aromatic hydrocarbons is imperfect due to the closing of the HOMO-LUMO gap and emerging radical character at the graphene edges, 43 Lazar et al. have shown that for seven small organic molecules, experimental adsorption enthalpies on graphene can be well predicted by ab initio adsorption enthalpies on coronene. 44 To ensure transferability from the finite-sized coronene model systems to the continuous graphene model, we fixed the charge of Ha atom type (C-H) and both charge and dipole parameters for the Ca0 atom type to be zero.
The total charge of the large model was balanced to zero in our parameterization. In keeping with our assumption that the nanoporous graphene should be overall neutral in the simulations, the simulation setup of various pores included a uniform neutralizing charge to all "pore" atoms with a magnitude that never exceeded 0.0001 elementary charge.
Following the generation and organization of ab initio data, the force field parameters were optimized using ForceBalance. 45 An initial set of bonded parameters and vdW parameters were adopted from the carbon and hydrogen parameters of benzene and methane in the AMOEBA09 force field for organic molecules. 24 For the electrostatic interactions, we followed the protocol utilizing the GDMA program, 46 described in the supporting information of the AMOEBA09 paper 24 to generate the initial guess for the charge, dipole, and quadrupole moments. The AMOEBA03 water and ion model parameters were employed in our parameterization and simulation. 22,47   All desalination simulations used a Langevin thermostat with a collision frequency of 0.1 ps −1 to maintain a temperature of 348.15 K. The particle-mesh Ewald (PME) method 22 with cut-off distance at 0.9 nm was adopted when evaluating non-bonded interactions. The systems are first equilibrated for 5 ns with step length of 1 fs; a Monte Carlo barostat with 1 atm pressure was used to adjust the size of the periodic box once every 25 MD steps. All three dimensions of the periodic box were allowed to change individually to relax internal tensions within the graphene model. To prevent ions from leaking during this equilibration period, a virtual blocking force is added that occupies the space of the pore and only repels the ions (Supporting Information Section S5).

Simulation Setup
Following equilibration, the volume is held constant and an external force is applied to the graphene wall to simulate the driving pressure, pushing water from the feed side to the permeate side. Simulation data is collected for 10 ns and repeated ten times from different equilibrated geometries to obtain averages and standard error estimates, for a total of 100 ns per system. Each of the seven Si-H and seven Si-OH pores were simulated at four driving pressures (100, 200, 500 and 1000 bar) for a total of 5.6 µs of simulation data in the main study. We also performed several simulations with the graphene removed and/or ions removed to address specific questions, described in later sections.

Trajectory Analysis
To gain further insights on the water and ion transfer process, we developed a trajectoryanalyzing tool called the "buffered-region method". Our analysis method first estimates the z-coordinate and thickness of the nanoporous graphene sheet. The z-coordinate is computed as the averaged geometric center of C atoms, and the thickness is estimated as the distance where the radial distribution function (RDF) between pore and water atoms reaches the chosen threshold of 0.3. As shown in Figure 5a, the membrane thickness turns out to be 0.28 nm from center to surface, though this is rather insensitive to the chosen RDF threshold.
The simulation box is divided into three sections: feed, middle, and permeate, and the events for each individual water molecule crossing the boundary of two sections are recorded. As illustrated in Figure 5b, a forward transfer event is counted when a water molecule crosses from the feed to middle region, followed by a crossing from the middle to permeate region; the reverse ordering is used to count backward transfer events. Any other crossing sequences, such as feed to middle then back to feed, are ignored. The total exchange rate is defined as the sum of forward and backward transfer rates, while the net transfer rate is the difference between forward and backward transfer rates. The analysis tool also allows us to access the characteristics of individual transfer events including transfer times and changes in the hydrogen bonding environment. Our analysis of the water transfer and ion rejection mechanism necessitate two complementary definitions of pore area. Based on the total water exchange, we calculate the effective pore area as: where A pore;eff is proportional to the total water exchange N pore , and A 0 /N 0 is a normalization factor calculated from water exchange across an unimpeded cross section of the bulk solution.
We also define the accessible pore area by integrating the region inside the pore with a nonzero water density; (Figure 6) in contrast to the former definition, the accessible pore area does not depend on dynamical properties such as the total exchange rate. In practice, the accessible pore area is converged to within 1% after collecting 100 ns of the simulation data. All molecular dynamics simulations were performed using the OpenMM package. 48 The three-dimensional rendering were generated using VMD. 49 Trajectory analyses were carried out in Python using the MDTraj library. 50 3 Results and Discussion

Optimized Force Field Parameters
The final force field with optimized parameters accurately reproduces the QM interaction energy data. For example, Figure 7 shows the SiH-Cl The partial charges and dipole moments of the optimized pore models are shown in

Water Transfer
The net water flux as a function of pore size and driving pressure is plotted in Figure 9a; an increase in net transfer proportional to applied pressure is observed as expected. The effect of functional groups is significant, as the flux for each size of the Si-OH pore is at most half of the corresponding Si-H pore.

Pore Area Comparison
SiH, Accessible SiH, Effective SiOH, Accessible SiOH, Effective SiOH, No Ions, Effective Figure 10: Comparison between accessible pore areas and effective pore areas; the latter is computed using total water exchange under 100 bar driving pressure.
The first hypothesis is tested directly by computing the accessible and effective pore areas. Figure 6 illustartes how the accessible pore area is computed for the SiH 2 -9 pore; the water density pattern is hexagonal, indicating that three SiH 2 groups interact less closely with water than the other six. The effective pore area is computed from the total water exchange as in Eq. 2, and the comparison of two area measurements is shown in Figure 10.
Despite the fact that water density is not equally distributed inside the pores, we found that the accessible pore area for each Si-H pore matches very well with the effective pore area. This agreement indicates that pore permeability is proportional to the accessible pore area, consistent with the macroscopic description of maximum flow rate in the limit of short pipe lengths predicted by Bernoulli's principle, Q max = A 2∆P/ρ. 51 For Si-OH pores, the accessible pore areas are smaller than Si-H by around 5-10Å 2 , but are not enough to rationalize their much smaller effective pore areas; this indicates the larger size of OH functional groups (Hypothesis #1) cannot fully account for the observed difference between the pores.
After carefully inspecting the simulated trajectories, we found that the Si(OH) 2 -9 pore is occupied by three or more Na + ions. (SI section 6) These ions become immobile as they balance the negative partial charges from oxygen atoms on the Si-OH edge groups, and as a result the water transfer through the pore is greatly slowed down. The average transfer time is 100 ps through the Si(OH) 2 -9 pore, compared to 10 ps through the SiH 2 -9 pore. To confirm the effect of ion blockage, we ran a separate set of simulations with the ions removed.
Without the ions blocking the pore, the water flow rates are much higher, and the resulting effective pore area for Si-OH closely matches the accessible pore area. This rules out the kinetic barrier hypothesis (#2) and confirms the ion blockage hypothesis (#3). The interaction between water and pore edges can be further understood by quantifying how the hydrogen bonding network changes as a water molecule passes through the pore. Figure 11 shows the number of hydrogen bonds (N hyd ) each water molecule is involved in, calculated using the hydrogen bond criterion of Baker and Hubbard. 52 We computed a reference value of < N hyd >= 3.6 for the average number of hydrogen bonds in bulk water simulated with the AMOEBA force field at 348 K; the KDE shows a clear peak at n=4 reflecting the underlying histogram. When the analysis is narrowed down to only the molecules undergoing transfer events, the distribution of N hyd is shifted towards smaller numbers. For the smallest and largest Si-H pore, < N hyd >= 1.7 and 3.2 respectively; in particular < N hyd >= 2.88 for SiH 2 -9. The reduced hydrogen bonding indicates that water molecules are partially desolvated as they are transferred through the pore, with < N hyd > decreasing by 0.5-2 depending on pore size. For the Si(OH) 2 -9 pore, the OH group pore edges can serve both as hydrogen bond donor and acceptor, and we expected the transferring water to make more hydrogen bonds than in the corresponding SiH 2 -9 pore; however, the distribution is highly similar to the SiH 2 -9 pore with < N hyd > = 2.83. The reduction in hydrogen bonding is again caused by the Na + ions filling the pore; when ions were removed from the system, the water molecules transferring through the Si-OH pores have < N hyd > = 3.76, even a little higher than bulk water. Additionally, we noticed that the population for zero hydrogen bonds is negligible, suggesting that the evaporation-condensation mechanism 17 is not playing a significant role in these simulations.
We investigated whether water transfer events are correlated in time by examining the time windows in which a water molecule occupies the middle region for all transfer events.
For the SiH 2 -9 pore we found that 70% of water transfer events overlap with other water transfer events in time. However, 70% of the entire trajectory contains at least one water transfer event, which means if the water transfer events are distributed randomly with no correlation, there is still a 70% chance that any one of them would overlap with another.
Therefore, we conclude that existing water transfer events did not have a significant effect on the possibility of the occurrence of another event.  Figure 12: Comparison of the Na + and Cl − rejection rates for Si-H and Si-OH pores, under 100 bar driving pressure (left) and 500 bar driving pressure (right). Negative rejection rates were rectified to zero. Extremely small number of water transfer prevents us from calculating a statistically meaningful rejection rate for Si(OH) 2 -6 pores.

Ion Rejection
The ion rejection rates are calculated by comparing the ion concentration of transferred water molecules and ions relative to the initial feed solution as: where R = 100% corresponds to zero ion leakage, and R = 0% means the transferred molecules have the same or higher ion concentration than in the feed; only transfer events in the forward direction are counted. We find significant differences in ion rejection between positively and negatively charged ions and a strong dependence on the functional groups on the pore edges. The Si-H pores tend to reject Na + ions with a higher percentage than the Cl − ions, as shown in Figure 12; by contrast, the Si-OH pores selectively blocks the Cl − ions and does not block Na + ions at all. These differences can be explained by two distinct mechanisms of ion rejection and the characteristics of Na + and Cl − hydration layers.  Figure S135) shows that a broad peak at low values of the acute angle vanishes for ions inside the pore, indicating that water oxygen atoms in the Cl − hydration layer can undergo significant conformational changes. Figure 13 shows the solvation structure of ions found inside a SiH-9 nanopore and provides intuition for how a more flexible hydration layer around Cl − enables it to permeate the pore more effectively.
(a) (b) Figure 14: Two snapshots of a Na + ion blocking a SiOH-9 nanopore. Left: The Na + hydration layer easily makes hydrogen bonds with the SiOH groups. Right: The SiOH groups can also replace water molecules in the Na + hydration layer.
By contrast, the Si-OH pores repel ions primarily by electrostatic interactions (also called the charge rejection mechanism); the OH groups on the pore edge create a negatively electrostatic potential that attracts cations while exclusively blocking anions. Figure 14 shows two snapshots of a Na + ion occupying the interior of a Si(OH) 2 -9 nanopore; Na + may retain its first hydration shell (left panel), or Si-OH groups may take the place of water molecules and coordinate to Na + directly (right panel). This is further substantiated by a large drop in the Na + -O coordination number from 5.7 in the bulk to 4.0 inside the Si-OH pore, which does not occur for Si-H pores (SI Table S16). The ability of Si-OH groups to substitute for water in the first hydration shell implies that the accessible pore area for Na + inside Si-OH pores is greater than in Si-H pores, and explains how Na + may pass through even the smallest Si(OH) 2 -6 pore with a diameter less than 6Å. By contrast, when a Cl − ion permeates the Si-OH pore (which occurs only rarely), we find there are no strong interactions with the pore edges; a snapshot of this system is shown in SI Figure S136.
Because the permeate solution cannot develop a macroscopic net charge, blocking one type of ion will effectively prevent the large-scale leakage of counter-ions as well; thus we should use the Na + rejection rates for the Si-H pores, and Cl − rejection rates for the Si-OH pores to characterize their overall salt rejection performance. The ion rejection rates for each pore size under different driving pressures are compared in the left and right panel of Figure 12; the ion rejection rates of Si-H pores are nearly 100% for small pores and starts to decrease for SiH 2 -9 and larger pores, in agreement with the size rejection mechanism.
On the other hand, Si-OH pores have nearly 100% ion rejection rates even for the relatively large Si(OH) 2 -11 whose diameter exceeds the hydration shell of Cl − . These results differ from previous studies 9,58 in which the authors concluded a lower salt rejection for hydroxylated (C-OH) pores than hydrogenated (C-H) pores; the discrepancy could be attributed to differences in the pore structures used in the simulations, levels of detail in the force fields, and/or simulation analysis approaches. Ion rejection also decreases with an increase of driving pressure; if we try to maximize the water flux while holding ion rejection at a constant value, we find that smaller pores with higher driving pressures have superior performance to larger pores with lower driving pressure. For example, water flux through SiH 2 -9 under 1000 bar (83 water molecules per ns) is higher than SiH 2 -11 under 200 bar (39 water molecules per ns), yet both have ion rejection rates above 90%.

Discussions
In our simulations, the largest Si-H and Si-OH pores with close to 100% ion rejection rates are SiH 2 -9 and Si(OH) 2 -11. Although Si(OH) 2 -11 has a larger accessible pore area, its average net water transfer rate (8.5 water molecules per ns) is smaller than SiH 2 -9 (12.2 water molecules per ns) under 200 bar of driving pressure; the second value is consistent with a previous study that estimates 10.9 water molecules per ns through C-H pore with 5.5 A diameter. 59 Although the Si-OH pores suffer from cations entering and blocking the pore, the potential transfer rate can be as high as 32.9 water molecules per ns if the ion blockage could be avoided. In experiments, ion blockage was suspected as the main reason causing the degrading of synthesized nanoporous graphene, but smaller ion concentrations were used (6 mM/mol) and the performance degradation happens on the time scale of 24 hours. Another possible route to reduce ion blockage is to reduce the negative charge density by replacing a fraction of OH groups with H or other less charged functional groups; ideally a good balance could be achieved that does not diminish the good Cl − rejection performance.
In the experimental study of Ref., 10 the authors concluded on a massive 3000 water molecules ns −1 pore −1 , more than two orders of magnitude higher than our simulated results.
The large difference might be explained by the large area of nanoporous graphene in direct contact with water. In the experimental setup, the nanoporous graphene prepared by oxygen plasma was suspended over a SiN wafer, in which a 5 µm diameter hole was made for water to pass through, and the hole area was considered as the total nanoporous graphene area that was operational in the desalination experiment. However, the area of the nanoporous graphene in direct contact with water is hundreds of times larger than the area of the 5 µm diameter hole in the SiN wafer. Therefore, if water could be transported laterally between the graphene and SiN wafer, the nanoporous graphene area responsible for the water transfer may be larger than the hole area. This hypothesis can be tested by expanding the 5 µm hole to 10 µm and even larger, and see if the water transfer rate approaches a limit instead of increasing proportionally to the hole area.
The same paper reports a separate experiment where the authors calculate the flux of forward osmosis (water entering the feed side) to be 70 g m −2 s −1 atm −1 . 10 Based on this number and assuming a pore density of 1 per 100 nm 2 , we calculated the initial water transfer rate under 50 atm of osmotic pressure to be 10 molecules per ns per pore. This is in close agreement with our simulated forward osmosis transfer rate with 1 mol/L NaCl solution and no driving pressure, which is 6 water molecules per ns per pore for the SiH-9 pore. The level of agreement between forward osmosis simulations and experiment suggests that the microscopic description of the pore and ion/water transfer is reasonable, and the remarkable desalination performance may originate from macroscopic features of the experimental configuration.
Another possibility for explaining the higher water transfer numbers in experiment is that our optimal pore size could be too small. We estimated that the optimal pores for high flux and salt rejection should have a diameter of about 6Å, which is consistent with earlier theoretical work 9 using non-polarizable force fields that predicted high flux and salt rejection with pores approximately 5.5Å in diameter). It is possible that long-range charge transfer in the graphene layer could allow nanopores to accumulate significant amounts of net charge, thereby allowing more ions to be rejected via the charge rejection mechanism (though Ref. 20 predicts an inverse effect when charges are too high).

Conclusions
In the present work, we investigated the molecular mechanism of water desalination using nanoporous graphene passivated by Si-H and Si-OH groups. Our simulation used a newly developed AMOEBA-type polarizable force field that reproduces high-level ab initio quantum chemistry data with good accuracy. Reverse osmosis simulations were carried out with various pore sizes and driving pressures to measure net water transfer and ion rejection rates; the effective pore area of Si-H pores based on water transfer rates is commensurate with the solvent-accessible pore area, while for Si-OH pores the water transfer is impeded due to Na + ions occupying the pore interior. Size-based and charge-based ion rejection mechanisms were discussed for the Si-H and Si-OH types of pores respectively. Our simulations and analyses revealed complex relationships between water, ions, and pore edges at the nanoporous graphene interface, suggesting that intermediate values of pore-ion interaction strengths may provide a route toward optimizing this system. While our models provide an accurate and physically realistic description of intermolecular interactions between pore edges and water/ions, they do not account for long-range charge transfer and local deformations of the graphene monolayer that may also impact desalination performance. Inclusion of these additional effects would further improve on the modeling of nanoporous graphene and further advance our understanding of this promising system.