Potential energy barrier for proton transfer in compressed benzoic acid

Benzoic acid (BA) is a model system for studying proton transfer (PT) reactions. The properties of solid BA subject to high pressure (exceeding 1 kbar = 0.1 GPa) are of particular interest due to the possibility of compression-tuning of the PT barrier. Here we present simulations aimed at evaluating the value of this barrier in solid BA in the 1 atm – 15 GPa pressure range. We find that pressure-induced shortening of O⋯O contacts within the BA dimers leads to a decrease in the PT barrier, and subsequent symmetrization of the hydrogen bond. However, this effect is obtained only after taking into account zero-point energy (ZPE) differences between BA tautomers and the transition state. The obtained results shed light on previous experiments on compressed benzoic acid, and indicate that a common scaling behavior with respect to the O⋯O distance might be applicable for hydrogen-bond symmetrization in both organic and inorganic systems.


Introduction
Proton transfer reactions lie at the heart of many chemical transformations, such as acid-base reactions, catalytic processes or hydrogen transport in the Earth's mantle. [1][2][3] Studying proton dynamics not only gives insight into these transformations, but also leads to a deeper understanding of the quantum effects that are strongly manifested during the PT process. [4][5][6][7][8][9] Benzoic acid, with its centrosymmetric dimers bound by strong hydrogen bonds, is considered a model system for studying proton transfer processes. 10,11 This molecule was studied with a multitude of experimental techniques, including X-ray and neutron scattering, [12][13][14][15] optical spectroscopy, [16][17][18][19][20] and nuclear magnetic resonance (NMR), [21][22][23][24][25] as well as ab initio modelling. [26][27][28][29][30][31][32][33][34][35][36] High pressure (exceeding 0.1 GPa ¼ 1 kbar) can greatly alter proton dynamics in solids. 37 Hence, experiments were conducted for solid benzoic acid at pressures reaching 18 GPa. [38][39][40][41][42][43] However, the challenging nature of these experiments led to ambiguity in the interpretation of the obtained data. In particular, some studies gave conicting values of phase transition pressures, 40,42 while others reported no phase change in compressed BA. 41,43 Furthermore, structural data obtained by Kang and co-workers suggested that the O/O distance within the BA dimers remains constant upon compression to 18 GPa, 43 while single-crystal X-ray measurements conducted by Cai and Katrusiak indicated drastic reduction of this contact with pressure. 41 It was also suggested that at 17 GPa symmetrization of hydrogen bonds occurs in solid BA, 43 but this conclusions was mostly based on Density Functional Theory (DFT) calculations utilizing the local-density approximation (LDA) which is an inadequate method for modelling hydrogen-bonded systems. At the same time theoretical studies on benzoic acid at high pressures are scarce. 44,45 With the aim of elucidating the problems mentioned above we present a computational study on how high pressure inuences the properties of solid benzoic acid. We model the changes in geometry of BA crystals upon compression up to 15 GPa with the use of the recently proposed SCAN functional, 46 which was shown to correctly reproduce the properties of liquid water, an archetypical hydrogen-bonded system. 47 We also asses the potential energy surface (PES) barrier towards the concerted double proton transfer within the benzoic acid dimers taking into account the differences in the ZPE energies between the initial and transition state. The obtained data sheds light on the experimental results and lead to a better understanding of the proton-transfer process in compressed solids.

Computational tools
DFT calculations utilizing periodic boundary conditions were conducted in VASP 6.2.1 (ref. 48 and 49) for the ambient pressure structure of benzoic acid of P2 1 /c symmetry containing four molecules (two dimers) in the unit cell. 15 A plane waves basis set with a 800 eV cut-off was used for the description of valence electrons, while the 1s electrons of C and O were modelled with a projector-augmented-wave potential. Full geometry optimization (only symmetry constraints were applied) of the lattice vectors and fractional atomic coordinates of crystals containing  BA tautomers was performed with the use of the regularized version of the SCAN functional (r 2 SCAN). 50 The convergence criterion for the electronic minimization was 10 À8 eV per cell. Sampling of the Brillouin zone was done through a Monkhorst-Pack mesh with a 2p Â 0.05Å À1 spacing. For an optimization at a given pressure its value was added to the diagonals of the stress tensor and the p$V term was added to the calculated electronic energy yielding the enthalpy of the system (not corrected for zero-point energy motion, ZPE).
The ZPE correction was evaluated from G-point vibrational frequencies obtained from the optimized structures with the use of r 2 SCAN and the nite-displacement method (0.005Å displacement), as implemented in VASP. These vibrations include both intermolecular modes (translations and librations of BA molecules), as well as intramolecular modes, including O-H stretching modes. For obtaining proton-transfer barriers single-point energy calculations on the r 2 SCAN-optimized structures were done with the use of the hybrid HSE06 functional. 51 These calculations utilized the same computational parameters as used in the r 2 SCAN calculations.

Results
At ambient and high pressure benzoic acid crystallizes in the P2 1 /c space group with two equivalent dimers present in the unit cell. 15 The two tautomeric forms of these centrosymmetric dimers (A and B) are shown in Fig. 1. They can convert to each other through a transfer of two protons. Due to crystal packing effects A and B tautomers are not energetically equivalent and at ambient conditions crystals containing the A form have an energy lower by 0.2-0.8 kJ mol À1 (throughout the text the energies are referenced to one dimer), as determined by various experimental techniques. 25 This is one of the smallest energy difference between tautomeric forms of solidcarboxylic acid dimerse.g. nearly 8 kJ mol À1 is found for ibuprofen. 52 It can be expected that the enthalpy difference between crystals containing A and B tautomers, hereinaer referred to as asymmetry enthalpy, will increase upon compression. To verify this hypothesis geometry optimization, utilizing the r 2 SCAN functional, was performed in the 1 atm -15 GPa pressure range for solids containing in their unit cells solely A or B dimers, and subsequently the enthalpy of the two crystals was compared. In agreement with experiment we nd that at ambient conditions the energy of the A tautomer is lower by 0.6 kJ mol À1 (0.4 kJ mol À1 aer including ZPE differences between the tautomers) than that of the B form. As can be seen in Fig. 2(a) compression up to 10 GPa increases the enthalpy difference between the two tautomers to 2.6 kJ mol À1 (1.6 kJ mol À1 aer including the ZPE correction). Above this pressure the ZPEcorrected asymmetry enthalpy starts to decrease due to the increased differences in the zero-point energy between the two tautomers. The lower ZPE of the B form stems from a more elongated O-H bond (1.04Å and 1.05Å for A and B, respectively at 10 GPa), and consequently lower O-H stretching frequencies (2309/2173 cm À1 for A g modes for A/B). The lengthening of the O-H bond in B as compared to A can be linked with the destabilizing effect of the interdimer O/H-C contact of the hydroxyl atom which is shorter in this form (see Fig. 1).
The calculations not only correctly reproduce the enthalpy difference between A and B tautomers, but also yield geometry parameters in good agreement with experiment. As shown in Table S1 in the ESI, † differences in lattice parameters and bond lengths do not exceed 2.2% when comparing with both ambient pressure, 15 and high pressure data, 39,41 with the exception of the X-ray data of ref. 43 for which much larger discrepancies are found. However, the structural parameters reported in that study yield unphysically constant O/O distances at high pressure, as shown by black stars in Fig. 2(b), hinting at a possible error in the experimental structure determination. We note that inclusion of dispersion corrections, in the form of the D3 correction, 54,55 does not lead to better agreement with experimental values.
Upon compression the O/O contact is predicted to shorten in both forms of BA, as seen in Fig. 2(b), although this decrease is less abrupt than in the experimental data. Shorter O/O distance lead to strengthening of the hydrogen bond, and consequently the decrease in the difference between the O-H bond and H/O distancesee Fig. 2(c). Nevertheless, even at 15 GPa the system is far from hydrogen-bond symmetrization, in contrast to the results of previous LDA calculations. 43 Pressure-induced strengthening of the hydrogen bond should lead to the decrease of the barrier for the double proton transfer within the benzoic acid dimers. In order to evaluate this effect a structure being a linear combination of A and B forms was constructed at each pressure point. The coefficients of this combination, both close to 0.5, where ne-tuned in such a way that the transferring protons in the resulting structure were position half-way between the O atoms. This structure was be used as a proxy for the transition state along the linear path (LP) trajectory for the concerted double proton transfer. This trajectory is characterized by minimal displacement of heavy The barrier along LP is larger than along the minimum energy path (MEP). However, it is postulated that the proton transfer trajectory follows a path of minimum action with a barrier height lying in between that of LP and MEP. 56 Therefore, the LP barrier can be treated as the upper bound, while the MEP barrier as a lower bound for the true PT barrier. In order to obtain information on the MEP barrier, a constrained geometry optimization of the LP transition state was performed at each pressure point. In the optimization, done with the use of the r 2 SCAN functional, the fractional positions of the acidic hydrogen atoms were xed while positions of other atoms, as well as the lattice parameters, were allowed to vary. Constraining the position of the transferring protons allowed for relaxing the structure without it falling back to one of the tautomers. The structure obtained aer this procedure was used as a proxy for the transition state along the MEP. As can be seen in Fig. 2(b), this structure is characterized by O/O contacts below 2.41Åconsiderably shorter than those found in the A and B tautomers, and much smaller that O/O distances observed so far for carboxylic acid dimers at ambient conditions. 57 The predicted shortening of oxygen distance in the MEP transition state of solid BA mirrors the result of DFT calculations performed for isolated dimers, 30,32,33 and highlights the fact that concerted proton transfer along the MEP path is connected with considerable movement of heavy atoms. In fact, comparison of the tautomeric structures and the MEP transition state indicated that the shortening of O/O contacts is a result of the movement of whole BA molecules.
By constructing the proxies of the LP and MEP transition state structures we are able to estimate the value of the barrier for a concerted transfer of both protons within the BA dimers, which is the dominating mechanism of proton transfer in carboxylic acid dimers. 58 The barrier was estimated at each pressure point by subtracting the enthalpy of the A tautomer from the enthalpy of the given transition state (LP or MEP)in analogy with the methodology used in previous studies on BA dimers. 30,32,33 The correct description of PT barriers by DFT methods requires the use of hybrid functionals which incorporate a portion of Hartree-Fock exchange. 59,60 However, the large computational cost of these methods precludes its use for geometry optimization and phonon frequency calculations. Therefore the values of the barriers were extracted from enthalpy differences obtained from single-point calculations (no geometry optimization) performed with the use of the hybrid HSE06 functional. 51 ZPE corrections were evaluated at the r 2 SCAN level. We note that the barriers extracted directly from r 2 SCAN calculations are substantially lower (by about 30%) than those obtained with HSE06.
The PT barriers obtained for solid BA with the methodology described above are presented in Fig. 3 as a function of pressure and the O/O distance in the A tautomer. The MEP barrier at 1 atm (4.5 kcal mol À1 ) is smaller than theoretical values previously reported for isolated dimers of this compound (6.5-7.4 kcal mol À1 ). 6,29,30,32,33 Inspection of Fig. 3(d) reveals that this can be traced back to the shorter O/O distance in the r 2 SCANoptimized structure, as the MEP barrier decreases upon compression due to the shortening of the O/O distance.
At ambient pressure the LP barrier (10.1 kcal mol À1 ) is substantially larger than MEP barrier, but due to a more pronounced dependence on the O/O distance it decreases much faster with increasing pressure, as shown in Fig. 3(a). Nevertheless, even at 15 GPa the barriers for both the LP (2.6 kcal mol À1 ) and MEP (1.5 kcal mol À1 ) trajectories are still substantially larger than thermal energy at room temperature (0.6 kcal mol À1 ). However, inclusion of ZPE corrections (calculated with the r 2 SCAN functional) leads to a substantial lowering of the barriers at all pressures (Fig. 3). This is due to the fact that the high-energy symmetric O-H stretching modes of A g symmetry have imaginary frequencies in the transition states and therefore do not contribute to the zero-point energy  43 of LP and MEP. The ZPE-corrected value of the MEP barrier is only 0.7 kcal mol À1 at 1 atm, and is predicted to vanish upon compression to 0.8 GPa.
Above this pressure, which corresponds to an O/O distance in the A tautomer of 2.55Å, benzoic acid crystals with a symmetric hydrogen bond have a lower ZPE-corrected enthalpy than either the A or B tautomers, which is reected by a negative value of the barrier in Fig. 3. For the LP trajectory the ZPEcorrected barrier (4.5 kcal mol À1 at 1 atm) becomes zero at 9.7 GPa which corresponds to an O/O distance in A of 2.46Å. We note that the range of O/O distances spanned by the values at which the MEP and LP barriers vanish are close to the 2.5Å value recently proposed as the borderline for ZPE overcoming the PT barrier height in hydrogen-bonded systems. 61 Moreover the distance at which the LP barrier becomes null (2.46Å) is close to the critical O/O distance (2.44Å) at which hydrogen-bond symmetrization is proposed for a number of inorganic systems. 62 This hints that the energetics of the PT process may be similar regardless of the environment surrounding the hydrogen bond, although more studies are needed to conrm this.
The current results indicate that pressure-induced symmetrization of the hydrogen bond should occur in solid BA, but retrieving this effect requires taking into account the quantum nature of the system. Calculations suggest that this process should occur between 0.8 and 9.7 GPa. Taking into account the 2.2% underestimation of the length of O/O contacts in both tautomers, leading to the underestimation of PT barriers, the experimental pressure at which the MEP and LP barriers vanish can be estimated as closer to 6 and 20 GPa, respectively. Nonetheless, even at lower pressure the barrier for the PT transition, which as mentioned earlier should lie in between the MEP and LP barriers, should become comparable to thermal energy. As a result, structural data obtained from experiments yielding time-averaged information (e.g. X-ray diffraction) should contain a considerable admixture from the PT transition state, particularly in the case of room-temperature experiments. This, together with the fact that the MEP transition state is predicted to exhibit very short O/O contacts, might explain the dramatic decrease of O/O distances in the room-temperature measurements. 41

Conclusions
Simulations presented in this work show how pressure can be used to tune the barrier of the proton transfer reaction in solid benzoic acid. It is shown for the rst time that compression induces a barrierless transfer of proton, and that the vanishing of the barrier is strongly connected with differences in the ZPE energy. This indicates that quantum effects retain their important role in PT processes even at large compression. The obtained results shed light on previous experiments on compressed benzoic acid, in particular the drastic reduction of the O/O distance with pressure. The current ndings are also of importance for other high-pressure studies on strongly hydrogen bonded systems, 62 most notably water. [63][64][65] In particular, it seems that the recently proposed critical O/O distance for hydrogen bond symmetrization (ref. 62) might be common also to organic compounds. Benzoic acid is of special interest in this context, as this distance can obtained for this compound at relatively low pressure (p < 20 GPa). Therefore new experimental investigation aimed at establishing the proton dynamics in this compound are required.

Conflicts of interest
There are no conicts to declare.