Deuteron

Oct 2, 2017 - These quantities strongly (orders of magnitude) depend on the energy below the barrier top, and the proton is more likely to transfer th...
2 downloads 11 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCC

Cite This: J. Phys. Chem. C XXXX, XXX, XXX-XXX

Relevance of the Nuclear Quantum Effects on the Proton/Deuteron Transmission through Hexagonal Boron Nitride and Graphene Monolayers Niranji Thilini Ekanayake,† Jingsong Huang,‡ Jacek Jakowski,‡ Bobby G. Sumpter,‡ and Sophya Garashchuk*,† †

Department of Chemistry & Biochemistry, University of South Carolina, Columbia, South Carolina 29208, United States Center for Nanophase Materials Sciences and Computational Sciences & Engineering Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States

‡

S Supporting Information *

ABSTRACT: According to recent experiments, atomically thin hexagonal boron nitride and graphene are permeable to protons and deuterons (and not to other atomic species), and the experimental estimates of the activation energy are lower than the theoretical values by about 0.5 eV for the isolated proton−membrane transfer model. Our analysis of the electronic potential energy surfaces along the normal to the transmission direction, obtained using correlated electronic structure methods, suggests that the aqueous environment is essential to stabilize the proton transmission, as opposed to the hydrogen atom. Therefore, the process is examined within a molecular model of H2O−H(D)+−material−H2O. Exact quantummechanical scattering calculations are performed to assess the relevance of the nuclear quantum effects, such as tunneling factors and the kinetic isotope effect (KIE). Deuteration is found to affect the thermal reaction rate constants (KIE of 3−4 for hexagonal boron nitride and 20−30 for the graphene) and to effectively lower the barriers to the proton transfer by 0.2 and 0.4 eV for the two membranes, respectively. This lowering effect is reduced for the deuteron by approximately a factor of 3. A more comprehensive description of the proton transmission is likely to require an extended explicit aqueous environment. the literature. The hydrogenated hBN was synthesized recently8 and studied theoretically.9,10 Though the barriers heights were not reported, one would expect a similar trend, i.e., a ratio of about 3−4 in the barrier heights for the hydrogen atom penetrating the hBN monolayer, compared to that for the proton, since hBN is isoelectronic to graphene. Returning to the proton transmission, subsequent experiments with deuterated species11 show that transmission of the protons is 10 times greater than that of the deuterons regardless of the membrane material. These findings are based on the conductance measurements for the device mentioned above1 and on the mass spectrometry of the H2/D2/HD gases in a different device.11 The gases form after the protons and deuterons permeate the Nafion-covered membrane, i.e., monoand bilayer hBN or a monolayer graphene, from a liquid cell to vacuum. The fraction of H/D in the gas phase, compared to the fraction of H/D in liquid, defines the isotope effect on the permeability of a membrane. The remarkable independence of the ratio (equal to 10) of the H/D transmission probabilities on the membrane material has been attributed11 to the difference in the zero-point energy of typical OH/OD bonds in the device

1. INTRODUCTION Two-dimensional (2D) atomically thin crystals, such as hexagonal boron nitride (hBN), graphene, and molybdenum disulfide, are of great interest to experimentalists and theorists as materials for the membranes, suitable for rapid and efficient sieving of gases and liquids and even for separation of hydrogen isotopes. In 2014, Hu and co-workers1 have demonstrated that monolayers of hBN and graphene and a bilayer of hBN are permeable to protons near room temperature. The core of the experimental devices consisted of 2D membranes covered with the proton-conducting polymer Nafion2 topped with the proton-injecting electrodes, enabling measurement of the proton current as a function of voltage. The dependence of the areal conductivity on the temperature, varied in the range of [275, 334] K, was fit with the Arrhenius-type curves from which the activation energies, Ea, were determined. The values of Ea for the proton transmission were estimated as 0.3 ± 0.02 and 0.78 ± 0.03 eV for the monolayers of hBN and graphene, respectively. These activation energies are noticeably lower than the barriers to proton transmission from theoretical studies, reporting values in the ranges of [0.7, 1.0] eV and [1.2, 1.5] eV for hBN and graphene single-layer membranes, respectively.1,3−6 The same materials at room temperature are impenetrable to the hydrogen atom. For example, the barrier for graphene is 4.2 eV,7 and even higher values are reported in © XXXX American Chemical Society

Received: August 15, 2017 Revised: October 1, 2017 Published: October 2, 2017 A

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Table 1. Theoretical Barriers to the Proton Transmission through the hBN and Graphene Monolayersa ref

hBN [eV]

graphene [eV]

method

1

0.3 ± 0.02 0.68

0.78 ± 0.03 1.26 (1.30−1.40) 1.17 1.48 1.56 1.5 1.41 (2.21) 1.6 1.37

Arrhenius fit CI-NEB (AIMD) AIMD PES PES CI-NEB PES (CI-NEB) AIMD PES

3 4 5 6 15 17 this work

1.02 0.91 0.9

0.77

theory

system

DFT DFT MP2 B3LYP DFT DFT DFT LRC-wPBEh

experiment periodic periodic 7 hexagons 19 hexagons periodic periodic periodic 19 hexagons

The energy values in parentheses have been obtained employing the methods in parentheses. “PES” indicates that the PES scans were used to determine the barrier heights. Consult the original references for full details.

a

environment at room temperature, which is the only temperature for which the results have been reported. This implies a purely classical transmission process for the proton. It is noted, however, that the transferring proton is likely to be a part of the hydrogen bonding network, an environment where the H/D nuclear quantum mechanical effects are known to play a role.12,13 Motivated by these intriguing experimental results and by the discrepancies with and within recent theoretical studies, herein we evaluate the nuclear quantum effect (NQE) on the proton/ deuteron transmission through a 2D membrane within a model that includes water. Unlike most prior theoretical studies, we examine the interaction of a proton with hBN and graphene monolayers beyond the proton−membrane adsorption region, in the presence of two water molecules, one on each side of the 2D membrane. Most of the results and discussion are focused on hBN, while comparison with graphene flakes treated at the same theoretical level as hBN is given as appropriate. The membrane is represented as a flat hydrogen-terminated “flake” of hBN, consisting of a few (up to 37) borazine rings. This setup enables us to benchmark the density functional theory (DFT) methods against more reliable correlated methods of electronic structure (ES). The proton and deuteron quantum dynamics restricted to one dimension is carried out using the DFT potential energy curves obtained for the 19-ring flakes of fixed geometry. The ES analysis carried out using Q-Chem 4.3.014 is described in section 2. It is followed by the analysis of the NQE on the membrane permeability, including the kinetic isotope effect (KIE) and isotope dependence of the activation energy based on quantum scattering on a barrier, presented in section 3. A more comprehensive description of the proton transmission is likely to require a chain of water molecules as suggested in ref 11. Our molecular model may be viewed as the last link in such a chain. Finally, the conclusions are drawn in section 4.

ring hBN flake adequately represents the membrane. The 19ring PES is used in the dynamics study of section 3. The experiments of Hu and co-workers1 indicate that the proton transmission is thermally activated with the activation energies of Ea = 0.3 ± 0.02 eV and Ea = 0.78 ± 0.03 eV for the hBN and graphene monolayers, respectively. These values are noticeably lower than the theoretically determined barriers.1,3−6,15−17 The barrier heights, summarized in Table 1, are given for the classical proton going through an infinite or finite planar monolayer. In the latter case (refs 4 and 5 and this work), the cluster size and the theory level are given as comments. Most of the theoretical values fall in the ranges of [0.7, 1.0] eV and [1.2, 1.6] eV for hBN and graphene, respectively. These electronic structure studies are mainly based on a model of an isolated proton approaching a 2D membrane along the normal axis, where the barrier to the proton transfer is defined with respect to the energy of a proton adsorbed on the membrane. An alternative is to define the barrier height with respect to the energy of the proton + hBN asymptotic channel, labeled channel A (H+ + hBN0) in Figure 1. However, the ES

2. ELECTRONIC STRUCTURE ANALYSIS OF THE PROTON TRANSMITTING THROUGH A SINGLE-LAYER HBN AND GRAPHENE In this section, we start with a justification and a description of a model describing the proton transmission through a 2D membrane in the presence of two water molecules, H2O-hBNH+-H2O rather than in gas phase. Then, using a minimal representation of hBN as a single borazine ring, we benchmark the ES methodology and select a practical DFT method for this model. Finally, we use the chosen ES method to examine larger models of hBN, i.e., 7, 19, and 37 rings, and show that the 19-

analysis of this asymptotic channel shows that the energy of removing the proton from the membrane is about 3.7 eV higher compared to the energy of removing the hydrogen atom (channel B: H0 + hBN+). Technically, this makes it difficult to follow the high-energy channel A with the state-of-the-art ab initio quantum chemistry methods, and even for the low-energy channel B, convergence to the correct electronic wave function at large proton-hBN distances is not straightforward. Figure 2a of ref 5 illustrates this challenge: the energy of the channel B at large proton−membrane separation is estimated from the triplet spin-state calculations, in which one electron is fully localized on the proton at all distances. DFT at the B3LYP/6-

Figure 1. Isolated proton transmission model. Separation of H+ beyond the chemisorption energy minimum leads to a high energy channel A, which is more than 3 eV higher than channel B, describing a neutral H atom and a charged hBN flake.

B

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C 31G(d,p) level has been employed. The authors comment on the difficulties of obtaining correct wave functions near and beyond the single/triplet crossing point at the borazine−proton distance of 1.4 Å. Therefore, we start the ES analysis by considering a small borazine−proton system to be able to test wave function methods describing the electron correlation. We find that post Hartree−Fock methods (MP2 and CCSD) are adequate at large (>2 Å) borazine−proton separations, while B3LYP is not, in the following sense. As seen in Figure 2 unlike the CCSD

Figure 3. Solvent effect on the PESs for the borazine−proton model. The charge on the transmitting hydrogen as a function of the borazine−proton distance is shown in the inset. The ES method is CCSD/6-31G** with the polarizable continuum model (PCM) of water as implemented in QChem.14 The PESs correlating with H+ as a fragment at large distances is labeled A. The PESs correlating with the neutral H is labeled B. The asymptotic energies are marked with symbols. The chemisorption minimum is taken as zero for all curves.

mechanism of the proton transfer suggested by experiments. To estimate the solvent effect on the channel B, which becomes a higher energy channel in the solvent, we computed the energy of the fragments within the PCM model. The PCM energies of H and H+ are EPCM = −0.49832Eh and EPCM = −0.19789Eh. H H+ The difference of the PCM energies of B3N3H6+ and B3N3H6 is EhBN+ − EhBN = 0.34517Eh. This makes solvated channel B (neutral H) higher in energy than channel A (H+) by 1.22 eV. The energy level diagram for the gas phase and solvated channels A and B is provided as Supplementary Information. With that, we construct a more realistic molecular model of the proton/deuteron transmission by explicitly including environmental water molecules to stabilize the experimentally relevant reactant channel. We add two water molecules with the proton transferring between the oxygen atoms through the membrane along the normal to hBN path:

Figure 2. Potential energy (a) and charge on the transferring proton (b) as a function of the borazine−H distance. The CCSD and B3LYP results are shown as solid lines and dashed lines, respectively. The calculations were performed for the optimized geometry of the borazine ring employing the 6-31G** basis.

method, the DFT wave function does not converge to the correct asymptotic form of channel B, which is manifested in the sloping tails in the DFT potential energy curve (similar to that in ref 5) and in the partial charge on H at large separations. Therefore, we use the CCSD theory to guide our choice of a DFT method, which will be practical for larger flakes of hBN. As seen in Figure 2, both CCSD and B3LYP methods give the barrier to chemisorption around 2.0 Å, which is near the intersection of PESs of channels A and B (see Figure 3a). Since the asymptotic energy of channel A is much higher (by 3.7 eV) in energy than that of channel B, the natural next step is to include the device environment in some way to lower the energy of the desired channel A. Although Nafion instead of aqueous electrolyte was used in the devices, it is known that Nafion’s functionality relies on percolating water channels spanning the polymer membrane. The proton conduction in Nafion mainly relies on the Grotthuss mechanism18 similar to that in bulk water, in which the proton jumps between water molecules are followed by local rotational rearrangements of water molecules.19 Therefore, it is reasonable to adopt a water environment in the modeling. A simple approach of including water is via the polarizable continuum model (PCM), as implemented in Q-Chem 4.3.0 with the CCSD/6-31G** ES method. Figure 3 shows that indeed aqueous environment stabilizes the proton+borazine wave function and lowers the energy of channel A by 3.0 eV; it also creates a barrier to adsorption of 0.1 eV (with respect to the asymptote) around 2.5 Å. However, within the PCM model the proton transmission barrier at the ring center is lower than the asymptotic energy by 0.7 eV; the barrier to adsorption of 0.1 eV is too small to be consistent with the barrier activation

H3O+−hBN−H 2O → H 2O−hBN−H3O+

as shown in Figure 4 for a single ring. The two water molecules are simple representations of the water channels on both sides

Figure 4. Proposed molecular model of the proton transfer through hBN stabilized by water molecules (red = oxygen, gray = hydrogen, blue = nitrogen, orange = boron). The oxygen to hBN distances are taken as equal. C

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 5. PESs, computed at CCSD/6-31G** level, as a function of the borazine−proton distance. (a) Curves obtained for the fixed hBN-oxygen distance R = [3.5, 8.0] Å as marked on the legend. (b) Oxygen-ring distance is fixed at R = 4.0 Å. The ZPE-corrected energies at the minima are indicated with dashes. (Only the ZPE along the proton transfer coordinate is considered.) The energies of the asymptotic channels are not to scale; their values are annotated.

of hBN flake at r = 0 and barrier II near r = 2.0 Å. The difference between them, Δbar = Ebar I − Ebar II, is about 0.5 eV. The electron density for the two barriers are shown in Figure 6a,c. The local minimum at approximately r = 1.0 Å describes

of the 2D materials. As said above, the proton conduction in the water channels of Nafion mainly relies on the Grotthuss mechanism as in bulk water, which can be essentially considered a lone-pair−lone-pair (LP-LP) relay of protons. Once the proton reaches the water−hBN interface, it resides on the interfacial water molecule’s LP orbital. Then, the proton is transported from water on the left side to hBN by a LP−π relay, because when the proton is given to the hBN sheet it is bound on the π-electron cloud of the hBN sheet. After the proton goes across the hBN sheet, it resides on the π-electron cloud on the other side of the hBN sheet. Then the proton is passed from hBN to the other interfacial water molecule by a π−LP relay. Afterward, the proton enters the water channel on the other side and is again conducted by a LP-LP relay mechanism. Here we only focus on the proton transmission between the two interfacial water molecules and across the hBN sheet. To scan the PES, we start with the geometry of the membrane optimized using a DFT method and fixed henceforth. Next, a water molecule and a hydronium ion are added; the two oxygens, the ring center M, and the transferring proton are aligned along the z-axis normal to the ring. The oxygen-ring distance, R, and the proton-ring distance, r, are fixed (r is fixed so that the proton is at the equilibrium distance from oxygen in a hydronium ion). Then, the positions of the four hydrogens bonded to oxygen are optimized to yield the lowest energy. The scan over the proton-hBN distance r gives a 1D PES shown in Figure 5a for several oxygen-hBN distances R = [3.5, 8.0] Å computed at the CCSD/6-31G** theory level. For R > 4.0 Å we observe plateau regions suggesting PES crossing (Figure 3). Examination of the Mulliken charges on the transferring proton demonstrates that for large R, even in the presence of water molecules, the electron jumps from hBN to the transferring proton to form a neutral H atom in the approximate range of r = [2.2, 6.1] Å and R = 8.0 Å which is effectively the same channel B as in Figure 2. The plateau regions of the curves correspond to the neutral hydrogen and a positively charged ring (channel B). They are adequately modeled at the CCSD and MP2 levels of theory but are not likely to describe a physically relevant situation. At a shorter oxygen-ring distance of R = 4.0 Å shown in Figure 5b, the PES exhibits two barriers: barrier I at the center

Figure 6. Electronic density for borazine-H distance r for R = 4.0 Å. The panels (a−d) correspond to the values r = {0.0, 1.0, 2.0, 3.0} Å, respectively; the ES method is MP2/cc-PVTZ. The distance of r = 2.0 Å corresponds to the curve crossing of PESs describing channels A and B in Figure 3.

the proton bound to the π-electron cloud of the borazine ring; the global minimum near r = 3.0 Å describes the proton bound to the electron lone pair of the water molecule, forming the hydronium + hBN configuration. The minima are labeled as min I and min II, respectively; the corresponding electron density is shown in Figures 6b,d. The adsorption barrier II is located at the charge hopping distance of approximately 2.0 Å, as in Figure 2. In further ES analysis, we use the oxygen-hBN distance of R = 4.0 Å, because this value is close to the hydronium−hBN contact distance of Rc ≈ 3.75 Å. Rc is estimated from the van der Waals radii20,21 as the sum of the radii for hydrogen, nitrogen, and the OH bond distance in the hydronium ion. Note, however, that in applications (for example PCM) the van der Waals radii are often scaled by 10−20% and that the relative heights of the barriers noticeably depend on the flake size and on the ES level. The 1D treatment also ignores the coupling to D

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C lateral modes of the flake, and there is no ZPE correction of the PES. Thus, our theoretical model and the parameter values are appropriated for the qualitative not the quantitative analysis. The PES for R = 4.0 Å and the energies of the asymptotic channels are shown in Figure 5b. For reference, the ZPEs of the global and local minima are estimated from the curvature of the PES for just the proton-transfer mode. The energy of the central barrier with respect to the adsorption minimum, bar I − min I ≈ 1.0 eV, is in the range of most of the prior studies, while the energy difference of the barriers bar I − bar II ≈ 0.5 eV is closer to the experimental value of Ea. The global energy minimum is about 4.0 eV lower than the central barrier, but we argue that this number maybe not directly relevant to the proton transfer process: most likely the proton approaches the membrane through the chain of water molecules. (If so, then the ZPE estimates for OH/OD bonds may also be different.) Further adding an implicit aqueous environment within PCM lowers the global minimum of the PES, while the shape of the PES between the two maxima is essentially unchanged. Therefore, we will focus on the potential energy curves within 2 eV from their respective maxima without implicit solvent. The ES calculations for the borazine−proton system have shown that MP2 and CCSD levels of theory give a physically correct wave functions in gas phase, while the Hartree−Fock and B3LYP do not. For the borazine−hydronium model, we have also checked the CCSD(T)/6-31G** ES method; we find that between MP2, CCSD and CCSD(T) the barrier heights with respect to the local minimum (min I in Figure 5) vary by less than 0.05 eV as listed in Table 2. The CCSD energies will be used for ES benchmarking in the remainder of this section.

distance r. This lowers the central barrier (bar I), but does not change the barrier to adsorption (bar II) and the overall features of the PES. Given the analysis above, we choose LRC-wPBEh/6-31G**, referred to as simply “PBE” henceforth for clarity, as the ES method for examination of larger flakes representing a 2D membrane. The PESs are computed for the fixed-geometry flakes consisting of B12N12H12 (seven hexagons), B27N27H18 (19 hexagons), and B48N48H24 (30 seven hexagons). As seen from Figure 7c, the size effect essentially disappears at the size of 19 rings, and even the 7-ring flake gives the PES parallel to those of larger flakes. The PBE energies at the key points of PES are summarized in Table 3. The PES obtained for 19 borazine rings with two water molecules at fixed geometry and hBN-oxygen distance of 4.0 Å is used to examine proton/deuteron transmission in one dimension employing the exact quantum dynamics approach as in section 3.

3. PROTON/DEUTERON DYNAMICS On the basis of the above PES analysis, in this section we are ready to examine the NQE on the proton/deuteron transmission through hBN and graphene monolayers and the isotope effect on the activation energy and the KIE. 3.1. Transmission through the Barrier. The exact fulldimensional QM studies involving large-amplitude motion of atoms are feasible for isolated chemical systems involving just a few atoms; gas-phase reactions of methane with hydrogen or halogen atoms are the state-of-the-art in this area.22,23 Beyond that, a popular approach of assessment of the transmission through a single barrier, including NQE, is achieved by developing 1D models along the (sometimes modified) reaction coordinate, within the parabolic approximation to the PES, such as the Wigner correction to account for shallow tunneling,24 or the quasiclassical approaches, such as the Wentzel−Kramers−Brillouin (WKB) approximation to account for deep tunneling.25 The merits and limitations of these semiclassical-type approaches are well-known and numerous sophisticated extensions have been developed (for example, a very accurate correction based on ring polymer molecular dynamics);26,27 their discussion is beyond the scope of this work. Let us note that in the current context the WKB approach has been used to estimate the proton/deuteron transmission through graphene monolayer in ref 17. Given the multiple barriers in our model PESs, to assess the magnitude of the H/D isotope effect on the membrane permeability we have performed exact QM dynamics calculations in one dimension, describing the proton transfer along the coordinate (denoted z, |z| = r) normal to the plane of hBN monolayer modeled by 19 borazine rings in the presence of two water molecules. Positions of all atoms except for the transferring proton are fixed. The PES is fitted to an analytical form given in the Supporting Information. The transmission through a barrier is characterized by (i) the energy-resolved probability, N(E), from which (ii) the thermal reaction rate constants, k(T), are obtained. N(E) is computed within the quantum scattering formalism.28,29 (The PES is set to zero for |z| > 3.0 Å.) The transmission probability is defined in terms of the energy eigenstates, equivalent to the plane waves as |z| → ∞ and satisfying the scattering boundary conditions. N(E) is computed from the time-correlation function30 of the reactant and product wavepackets, resolved and normalized in the energy domain. The scattering calculations are performed on the equidistant grid using the

Table 2. Barrier Heights for the Borazine−Hydronium Model Obtained with the Post-HF Methodsa

a

method

bar I − min I

bar II − min I

MP2 CCSD CCSD(T)

1.25 eV 1.26 eV 1.29 eV

0.80 eV 0.83 eV 0.78 eV

The basis is 6-31G**.

The wave function methods such as MP2 and CCSD generally show stronger basis set dependence on molecular geometry than the DFT methods. Thus, to confirm that the entire PES is obtained with comparable accuracy, we have computed the energies at four key points, r = {0, 1, 2, 3} Å, using 6-31G** and larger basis sets, namely, cc-pVTZ and 6311++G(2df,2pd) (not shown) for the seven ring membrane. As seen from Figure 7a, the resulting curves are close to being parallel to each other between the barriers, even though they are shifted with respect to the global minimum (hydronium + hBN) near r = 3.0 Å, which supports our choice of the CCSD/ cc-pVTZ for benchmarking. To describe the proton transmission within a larger, more realistic representation of the 2D membrane, we need an ES method that is more practical than the CCSD theory. Therefore, we have obtained PESs using several DFT methods. The comparison of these PESs to the CCSD curve in Figure 7b shows that unlike B3LYP the long-range corrected (LRC) PBE method, namely, LRC-wPBEh/6-31G**, tracks the CCSD curve rather well. The dispersion correction does not improve the agreement between the CCSD and both DFT methods. The PBE curve is also shown for the relaxed borazine ring, i.e., the geometry of the ring is optimized for each proton-ring E

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 7. Dependence of the PES for the borazine-hydronium model on (a) the basis set, (b) theory level, and (c) flake size. The oxygen−borazine distance is fixed at R = 4 Å. The energy of the global minimum at r = 3.0 Å (min II) is taken as the zero on the vertical axis. The ES methods are marked on the legend. In (b), the extension “-d” stands for the dispersion correction and LRC-wPBEh is labeled “PBE” for short; in (c), the ES method is LRC-wPBEh/6-31G**.

for the energy-resolved probabilities. The wavepacket correlation function C(t), i.e., the overlap of the reactant/product wavepackets

Table 3. Selected Electronic Energies (eV) for the H2O− hBN−H3+O Modela hBN flake energy [eV] min I bar I bar II bar I −min I

B3N3H6 2.37 3.52 3.03 1.15

(2.22) (3.17) (3.01) (0.95)

B12N12H12

B27N27H18

B48N48H24

2.06 2.91 2.88 0.85

1.94 2.81 2.79 0.87

1.90 2.77 2.77 0.87

̂

C(t ) = ⟨ψP|e−iHt / ℏ|ψR ⟩

is computed as a function of time. The scattering matrix elements, SRP, are obtained by Fouriertransforming and normalizing C(t)

a

The ES method is PBE/6-31G**. The oxygen-hBN separation is fixed at R = 4.0 Å. The H-hBN distances of r = {0, 1, 2, 3} Å correspond to the geometries bar I, min I, bar II, and min II, respectively. The values in parentheses are given for the relaxed geometry of the ring. Otherwise the geometry is fixed as that of an isolated hBN as described in text.

SRP(E) =

(3)

(4)

Assuming that the reactant/product Gaussian functions consist of just the incoming/outgoing plane waves, the normalization functions are obtained analytically. ηR/P(E) =

1/4

exp( −α(z ∓ z 0)2 + ip0 (z ∓ z 0))

∞

∫−∞ C(t )eiEt /ℏ dt

N (E) = |SRP(E)|2

split-operator/fast Fourier transform. A parabolic absorbing potential is used to prevent reflection off the grid boundaries at long times. The parameters pertaining to the dynamics calculation are listed in Table S2. The reactant/product wavepackets ψR/P are taken as Gaussian functions, ⎛ 2α ⎞ ⎜ ⎟ ⎝π ⎠

(2π ℏ)−1 ηP*(E)ηR (E)

The transmission probability is simply

31,32

ψR/P(z) =

(2)

m 2π ℏk

∞

2 2

∫−∞ e−ikzψR/P(z) dz , E = ℏ2mk

(5)

where m is the mass of the particle. This theoretical framework (including the PES construction) is used to simulate proton/ deuteron transmission through the 19-ring flake of hBN and, for comparison, of graphene. The PESs for the 19-ring flake of hBN and graphene (for the oxygen-membrane distance R = 4.0 Å) are shown in Figure 8a. Note that the maximum of the graphene PES is 0.5 eV higher

(1)

localized to the right/left of the barrier. The initial momentum p0 is chosen to have appreciable transmission through the barrier. The width parameter α is chosen to reproduce the vibrational eigenstate of the OH bond, but it is not necessary F

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 8. (a) 1D electronic PES as a function of the proton−membrane separation. The membrane is modeled with 19 rings. The PES for H2O/ hBN(graphene)/H3O+ at fixed geometry is given for the oxygen-membrane distance of R = 4 Å, and the global minimum of this PES is taken as zero. These PESs are labeled as “graphene” and “hBN”. The PESs for 19-ring hBN(graphene)−proton (no water molecules), shown with symbols and labeled with the superscript “p”, are shifted to align the central barriers of the respective curves. Note the overlap of PESs for the models with and without water molecules in the central barrier region. (b) Energy-resolved transmission probabilities for hBN (black) and graphene (green). Dashes indicate the curves for deuterium.

their ratio for the hydrogen and deuterium, KIE = kH/kD. Formally, the reaction rate constant is defined in terms of N(E) as33

than the maximum of the hBN PES which is consistent with the difference in the activation energies determined experimentally.1 Within our molecular model, we also see the difference in the shape of the barrier: For graphene there is an obvious global maximum at r = 0 (bar I), while for hBN the PES maxima are of nearly equal height. In the same figure, the curves labeled with the superscript p are the PESs for the membrane-proton without the water molecules, aligned with the membrane−hydronium curves at r = 0. For r < 0.75 Å, once the proton strongly interacts with the membrane, the PESs with and without water molecules match very well as expected: Interaction of the proton with the membrane at short separations should not depend on the environment. The energy-resolved N(E) for the proton and deuteron is shown in Figure 8b; the triple barrier for hBN results in the resonant features of N(E), absent in the case of graphene monolayer. The probabilities for the proton and deuteron and their ratio follow the expected trend. These quantities strongly (orders of magnitude) depend on the energy below the barrier top, and the proton is more likely to transfer through the membrane than the deuteron. When E is comparable to the barrier height, shallow tunneling and above-the-barrier reflection result in the nonmonotonic behavior of N(E), as seen even on the log-scale in panel Figure 8b for the hBN system. Another interesting feature of the hBN probabilities are resonant peaks for E = [2.8, 2.9] eV. N(E) in this range, however, is low, and these subtle features are typically suppressed in the thermal averaging, yielding the reaction rate constants. Overall, comparing N(E) for the hBN and graphene for the PESs of Figure 8, i.e., 19-ring flakes with oxygen distance fixed at R = 4 Å, the below-the-barrier transmission is orders of magnitude lower for hBN due to the triple barrier rather than a single barrier. Therefore, we expect a more classical transmission and a smaller proton/deuteron isotope effect for hBN membrane compared to that for graphene after N(E) is thermally averaged. 3.2. Thermal Reaction Rate Constants and the Isotope Effect. We now turn to calculation of the thermal reaction rate constants k(T) and the kinetic isotope effect (KIE) defined as

k q (T ) =

1 2πQ (T )

∫0

∞

N (E)e−βE dE , β = (kBT )−1

(6)

where β is the inverse temperature variable and kB is the Boltzmann constant. Q(T) is the partition function of the reactants, which in 1D is simply that of translation Q (T ) ≈

2πmkBT

(7)

Denoting the barrier height as Vb, and treating the particle as classical, i.e., N(E) = 1 for E > Vb and zero otherwise), the rate constant becomes kcl(T ) =

kBT −βVb e 2πQ

(8)

The ratio of the quantum and classical rate constants defines the tunneling factor κ κ (T ) =

k q(t ) kcl(T )

(9)

In practice, the experimental rate constants for a surprising variety of reactions are well described by the empirical Arrhenius equation k A(T ) = A e−βEa

(10)

where Ea is the activation energy. The exponential term describes the dominant temperature dependence of k(T), and the pre-exponential factor A describes an additional empirical dependence due to the collision frequency orientation; the quantum tunneling correction can be introduced as a prefactor, such as those in refs 34 and 35 as well. Given the similarity of the classical and Arrhenius expressions for k(T) (eqs 8 and 10), Ea is often equated with Vb within the approximation of temperature independent pre-exponential factors as done in ref.11 For a limited range of experimental temperature, Ea can be defined by the slope of the Arrhenius plot of ln(k(T)) on the G

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C vertical axes vs β on the horizontal axes. Fitting the experimental Arrhenius plot to a straight line over the temperature range of [275, 334] K, gave the activation energies of 0.3 ± 0.02 and 0.78 ± 0.03 eV for the proton transmission through the hBN and graphene monolayers.1 These experimental Ea are lower than the theoretical Vb by 0.5 eV or more as discussed in section 2. While the ab initio determination of the reaction barriers is, generally, prone to large errors, our ES benchmarking and the level of agreement with and across prior theoretical studies1,3−6,15−17 give credibility to the computed barrier heights Vb with respect to the chemisorption minimum. Therefore, let us consider that Ea is not equal to Vb (and that there is no straightforward relation between the two)36 and define Ea from the quantum N(E) computed for the 1D models. The general definition of the activation energy EA (β) = −

d ln k(β) dβ

Figure 9. Scattering on a barrier in 1D for the proton and deuteron: (a) rate constant as a function of inverse temperature; (b) activation energy obtained from the rate constant computed for quantum H, quantum D, and classical particle, using A = kBT, shown for experimental range of temperatures. T = 300 K is marked with the vertical dash.

(11)

is used to theoretically estimate the difference between Ea and Vb, using the formal definition of k(T) given by eq 6. The collision kinetics arguments (Q ≈ mβ is defined by eq 7), and the transition state theory arguments (Q ≈ β) both give Ea, which is lower than Vb and slowly depends on temperature over the experimentally relevant range (T = [275, 334] K) which in the Arrhenius plot unit of 1000/K is equivalent to the range of [3.0, 3.65]. The prefactors cancel in the KIE and κ, and the values of the activation energy at room temperature, computed with two forms of the pre-exponential factor, are so close to each other that we give only the data for the latter form. The room-temperature KIE and the activation energies were determined for three pairs of potentials computed for 19-ring flakes of hBN and graphene: (i) for the membrane−proton model, (ii) for the membrane−hydronium model with oxygenmembrane distance of R = 4.0 Å (shown in Figure 8), and (iii) for the membrane−hydronium model with oxygen−membrane distance of R = 3.7 Å. In each case C(t), N(E), and k(T) were computed for the quantum proton and deuteron. The reaction rate constants are also determined using classical N(E). For illustration, the quantum and classical k(T) for the graphene− hydronium PES at R = 4.0 Å is shown in Figure 9a as a function of 1000/K for the range of T = [250, 500] K. The reaction rate constant for quantum hydrogen is higher than that for deuterium, with the latter still higher than the classical value. The tunneling coefficients κ and KIE at 300 K, listed in Table 4, consistently show that the NQE is higher for hydrogen than for deuterium and that KIE is higher for graphene than for hBN membranes. In the case of graphene, the KIE is on the order of 20−30, and in case of hBN, it is on the order of 1.5−4 depending on the particular PESs. Comparing the results for the PESs (i) and (iii), we conclude that the KIE correlates with the shape of the central barrier. This barrier is effectively blocked by the adsorption barrier for the hBN PES (ii) defined for R = 0.0 Å, which explains the low value of the KIE in this case. In a realistic dynamics calculation the range of R should be sampled, similar for example to a proton transfer study in SLO1,37,38 and the fact that the KIE obtained for the model PESs here is of the same order of magnitude indicates that our model is reasonable. We plan to perform a direct dynamics study in the near future.

Having computed k(T) we estimate the activation energies using a formal definition shown in eq 11. Ea is computed for both isotopes and shown along with the barrier height in Figure 9b for experimentally relevant range. We observe that Ea weakly depends on temperature, and also on the isotope. The NQE lowers Ea compared to the barrier height by ∼0.4 eV for the proton−graphene models and by 0.15−0.2 eV for the proton− hBN models. Thus, this effect could be used to explain the difference between the experimental Ea and theoretical barrier heights. The trend for graphene is consistent with that of ref 17 where the Ea reduction compared to the barrier height has been obtained using 1D WKB analysis, free energy thermodynamic integration on a 3D PES and the path-integral-based centroiddensity quantum transition state theory. (The Ea reduction is 1.2 eV for the latter method.) Given the agreement of our results with the WKB results, reported in ref 17 for just graphene, we conclude that the effect of Ea reduction is mostly due to the properties of the central barrier. Our values for the Ea reduction, given in the columns Δ = V0 − Ea, are consistent among the PESs (with the exception of hBN at R = 4 Å discussed above). The values Δcl are computed using the classical rate constants and are zeros within two decimal places. The relevant energy barriers are listed as central columns of Table 4. The role of the NQE and of the non-Arrhenius behavior of the reaction rate constants at room temperature could be tested by performing the reaction rate measurements for the protonated and deuterated systems.

4. SUMMARY To summarize, motivated by recent experiments1,11 we have examined NQE including KIE on the proton/deuteron transmission through the hBN and graphene monolayers within the H2O−membrane−H3O+ model which stabilizes the H+ + hBN channel relative to the other one, H0 + hBN+, and allows one to explore the PES beyond the physisorption minimum, which was a limitation of the previous ES analysis. This model is relevant to the experimental set up, where the membranes are covered with the Nafion layer. The proton conduction in the water channels of Nafion mainly relies on a lone-pair−lone-pair (LP-LP) relay. Once the proton reaches the water−hBN interface, it resides on the interfacial water H

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Table 4. Transmission through the Barrier in 1D at T = 300 Ka V0 − Ea [eV] system

ΔHa

ΔDa

hBN graphene

−0.18 −0.41

−0.06 −0.12

hBN graphene

−0.01 −0.44

0.01 −0.14

hBN graphene

−0.15 −0.37

−0.05 −0.11

PES [eV] Δcla

V0

Vads

(i) Membrane−Proton Model 0.00 0.77 0.00 1.37 (ii) Membrane−Hydronium Model R = 4.0 Å 0.00 2.92 0.88 0.00 3.41 1.58 (iii) Membrane−Hydronium Model R = 3.7 Å 0.00 2.88 0.87 0.00 3.47 1.53

NQE κH

κD

KIE

9.3 24.7

2.4 4.8

3.9 25

0.0 0.74

0.54 210

0.34 5.9

1.6 35

0.40 1.11

7.8 93

2.4 5.0

3.2 19

Vbar

Δa = Vb−Ea; V0, Vads, and Vbar are the central barrier values relative to the global minimum, to the adsorption minimum and to the barrier II, respectively. κH = kH/kcl is the tunneling coefficient (defined analogously for D), KIE = kH/kD.

a

molecule’s LP orbital. After testing several wave function and DFT electronic structure methods on a borazine−hydronium model, we conclude that correlated methods (MP2 or better) are needed to benchmark more practical DFT methods. We found that LRC-wPBEh/6-31G** method reproduces the features of the PES within 2 eV of the barrier to the proton transfer and used it to model larger (up to 37 rings) flakes of the material. The PES along the proton transfer coordinate was deemed converged with respect to flake size at 19 rings. We have studied NQE on the proton/deuteron transmission within this 1D model using exact quantum scattering approach and found that the nuclear quantum effects influence the thermal reaction rates and lower the activation energy by 0.4(0.1) eV for graphene and by 0.18(0.06) eV for hBN models, where the values in parentheses refer to the deuteron case. This effect likely contributes to the discrepancy between Ea determined under the Arrhenius-behavior assumption experimentally1 and theoretical barrier heights: Ea was found to be about 0.5 eV lower than the theoretical estimates of the barrier heights with respect to the adsorption minima. The isotope dependence of the Ea, associated with the NQE on the proton/deuteron transmission through the barrier, could be tested in experimental measurements of the reaction rate constants for the protonated and deuterated systems. Our estimates for the room temperature KIE are 3−4 for hBN and 20−30 for graphene, which is different from the experimental values of 10 for both systems. We note, however, that 1D models tend to exaggerate quantum effects and that dynamics studies with the proton/deuteron as 3D particles and flexible membranes employing on-the-fly quantum-trajectory−ES approach39 are planned for the future.

■

Bobby G. Sumpter: 0000-0001-6341-0355 Sophya Garashchuk: 0000-0003-2452-7379 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This material is based upon work supported by the National Science Foundation under Grant No. CHE-1056188 and CHE1565985 and by an ASPIRE grant from the Office of the Vice President for Research at the University of South Carolina. Part of the work was conducted at the Center for Nanophase Materials Sciences, a U.S. Department of Energy Office of Science User Facility. The research used an XSEDE allocation TG-DMR110037, resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC0205CH11231 and the USC HPC cluster, funded by the National Science Foundation under Grant No. CHE-1048629. We thank Professor Vitaly Rassolov for many helpful discussions. S.G. acknowledges Nanomaterials Theory Institute of the Center for Nanophase Materials Sciences at ORNL for hosting her as a short-term visitor, while on sabbatical leave from the University of South Carolina.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b08152. Parametrization of the potential energy surfaces and details of the time-dependent dynamics of section 3 (PDF)

■

REFERENCES

(1) Hu, S.; Lozada-Hidalgo, M.; Wang, F. C.; Mishchenko, A.; Schedin, F.; Nair, R. R.; Hill, E. W.; Boukhvalov, D. W.; Katsnelson, M. I.; Dryfe, R. A. W.; Grigorieva, I. V.; Wu, H. A.; Geim, A. K. Proton Transport Through One-Atom-Thick Crystals. Nature 2014, 516, 227−230. (2) Mauritz, K. A.; Moore, R. B. State of Understanding of Nafion. Chem. Rev. 2004, 104, 4535−4585. (3) Wang, W. L.; Kaxiras, E. Graphene Hydrate: Theoretical Prediction of a New Insulating Form of Graphene. New J. Phys. 2010, 12, 125012. (4) Zhang, Q.; Ju, M.; Chen, L.; Zeng, X. C. Differential Permeability of Proton Isotopes Through Graphene and Graphene Analogue Monolayer. J. Phys. Chem. Lett. 2016, 7, 3395−3400. (5) Seel, M.; Pandey, R. Proton and Hydrogen Transport Through Two-Dimensional Monolayers. 2D Mater. 2016, 3, 025004. (6) Kroes, J. M. H.; Fasolino, A.; Katsnelson, M. I. Density Functional Based Simulations of Proton Permeation of Graphene and Hexagonal Boron Nitride. Phys. Chem. Chem. Phys. 2017, 19, 5813− 5817. (7) Tsetseris, L.; Pantelides, S. T. Graphene: An Impermeable or Selectively Permeable Membrane for Atomic Species? Carbon 2014, 67, 58−63.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jingsong Huang: 0000-0001-8993-2506 Jacek Jakowski: 0000-0003-4906-3574 I

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (8) Spath, F.; Gebhardt, J.; Dull, F.; Bauer, U.; Bachmann, P.; Gleichweit, C.; Gorling, A.; Steinruck, H.; Papp, C. Hydrogenation and Hydrogen Intercalation of Hexagonal Boron Nitride on Ni(1 1 1): Reactivity and Electronic Structure. 2D Mater. 2017, 4, 035026. (9) Tang, S.; Cao, Z. Structural and Electronic Properties of the Fully Hydrogenated Boron Nitride Sheets and Nanoribbons: Insight from First-Principles Calculations. Chem. Phys. Lett. 2010, 488, 67−72. (10) Zhou, J.; Wang, Q.; Sun, Q.; Jena, P. Electronic and Magnetic Properties of a BN Sheet Decorated with Hydrogen and Fluorine. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 085442. (11) Lozada-Hidalgo, M.; Hu, S.; Marshall, O.; Mishchenko, A.; Grigorenko, A. N.; Dryfe, R. A. W.; Radha, B.; Grigorieva, I. V.; Geim, A. K. Sieving Hydrogen Isotopes Through Two-Dimensional Crystals. Science 2016, 351, 68−70. (12) Li, X.-Z.; Walker, B.; Michaelides, A. Quantum Nature of the Hydrogen Bond. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 6369−6373. (13) Ceriotti, M.; Cuny, J.; Parrinello, M.; Manolopoulos, E. Nuclear Quantum Effects and Hydrogen Bond Fluctuations in Water. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 15591−15596. (14) Shao, Y.; Rassolov, V. A.; Slipchenko, L. V.; Subotnik, J. E.; Van Voorhis, T.; Herbert, J. M.; Krylov, A. I.; Gill, P. M. W.; Head-Gordon, M. Advances in Molecular Quantum Chemistry Contained in the QChem 4 Program Package. Mol. Phys. 2015, 113, 184−215. (15) Miao, M.; Nardelli, M. B.; Wang, Q.; Liu, Y. First Principles Study of the Permeability of Graphene to Hydrogen Atoms. Phys. Chem. Chem. Phys. 2013, 15, 16132−16137. (16) Xin, Y.-B.; Hu, Q.; Niu, D.-H.; Zheng, X.-H.; Shi, H.-L.; Wang, M.; Xiao, Z.-S.; Huang, A.-P.; Zhang, Z.-B. Research Progress of Hydrogen Tunneling in Two-Dimensional Materials. Acta. Phys. Sin. 2017, 66, 56601. (17) Poltavsky, I.; Zheng, L.; Mortazavi, M.; Tkatchenko, A. Quantum Tunneling of Thermal Protons Through Pristine Graphene. 2017, arxiv.org:1605.06341v2. arXiv.org e-Print archive. https://arxiv. org/abs/1605.06341v2. (18) Cukierman, S. Et Tu Grotthuss! and Other Unfinished Stories. Biochim. Biophys. Acta, Bioenerg. 2006, 1757, 876−885. (19) Devanathan, R. In Energy Production and Storage - Inorganic Chemical Strategies for a Warming World; Crabtree, R. H., Ed.; Wiley: Chichester, U.K., 2010; pp 89−100. (20) Bondi, A. Van Der Waals Volumes + Radii. J. Phys. Chem. 1964, 68, 441−451. (21) Rowland, R. S.; Taylor, R. Intermolecular Nonbonded Contact Distances in Organic Crystal Structures: Comparison with Distances Expected from Van Der Waals Radii. J. Phys. Chem. 1996, 100, 7384− 7391. (22) Yang, M.-Y.; Yang, C.-L. Isotope Effects of Quantum Dynamics for the Cl+CH4/CD4 Reaction. J. Mol. Struct.: THEOCHEM 2010, 945, 23−26. (23) Schiffel, G.; Manthe, U. Quantum Dynamics of the H+CH4 → H2+CH3 Reaction in Curvilinear Coordinates: Full-Dimensional and Reduced Dimensional Calculations of Reaction Rates. J. Chem. Phys. 2010, 132, 084103. (24) Wigner, E. P. On the Quantum Correction for Thermodynamic Equilibrium. Phys. Rev. 1932, 40, 749−759. (25) Landau, L. D.; Lifshitz, E. M. Quantum Mechanics; ButterworthHeinemann: Oxford, U.K., 1999. (26) Zhang, Y.; Rommel, J. B.; Cvitas, M. T.; Althorpe, S. C. ShallowTunnelling Correction Factor for Use with Wigner-Eyring TransitionState Theory. Phys. Chem. Chem. Phys. 2014, 16, 24292−24300. (27) Craig, I. R.; Manolopoulos, D. E. Chemical Reaction Rates from Ring Polymer Molecular Dynamics. J. Chem. Phys. 2005, 122, 084106. (28) Newton, R. G. Scattering Theory of Waves and Particles; SpringerVerlag: New York, 1982. (29) Tannor, D. J. Introduction to Quantum Mechanics: A TimeDependent Perspective; University Science , Books: Sausalito, CA, 2006. (30) Tannor, D. J.; Weeks, D. E. Wave Packet Correlation-Function Formulation of Scattering Theory - the Quantum Analog of Classical S-Matrix Theory. J. Chem. Phys. 1993, 98, 3884−3893.

(31) Kosloff, R. Time-Dependent Quantum-Mechanical Methods for Molecular Dynamics. J. Phys. Chem. 1988, 92, 2087−2100. (32) Feit, M. D.; Fleck, J. A., Jr.; Steiger, A. Solution of the Schrödinger Equation by a Spectral Method. J. Comput. Phys. 1982, 47, 412−433. (33) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. Quantum Mechanical Rate Constants for Bimolecular Reactions. J. Chem. Phys. 1983, 79, 4889−4898. (34) Garrett, B. C.; Truhlar, D. G. Accuracy of Tunneling Corrections to Transition State Theory for Thermal Rate Constants of Atom Transfer Reactions. J. Phys. Chem. 1979, 83, 200−203. (35) Skodje, R. T.; Truhlar, D. G. Parabolic Tunneling Calculations. J. Phys. Chem. 1981, 85, 624−628. (36) Hand, M. R.; Rodriquez, C. F.; Williams, I. H.; Balint-Kurti, G. G. Theoretical Estimation of the Activation Energy for the Reaction HO + H2O → H2O + OH: Importance of Tunneling. J. Phys. Chem. A 1998, 102, 5958−5966. (37) Iyengar, S. S.; Sumner, I.; Jakowski, J. Hydrogen Tunneling in an Enzyme Active Site: A Quantum Wavepacket Dynamical Perspective. J. Phys. Chem. B 2008, 112, 7601−7613. (38) Mazzuca, J.; Garashchuk, S.; Jakowski, J. Description of Proton Transfer in Soybean Lipoxygenase-1 Employing Approximate Quantum Trajectory Dynamics. Chem. Phys. Lett. 2012, 542, 153−158. (39) Garashchuk, S.; Jakowski, J.; Wang, L.; Sumpter, B. G. Quantum Trajectory-Electronic Structure Approach for Exploring Nuclear Effects in the Dynamics of Nanomaterials. J. Chem. Theory Comput. 2013, 9, 5221−5235.

J

DOI: 10.1021/acs.jpcc.7b08152 J. Phys. Chem. C XXXX, XXX, XXX−XXX