Density Functional Based Simulations of Proton Permeation of Graphene and Hexagonal Boron Nitride

TL;DR

This study uses density functional theory to investigate proton permeation through graphene and hexagonal boron nitride, analyzing various factors affecting permeation barriers and charge transfer, but finds theoretical results do not fully match experimental observations.

Contribution

It provides a detailed computational analysis of proton permeation mechanisms and factors influencing barriers in graphene and h-BN, highlighting discrepancies with experimental data.

Findings

Permeation barriers are influenced by structural and environmental factors.

Charge transfer from membrane to proton is significant but does not resolve theory-experiment discrepancies.

None of the considered effects fully explain the observed proton permeation rates.

Abstract

Using density functional theory, we study proton permeation through graphene and hexagonal boron nitride. We consider several factors influencing the barriers for permeation, including structural optimization, the role of the solvent, surface curvature and proton transport through hydrogenated samples. Furthermore, we discuss the ground state charge transfer from the membrane to the proton and the strong tendency for bond formation. If the process is assumed to be slow we find that none of these effects lead to a satisfactory answer to the observed discrepancies between theory and experiment.

Figures

b) which brings the atoms closer in height to follow the (moving) SMD equilibrium spring position. This shows that it is favourable to bend the substrate rather than to directly approach the surface. The resulting work done is shown in Fig. [4(f)](#S4.F4.sf6). The barrier estimated from these curves contains two separate effects. The first involves the removal of the proton from the liquid and the second involves passing through the membrane. These steered dynamics simulations are equivalent to the Tomlinson model (see e.g. Ref. [24](#bib.bib24) and references therein) used in the field of fri

Equations

$$F^{\mathrm{guide}}=-\kappa {\left[{\xi }^{\mathrm{pr}}(t)-{\xi }^{\mathrm{pr}}(0)-v_{\xi }t\right]}^2,$$

$$W(\xi )=\kappa v{\int }_0^{t}d{t}^{\prime }\left[\xi -{\xi }^{\mathrm{target}}(t^{\prime })\right]d{t}^{\prime }.$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: Radboud University, Institute for Molecules and Materials, Heijendaalseweg 135, 6525 AJ Nijmegen, The Netherlands.

1 Introduction

Graphene and hexagonal boron nitride (h-BN) are chemically inert and strongly hydrophobic materials. As a result, it has generally been expected that no ions could penetrate their dense electron clouds (Fig. 1). However, recently it was nevertheless observed 1 that it is possible for protons to do precisely that.

In the experiments by Hu et al. 1, hydrogen was split into electrons and protons, with protons passing through the membrane and electrons passing through an electrical wire under the influence of an applied bias voltage. These protons and electrons then recombine into hydrogen again on the other side. From the resulting IV-curve the conductivity could then be extracted. The temperature dependence of the conductivity was used to determine the energy barrier for permeation. These experiments are characterized by a highly reproducible and linear IV-behaviour, indicative of an independence of specific details of the sample such as defects.

Hu et al. found proton permeation to be easiest for monolayer h-BN but also observable for bilayer h-BN and monolayer graphene 1. Notably, no measurable conductance was observed for bilayer graphene, which was understood to be a result of its AB layered stacking. From the near-perfect Boltzmann behaviour of the conductance as a function of temperature the associated barriers were estimated to be 0.3, 0.6 and 0.8 eV for monolayer h-BN, bilayer h-BN and monolayer graphene respectively. At a typical bias voltage of 0.1 V, the resulting proton current through h-BN corresponds to roughly 600 protons/s/nm2. We therefore consider this process to be slow in comparison to typical ionic timescales and attempt rates.

More recently it has also been observed that the permeation barrier depends on the ionic mass 2, making these materials naturally selective membranes for separating protons and deuterons. This was explained 2 as originating from the isotope effect, i.e. the difference in zero-point energies of the proton and deuteron. The zero-point energy, $\Delta E^{0}$, comes into play only for the initial state, where the proton is bound to an oxygen atom. In water, $\Delta E^{0}$ is principally determined by the OH- stretching frequency, which is estimated to be $\Delta E^{0}=0.20$ eV for protons and $\Delta E^{0}=0.14$ eV for deuterons based on literature values 3, 2. In DFT, the barrier ($\Delta E$) is calculated as difference between the transition state and the initial state (without zero-point energy), whereas experimentally, the barrier is calculated between the transition state and the initial state which is raised by its zero-point energy with respect to the value given by DFT. Therefore, the barrier calculated within DFT is larger than the experimental barrier which naturally includes the zero-point motion, namely $\Delta E^{\textrm{exp}}=\Delta E^{\textrm{DFT}}-\Delta E^{0}.$

This selective permeation of protons versus other ions as well as the separation of protons and deuterons, has multiple direct technological applications making a fundamental understanding of this process important (see also Refs. 1, 2 and references therein). The most obvious application lies perhaps in proton exchange membrane fuel cells, which can be scaled down in size by the use of such nano-sized membranes. Another application may be found in the separation of protons and deuterons for example in the production of heavy water. This is typically a complicated and costly process which may be greatly simplified by the use of naturally semi-permeable membranes.

On the theoretical side, several efforts have been undertaken to understand the proton permeation of 2D membranes 1, 4, 5, 6, most of them using density functional theory (DFT). However, significant discrepancies remain between the barriers found in experiment and those calculated. For example, for graphene, calculated values are 1.17 eV 4, 1.41 or 2.21 eV 5 (the latter after relaxation), 1.25–1.40 eV 1 or 1.56 eV 7, all of which overestimate the experimental value of 1.0 eV 2. More severe discrepancies arise when attempting to improve upon these calculations as discussed in the following sections.

One may question the accuracy of the applied theoretical methods for the problem at hand and whether DFT based calculations are indeed suitable to describe the permeation of a proton through a highly polarizable membrane or the quantum nature of the ions or excited states should be taken into account. Classical molecular dynamics simulations 8 using the REBO potential 9 found much higher barriers for hydrogen permeation (11 eV for REBO compared to 2.5 eV for DFT) illustrating the need to go beyond classical force fields for quantitative estimates of the barriers. Recent calculations also aimed to take into account the quantum nature of the ions by means of path integral molecular dynamics 10. These calculations predict lower barriers as result of quantum tunneling but at the same time predict a much larger isotope effects ($\sim$1 eV instead of 0.06 eV) than observed experimentally. Therefore, despite the claims made in Ref. 10, the puzzle of low barriers is still unsolved.

Aside from these theoretical considerations, many of the details of the experimental setup may prove essential. The theoretical difficulty to reproduce these barriers stems at least partially from the difficulty in identifying the relevant experimental parameters. For example, the role of the solvent (Nafion) in which the proton diffuses is largely unknown. The effect of temperature, surface curvature or applied electric field may also prove crucial. Finally, it can also not be excluded that (ad-)atoms or other defects are present that may modify the barrier as suggested by experiments with varying pH concentrations 11, 12.

2 Methods

Here we consider various factors relevant to the calculated proton permeation barriers within the framework of DFT. Our calculations were done using the non-empirical PBE exchange-correlation functional 13, 14 implemented in the CP2K code 15. To account for van-der-Waals interactions, especially relevant for the simulations of the membrane in water, we use the Grimme-D3 dispersion correction term 16. Unless mentioned, for graphene (h-BN) we use a model consisting of 6x12 orthorhombic cells, corresponding to 288C (144B and 144N) atoms and a supercell of 25.6$\times$29.5 Å2 (26.2$\times$30.3 Å2) with periodic boundary conditions. For the perpendicular cell size we use 15 Å. The lattice constant of graphene and h-BN are fixed such that $r_{\textrm{CC}}=1.421$ Å and $r_{\textrm{BN}}=1.457$ Å when flat. We considered non-spin polarized calculations at the $\Gamma$-point of the Brioullin zone , used a Fermi-Dirac electronic smearing with a width of 300 K. The Quickstep method is employed, with wave functions expanded onto a localized double-$\zeta$-valence-polarized basis set and the electronic density expanded onto a plane-wave basis set with a kinetic energy cutoff of 500 Ry. Goedecker-Teter-Hutter pseudopotentials 17 are used to describe the interactions with the core electrons.

For the optimization of the barrier we used climbing-image nudged elastic band (CI-NEB 18) method with 18 images, of which one was the climbing image.

We also performed Born-Opperheimer molecular dynamics (MD) simulations of a proton, in water solution, permeating through graphene or h-BN. For this we used a velocity-Verlet integration timestep of 0.5 fs. These simulations were done using using a periodic model of graphene or h-BN consisting of 60 atoms with a solvent consisting of 88 water molecules in an orthorhombic super cell, corresponding to a 1 kg/L density of water when taking into account a 3.5 Å vacuum due to the hydrophobic nature of the membrane. The cell sizes are 12.8$\times$12.3$\times$20.2 Å3 for graphene and 13.2$\times$12.7$\times$19.4 Å3 for h-BN. In order to have well converged wavefunctions, a 1 ps NVE benchmark simulation was performed to determine the necessary convergence criteria to have total energy drifts of at most 1 meV/ps, which was found to be $\epsilon^{\textrm{SCF}}=10^{-7}$ Ha. For temperature control we used the coloured noise NVT thermostat 19 with a coupling constant of 1 ps and a constant temperature of 325 K. After an equilibration period of 30 ps, a proton is pulled through the membrane by means of a steered molecular dynamics (SMD) simulation in 2 ps.

To compute the free energy barrier at room temperature with SMD, a guiding spring-force is added between the proton and a (moving) target point using a collective variable, $\xi$, equal to the height difference between the proton and an atom in the membrane

[TABLE]

with spring constant $\kappa=10$ eV/Å, where $v_{\xi}=5$ Å/ps is the target speed, chosen such that during the simulation time the proton will pass through the membrane. The work done then gives an estimate of the barrier and is computed as

[TABLE]

These simulations were done using the PLUMED package 20, integrated with CP2K.

3 Static energy barriers

In the simplest calculations 1, 4, 5, 6, the proton passes through the center of a hexagon with the lattice remaining fixed and the proton permeation path being a straight line perpendicular to the material. The energy barriers found in this way are about 0.9 or 1.5 eV for a proton passing through h-BN or graphene respectively as shown in Fig. 2a. These values are higher than the experimentally observed barriers (0.5 and 1.0 eV), but appear to be in reasonable agreement. Tensile strain may further reduce the calculated barriers slightly but, as shown in Fig. 2a, even 2% of tensile strain reduces the calculated barriers by only about 0.1 eV. We verified that the perpendicular cell size is sufficient, increasing it from 15 to 30 Å, the barrier changes less than 2 meV.

The interpretation of the proton by itself in DFT calculations is however not obvious. This is because, in the presence of a highly polarizable membrane such as graphene or hBN, an electron will be taken from the membrane leaving a charged membrane with neutral hydrogen in the ground state as illustrated in Fig. 2b and described in more detail in the supplementary information (Fig. S1, S2 and Table S1). As a result, at larger proton-membrane separation distances (above $\sim$1 Å), the energy curves depend strongly on the model size as shown in the right panel of Fig. 2a but flatten out with increasing sample size.

If we indeed assume this process to be slow, relaxation of the nuclear coordinates should be possible. When such a relaxation is done using a quasi-Newton method 21, the initial state is found to be unstable and the proton will chemisorb on the surface. If we then compute the barrier by means of the nudged elastic band (NEB) method with one climbing image (CI-NEB 18) the barrier (for h-BN) increases from 0.5 to 3.0 eV as shown in Fig. 3, a value similar to the barrier for permeation of a hydrogen atom 22. Thus, contrary to the intuition that optimization of the path should lead to a lower barrier, the barrier increases as a result of relaxation because the initial state is in fact unstable.

4 Effect of the solvent

In order to consider the role of the solvent we performed room-temperature Born-Oppenheimer ab initio molecular dynamics simulations of the membrane plus the explicit solvent with one additional solved proton. An example of the proton permeation process in the presence of water for h-BN is shown in Fig. 4. In this case we consider proton-membrane separation distances larger than 1 Å because the proton is stabilized by the presence of the water as demonstrated by the Mulliken population analysis 23 in Table S2.

After equilibration, the hydrogen nucleus closest to the membrane, together with the nearest atom in the lattice (B or N) are chosen for the SMD forcing. The SMD spring force is applied to both these atoms, forcing a change in z-coordinate difference (perpendicular to the membrane) from its initial value to a value of opposite sign and large enough that the SMD work is expected to flatten out. We note that, because the excess charge is spread over different nuclei in the liquid (see Fig. S3 and Table S2), the identification of a single proton is not possible and our selected H nucleus is in fact close to neutral initially. Rather than passing directly through the membrane, we observe the selected atom to move sideways away from its initial position. Simultaneously, a curvature is induced in the membrane (Fig. 4b) which brings the atoms closer in height to follow the (moving) SMD equilibrium spring position. This shows that it is favourable to bend the substrate rather than to directly approach the surface. The resulting work done is shown in Fig. 4(f). The barrier estimated from these curves contains two separate effects. The first involves the removal of the proton from the liquid and the second involves passing through the membrane. These steered dynamics simulations are equivalent to the Tomlinson model (see e.g. Ref. 24 and references therein) used in the field of friction at the atomic scale. The top of the barrier relevant to the experimental setup is therefore defined by the moment of slipping, whereas the starting configuration is given by the moment of OH bond breaking. As can be seen from Fig.4(f) and in more detail in Fig. S4, these barriers are then of the order of 3–5 eV, in qualitative agreement the nudged elastic band calculations of Ref. 11. After the reformation of an OH bond on the opposite side of the membrane the SMD work done indeed tends to flatten out as expected. Compared to the static calculations, rather than decreasing the barrier, barriers are thus increased because of the stabilizing effect the water has on the proton, lowering the energy of the initial state. We thus conclude that the inclusion of water in our model (one proton for 87 water molecules, corresponding to a pH close to zero) cannot explain the observed discrepancy between theory and experiment.

5 Curvature

We also considered the effect of curvature. For this we return to the case of static permeation, i.e. without optimization of the path. We consider two different cases of curvature, namely (i) a (10,0)–carbon nanotube (CNT) compared to graphene, and (ii) a buckled h-BN surface, where all B (N) atoms are displaced up (down) by 0.1 Å. The latter configuration could for example be the result of the ionic nature of h-BN causing B and N atoms to respond oppositely to the applied electric field. For (i) we use a model consisting of 360C atoms (9x10 orthorhombic unit cells) and recompute the energy curve with the same amount of atoms for graphene, while for (ii) we keep the usual 6x12 orthorhombic cells. The resulting barriers are shown in Fig. 5. While the barriers indeed decrease when considering permeation from the side favourable by curvature, the net effect is small compared to the total barrier, $\sim$0.2 eV in the case of a buckled h-BN sheet and less than 0.1 eV in the case of the CNT. These modifications therefore remain insufficient to explain the discrepancies between theory and experiment despite the relatively high curvature induced in these model system. Moreover, if it is indeed the case that the process is slow, we should rather expect it to chemisorb as discussed before and therefore only increase the barrier further.

6 Hydrogenation

In Ref. 22 we considered in detail the hydrogenation of h-BN. We also considered the permeation of a proton through these hydrogenated surfaces. However, as in the case of static permeation in vacuum, we find the initial state with a proton 3 Å above (measured from the center of mass) a group of chemisorbed hydrogen atoms to be unstable. In this case however, the intermediate state is not a chemisorbed one, but rather the formation of molecular hydrogen (H2) as shown in Fig. 6. We find similar formations of molecular hydrogen in other structures considered in Ref. 22. As a result, these intermediate stable states again lead to a increase in barrier when optimization of the path is considered.

7 Conclusion

We considered several effects on calculated barrier heights for proton permeation through h-BN and graphene. Despite our efforts to include more features of the environment we systematically find an increase of the permeation energy barriers rather than the desired decrease needed to match experimental observations. Several complications are discussed such as the high polarizability of the membranes and the existence of stable intermediates. These naturally arising complications have not been considered in previous work and we therefore found it necessary to comment explicitly on them in this work, despite not providing a definitive explanation for the proton permeation barrier heights observed experimentally. We therefore believe this problem may lie outside the scope of our DFT-PBE approach if the experimental membranes are indeed assumed to be defect-free and the process is slow compared to atomistic timescales.

The research leading to these results received funding from the Foundation for Fundamental Research on Matter (FOM), part of the Netherlands Organisation for Scientific Research (NWO). The work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 696656 – GrapheneCore1. We thank A.K. Geim, W. Brandeburgo and E.J. Meijer for useful discussions.