In Silico Demonstration of Fast Anhydrous Proton Conduction on Graphanol

Development of new materials capable of conducting protons in the absence of water is crucial for improving the performance, reducing the cost, and extending the operating conditions for proton exchange membrane fuel cells. We present detailed atomistic simulations showing that graphanol (hydroxylated graphane) will conduct protons anhydrously with very low diffusion barriers. We developed a deep learning potential (DP) for graphanol that has near-density functional theory accuracy but requires a very small fraction of the computational cost. We used our DP to calculate proton self-diffusion coefficients as a function of temperature, to estimate the overall barrier to proton diffusion, and to characterize the impact of thermal fluctuations as a function of system size. We propose and test a detailed mechanism for proton conduction on the surface of graphanol. We show that protons can rapidly hop along Grotthuss chains containing several hydroxyl groups aligned such that hydrogen bonds allow for conduction of protons forward and backward along the chain without hydroxyl group rotation. Long-range proton transport only takes place as new Grotthuss chains are formed by rotation of one or more hydroxyl groups in the chain. Thus, the overall diffusion barrier consists of a convolution of the intrinsic proton hopping barrier and the intrinsic hydroxyl rotation barrier. Our results provide a set of design rules for developing new anhydrous proton conducting membranes with even lower diffusion barriers.


DFT Method Details
Density functional theory (DFT) data were generated using the Vienna ab initio simulation package (VASP). We used the projected augmented-wave (PAW) method 1 to describe electron-ion interactions and set a plane-wave cutoff of 520 eV. The generalized gradient approximation (GGA) exchange-correlation functional of Perdew-Burke-Ernzerhof (PBE) 2,3 was used. The Brillouin zone sampling was performed using a Monkhorst-Pack k-point grid size of 3 × 3 × 1 for the 24C cells. We determined this to be the optimal grid size as shown below. The ionic positions were relaxed until the forces were lower than the tolerance of 10 −3 eV/Å for optimization calculations. We used a vacuum spacing of 20Å above the surface of graphanol to mitigate interactions between layers of the material under periodic boundary conditions. A u24C cell was transformed into a c24C cell by the introduction of an additional proton to the hydroxylated surface, giving a system with a +1e charge. Calculations on charged systems under periodic boundary conditions are commonly handled with a homogeneous neutralizing background jellium charge. This background charge was shown to have a very small influence on proton dynamics. 4 Training data were generated using these DFT settings for the u24C, c24C and u96C systems. We employed both DFT-MD and single-point VASP calculations generate a suf-S2 ficiently varied training dataset. DFT-MD simulations in the N V T ensemble were carried out using a Nosé-Hoover thermostat. 5 The frequency of temperature oscillations were set to 40 time steps (SMASS=0).

k-point Mesh Optimization
We determined the optimal number of k-points required to converge the total energy of graphanol by running multiple single-point VASP calculations. We considered two sets of kpoints, only one of which included the Γ point. Each set comprised ten different calculations in which the number of k-points was progressively increased by varying the VASP KSPAC-ING tag from 0.1 to 1.0 in increments of 0.1. KSPACING is defined as the smallest allowed spacing between k points and takes the units ofÅ −1 . The results are shown in Figure S2.
We chose a KSPACING of 0.4 with the gamma point not included to be the best option. This is not necessarily the optimal case, because the energy does not vary monotonically with the number of k-points. However, a spacing 0.4 seems to be accurate enough compared to the fully k-point converged calculation. This spacing corresponds to a Monkhorst-Pack k-point grid size of 3 × 3 × 1.

Center of Excess Charge (CEC)
Proton transport (PT) was evaluated by tracking the center of excess charge (CEC) in charged graphanol. The CEC of a configuration is the position of the oxygen atom that contains a proton, thus carrying an excess charge. The CEC is defined as the oxygen atom that has two H nearest neighbors (cutoff radius of 1.2Å). It is possible, because of fluctuations in the O−H bond distances, to encounter situations where there are multiple O atoms having two nearest H neighbors, even when there is only 1 proton in the system. When this arises and when we have no history of where the CEC was at the previous time step (such as when starting the CEC calculation after equilibration), then we assign the CEC by comparing the Figure S2: Convergence of energy with the spacing between k-points. KSPACING was varied from 0.2 to 1.0 with (orange circles) and without (black stars) Γ point. were computed from multiple independent production runs using Einstein's relation, where t is the time, d = 2 is the dimensionality of the system since the diffusion is 2-D, r CEC is the position of the O atom that is labeled as the CEC at time t. It is important to note that these graphanol systems have only a single added charge and thus will only have a single CEC position at any given time. Running multiple production simulations was required to achieve improved statistics in our results. The uncertainties in D s are reported as twice the standard deviation of the independent values. The calculated D s at different temperatures were used to measure the PT diffusion energy barrier (E A ) using the Arrhenius equation, Other

Projected Phonon Density of States
Projected phonon density of states for the primitive cell of 2-D single-sided graphanol (u2C) were computed with DP and compared with VASP, as shown in Figure S4.

Mean Square Displacement Plots
We have computed the mean square displacements of the CEC, which is the ensemble averaged term in Eq. S1, M SD(t) = ⟨ |r CEC (t) − r CEC (0)| 2 ⟩. We used multiple time origins S6 in these calculations. We see from Eq. S1 that a plot of M SD(t)/t must approximately be a straight line with zero slope at long times if the system obeys the Einstein diffusivity relation. If M SD(t)/t consistently decreases with time then the system is sub-diffusive (as is seen for systems that obey single-file diffusion) and if it increases with time the system is super-diffusive (e.g., ballistic diffusion). Plots of the M SD(t)/(4t) are provided in Figure S5.
The factor of 4 is 2d, where d = 2 is the dimesionality of the problem, so that M SD(t)/(4t) gives  Table 1.

Proton Hopping Energy
The proton hopping energy, E hop , is defined as the intrinsic activation energy for a single proton to hop from one O atom to a neighboring O atom. E hop is expected to be significantly smaller than the overall PT barrier E A . We used our final DP to estimate the value of E hop .
Two images from a c24C DP-MD simulation were taken and geometrically relaxed using the DP. These two images constitute the start and end state of a proton hopping event, as shown in Figure S7. We then performed a geometric interpolation to identify 11 intermediate states between the start and end images. Single-point energies at T = 0 K were calculated for start, end and these 11 intermediate images. The energy required for a proton to move S9 from the start state to the maximum energy transition state is 12.25 meV. The barrier for the reverse transition is 13.64 meV. The average of the two barriers is 12.94 meV, which we take as an estimate of E hop . This energy is about one order of magnitude smaller than the overall PT barrier, E A . We note that our calculated value of E hop is similar to the reported barrier for proton hopping in water, computed from potential of mean force. 8 Figure S7: Estimate of the energy barrier for a single proton hopping event (E hop ) using the DP. All energies are estimated by performing single point calculations at T = 0 K. Visualization of the proton hopping event in shown at the bottom. The location of the proton changes from the "Start" state to the "End" state. A blue dashed box is added to mark the position where the proton hopping is occurring.

Hydroxyl Group Rotation Energy
Calculating the activation energy for hydroxyl group rotation in graphanol, E rot , is challenging because the rotation of each hydroxyl group is coupled to the rotation of other groups.
We therefore calculated rotational diffusion coefficients, D r , 9 as a function of temperature as a way to estimate E rot . Angular mean-squared displacements (AMSD) ⟨ϕ 2 ⟩ were first computed by tracking the rotation of dihedrals of all the OH groups from DP-MD simulations.

S10
The C-C-O-H chain forms the dihedral ϕ i for the OH group i, as shown in Figure S8a. The rotational diffusion coefficient D r was calculated using the Einstein relation, where (∆ϕ) 2 is the time averaged AMSD and t is the simulation time. We used the DP to run several 2.5 ns N V T -MD simulations at three different temperatures T = 800, 900 and 1000 K. We ran 20 independent runs for each temperature to obtain improved statistics. We then estimated E rot by fitting the D r to the Arrhenius equation, The fitting gives E rot = 293± 45 meV, which is higher than our estimates for E A and significantly higher than our estimate for E hop .  Table S1: Membrane types, labels used in Figure 3c, reinforcements material used, and the activation energy for PT in meV.

Data Generation, DP Training and Validation
A schematic of the overall data generation, DP training and validation is shown in Figure 6 of the main text. The initial data set was generated to target essential aspects of the which were validated with DFT. This procedure determines the quality of the "best" DP for that given active learning iteration.
Another key part to the active learning process is the exploration step (MD exploration in Figure 6), which involves sampling of the configuration space and transferring these configurations to the labeling step. The error indicator from the exploration step was used as a metric to test the quality of DPs for that iteration. A successful iteration was one where the "best" DP passed the phonon test, bond energy tests, and achieved a low error (ϵ) from times. We ensured that no two trace-back sequences overlapped at a given time, to avoiding over-counting. The time-distribution histogram obtained from this analysis is plotted in Figure 5, which contains all trace-back sequences (of various lengths). A log transform of this plot ( Figure S9a). exhibits a normal distribution. This was verified using the QQ plot and the D'Agostino and Pearson's statistical tests. 21,22 According to our analysis, the mean proton trace-back time is 0.24 ps. A chain-length resolved histogram is plotted in Figure S9bd. Trace-back sequences involving three unique O atoms dominate the distribution, whereas sequences of increasing length appear less and less frequently. This finding suggests that it is unlikely that a GC much longer than 6 survives for sufficiently long times for trace-back to occur, or that GCs much longer than 6 are unlikely to occur. Our analysis cannot distinguish between these two possibilities.