Methyl Dynamics Flattens Barrier to Proton Transfer in Crystalline

Feb 2, 2012 - National Institute of Chemistry, Hajdrihova 19, SI-1000 Ljubljana, Slovenia ... University of Sydney, Sydney, New South Wales 2006, Aust...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCA

Methyl Dynamics Flattens Barrier to Proton Transfer in Crystalline Tetraacetylethane Gordon J. Kearley,*,† Jernej Stare,‡ Ramzi Kutteh,†,§ Luke L. Daemen,∥ Monika A. Hartl,∥ and Juergen Eckert⊥ †

The Bragg Institute, Australian Nuclear Science and Technology Organisation, Lucas Heights, 2234 New South Wales, Australia National Institute of Chemistry, Hajdrihova 19, SI-1000 Ljubljana, Slovenia § School of Physics, University of Sydney, Sydney, New South Wales 2006, Australia ∥ Los Alamos Neutron Science Centre LANSCE, Los Alamos, New Mexico 87545, United States ⊥ Department of Chemistry, University of South Florida, 4202 East Fowler Avenue, Tampa, Florida 33620, United States ‡

ABSTRACT: We analyze the interplay between proton transfer in the hydrogen-bond bridge, O···H···O, and lattice dynamics in the model system tetraacetylethane (TAE) (CH3CO)2CHCH(COCH3)2 using density functional theory. Lattice dynamics calculations and molecular dynamics simulations are validated against neutron scattering data. Hindrance to the cooperative reorientation of neighboring methyl groups at low temperatures gives a preferred O atom for the bridging proton. The amplitude of methyl torsions becomes larger with increasing temperature, so that the free-energy minimum for the proton becomes flat over 0.2 Å. For the isolated molecule, however, we show an almost temperature-independent symmetric double-well potential persists. This difference arises from the much higher barriers to methyl torsion in the crystal that make the region of torsional phase space that is most crucial for symmetrization poorly accessible. Consequently, the proton-transfer potential remains asymmetric though flat at the base, even at room temperature in the solid.

1. INTRODUCTION Short strong hydrogen bonds in which a proton bridges donor and acceptor atoms have been a continuing topic of research as model systems for proton transfer,1−7 with particular interest in their structure and dynamics. Here, we present an accurate description of the dynamics of proton transfer in one such system, namely, tetraacetylethane (TAE) (CH3CO)2CH CH(COCH3)2, by bringing together time-average structure and thermal parameters from neutron-diffraction crystallography,6 vibrational spectroscopy from inelastic neutron scattering (INS), and density functional theory (DFT) calculations. In addition to its importance noted above, we have chosen to study TAE, depicted in Figure 1, also because the potential-energy (PE) surface of its bridging proton has important components from intramolecular and intermolecular interactions with methyl groups. Consequently, using ab initio molecular dynamics (MD) simulations it should be straightforward to separate and analyze these interactions. The relevant features common to all of the aforementioned model systems are the vibrational modes of the proton, how these are coupled to other coordinates in the system, and how this coupling evolves with temperature. The general picture that emerges from past studies on these mainly nonsymmetric systems is that as the temperature increases the time-average position of the proton moves away from the donor toward a midposition and may even be transferred permanently (on time © 2012 American Chemical Society

Figure 1. Schematic illustration of the TAE molecule in the lowest energy conformation (isolated molecule).

average) to the acceptor.4 Similarly, crystallography has shown that in TAE, despite the symmetry of the molecule, there is a tendency for the bridging proton to be closer to one of the O atoms even at room temperature,6,7 where we henceforth denote the O atoms as O1 and O2, with the preferred O atom for the bridging proton as O1. On increasing the temperature Received: October 24, 2011 Revised: January 31, 2012 Published: February 2, 2012 2283

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

2. EXPERIMENTAL AND COMPUTATIONAL DETAILS The INS spectrum was collected using the Filter Difference Spectrometer (FDS)11 at the Lujan Centre, Los Alamos National Laboratory. Samples were contained in thin-walled aluminum containers, with temperature control being achieved using a standard closed cycle helium cryostat. Data correction and conversion to the energy scale were made using locally written algorithms. The five CIF files containing the neutron diffraction crystal structures6 of TAE at temperatures of 20, 40, 70, 110, and 298 K were obtained from the Cambridge Crystallographic Data Centre. The crystal structure of TAE has an orthorhombic crystal system with the Pbcn space group. The unit cell contains 112 atoms, where a TAE molecule has 28 atoms as shown in Figure 1. The lattice parameters at the above five temperatures are given by Piccoli et al.6 DFT calculations on the isolated TAE molecule were carried out using [DMol3]12,13 and CPMD v.3.13.14 The DMol3 calculations were performed with an all-electron treatment using the generalized gradient approximation (GGA) with the Perdew−Wang exchange-correlation functional (PW91).15 Calculations on the periodic TAE systems were performed using CPMD v.3.13 and the plane wave DFT code VASP.16−20 CPMD calculations of methyl torsion potentials of isolated and crystalline TAE were carried out with the BLYP functional21,22 coupled with a plane-wave basis set with a kinetic energy cutoff of 70 Ry; the core electrons were approximated by atomic pseudopotentials of Trouiller and Martins.23 While the simulation box of solid TAE was constructed according to the crystallographic data, a cubic simulation box of 10 Å was used for the isolated molecule and periodic boundary conditions were not employed. For the complex low-frequency dynamics we found that VASP gave results that were in better agreement with experiment than CPMD. For stepwise calculation of barriers there is little difference, but CPMD offers better control over geometry constraints. Three types of calculations on the periodic TAE systems were performed with VASP: geometry optimization, phonon frequency calculation, and MD simulation. In all VASP calculations we adopted the projector-augmented wave (PAW) potential24 and the Perdew−Burke−Ernzerhof25 (PBE) exchange-correlation functional. For accurate calculations a 4 × 4 × 3 Γ-centered Monkhorst-Pack26 k-point mesh and 400 eV energy cutoff were used. For the MD simulations, 14 k points were used for the simulation at 40 K for generating the INS spectrum. Only slight differences were found when 8 k points were used, so MD simulations at all other temperatures used 8 k points for computational efficiency. For all MD simulations the energy cutoff was reduced to 300 eV and a time step of 1 fs was used. 2.1. Geometry Optimization. Geometry optimization was performed on the 20, 40, 70, 110, and 298 K neutron diffraction structures of TAE using VASP with the cell volume and shape held constant at the experimental values but allowing the ions to relax to their minimum energy positions. The values of the cell constants are known from experiment with considerable accuracy, and holding these fixed overcomes the poor treatment of dispersion interactions in the DFT methods that typically leads to artificial cell expansion. We ensured that the full symmetry of the TAE crystal structure was recognized by VASP and maintained during optimization. The 20 K optimized structure formed the basis for the vibrational

further, the proton is found to migrate toward a position that is closer to the midpoint between the donor and the acceptor. It is generally accepted that in nonsymmetric systems neither lattice expansion nor formal (Donor···H) stretching vibrations plays a significant role in this proton behavior and transfer but that they are due instead to specific low-frequency phonons that modulate the geometry of the H-bonded system.2,3 In TAE however we will show that this general proton behavior and transfer result from a rather different mean field effect. The Hbond system in TAE is symmetric and sufficiently simple that we can be reasonably confident that this effect involves the methyl groups, making the analysis and understanding of proton-transfer dynamics more straightforward in TAE, as mentioned before. This situation is analogous to an earlier study of acetylacetone8 (AC) in which the barrier to methyl rotation is rather low so that methyl rotation is much more rapid than proton transfer. In the present case however the barrier to methyl rotation is considerably higher, giving the possibility of either assistance or hindrance of the protontransfer process. INS is particularly sensitive to H-atom motion of the Hbonding structure and of the methyl groups. Because it is straightforward to calculate the INS spectrum from either vibrational frequencies with their displacement vectors or atomic MD simulation trajectories, we are able to validate our DFT calculations not only against the measured crystallographic data but also against the measured INS spectrum. Essentially, we are validating calculated structure and dynamics directly against measured counterparts. After describing in section 2 the details of our experimental and computational approaches to studying the TAE system, such a validation is given in section 3. Following validation our task is then to understand the correlations between the movements of different parts of the simulated system and derive components of the PE surface of the bridging H atom in section 4. The current study differs from most reported work because it does not single out MD frames for which vibrational coordinates are then analyzed3,9,10 but analyzes the PE surfaces based on all MD frames and correlations between coordinates. The dynamics are determined by the time-dependent PE surface, which can be characterized in a meaningful way by use of MD simulation. Effective one-dimensional (1D) PEs are easily obtained by mapping the energy against some chosen coordinate (or combination of coordinates), but these may provide a poor estimate of the PE barrier because they miss dynamic effects, which are strong in the present case. Instead, we derive the free-energy (FE) surface from the time-average nuclear density distributions by assuming a Boltzmann distribution at the MD simulation temperature. The FE surface is closely related to the PE surface but has the advantage that it takes account of volume and entropy terms that are of particular importance here. By using calculations from fully periodic systems representing the crystal and from an isolated TAE molecule we are able to highlight the intermolecular interactions. By understanding these we can offer a consistent picture of how, on one hand, phonons may assist proton transfer and, on the other, stabilize a particular geometry leading to a polaron effect that tends to localize the proton at one end on the bridge, especially at lower temperatures. 2284

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

transformed to obtain the scattering function, S(q,ω). In-house software was used to ensure the appropriate cut was taken through this surface to coincide with the experimental INS data, and this was then convoluted with the instrumental resolution function to obtain the calculated INS spectrum, which is compared with the observed spectrum in Figure 3. This will be discussed further in section 3.3.

calculations described below (see Figure 2), while the remaining optimized structures were used as the starting points for the MD simulations.

Figure 2. Comparison of the observed INS spectrum (solid line) of TAE at 20 K with the results of DFT lattice dynamics calculation using the harmonic approximation (dashed line). (Inset) Low-frequency region where the methyl torsions arise.

2.2. Phonon and INS Calculations. We performed phonon calculations using the code PHONON27,28 in which a supercell is built from the optimized unit-cell structure and the nonzero force constants are determined by successive small displacements of the crystallographically inequivalent atoms. Forces in these displaced structures were calculated with VASP and used to generate the dynamical matrix in PHONON, which was diagonalized for a large number of random points in the reciprocal lattice to obtain the vibrational density of states. By taking account of the neutron scattering cross-section of the different atoms and the experimental conditions used for measuring the INS spectrum it is now straightforward to calculate the expected INS spectrum. This spectrum was compared with the measured spectrum (Figure 2) to assess the validity of the DFT calculations and is discussed further in section 3.2. 2.3. MD Simulations and INS Calculations. We performed microcanonical ensemble molecular dynamics simulations on the 40, 70, 110, and 298 K structures of TAE, optimized as described above. Starting with the TAE structure at temperature T, the following procedure was followed: (1) 4 ps of equilibration in the isokinetic ensemble (i.e., velocity scaling) with the temperature fixed to T; (2) 4 ps of equilibration in the microcanonical ensemble; and (3) at least 16 ps of production in the microcanonical ensemble with the thermodynamic temperature fluctuating about T. The dynamical structure factor, S(q,ω), was calculated from the atomic trajectories of the MD simulations using the code NMoldyn.29 First, the incoherent intermediate scattering function, F(q,t) was calculated on a regular grid in both momentum transfer q and time t as F (q , t ) =

Figure 3. Comparison of the observed INS spectrum (solid line) of TAE at 20 K with the spectrum derived from the DFT/MD calculation at 40 K. (Inset) Low-frequency region where the errors due to limitations on cell size and simulation time have the largest effect.

3. VALIDATION OF CALCULATIONS WITH EXPERIMENT Before drawing any conclusions based on our computational TAE results, we validate the TAE model and parameters by comparing with experiment the results of the geometry optimization, vibrational calculations, and molecular dynamics simulations. 3.1. Structure. Experimentally determined structures at 40, 70, 110, and 298 K were optimized, partly to provide a starting model for the MD simulations but also to ensure reasonable agreement between the experimental and the DFT energyminimized structures. This was achieved at all temperatures, and in particular, the root-mean-square difference between the starting 20 K neutron diffraction structure of TAE and the final DFT structure optimized as described above is only 0.038 Å. The same temperatures were chosen for the MD simulations so that the experimental unit-cell parameters could be used. 3.2. Vibrational Spectra. Structural calculations give information on the potential-energy minima where the atoms are situated at the lowest temperatures, but the main interest in the current work is dynamics, which probes the curvature of the potential-energy surface around these minima. Intermolecular interactions in the solid will lead to dispersion in the vibrational dynamics so that the measured INS spectrum will be the envelope of the vibrational density of states (VDOS) rather than the more familiar sharply resolved spectrum of the vibrations at the Brillouin zone center seen in IR and Raman spectoscopies. Inclusion of dispersion is crucial to the present study, and we calculate the INS spectrum from the VDOS (using PHONON) and compare this with the experimental INS spectrum. A comparison of the spectra shown in Figure 2 reveals that both the calculated vibrational frequencies and the amplitudes (closely related to intensities) provide an INS

∑ bi⟨exp[−iq·R i(0)]exp[iq·R i(t )]⟩

where R(0) is the position of the atom at t = 0, b is the total scattering cross-section (using the incoherent approximation), and the summation is over the atoms of interest in the calculation. In practice, b for the H-atom nucleus is so much larger than for any other nucleus in the system that only H atoms need be considered. Second, F(q,t) was then Fourier 2285

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

profile which resembles that measured. The agreement is good enough for us to make reasonable assignments for most of the observed peaks to the atomic displacements. The regions of disagreement between these profiles, particularly around 1000− 1200 cm−1, reflect the deficiencies of the harmonic approximation for the H-bond dynamics in TAE . Although we only compare the calculated vibrational density of states with the experimental spectrum, the PHONON calculations also provide the dispersion curves for each of the phonon branches. Although these should be used with caution in the spectral region alluded to above, they provide an estimate of the contribution of each atom to these branches and the participation of the H-bonding proton in intermolecular interactions that are evident as dispersion across the Brillouin zone. 3.3. Molecular Dynamics. The structural and vibrational calculations detailed above can be carried out at high accuracy because the total number of energy and force calculations is fairly modest. However, these calculations are limited to 0 K and the harmonic approximation, which actually works reasonably well except in the 1000−1200 cm−1 region. In order to overcome these limitations and explore the temperature regions in which proton transfer arises, we need to use MD simulations. The INS spectrum calculated from the MD simulation at the lowest temperature (40 K) that will reasonably populate the vibrational levels of interest is illustrated in Figure 3. Agreement between experimental and calculated spectra is improved compared with Figure 2, indicating the importance of anharmonicity. The remaining difference between the observed and the calculated spectra arises from our limited system (a single unit cell) and the neglect of quantum aspects of the proton dynamics, both of which would require considerably more computing resources that we have currently available. From the general agreement in Figure 3 we consider that our MD simulations reproduce the most essential elements of the dynamics will also realistically reproduce the temperature dependence of these dynamics.

Figure 4. VDOS of the bridging proton (dashed line) and methyl H atoms (solid line) calculated with PHONON. Three fundamental displacements of the proton are approximately along the O−H···O (O−H stretch, 1940 cm−1) perpendicular to that direction and in the C−O···O−C plane (COH in-plane angle bend, 1600 cm−1) and perpendicular to both of these directions (out-of-plane COH deformation, 1180 cm−1). (Inset) Low-frequency region dominated by methyl torsions with little presence of bridging proton motions.

made by visual inspection of the corresponding MD atomic displacements. The much weaker low-frequency features in the inset of Figure 4 correspond to other modes that cause relatively small displacements of the bridging H atom. These modes will be discussed (in section 4.2), particularly in relation to the dominant methyl torsions evident in the inset of Figure 4. Surprisingly high frequencies for the O−H stretching mode (here 1940 cm−1), unobserved in experimental spectra, have also been calculated in other H-bonded systems;2 these high values are a consequence of the harmonic approximation inherent in the lattice dynamics approach. The MD simulation is not subject to the harmonic approximation, and the O−H stretching mode deduced from these simulations will be discussed in section 4.2. We mention that a third approach to this problem, intermediate in accuracy, is through solution of the Schrödinger equation on a multidimensional PE surface.9 Figure 5 shows dispersion curves for the three vibrational modes of the bridging proton. A relatively large spread in the

4. RESULTS AND DISCUSSION In section 4.1 we discuss the fundamental vibrational modes of the bridging H atom, based on its VDOS and the dispersion curves of its modes, computed with the method of small displacements. In section 4.2 we address the question of possible correlation between methyl torsions and protontransfer dynamics in both crystalline TAE as well as the isolated TAE molecule. The most important free-energy profiles and surfaces for proton transfer are analyzed in section 4.3, and the interplay between methyl torsion and the symmetry of the H bond is described in section 4.4. Finally, in section 4.5 we look in more detail at the overall interaction between the H-bond structure and the methyl torsions and how the net effect of these decreases with increasing temperature. 4.1. Vibrational Characteristics of the H-Bonded Bridge. Referring to Figure 1, the bridging H atom has three degrees of freedom that are commonly used to classify hydrogen-bond vibrational dynamics: O−H stretch, COH inplane angle bend, and torsion about the C−O bond (or out-ofplane COH deformation). The vibrational frequencies of these degrees of freedom are 1940, 1600, and 1180 cm −1, respectively, as evidenced by the prominent features of the VDOS of the bridging H atom in Figure 4, calculated with the lattice dynamics program PHONON (section 2.2). Attributions of the vibrational frequencies to degrees of freedom are easily

Figure 5. Dispersion of the fundamental vibrations involving the bridging proton, computed with PHONON27 Note the comparatively large dispersion of the in-plane C−O−H bending mode at 1600 cm−1 compared with the much flatter dispersion of the other modes.

frequency of the COH in-plane deformation mode at 1600 cm−1 is observed, giving rise to the comparatively large width of the peak at 1600 cm−1 in the VDOS of Figure 4. This spread is unusually large for formally internal molecular vibrations at 2286

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

determined by very small atomic displacements cannot explore the H-transfer process, whereas the MD even at 70 K shows significant proton transfer, and the simulation can be used to analyze that process. In particular, from the trajectory of the MD simulation at 110 K many proton transfers occur, but the maximum torsional amplitude of the methyl groups is less than 30° on either side of the mean. Even in the MD simulation at 298 K, where the torsional amplitudes are larger and methyl reorientations do occur, there is still no significant correlation between proton transfer and methyl reorientation.. Fourier transforming the time evolution of the O···H···O antisymmetric stretching coordinate and the methyl torsional coordinate (φ1,φ2) from the MD simulation at 70 K gives their frequency spectra in Figure 7, showing methyl torsions between 100 and

such high frequency, reflecting strong intermolecular coupling of this mode and the importance of H-atom displacements perpendicular to the O···O vector of the H-bonding bridge. The many branches at around 1600 cm−1 in Figure 5 separate into two groups (in phase and out of phase for the intramolecular H-bond bridges), with a corresponding splitting in the VDOS at this energy (Figure 4). In contrast, the flat dispersion curves of the stretching and torsional modes show weak intermolecular coupling. 4.2. Are the Methyl Torsions Coupled to the ProtonTransfer Dynamics? The answer depends on whether the system in question is the free TAE molecule or crystalline TAE. In the free molecule the answer is yes, as illustrated by the 1-D proton-transfer PE curves in Figure 6, which were obtained by

Figure 7. Fourier transforms of the antisymmetric O−H stretch coordinate (blue) and methyl torsion coordinate (black) from the MD simulation of crystalline TAE at 70 K. There is little overlap between these two except around 200 cm−1, which corresponds to phonon assistance of the proton transfer.

Figure 6. BLYP/6-31G(d,p) relaxed proton potentials of a free TAE molecule, calculated at different torsional orientations of the methyl groups, φ1 and φ2 (see Figure 1). Color/curve correspondence is as follows: red, φ1,φ2 optimized; blue, φ1,φ2 fixed to 0°; magenta, φ1,φ2 fixed to 180°; green, φ1,φ2 fixed to 0° and 180°; cyan, φ1,φ2 fixed to the fully optimized values in the crystal field 18.5° and 172.0°, respectively.

200 cm−1, in agreement with the VDOS in Figure 4. Most of the O···H···O antisymmetric stretching magnitude is between 200 and 500 cm−1. There is some coupling between the two coordinates, as evidenced by their weak spectral amplitudes over the entire region, but the intersection is small. In contrast to the isolated TAE molecule system, there is no significant correlation between methyl-group torsions and proton transfer in the crystalline system. The barriers to methyl rotation are clearly much higher in the crystal, and this drastically reduces correlation to the proton-transfer coordinate. Finally, Figure 7 shows a broad feature of significant intensity for the antisymmetric O−H stretch coordinate centered around 1600 cm−1, in reasonable agreement with the expected frequency for this mode and in contrast to the unexpectedly high value (1940 cm−1) in Figure 4, which is due to the limitations of the harmonic approximation, as noted in section 4.1. 4.3. Proton-Transfer Dynamics. Analysis of the bridging H-atom MD trajectory for the transfer from one side of the O···O vector to the other is straightforward. However, many other degrees of freedom are involved in the proton transfer, and these are difficult to quantify in a simple way due to their complex temporal variations. A convenient means of including at least the resultant of these variations is by calculating the potential of mean force (PMF) from the time-average density distribution of the appropriate coordinates. The PMF is equivalent to the FE (free energy) that we can obtain from the MD simulation at temperature T for any coordinate r as FEr

mapping the energy as a function of proton position. These curves show clearly that the deepest symmetric double well (highest barrier) corresponds to proton transfer with the torsional angles φ1, φ2 of the methyl groups (see Figure 1) optimized at each step (red line). Constraining the methyl groups to be symmetrically related by having a C−H that is either cis (φ1, φ2 = 180°) or trans (φ1, φ2 = 0°) to the CO, also results in double-well potentials, albeit shallower. For the asymmetric cis/trans arrangement (φ1 = 0°, φ2 = 180°), the deeper minimum corresponds to the proton being closer to the trans conformation, which as we will show accounts for the crystallographic inequivalence of the donor and acceptor atoms. There is some correlation between the orientations of methylgroup pairs and the proton position in the isolated TAE molecule. Methyl groups on either side of a H-bond bridge tend to be staggered with respect to each other. The proton spends more time closer to the C−O staggered with respect to its methyl group. In the TAE crystal however the answer to the question in the title of this section is somewhat different. The VDOS of the bridging H atom and the methyl H atoms in the crystal, computed with PHONON27 and compared in Figure 4, show clearly that there is only slight displacement of the bridging H atom during the methyl librations between 100 and 200 cm−1. Clearly, a vibrational analysis based on force constants 2287

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

= −kT ln(Pr/Pmax), where Pr is the time-average population of r and Pmax the maximum value of Pr. The most condensed representation of the proton-transfer FE surface is a one-dimensional FE of the proton at distances from the average of the two O···H distances, with all calculated at each MD frame. This FE is illustrated in Figure 8 for MD of

Figure 8. FE from the occupation frequency of the bridging proton at distances from the average of the two r(O−H) lengths at each frame of the MD simulation. Note the larger barrier and greater distance for proton transfer in the isolated TAE molecule (black and inset). Right and left sides of the graph correspond to O1 and O2, respectively.

the crystalline TAE at four temperatures and the isolated TAE molecule at 120 K. For the crystalline TAE system there is at all temperatures studied a single PE minimum near 0.07 Å, on the O1 side of the O2···H···O1 bridge, with no symmetrically placed minimum on the O2 side of the bridge. There is a marked flattening of the FE with increasing temperature with an almost flat-bottomed minimum at 298 K, stretching over ∼0.1 Å on either side of the average position. For the isolated TAE molecule there are two almost symmetrically placed minima, 0.18 Å away from the average position and separated by a barrier of 0.3 kcal/mol (150 K). Although all degrees of freedom of the FE surface cannot be displayed simultaneously, we can gain insight by considering the most important dimensions in pairs. The main limitations of 1-D FE described above can be overcome by the FE function of the two O···H distances, as illustrated in Figure 9a and 9b for the MD simulation of crystalline TAE at 298 K (and 40 K) and the isolated TAE molecule at 180 K, respectively. The most obvious feature of these figures is the curvature that arises because the well tends to be deeper as the O1···H and O2···H distances become equal (at ∼1.25 Å) and their sum becomes a minimum. It is interesting to note that this midpoint (1.25 Å) is the same for the crystal and the isolated molecule, but for the latter this is a saddle point of height 0.3 kcal/mol, while for the crystal at 298 K this point is part of a long, wide trough. The width of this trough indicates that interactions that are not parallel to the O···O vector play an important role in the phonon assistance, this being related to the high dispersion of the COH in-plane deformation mode (Figure 5) discussed in section 4.1. For the crystal at 40 K (inset in Figure 9a) there is a single asymmetric minimum that forms as the polaron effect increases. By plotting the potential as a function of how far the proton is from the average of the two O···H distances (for each time step) against the O···O separation, we retrieve the classic

Figure 9. (a) Contours of FE surfaces (kcal/mol) determined from the occupation frequency at pairs of distance (of proton from O1 and O2) for the MD simulation of crystalline TAE at 298 K. Curvature reflects the overall trajectory of the proton between the two O atoms. (Inset) Surface for the MD simulation of the TAE crystal at 40 K, with axes covering the same range as the main axes. (b) Contours of the FE surface (kcal/mol) as for part a but from MD simulation of the isolated TAE molecule at 180 K.

picture of the hydrogen-bond potential (Figure 10a and 10b). This is more obvious for the isolated molecule (Figure 10b), where it can be seen that as the O···O distance increases from 2.4 Å (the saddle point) the two symmetrically placed minima separate. For the crystal at 298 K (Figure 10a) the situation is rather different, with a well that is 0.1 Å wide in the O···O coordinate (energy < 0.34 kcal/mol) and 0.4 Å along the approximate O···H coordinate. Finally, as mentioned at the start of this section, several other coordinates are involved in the proton-transfer process, and rather than examining their FE surfaces, we will simply give the most important correlation coefficients at 298 K. The PE surface of the proton has an important contribution from the C−O···H angle deformation, with a correlation coefficient of −0.30 between this angle and the corresponding O···H distance, reflecting the tendency for the proton to move around the O atom when the approach becomes close. This behavior also contributes to the curved shapes seen in Figures 9 and 10. Not surprisingly, there is a strong correlation, −0.56, between the O···H separation and the corresponding C−O bond length, which agrees well with the time-average results of the crystallographic study.6 This correlation is similar to that between the O···O and the O···H distances (0.51), from which it follows that the C−O and O···O coordinates are equally important in determining the net PE barrier for proton transfer. 4.4. Methyl Torsions and Symmetry of the Hydrogen Bond. To critically examine the conditions for a symmetric 2288

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

Figure 11. BLYP/PW torsional potentials of the methyl groups of TAE at the donor and acceptor sides, respectively, with the solid line corresponding to crystalline TAE and the dashed line to the isolated TAE molecule.

Figure 10. (a) Contours of the FE surface (kcal/mol) from the O···H and O···O distances from the MD simulation of the TAE crystal at 298 K. (b) Contours of the FE surface (kcal/mol) as for part a but from MD simulation of the isolated TAE molecule at 180 K.

significant quantum rotational tunneling of the methyl groups. Consequently, the H-bonded bridge in AC contributes to the rotational barrier of the methyl groups, but the rate of their rotation is much higher than that of proton transfer, and again dynamical coupling is absent. We can estimate the height of the methyl rotation barrier in TAE from the observed methyl libration features in the INS at 130 and 180 cm−1 (inset in Figure 2) and in conjunction with the Schrödinger equation. These estimated heights are 1.24 and 2.25 kcal/mol, respectively, and given that such low-energy vibrations are unlikely to be pure torsions, these experimentally derived barriers are in rather good agreement with those in Figure 11. 4.5. Proton-Transfer Mechanism in TAE. We examined the correlations in crystalline TAE between proton displacements along the O···O vector and various changes in internal and intermolecular coordinates of neighboring molecules, but the highest correlation coefficient we found is only −0.1, which is between the average separations of the four nearest methyl H atoms from the O2 atom and the corresponding O1···H separation. Such low correlations are consistent with the frequency of the proton transfer being considerably higher than that of the methyl libration and much higher than that of the methyl rotational jump. Therefore, the hydrogen-bond bridge experiences the mean field of the neighboring methyl groups, which becomes smoother as the librational amplitudes of the methyl groups increase at higher temperatures. In classical terms the interaction between the proton and its immediate environment of methyl groups is part Coulomb and part van der Waals. The net interaction is approximated in Figure 12 as between the sum of 1/r2 for all H atoms on the nearest methyl groups and the O1 atom of the O1···H···O2 bridging structure. We choose all of the methyl H atoms so that the effects of methyl rotation at higher temperature are included, and we use the O2 atom rather then the bridging H atom because it

proton potential in Figure 8, we calculated the torsional potentials of methyl groups for both the isolated molecule and crystalline TAE (Figure 11). From the relative amplitudes in Figure 11 it is clear that the crystal field interactions substantially affect the torsional barriers. Although individually weak, the C−H···O interactions in the crystal provide an increase of the torsional barrier from 0.4 to 1.4 kcal/mol at the donor side and to 1.8 kcal/mol at the acceptor side, resulting in the following profound effect on the torsional dynamics. As the methyl groups are anchored by the C−H···O interactions and the barriers are greater than kBT even at 298 K substantial portions of the torsional phase space are far less accessible than others, in particular those portions essential for symmetrization of the proton potential in Figure 8. Consequently, in the course of MD simulation the effective proton potential remains asymmetric and the proton exhibits an affinity to be nearer to O1 than O2. As demonstrated in section 4.2, no correlation between proton-transfer dynamics and methyl torsions could be observed from the MD trajectories of crystalline TAE. Correlation is absent because the methyl torsion barriers are high while the proton-transfer barrier is low, resulting in several proton-transfer events in the course of the MD simulation but few substantial torsional rearrangements of methyl groups. Nevertheless, the methyl groups dominate the proton-transfer potential in the TAE crystal, both intramolecular and intermolecular, but their slow rotation rate leads to little dynamical coupling to the proton transfer, a situation similar but somewhat the reverse of that found in the AC8 system. The low barrier to methyl rotation (∼0.46 kcal/mol), even in the crystal (similar to that in the gas-phase barrier), allows 2289

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

neighboring methyl groups are “blocked”, preventing the proton-transfer free energy from being symmetrized even at room temperature. Longer range phonons that include methyl libration components provide phonon assistance via a meanfield smoothing of proton-transfer free-energy effect, and it would be interesting to study the rotational dynamics of the methyl groups in more detail in other systems with a view to long-range proton transfer that could lead to high proton conductivity in an organic system. The term “phonon assistance” is perhaps inaccurate in the present case, as the phonons do not drive the proton transfer, but the effect is more analogous to negative thermal expansion where the mean field is modified as the amplitude of librational modes increases with temperature.30 Our neglect of nuclear quantum effects in calculations may seriously affect the computed frequency of proton dynamics. Nevertheless, the FE at low temperature should still be a reasonable approximation to the real PE, although the energy range over which we determine this will be smaller than that explored by the quantum system. This has been established using displacements of classical methyl group H atoms to produce calculated rotational PE functions that agree well with those measured by quantum rotational tunnelling.31 In the TAE crystal the proton-transfer potential energy surface is largely determined by the methyl groups, and their rotational potential is high enough that a classical treatment of their librational dynamics is appropriate. In sections 1 and 4.4 we contrasted TAE with AC because in both systems dynamical coupling of the methyl rotation to the proton-transfer process might occur. In AC this coupling was known8 to be absent because the crystal field is too weak and methyl rotation too rapid. Here we showed that in TAE the crystal field is too strong so that methyl rotation is too slow. This result implies that “ideal” systems could be engineered in which the methyl rotation rate closely matches that of proton transfer, potentially providing a route to efficient proton transport in the solid state.

Figure 12. Temperature dependence of an effective interaction distribution between the closest methyl H atoms and O2. This moves to longer distances with increasing temperature, showing how the potential becomes flat and symmetric for the H atom at 298 K.

represents an approximately stationary point. Figure 12 shows the gradual shift of this interaction to longer distances as the temperature increases from 40 to 110 K, which arises from the increasing amplitude of the methyl librations. The presence of methyl rotation at 298 K has a more pronounced effect giving a much wider range of distances. The nonuniform distribution of methyl-group orientations arises from the weak C−H···O interactions in the crystal. Clearly, the barrier provides unfavorable torsion angles that prevent symmetrization of the H-transfer potential.

5. CONCLUSION In hydrogen-bonded systems a significant part of the potentialenergy well for the proton comes from relaxation of the environment around the proton, this being easy to see for TAE as the methyl-group orientations for the isolated molecule in Figure 6. Perhaps surprisingly, the resulting more complex interactions with methyl groups on neighboring molecules in crystalline TAE reduces this polaron effect, notwithstanding the higher rotational potential for the methyl groups in the crystal. This leaves a rather stable but weak polaron that allows a proton-transfer rate that is too high for the polaron to follow and hence the natural asymmetry in the system that was first found crystallographically. As the librational amplitude of the methyl groups increase with temperature the polaron is weakened further, but its effects are still discernible at room temperature. We analyzed the results of MD simulations and lattice dynamics calculations of crystalline TAE with corresponding MD simulations for the isolated TAE molecule using DFT in considerable detail. This level of detail is justified by the good agreement we obtained between the results of the calculations and experimental crystallographic and spectroscopic results, as shown in Figures 2 and 3. Using a thermodynamic approach to describing the energy surfaces we were able to discern the thermal effects of various coordinates on the FE (see Figures 8−10) and hence on the proton-transfer process. Simple models that use displacement of the proton along the O···O vector may be able to discern the polaron effect described above, but it would be very difficult for any technique that does not invoke lattice dynamics in some way to find either phonon assistance or phonon hindrance. The experimentally observed asymmetry of the H-bond bridge in crystalline TAE arises because the reorientations of

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Prof. Dusan Hadzi, National Institute of Chemistry, Ljubljana, Slovenia, for initializing our interest in TAE and many stimulating discussions. Financial support from the Slovenian Research Agency (program no. P1-0012 and project no. J1-2014) is gratefully acknowledged.

■

REFERENCES

(1) Dopieralski, P. D.; Latajka, Z.; Olovsson, I. J. Chem. Theory Comput. 2010, 6, 1455−1461. (2) Fontaine-Vive, F.; Johnson, M. R.; Kearley, G. J.; Cowan, J. A.; Howard, J. A. K.; Parker, S. F. J. Chem. Phys. 2006, 124. (3) Fontaine-Vive, F.; Johnson, M. R.; Kearley, G. J.; Howard, J. A. K.; Parker, S. F. J. Am. Chem. Soc. 2006, 128, 2963−2969. (4) Martins, D. M. S.; Middlemiss, D. S.; Pulham, C. R.; Wilson, C. C.; Weller, M. T.; Henry, P. F.; Shankland, N.; Shankland, K.; Marshall, W. G.; Ibberson, R. M.; Knight, K.; Moggach, S.; Brunelli, M.; Morrison, C. A. J. Am. Chem. Soc. 2009, 131, 3884−3893. 2290

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291

The Journal of Physical Chemistry A

Article

(5) Morrison, C. A.; Siddick, M. M.; Camp, P. J.; Wilson, C. C. J. Am. Chem. Soc. 2005, 127, 4042−4048. (6) Piccoli, P. M. B.; Koetzle, T. F.; Schultz, A. J.; Zhurova, E. A.; Stare, J.; Pinkerton, A. A.; Eckert, J.; Hadzi, D. J. Phys. Chem. A 2008, 112, 6667−6677. (7) Koetzle, T. F.; Piccoli, P. M. B.; Schultz, A. J. Nucl. Instrum. Methods Phys. Res., Sect. A 2009, 600, 260−262. (8) Johnson, M. R.; Jones, N. H.; Geis, A.; Horsewill, A. J.; Trommsdorff, H. P. J. Chem. Phys. 2002, 116, 5694−5700. (9) Pirc, G.; Stare, J.; Mavri, J. J. Chem. Phys. 2010, 132, 224506− 224501. (10) Stare, J.; Panek, J.; Eckert, J.; Grdadolnik, J.; Mavri, J.; Hadzi, D. J. Phys. Chem. A 2008, 112, 1576−1586. (11) Daemen, L.: http://lansce.lanl.gov/lujan/instruments/FDS/. (12) Delley, B. J. Chem. Phys. 1990, 92, 508−517. (13) Delley, B. J. Chem. Phys. 2000, 113, 7756−7764. (14) CPMD, V3.13; Copyright IBM Corp.: 1990−2008, Copyright MPI fuer Festkoerperforschung Stuttgart: 1997−2001. (15) Perdew, J. P.; Burke, K.; Wang, Y. Phys. Rev. B 1996, 54, 16533− 16539. (16) Kresse, G.; Furthmuller, J. Phys. Rev. B 1996, 54, 11169−11186. (17) Kresse, G.; Furthmuller, J. Phys. Rev. B 1996, 54, 11169−11186. (18) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15−50. (19) Kresse, G.; Hafner, J. Phys. Rev. B 1993, 47, 558−561. (20) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758−1775. (21) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (22) Lee, C. T.; Wang, W. T.; Parr, R. G. Phys. Rev. B 1988, 37, 785− 789. (23) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 1993−2006. (24) Blochl, P. E. Phys. Rev. B 1994, 50, 17953−17979. (25) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (26) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13. (27) Parlinski, K. PHONON software; Cracow, 2010. http://www. computingformaterials.com/. (28) Parlinski, K. In Neutrons and Numerical Methods - N2M: Grenoble, France; Johnson, M. R., Kearley, G. J., Buttner, H. G., Eds.; Proceedings of the American Institute of Physics; AIP Press: Woodbury, NY, 1999; Vol. 479, pp 121−130. (29) Kneller, G. R.; Keiner, V.; Kneller, M.; Schiller, M. Comput. Phys. Commun. 1995, 91, 191−214. (30) Peterson, V. K.; Kearley, G. J.; Wu, Y.; Ramirez-Cuesta, A. J.; Kemner, E.; Kepert, C. J. Angew. Chem., Int. Ed. 2010, 49, 585−588. (31) Johnson, M. R.; Kearley, G. J. Annu. Rev. Phys. Chem. 2000, 51, 297−321.

2291

dx.doi.org/10.1021/jp210212q | J. Phys. Chem. A 2012, 116, 2283−2291