PCCP View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 24099
View Journal | View Issue
Ab initio metadynamics study on hydronium ion dynamics at acid-functionalized interfaces: effect of surface group density Swati Vartak, Anatoly Golovnev, Ata Roudgar and Michael Eikerling* This article presents an ab initio metadynamics study of elementary hydronium ion transitions at dense arrays of surface groups with sulfonic acid head groups. Calculations simulate minimally hydrated conditions of the interfacial ionic system. The specific focus is on the influence of the surface group density on hydronium ion transport. Results reveal a high sensitivity of the activation free energy of hydronium translocations to the surface group density. A spontaneous concerted transition with low activation barrier is found at a surface
Received 5th July 2014, Accepted 25th September 2014
group separation of 6.8 Å. When hydroniums translocate concertedly, the activation barrier of the transition
DOI: 10.1039/c4cp02937b
interaction constants of hydronium ions and anionic surface groups as well as the surface group flexibility
drops by more than a factor of two to the value of 0.25 eV. An approach is presented to determine from the analysis of frequency spectra. These properties are discussed in the context of a recently developed
www.rsc.org/pccp
soliton theory of interfacial proton transport.
1. Introduction Proton transport (PT) is a key phenomenon in acid–base chemistry, energy transduction in biology, and electrochemical energy conversion in polymer electrolyte fuel cells. The general interest in PT has stimulated a plethora of studies on PT in bulk water as well as at charged interfaces.1–8 PT in bulk water proceeds via structural diffusion of a protonic defect. Hydrogen bond breaking in the second solvation shell of an hydronium ion and reorganization of the surrounding water molecules contrives proton transport events in free bulk water.3,9 Bulk water-like proton transport is also the main contribution to the proton conductivity of well-hydrated polymer electrolyte membranes (PEMs). In the absence of water with bulk-like dynamic properties, proton transport inside of acid-lined nanochannels in PEMs or at acid-functionalized interfaces cannot proceed through the same mechanism as in free bulk water. Under low hydration conditions surface PT must make a significant contribution. Surface proton diffusion has been studied extensively in biologically relevant systems.8,10–15 It is believed that protons travel along a surface as a part of hydronium ions. Therefore, proton transport should be understood in the context of hydronium ion transport. Experiments at Langmuir monolayers of protogenic surface groups (SGs) enabled the control of the monolayer composition, surface pressure, and SG density. These experiments have revealed Department of Chemistry, Simon Fraser University, 8888 University Drive, Burnaby, BC, Canada V5A1S6. E-mail: [email protected]
This journal is © the Owner Societies 2014
that the rate of lateral proton transport depends on the nature of ionic groups at the interface as well as the distance between these groups. Owing to the different compositions of monolayer systems, varying pKa of ionic groups as well as experimental factors such as convective disturbances, a range of values of proton diffusion coefficients (105 to 106 cm2 s1) have been reported in the literature.11,12,16,17 In spite of some controversy surrounding these experimental results, it is clear that hydrogen bonded ionic groups at the water/membrane interface could enable a highly efficient mechanism of surface transport of protons. In hydrogen-bonded media, profiles of the proton free energy as a function of the proton coordinate, which is defined relative to equilibrium positions of donor and acceptor groups, exhibit a double well character.18 For PT in bulk water, donor and acceptor water molecules form an intermediate Zundel-ion state, with the transferring proton delocalized between two neighboring water molecules. At the hydrogen-bond distance that corresponds to the formation of the Zundel-ion, the double well potential is symmetrical and the barrier in the proton free energy vanishes. For proton transfer at an acid-functionalized interface, the potential barrier for PT is significantly higher, being related to the strong localization of water molecules around ionic SGs. It is however expected that decreasing the distance between ionic groups can lower the proton transfer barrier markedly. Upon densification of SGs, a critical SG separation has been found in experiment, at which a transition to highly efficient surface PT was observed. The experimental value of the critical SG separation is B7 Å.11,15,19
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24099
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
Paper
Driven by the interest in biochemical phenomena, experimental monolayer studies have been conducted mainly with surface groups (SGs) that contain carboxylic acid head groups. Recently, Fujita et al.20 have measured proton conductivity of mesoporous silica films functionalized with sulfonic acid groups. At higher interfacial acid density than that of Nafion, the proton conductivity was measured as 5.4 103 S cm1 at a relative humidity of 20% and the activation energy for PT was determined as 0.26 eV. Ilhan and Spohr carried out ab initio molecular dynamics simulations on proton networks in polymer electrolyte pores modeled using CF3CF3 and CF3SO3H groups.21 They compared proton and water distribution in two pore channels at different hydration levels. They observed high proton mobility in the case of a narrower pore, in which the distance between acid head groups was between 6 and 8.5 Å. Despite a number of studies at interfaces comprised of carboxylic and sulfonic acid head groups, the mechanism of proton transport is not understood yet. Theoretical approaches such as soliton models have been employed to explain a concerted motion of protons along the surface of a one-dimensional hydrogen-bonded system, resulting in high surface proton mobility.22–24 Water molecules bind strongly to the surface, forming hydronium ions whose motion along the surface results in proton transport. Dislocation of one hydronium causes dislocation of a neighboring hydronium ion. Solitary waves of such dislocations travel through the system, leaving a trace, because each translocation changes the local configuration of SGs. Therefore, in a one-dimensional system, the passing of a protonic soliton should be followed by a reorientation of SGs so that another solitary wave can pass through the system; such a requirement results in low proton conductivity. As a consequence of this restriction, the onedimensional soliton approach was essentially abandoned. However, the problem can be overcome by considering a two-dimensional system25 since it allows multiple paths for a soliton to pass through the system and restore the initial configuration. Our research on this topic focuses on ordered interfacial systems of SGs with sulfonic acid head groups. Previous ab initio simulation studies on this system have explored the impact of molecular structure and packing density of SGs on spontaneous ordering, acid dissociation, water binding, and hydronium ion dynamics at interfaces.26–29 The main finding of structural
PCCP
studies was a transition to a condensed surface state at a critical value of the SG density. In the condensed state, completely dissociated sulfonate anions and hydronium ions form a hexagonally ordered ionic crystal in 2D. The saturated number of H bonds renders the interface superhydrophobic. In ref. 29, we have employed the metadynamics method to study interfacial hydronium transitions. Proton transitions were evaluated at an SG spacing of 6.6 Å, which corresponds to an area of 37.5 Å2 per SG. We observed a strong energetic preference for a highly concerted proton transfer. Based on these findings, we have developed a soliton theory of concerted interfacial proton motion in this system.25,30 All of these studies on interfacial hydronium ion dynamics have been performed at a uniform surface group density, which corresponds to the critical value for the transition to the superhydrophobic ordered surface state. At interfaces in real PEMs, acid-terminated sidechains are not uniformly grafted on the pore walls formed by ionomer aggregates. Moreover, the flexibility of sidechains contributes to fluctuations in sidechain separations. It is expected that the interfacial dynamics of hydronium ions and SG will be influenced by these fluctuations in SG density. In this work, we have extended the calculations to different SG densities to understand the impact of density variations on mechanisms and energetics of interfacial hydronium ion motions. Section 2 presents the model system and reviews the main results obtained in ref. 29, including the importance of the SG spacing (dCC) of 6.6 Å. Section 3 provides computational details of the simulations. Section 4 presents the analysis of hydronium ion transitions in terms of free energy profiles, hydrogen bond strength, and elastic constants of the coupled system of hydronium ions and anionic SGs found from the frequency spectra.
2. Model system The model system is depicted in Fig. 1. It consists of a densely packed array of CF3SO3H (triflic acid) surface groups. Terminal C atoms of all SGs are fixed at hexagonal lattice positions. This is the system that was evaluated in previous studies of stable interfacial configurations and wettability.26,28 The shortest
Fig. 1 Model system consisting of 2D arrays of completely ionized CF3SO3H surface groups (SG) with one water molecule added per SG. F, S, O, H and C atoms are shown by yellow, red, white, orange and blue circles, respectively. An upright structure is shown on the left and a tilted structure is shown on the right side of the figure. The dotted lines indicate the unit cell used in the simulations in which C atoms are fixed at their positions. Donor, acceptor and spectator SG are labeled with D, A and S respectively and the distances that define the collective variable for the lateral shift of an H3O+ ion are shown, CV = d12 d23. On the left, red arrows indicate hydronium ion transitions that will transform the tilted structure into the upright structure.
24100 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
Paper
3. Computational details
concerted proton translocations, we considered a unit cell consisting of 12 SGs and 12 water molecules with periodic boundary conditions as indicated by dashed lines in Fig. 1. The unit cell area for SG separations of 6.4, 6.6 and 6.8 Å was 424, 456 and 482 Å2 respectively. Ab initio metadynamics utilizes a repulsive history-dependent potential in the form of Gaussian functions to achieve an efficient sampling of the configuration space. The method allows describing the system by a small number of collective variables (CV); after numerous trials, we concluded to use one CV. During simulations, the reduced phase space of the CV is explored and Gaussian potentials are added at points visited by the system. The accumulation of Gaussian functions discourages the system from revisiting points already visited. The calculation is completed when the system reaches a diffusive state, in which the addition of Gaussian functions has evened out the free energy surface. The ensemble of Gaussian functions is then subtracted revealing a positive image of the free energy landscape. The sampling efficiency and the precision of the simulation depend on the height and width of the Gaussian functions and hence these parameters must be optimized. Adding Gaussian functions to the system adds energy to the degrees of freedom associated with CVs. This may lead to a significant inhomogeneity in the temperature distribution of the system, and may also affect the dynamics. To address this issue, a harmonic coupling between the real and fictitious CV is introduced in the metadynamics Lagrangian, which ensures uniform sampling of configurations in well regions of the free energy landscape.31 The CP2K package uses a mixed Gaussian and plane wave basis set.34 Valence electrons are represented by double-z augmented Gaussian basis sets (DZVP-MOLOPT).35 The pseudopotential of Goedecker, Teter and Hutter (GTH) represents core electrons.36 We used an energy cut-off of plane-wave expansions of 300 Ry (1 Ry E 13.6 eV). Exchange and correlation energies were computed within the GGA approximation using the BLYP functional.37 For every time step, the electronic structure was relaxed to an accuracy of 107 Hartree (1 Hartree E 27.2 eV). Optimized structures were thermalized for B3 ps. The temperature was set to 300 K using a Nose–Hoover thermostat and the timestep was set to 0.5 fs. In the metadynamics Lagrangian, we used a coupling constant k = 0.4 au (1 au E 1.55 103 N m1) and a fictitious particle mass M = 75 au (1 au E 1.67 1027 kg). Gaussian potentials with height h = 0.013 eV and width Dd = 0.014, 0.02, 0.016 Å for dCC = 6.4, 6.6 and 6.8 Å, respectively, were used. The height of Gaussians was chosen such that it is large enough to optimize the computational cost and small enough to conserve the total energy of the system and the width was calculated as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dd ¼ hs2 i hsi2 ; (1)
We employ the ab initio metadynamics method of Laio and Parrinello31,32 to simulate free energy surfaces and proton transition pathways in the highly charged interfacial system.29 Metadynamics at the DFT level is incorporated in the CP2K package.33 In order to explore the possibility of spontaneous
where s is the value of the collective variable. The CV was defined as the lateral shift of the transferring H3O+ ion, dCV = d12 d23, as illustrated in Fig. 1. Addition of consecutive Gaussians was done when the displacement of CV relative to the previous state reached 3Dd/2.
possible surface group is considered for computational convenience. In previous work, longer surface groups have been evaluated as well.28 These studies have revealed that a restriction to triflic acid does not pose any major restriction in terms of interfacial structure and dynamics. The model system is evaluated in the limit of high SG density. The interfacial area per surface group is in the range of 35–40 Å2 corresponding to nearest neighbor SG separation distances of 6.4 to 6.8 Å. A vacuum layer of 15 Å in height was introduced in the z-direction. One water molecule is added per group, representing conditions of minimal hydration. In spite of the small amount of water added, anionic head groups are found in a completely dissociated state owing to the high SG density. In ref. 26 it was found that at such SG separations there are two stable structures, upright and tilted, as shown in Fig. 1. Transitions between these structures occur by elementary hydronium ion transitions, as illustrated in Fig. 1 with an arrow. A sequence of such elementary transitions leads to hydronium transport along the surface. In ref. 25 we discussed what particular sequence of elementary transitions would facilitate hydronium transport most efficiently. The present paper focuses on the energetics of elementary interfacial proton transitions. We started this consideration in ref. 29, where a single transition of a hydronium ion was considered at a SG separation of 6.6 A. This separation corresponds to the case, where the energies of the upright and tilted structures coincide. Therefore, this case might be most favorable for the hydronium ion mobility because back and forth transitions between the structures are isoenergetic. On the other hand, in real systems SGs are never grafted precisely at regular lattice positions and as a result the distance between different SGs varies. This article explores the influence of a SG separation change on the energetics of a hydronium ion hop. We investigate what happens when the energy of one of the structures is larger. At 6.4 Å, the energy of the tilted structure is larger than upright structure by 0.7 eV. At 6.8 Å, the energy of the upright structure is larger than the tilted structure by 0.8 eV. At SG separations further away from the critical value, the energy difference between upright and tilted structures increases substantially, thereby rendering transformations between the two structures energetically very expensive. Also, the interfacial array will become strongly hydrophilic at larger SG separations. Mechanisms of surface hydronium ion transport, as considered in our work, will be stifled by surface to bulk transfer of protons followed by bulk-like proton motion. For the dependence of energies of the interfacial structures on the SG separation see ref. 26.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24101
View Article Online
Paper
PCCP
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
We obtained frequency spectra by taking the Fourier transform of the velocity autocorrelation function, gvv(t), calculated from the metadynamics trajectory. The function gvv(t) was calculated as an average of scalar products of atomic velocities, [next-line-equation], calculated as gvv ðtÞ ¼
hvi ð0Þ vi ðtÞi ; hvi ð0Þ vi ð0Þi
(2)
where vi(0) is the initial velocity of atom i and vi(t) is its velocity at time t. The average was calculated over various metadynamics steps chosen for vi(0) and vi(t). A Fourier transformation of this function was then performed to obtain the frequency spectrum, ð1 IðoÞ ¼ eiot gvv ðtÞdt (3) 0
The function I(o) gives information about specific frequency modes. The spectra were obtained for the motion of all atoms as well as specific modes, namely, stretching of O–H bonds and tilting of SGs.
4. Results and discussion In ref. 29, we applied the metadynamics method to study two cases: (1) collective transfer of hydronium ions and (2) disconcerted transfer of hydronium ions. In the collective mechanism, a transition of one hydronium ion from a donor to an acceptor SG triggers the transition of the neighboring hydronium ion from an acceptor SG to the next acceptor SG and so forth. The number of hydrogen bonds per SG is conserved in this mechanism. The disconcerted mechanism, on the other hand, creates a hydrogen bond defect with the formation of 4 hydrogen bonds at the acceptor SG and only 2 hydrogen bonds at the donor SG. In these calculations, we obtained barrier energies of B0.3 eV for the collective transfer and B0.6 eV for a disconcerted transfer of a single hydronium ion. These calculations showed that the collective mechanism of hydronium ion transfer is energetically favorable over the disconcerted mechanism.
A sequence of such hydronium relocation steps can lead to long-range transport of protons in the system. In the work presented here we applied the metadynamics method with a single CV to study hydronium ion transfer at dCC = 6.4, 6.8 Å. At the shortest SG separation distance of 6.4 Å, the upright structure is lower in energy than the tilted structure by 0.7 eV and at 6.8 Å, the tilted structure is more stable than the upright structure by 0.8 eV, which facilitates the proton transition. The metadynamics results at these two SG separations are compared with those obtained previously at dCC = 6.6 Å. Hydronium ion transitions were completed in 100, 115 and 40 ps at dCC = 6.4, 6.6 and 6.8 Å, respectively. The reconstructed Helmholtz free energy profiles for dCC = 6.4 and 6.6 are shown in Fig. 2 together with snapshots along the trajectory. The activation energy was calculated as 0.65, 0.6 and 0.25 eV, respectively, for dCC = 6.4, 6.6 and 6.8 Å. The hydronium ion transitions at dCC = 6.4 and 6.6 Å, shown in Fig. 2, give rise to the formation of a hydrogen bond defect, where the donor SG becomes deficient in hydrogen bonds and the acceptor has 4 hydrogen bonds. The structures formed along the trajectories are similar for the SG spacings of 6.4 and 6.6 Å. The free energy profile for dCC = 6.8 Å is shown in Fig. 3. At this SG spacing we observed a spontaneous transition of another hydronium ion. The hydronium ion labeled as 3 (see Fig. 3) underwent flipping following the transition of hydronium 1. However, it reverted back and hydronium 2 underwent flipping. The transition of hydronium 2 was accompanied by reorientation of the donor SG. The O–O distances between the oxygen nucleus of the transferring hydronium ion (OT) and oxygen nuclei of spectator (OS), donor (OD) and acceptor (OA) SGs, respectively, are shown in Fig. 4 at SG separations of 6.4, 6.6 and 6.8 Å. Red and black curves correspond to the OT–OS distance, the blue curve corresponds to the OT–OD distance and the green curve shows the OT–OA distance. The OT–OS distances remain in the range of 2.5 to 2.7 Å. As expected, the OT–OD distance increases and the OT–OA distance decreases.
Fig. 2 Reconstructed free energy surfaces (FES) for local defect-type H3O+ ion transitions at SG spacing of 6.4 Å (left panel) and 6.6 Å (right panel). Snapshots show structures at points A, B, C, D marked on the FES. Since similar structures are formed at both the SG densities, only one set of structures at the SG separation of 6.4 Å is shown.
24102 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
Paper
Fig. 3 Reconstructed free energy surfaces (FES) for H3O+ ion transitions at SG spacing of 6.8 Å. Snapshots show structures at points A, B, C, D marked on the FES. The transition of ion 1 is followed by the transition of ion 3, which then reverts back and ion 3 undergoes flipping. This concerted transition of hydronium ions 1 and 2 conserves the number of hydrogen bonds per SG.
Fig. 4 O–O distances as a function of simulation time for SG separations of 6.4, 6.6 and 6.8 Å. Black and red curves represent the distance between O atom of spectator SGs and O atom of a hydronium ion undergoing transition. The blue curve represents the distance between O atom of a donor SG and hydronium O and the green curve represents the distance between an acceptor SG and the transferring hydronium ion. Values are averaged over intervals of 1000 time steps to reduce the noise.
With only a slight decrease in SG density, the hydronium ion transport mechanism changes from disconcerted to collective with significantly lower activation energy. This effect is believed to be a consequence of the fact that at dCC = 6.8 Å the final tilted structure is energetically preferred over the initial upright structure. To analyze changes in dynamic mechanism with SG density, we calculated frequency spectra of the system, obtained from Fourier transformation of the velocity autocorrelation function, as depicted in Fig. 5. In order to be able to assign the frequency
peaks, we calculated separately the stretching frequencies of hydronium OH bonds at different SG separations, shown in Fig. 6, as well as tilting frequencies of SGs, shown in Fig. 8. We found 3 O–H stretching frequency bands in the 100–4000 cm1 range: the most intense band occurs in the 130–1330 cm1 range with a peak at around 675 cm1, a lower intensity band lies between 1450–2260 cm1 with a peak at around 1720 cm1 and the lowest intensity band is observed in the 2450–3600 cm1 range. The frequency spectra for all three SG separations are
Fig. 5 Full frequency spectra at different SG separations obtained from Fourier transform of the velocity autocorrelation function.
Fig. 6 Stretching frequency spectra of hydronium O–H bonds at different SG separations.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24103
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
Paper
PCCP
overlapping almost completely indicating that the variation in SG density exerts only a small effect on the stretching frequencies of O–H bonds. This is an important finding in this work since it shows that in this range of SG separations, the length and strength of hydrogen bonds is essentially unaltered. Buzzoni et al.38 studied IR and Raman spectra of water molecules interacting with Nafion. In order to understand the basic features of the IR spectra and assign absorption peaks, they proposed an empirical model for a solvated Zundel ion. They categorized hydrogen bonds into ‘internal’, ‘external’ and ‘solvated’ bonds based on O–O distance. For internal hydrogen bonds with O–O distance in the range of 2.45–2.57 Å, OH–O stretching frequencies are expected to be less than 2000 cm1. The external hydrogen bonds with O–O distance in the range of 2.6–2.7 Å are expected to show absorption within the range from 2500 to 3200 cm1. Finally for solvated hydrogen bonds, for which the O–O distance is larger than 2.7 Å, the frequencies are expected in the range of 3200 to 3600 cm1. We have observed the peaks roughly in similar frequency ranges. Comparing the relative intensities of the peaks we can say that the O–O distances mostly fall in the range of ‘internal’ H-bonds. To understand the effect of bond length fluctuations on stretching frequencies, we calculated the average O–O distances of different hydrogen bonds and the range of fluctuations around the average value. We divided hydrogen bonds into 3 categories: (1) donor hydrogen bonds, (2) hydrogen bonds of spectator SGs directly bonded to the transferring hydronium ion, and (3) hydrogen bonds of SGs farther away from the transferring hydronium ion. The range of average O–O distances of various H-bonds is given in Table 1. The time averaged O–O distances of different H-bonds lie in the range of 2.57–2.7 Å, which corresponds to a 2.5% variance. However, this dispersion might not be relevant because time-fluctuations of the O–O distance of any single H-bond is around 5%. All O–O distances in the system remain fairly constant upon densification of the SG array. This is achieved by a relaxation phenomenon where hydronium ions shift perpendicular to the plane spanned by SGs. Fluctuations in the height of the hydronium oxygen with respect to the plane drawn through S atoms of SGs is shown in Fig. 7. It shows that at shorter SG separation, the hydronium oxygen remains predominantly above the plane, whereas at longer SG separation, hydronium ions descend closer to the plane. We also calculated SG frequency spectra at 3 SG separations as depicted in Fig. 8. As shown in the figure, the SG separation
Fig. 7 Fluctuations in the height of O atom of a spectator hydronium ion. Values are averaged over intervals of 1000 time steps to reduce the noise.
Fig. 8 SG frequency spectra calculated at 3 SG separations.
in this small range does not alter the frequencies of the peaks. To interpret the spectra, let us consider small oscillations of SGs and hydronium ions about their equilibrium positions. As discussed in ref. 25, we consider only pair interactions in the nearest neighbor approximation. In general, in a two-component system, there are three distinct pairs: hydronium ion–hydronium ion, SG–SG, and hydronium ion–SG. The corresponding elastic constants are KHH, KSS, and KHS. These constants appear explicitly in theoretical models, see ref. 23–25. Additionally, SGs are fixed at the basal plane, which gives rise to the force returning SGs to their equilibrium position, described by an elastic constant KF. For such a system a standard approach for diatomic crystal lattices can be employed, with an additional term describing the flexibility of SGs. A schematic of the model is shown in Fig. 9.
Table 1 Influence of the hydronium ion transition on bond lengths of surrounding groups. The ranges of average O–O distances at different SG densities are given in the table. The nearest spectator groups are the SGs directly bonded to the transferring hydronium ion; the SGs further away from the transferring hydronium ion are designated as far spectators in the table
dCC (Å)
Range of average O–O distances (Å)
6.4
6.6
6.8
Donor Nearest spectators Far spectators
2.59 2.61–2.62 2.60–2.62
2.57 2.64–2.65 2.61–2.62
2.64 2.64–2.70 2.61–2.64
24104 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
Fig. 9 A schematic representation of SGs (black circles) connected to hydronium ions (red circles) by springs with spring constant KHS. The coupling constant between SGs is KSS and that between hydronium ions is KHH. Finally, KF describes the flexibility of SGs, which tends to return them to their equilibrium positions.
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
Paper
For simplicity, we will consider only a one-dimensional string of masses connected with springs. Such a 1D system captures essential features of the 2D lattice due to the hexagonal symmetry of the lattice. Because we consider only small oscillations of hydronium ions, each hydronium and each SG is subjected to a circularly symmetric potential, which is a function of the displacement from the equilibrium position. In 1D approximation, the displacement of the j-th atom is u2j = BH exp(i(2jkdCC ot)),
(4)
u2j+1 = BS exp[i((2j + 1)kdCC ot)],
(5)
where BH and BS are unknown amplitudes and k = 2p/l is a wave vector. Since the problem is 1D, we do not use vector notations. The value of l is defined by the boundary conditions in our simulations, l = 6 dCC/O3. According to Hook’s law, the force acting on the j-th atom is F2j = KHS(u2j+1 + u2j1 2u2j),
(6)
F2j+1 = KHS(u2j+2 + u2j 2u2j+1) KFu2j+1.
(7)
Inserting these equations into Newton’s second law and eliminating BH and BS, we get 1 1 KF þ þ o1;22 ¼ KHS mH mS 2mS KHS sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 1 1 KF 2KF 4 sin2 kdCC KHS þ þ mH mS 2mS KHS mH mS KHS mH mS (8) where mH and mS represent the masses of hydronium ions and SGs, respectively. The signs correspond to optical and acoustic modes. The third peak in Fig. 8 corresponds to SG–SG interactions. It can be modeled by a string of particles with the same mass connected by springs. The mathematical formalism is similar to the one shown above. The resulting frequency is 1 kdCC o32 ¼ þ KF : 4KSS sin2 (9) mS 2 The argument of the sin is different from the one in eqn (8) since the distance between SGs is twice as large as that between SGs and hydronium ions. To assign the frequencies o1,2,3 to the peaks in Fig. 8, one has to use the corresponding spectrum of hydronium ion oscillations. The frequencies o1,2 should be the same for SG and hydronium ions. The third peak, corresponding to oscillations in the hydronium subsystem only, is 1 kdCC 0 o3 2 ¼ ; (10) 4KHH sin2 mH 2 which is analogous to eqn (9). Unfortunately, it is impossible to single out the peaks in the hydronium spectrum because peaks are broad and merged together as shown in Fig. 10. However, it is clear that there are no peaks at frequencies above 700 cm1. Therefore, the peak observed at 1000 cm1 in Fig. 8 corresponds to o3. Hence, the leftmost peak is an acoustic mode, and the middle peak is the optical mode.
This journal is © the Owner Societies 2014
Fig. 10 Frequency spectrum of a hydronium ion.
The broadness of the peaks for hydronium ions, in contrast to well-distinguished peaks for SGs, can be explained based on the mechanism of energy loss suggested in ref. 30. When a hydronium ion or an SG shifts, it pulls on its nearest neighbors in transverse direction, which transfers energy to the remaining lattice. A heavy SG can easily pull a hydronium ion, whereas it is difficult for a light hydronium ion to pull a SG, which is fixed at the basal plane. Therefore, hydronium ions lose energy faster than SGs, which broadens the hydronium ion peaks. Since we can extract values of the characteristic mode frequencies from Fig. 8 and 10, eqn (8) and (9) can be solved for the elastic constants. Using o1 = 200 cm1, o2 = 625 cm1, o3 = 995 cm1, o30 = 700 cm1, mH = 19 amu, mS = 80 amu, k = 2.75 109 m1 and dCC = 6.6 1010 m, we obtain KF = 1808.68 N m1, KHS = 22.58 N m1 and KSS = 1165.02 N m1, KHH = 222.6 N m1. The comparatively small value of KHS can be explained based on the relaxation mode of hydronium ions relative to the plane spanned by SGs, discussed above. To accommodate a change of dCC, hydronium ions tend to move away from the basal plane when SGs are squeezed together and approach the basal plane on increasing the SG separation. This relaxation of hydronium ions perpendicular to the plane of SGs keeps the H-bond length relatively constant and softens the coupling between hydronium ions and SGs. The strength of coupling between hydronium ions plays a pivotal role for the collectiveness of the hydronium ion transitions. In ref. 25 and 30, we have developed a soliton theory of collective hydronium ion motion at the minimally hydrated interface of triflic acid SGs. In the fully dissociated condensed state of the interfacial system,26 the approach distinguishes two sublattices: the sublattice of sulfonate anions that fluctuate about their equilibrium positions and the sublattice of mobile hydronium ions. The two sublattices are coupled by strong hydrogen bonds. In the genuine soliton Hamiltonian,22–24 discussed in ref. 25, harmonic coupling terms represent interactions in the system of SGs and hydronium ions. To simplify the mathematical formalism, interaction terms for SG–SG coupling and SG–hydronium ion coupling are replaced by an effective potential, whereas hydronium–hydronium coupling is still treated explicitly in harmonic approximation. The resulting Hamiltonian for single file soliton motion, illustrated in Fig. 11, is X m k H¼ (11) u_ i2 þ ðuiþ1 ui Þ2 þV ðui Þ 2 2 i
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24105
View Article Online
Paper
PCCP
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
system. In the future, further theoretical studies and molecular simulations will be performed specifically to understand mechanisms of soliton creation and the energetic cost associated with it.
5. Conclusions Fig. 11 Lattice configurations of the 2D two-component lattice of H3O+ and SO3 ions. The native 2D crystal structure in (a) is shown at the critical SG density of dCC E 6.6 Å. In the introduced symbolic representation, blue triangles indicate hydronium positions and red dots represent positions of sulfonate anions. The structure for a single file hydronium ion transfer, which conserves the number of H-bonds per SG, is shown in (b).
where ui is the displacement of the ith hydronium ion, m is the mass of a hydronium ion, k is the coupling constant between neighboring hydronium ions and V(ui) is the effective potential. The Hamiltonian equations of motion of this system were solved analytically in ref. 25 for different potential functions, V(ui), in the continuum limit, assuming that the soliton width is much larger than dCC. The size of the soliton is given by sffiffiffiffiffiffiffi ka2 N0 ¼ (12) 2e It thus depends on hydronium–hydronium coupling constant, k, the hydronium relocation distance, a and the depth of the potential well, e. The effective potential V(ui) in the soliton Hamiltonian can be parameterized using free energy curves obtained from Metadynamics simulations, as reported in ref. 29. The main parameters, which define this potential and thus the static properties of solitons, are the hydronium relocation distance, a, and the depth of the potential well, e. The main parameter for the soliton dynamics is the hydronium–hydronium coupling constant, k. Using estimates for the coupling constants and scaling arguments, an explicit expression for the soliton mobility was derived in ref. 30 and a rough estimate of the mobility value was obtained. Except for the coupling constant k, values for other parameters were obtained from Metadynamics simulations. The value of the coupling constant was estimated to be B2000 N m1, based on eqn (12) and assuming that the soliton consists of 50 hydronium ions. The coupling constant KHH reported here is derived from frequency spectra, which consider only small displacements of hydronium ions, treated in harmonic approximation, whereas for soliton-like transport the coupling constant k is associated with large displacement of hydronium ions. In the limit of weak coupling, KHH would be much smaller than k and hence KHH would not provide a reasonable estimate of the soliton size. However, in the case of strong coupling, KHH could be used in eqn (12) to estimate the size of solitons. Thus, based on the KHH value obtained here, the static soliton would involve 20–25 hydronium ions. This size of soliton justifies a posteriori the continuum approximation used in the soliton theory in ref. 25 and 30. It implies a strong coupling of hydronium ions in our
24106 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
We have performed ab initio metadynamics simulation to determine an activation barrier for the transition of hydronium ions at different SG spacing. We observed that at a SG separation of 6.4 and 6.6 Å, the energy barrier for the hydronium ion transition is high and gives rise to the formation of a hydrogen bond defect whereas at a SG separation of 6.8 Å, a spontaneous transition of a second hydronium ion remedies the hydrogen bond defect and conserves the number of hydrogen bonds in the unit cell. The activation energy in this case (0.25 eV) was found to be much lower than that obtained for the other transitions that result in the formation of hydrogen bond defects (Z0.6 eV). In spite of this change in activation energy, the stretching frequency spectra showed that in this range of SG separations, hydrogen bond strength remains essentially unaltered. Coupling constants between surface groups and hydronium ions were calculated from frequency spectra. They represent the coupling strength of the system for small fluctuations relative to the equilibrium configuration. The coupling constant between neighboring hydronium ions has been employed to estimate the size of collective, soliton-like excitations of hydronium ions. We determined that a soliton include about 20–25 hydronium ions, which would be characteristic of a strongly coupled system. We believe that the insensitivity of frequency spectra, hydrogen bond strength and coupling constants to the density of surface groups is owed to a peculiar relaxation motion of hydronium ions perpendicular to the plane of SGs.
References 1 O. Markovitch and N. Agmon, J. Phys. Chem. A, 2007, 111, 2253. 2 O. Markovitch, H. Chen, S. Izvekov, F. Parsani, G. A. Voth and N. Agmon, J. Phys. Chem. B, 2008, 112, 9456. 3 D. Marx, M. E. Tuckerman, J. Hutter and M. Parrinello, Nature, 1999, 397, 601. 4 M. E. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Phys. Chem., 1995, 99, 5749. 5 M. E. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Chem. Phys., 1995, 103, 150. 6 A. Y. Mulkidjanian, J. Heberle and D. A. Cherepanov, Biochim. Biophys. Acta, 2006, 1757, 913. 7 H. Morgan, D. M. Taylor and O. N. Oliviera Jr., Biochim. Biophys. Acta, 1991, 1062, 149. 8 B. Gabriel and J. Teissie, Proc. Natl. Acad. Sci. U. S. A., 1996, 93, 14521. 9 D. Marx, M. E. Tuckermann and M. J. Parrinello, J. Phys.: Condens. Matter, 2000, 12, A153.
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
10 J. Heberle, G. Riesle, D. Thiedemann, D. Oesterhelt and N. A. Dencher, Nature, 1994, 370, 379. 11 J. Zhang and P. R. Unwin, J. Am. Chem. Soc., 2002, 124, 2379. 12 S. Serowy, S. M. Saparov, Y. N. Antonenko, W. Kozlovsky, V. Hagen and P. Pohl, Biophys. J., 2003, 84, 1031. 13 I. Sakurai and Y. Kawamura, Biochim. Biophys. Acta, 1987, 904, 405. 14 K. Leberle, I. Kempf and G. Zundel, Biophys. J., 1989, 55, 637. 15 V. B. P. Leite, A. Cavalli and O. N. Oliviera Jr, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 6835. 16 J. Teissie, M. Prats, S. Soucaille and J. F. Tocanne, Proc. Natl. Acad. Sci. U. S. A., 1985, 82, 3217. 17 A. Polle and W. Junge, Biophys. J., 1989, 56, 27. 18 L. I. Krishtalik, Biochim. Biophys. Acta, 2000, 1458, 6. 19 C. J. Slevin and P. R. Unwin, J. Am. Chem. Soc., 2000, 122, 2507. 20 S. Fujita, A. Koiwai, M. Kawasumi and S. Inagaki, Chem. Mater., 2013, 25, 1584. 21 M. A. Ilhan and E. Spohr, J. Phys.: Condens. Matter, 2011, 23, 234104. 22 S. Pnevmatikos, Phys. Rev. Lett., 1988, 60, 1534. 23 A. V. Zolotaryuk, S. Pnevmatikos and A. V. Savin, Phys. Rev. Lett., 1991, 67, 707. 24 M. Peyrard, S. Pnevmatikos and N. Flytzanis, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 903.
This journal is © the Owner Societies 2014
Paper
25 A. Golovnev and M. Eikerling, J. Phys.: Condens. Matter, 2013, 25, 045010. 26 A. Roudgar, S. P. Narasimachary and M. Eikerling, J. Phys. Chem. B, 2006, 110, 20469. 27 A. Roudgar, S. P. Narasimachary and M. Eikerling, Chem. Phys. Lett., 2008, 457, 337. 28 S. P. Narasimachary, A. Roudgar and M. Eikerling, Electrochim. Acta, 2008, 53, 6920. 29 S. Vartak, A. Roudgar, A. Golovnev and M. Eikerling, J. Phys. Chem. B, 2013, 117, 583. 30 A. Golovnev and M. Eikerling, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 062908. 31 A. Laio and F. L. Gervasio, Rep. Prog. Phys., 2008, 71, 126601. 32 A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562. 33 http://www.cp2k.org. 34 J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167(2), 103. 35 J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105. 36 S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703. 37 A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098. 38 R. Buzzoni, S. Bordiga, G. Ricchiardi, G. Spoto and A. Zecchina, J. Phys. Chem., 1995, 99, 11937.
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24107
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 24099
View Journal | View Issue
Ab initio metadynamics study on hydronium ion dynamics at acid-functionalized interfaces: effect of surface group density Swati Vartak, Anatoly Golovnev, Ata Roudgar and Michael Eikerling* This article presents an ab initio metadynamics study of elementary hydronium ion transitions at dense arrays of surface groups with sulfonic acid head groups. Calculations simulate minimally hydrated conditions of the interfacial ionic system. The specific focus is on the influence of the surface group density on hydronium ion transport. Results reveal a high sensitivity of the activation free energy of hydronium translocations to the surface group density. A spontaneous concerted transition with low activation barrier is found at a surface
Received 5th July 2014, Accepted 25th September 2014
group separation of 6.8 Å. When hydroniums translocate concertedly, the activation barrier of the transition
DOI: 10.1039/c4cp02937b
interaction constants of hydronium ions and anionic surface groups as well as the surface group flexibility
drops by more than a factor of two to the value of 0.25 eV. An approach is presented to determine from the analysis of frequency spectra. These properties are discussed in the context of a recently developed
www.rsc.org/pccp
soliton theory of interfacial proton transport.
1. Introduction Proton transport (PT) is a key phenomenon in acid–base chemistry, energy transduction in biology, and electrochemical energy conversion in polymer electrolyte fuel cells. The general interest in PT has stimulated a plethora of studies on PT in bulk water as well as at charged interfaces.1–8 PT in bulk water proceeds via structural diffusion of a protonic defect. Hydrogen bond breaking in the second solvation shell of an hydronium ion and reorganization of the surrounding water molecules contrives proton transport events in free bulk water.3,9 Bulk water-like proton transport is also the main contribution to the proton conductivity of well-hydrated polymer electrolyte membranes (PEMs). In the absence of water with bulk-like dynamic properties, proton transport inside of acid-lined nanochannels in PEMs or at acid-functionalized interfaces cannot proceed through the same mechanism as in free bulk water. Under low hydration conditions surface PT must make a significant contribution. Surface proton diffusion has been studied extensively in biologically relevant systems.8,10–15 It is believed that protons travel along a surface as a part of hydronium ions. Therefore, proton transport should be understood in the context of hydronium ion transport. Experiments at Langmuir monolayers of protogenic surface groups (SGs) enabled the control of the monolayer composition, surface pressure, and SG density. These experiments have revealed Department of Chemistry, Simon Fraser University, 8888 University Drive, Burnaby, BC, Canada V5A1S6. E-mail: [email protected]
This journal is © the Owner Societies 2014
that the rate of lateral proton transport depends on the nature of ionic groups at the interface as well as the distance between these groups. Owing to the different compositions of monolayer systems, varying pKa of ionic groups as well as experimental factors such as convective disturbances, a range of values of proton diffusion coefficients (105 to 106 cm2 s1) have been reported in the literature.11,12,16,17 In spite of some controversy surrounding these experimental results, it is clear that hydrogen bonded ionic groups at the water/membrane interface could enable a highly efficient mechanism of surface transport of protons. In hydrogen-bonded media, profiles of the proton free energy as a function of the proton coordinate, which is defined relative to equilibrium positions of donor and acceptor groups, exhibit a double well character.18 For PT in bulk water, donor and acceptor water molecules form an intermediate Zundel-ion state, with the transferring proton delocalized between two neighboring water molecules. At the hydrogen-bond distance that corresponds to the formation of the Zundel-ion, the double well potential is symmetrical and the barrier in the proton free energy vanishes. For proton transfer at an acid-functionalized interface, the potential barrier for PT is significantly higher, being related to the strong localization of water molecules around ionic SGs. It is however expected that decreasing the distance between ionic groups can lower the proton transfer barrier markedly. Upon densification of SGs, a critical SG separation has been found in experiment, at which a transition to highly efficient surface PT was observed. The experimental value of the critical SG separation is B7 Å.11,15,19
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24099
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
Paper
Driven by the interest in biochemical phenomena, experimental monolayer studies have been conducted mainly with surface groups (SGs) that contain carboxylic acid head groups. Recently, Fujita et al.20 have measured proton conductivity of mesoporous silica films functionalized with sulfonic acid groups. At higher interfacial acid density than that of Nafion, the proton conductivity was measured as 5.4 103 S cm1 at a relative humidity of 20% and the activation energy for PT was determined as 0.26 eV. Ilhan and Spohr carried out ab initio molecular dynamics simulations on proton networks in polymer electrolyte pores modeled using CF3CF3 and CF3SO3H groups.21 They compared proton and water distribution in two pore channels at different hydration levels. They observed high proton mobility in the case of a narrower pore, in which the distance between acid head groups was between 6 and 8.5 Å. Despite a number of studies at interfaces comprised of carboxylic and sulfonic acid head groups, the mechanism of proton transport is not understood yet. Theoretical approaches such as soliton models have been employed to explain a concerted motion of protons along the surface of a one-dimensional hydrogen-bonded system, resulting in high surface proton mobility.22–24 Water molecules bind strongly to the surface, forming hydronium ions whose motion along the surface results in proton transport. Dislocation of one hydronium causes dislocation of a neighboring hydronium ion. Solitary waves of such dislocations travel through the system, leaving a trace, because each translocation changes the local configuration of SGs. Therefore, in a one-dimensional system, the passing of a protonic soliton should be followed by a reorientation of SGs so that another solitary wave can pass through the system; such a requirement results in low proton conductivity. As a consequence of this restriction, the onedimensional soliton approach was essentially abandoned. However, the problem can be overcome by considering a two-dimensional system25 since it allows multiple paths for a soliton to pass through the system and restore the initial configuration. Our research on this topic focuses on ordered interfacial systems of SGs with sulfonic acid head groups. Previous ab initio simulation studies on this system have explored the impact of molecular structure and packing density of SGs on spontaneous ordering, acid dissociation, water binding, and hydronium ion dynamics at interfaces.26–29 The main finding of structural
PCCP
studies was a transition to a condensed surface state at a critical value of the SG density. In the condensed state, completely dissociated sulfonate anions and hydronium ions form a hexagonally ordered ionic crystal in 2D. The saturated number of H bonds renders the interface superhydrophobic. In ref. 29, we have employed the metadynamics method to study interfacial hydronium transitions. Proton transitions were evaluated at an SG spacing of 6.6 Å, which corresponds to an area of 37.5 Å2 per SG. We observed a strong energetic preference for a highly concerted proton transfer. Based on these findings, we have developed a soliton theory of concerted interfacial proton motion in this system.25,30 All of these studies on interfacial hydronium ion dynamics have been performed at a uniform surface group density, which corresponds to the critical value for the transition to the superhydrophobic ordered surface state. At interfaces in real PEMs, acid-terminated sidechains are not uniformly grafted on the pore walls formed by ionomer aggregates. Moreover, the flexibility of sidechains contributes to fluctuations in sidechain separations. It is expected that the interfacial dynamics of hydronium ions and SG will be influenced by these fluctuations in SG density. In this work, we have extended the calculations to different SG densities to understand the impact of density variations on mechanisms and energetics of interfacial hydronium ion motions. Section 2 presents the model system and reviews the main results obtained in ref. 29, including the importance of the SG spacing (dCC) of 6.6 Å. Section 3 provides computational details of the simulations. Section 4 presents the analysis of hydronium ion transitions in terms of free energy profiles, hydrogen bond strength, and elastic constants of the coupled system of hydronium ions and anionic SGs found from the frequency spectra.
2. Model system The model system is depicted in Fig. 1. It consists of a densely packed array of CF3SO3H (triflic acid) surface groups. Terminal C atoms of all SGs are fixed at hexagonal lattice positions. This is the system that was evaluated in previous studies of stable interfacial configurations and wettability.26,28 The shortest
Fig. 1 Model system consisting of 2D arrays of completely ionized CF3SO3H surface groups (SG) with one water molecule added per SG. F, S, O, H and C atoms are shown by yellow, red, white, orange and blue circles, respectively. An upright structure is shown on the left and a tilted structure is shown on the right side of the figure. The dotted lines indicate the unit cell used in the simulations in which C atoms are fixed at their positions. Donor, acceptor and spectator SG are labeled with D, A and S respectively and the distances that define the collective variable for the lateral shift of an H3O+ ion are shown, CV = d12 d23. On the left, red arrows indicate hydronium ion transitions that will transform the tilted structure into the upright structure.
24100 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
Paper
3. Computational details
concerted proton translocations, we considered a unit cell consisting of 12 SGs and 12 water molecules with periodic boundary conditions as indicated by dashed lines in Fig. 1. The unit cell area for SG separations of 6.4, 6.6 and 6.8 Å was 424, 456 and 482 Å2 respectively. Ab initio metadynamics utilizes a repulsive history-dependent potential in the form of Gaussian functions to achieve an efficient sampling of the configuration space. The method allows describing the system by a small number of collective variables (CV); after numerous trials, we concluded to use one CV. During simulations, the reduced phase space of the CV is explored and Gaussian potentials are added at points visited by the system. The accumulation of Gaussian functions discourages the system from revisiting points already visited. The calculation is completed when the system reaches a diffusive state, in which the addition of Gaussian functions has evened out the free energy surface. The ensemble of Gaussian functions is then subtracted revealing a positive image of the free energy landscape. The sampling efficiency and the precision of the simulation depend on the height and width of the Gaussian functions and hence these parameters must be optimized. Adding Gaussian functions to the system adds energy to the degrees of freedom associated with CVs. This may lead to a significant inhomogeneity in the temperature distribution of the system, and may also affect the dynamics. To address this issue, a harmonic coupling between the real and fictitious CV is introduced in the metadynamics Lagrangian, which ensures uniform sampling of configurations in well regions of the free energy landscape.31 The CP2K package uses a mixed Gaussian and plane wave basis set.34 Valence electrons are represented by double-z augmented Gaussian basis sets (DZVP-MOLOPT).35 The pseudopotential of Goedecker, Teter and Hutter (GTH) represents core electrons.36 We used an energy cut-off of plane-wave expansions of 300 Ry (1 Ry E 13.6 eV). Exchange and correlation energies were computed within the GGA approximation using the BLYP functional.37 For every time step, the electronic structure was relaxed to an accuracy of 107 Hartree (1 Hartree E 27.2 eV). Optimized structures were thermalized for B3 ps. The temperature was set to 300 K using a Nose–Hoover thermostat and the timestep was set to 0.5 fs. In the metadynamics Lagrangian, we used a coupling constant k = 0.4 au (1 au E 1.55 103 N m1) and a fictitious particle mass M = 75 au (1 au E 1.67 1027 kg). Gaussian potentials with height h = 0.013 eV and width Dd = 0.014, 0.02, 0.016 Å for dCC = 6.4, 6.6 and 6.8 Å, respectively, were used. The height of Gaussians was chosen such that it is large enough to optimize the computational cost and small enough to conserve the total energy of the system and the width was calculated as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dd ¼ hs2 i hsi2 ; (1)
We employ the ab initio metadynamics method of Laio and Parrinello31,32 to simulate free energy surfaces and proton transition pathways in the highly charged interfacial system.29 Metadynamics at the DFT level is incorporated in the CP2K package.33 In order to explore the possibility of spontaneous
where s is the value of the collective variable. The CV was defined as the lateral shift of the transferring H3O+ ion, dCV = d12 d23, as illustrated in Fig. 1. Addition of consecutive Gaussians was done when the displacement of CV relative to the previous state reached 3Dd/2.
possible surface group is considered for computational convenience. In previous work, longer surface groups have been evaluated as well.28 These studies have revealed that a restriction to triflic acid does not pose any major restriction in terms of interfacial structure and dynamics. The model system is evaluated in the limit of high SG density. The interfacial area per surface group is in the range of 35–40 Å2 corresponding to nearest neighbor SG separation distances of 6.4 to 6.8 Å. A vacuum layer of 15 Å in height was introduced in the z-direction. One water molecule is added per group, representing conditions of minimal hydration. In spite of the small amount of water added, anionic head groups are found in a completely dissociated state owing to the high SG density. In ref. 26 it was found that at such SG separations there are two stable structures, upright and tilted, as shown in Fig. 1. Transitions between these structures occur by elementary hydronium ion transitions, as illustrated in Fig. 1 with an arrow. A sequence of such elementary transitions leads to hydronium transport along the surface. In ref. 25 we discussed what particular sequence of elementary transitions would facilitate hydronium transport most efficiently. The present paper focuses on the energetics of elementary interfacial proton transitions. We started this consideration in ref. 29, where a single transition of a hydronium ion was considered at a SG separation of 6.6 A. This separation corresponds to the case, where the energies of the upright and tilted structures coincide. Therefore, this case might be most favorable for the hydronium ion mobility because back and forth transitions between the structures are isoenergetic. On the other hand, in real systems SGs are never grafted precisely at regular lattice positions and as a result the distance between different SGs varies. This article explores the influence of a SG separation change on the energetics of a hydronium ion hop. We investigate what happens when the energy of one of the structures is larger. At 6.4 Å, the energy of the tilted structure is larger than upright structure by 0.7 eV. At 6.8 Å, the energy of the upright structure is larger than the tilted structure by 0.8 eV. At SG separations further away from the critical value, the energy difference between upright and tilted structures increases substantially, thereby rendering transformations between the two structures energetically very expensive. Also, the interfacial array will become strongly hydrophilic at larger SG separations. Mechanisms of surface hydronium ion transport, as considered in our work, will be stifled by surface to bulk transfer of protons followed by bulk-like proton motion. For the dependence of energies of the interfacial structures on the SG separation see ref. 26.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24101
View Article Online
Paper
PCCP
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
We obtained frequency spectra by taking the Fourier transform of the velocity autocorrelation function, gvv(t), calculated from the metadynamics trajectory. The function gvv(t) was calculated as an average of scalar products of atomic velocities, [next-line-equation], calculated as gvv ðtÞ ¼
hvi ð0Þ vi ðtÞi ; hvi ð0Þ vi ð0Þi
(2)
where vi(0) is the initial velocity of atom i and vi(t) is its velocity at time t. The average was calculated over various metadynamics steps chosen for vi(0) and vi(t). A Fourier transformation of this function was then performed to obtain the frequency spectrum, ð1 IðoÞ ¼ eiot gvv ðtÞdt (3) 0
The function I(o) gives information about specific frequency modes. The spectra were obtained for the motion of all atoms as well as specific modes, namely, stretching of O–H bonds and tilting of SGs.
4. Results and discussion In ref. 29, we applied the metadynamics method to study two cases: (1) collective transfer of hydronium ions and (2) disconcerted transfer of hydronium ions. In the collective mechanism, a transition of one hydronium ion from a donor to an acceptor SG triggers the transition of the neighboring hydronium ion from an acceptor SG to the next acceptor SG and so forth. The number of hydrogen bonds per SG is conserved in this mechanism. The disconcerted mechanism, on the other hand, creates a hydrogen bond defect with the formation of 4 hydrogen bonds at the acceptor SG and only 2 hydrogen bonds at the donor SG. In these calculations, we obtained barrier energies of B0.3 eV for the collective transfer and B0.6 eV for a disconcerted transfer of a single hydronium ion. These calculations showed that the collective mechanism of hydronium ion transfer is energetically favorable over the disconcerted mechanism.
A sequence of such hydronium relocation steps can lead to long-range transport of protons in the system. In the work presented here we applied the metadynamics method with a single CV to study hydronium ion transfer at dCC = 6.4, 6.8 Å. At the shortest SG separation distance of 6.4 Å, the upright structure is lower in energy than the tilted structure by 0.7 eV and at 6.8 Å, the tilted structure is more stable than the upright structure by 0.8 eV, which facilitates the proton transition. The metadynamics results at these two SG separations are compared with those obtained previously at dCC = 6.6 Å. Hydronium ion transitions were completed in 100, 115 and 40 ps at dCC = 6.4, 6.6 and 6.8 Å, respectively. The reconstructed Helmholtz free energy profiles for dCC = 6.4 and 6.6 are shown in Fig. 2 together with snapshots along the trajectory. The activation energy was calculated as 0.65, 0.6 and 0.25 eV, respectively, for dCC = 6.4, 6.6 and 6.8 Å. The hydronium ion transitions at dCC = 6.4 and 6.6 Å, shown in Fig. 2, give rise to the formation of a hydrogen bond defect, where the donor SG becomes deficient in hydrogen bonds and the acceptor has 4 hydrogen bonds. The structures formed along the trajectories are similar for the SG spacings of 6.4 and 6.6 Å. The free energy profile for dCC = 6.8 Å is shown in Fig. 3. At this SG spacing we observed a spontaneous transition of another hydronium ion. The hydronium ion labeled as 3 (see Fig. 3) underwent flipping following the transition of hydronium 1. However, it reverted back and hydronium 2 underwent flipping. The transition of hydronium 2 was accompanied by reorientation of the donor SG. The O–O distances between the oxygen nucleus of the transferring hydronium ion (OT) and oxygen nuclei of spectator (OS), donor (OD) and acceptor (OA) SGs, respectively, are shown in Fig. 4 at SG separations of 6.4, 6.6 and 6.8 Å. Red and black curves correspond to the OT–OS distance, the blue curve corresponds to the OT–OD distance and the green curve shows the OT–OA distance. The OT–OS distances remain in the range of 2.5 to 2.7 Å. As expected, the OT–OD distance increases and the OT–OA distance decreases.
Fig. 2 Reconstructed free energy surfaces (FES) for local defect-type H3O+ ion transitions at SG spacing of 6.4 Å (left panel) and 6.6 Å (right panel). Snapshots show structures at points A, B, C, D marked on the FES. Since similar structures are formed at both the SG densities, only one set of structures at the SG separation of 6.4 Å is shown.
24102 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
Paper
Fig. 3 Reconstructed free energy surfaces (FES) for H3O+ ion transitions at SG spacing of 6.8 Å. Snapshots show structures at points A, B, C, D marked on the FES. The transition of ion 1 is followed by the transition of ion 3, which then reverts back and ion 3 undergoes flipping. This concerted transition of hydronium ions 1 and 2 conserves the number of hydrogen bonds per SG.
Fig. 4 O–O distances as a function of simulation time for SG separations of 6.4, 6.6 and 6.8 Å. Black and red curves represent the distance between O atom of spectator SGs and O atom of a hydronium ion undergoing transition. The blue curve represents the distance between O atom of a donor SG and hydronium O and the green curve represents the distance between an acceptor SG and the transferring hydronium ion. Values are averaged over intervals of 1000 time steps to reduce the noise.
With only a slight decrease in SG density, the hydronium ion transport mechanism changes from disconcerted to collective with significantly lower activation energy. This effect is believed to be a consequence of the fact that at dCC = 6.8 Å the final tilted structure is energetically preferred over the initial upright structure. To analyze changes in dynamic mechanism with SG density, we calculated frequency spectra of the system, obtained from Fourier transformation of the velocity autocorrelation function, as depicted in Fig. 5. In order to be able to assign the frequency
peaks, we calculated separately the stretching frequencies of hydronium OH bonds at different SG separations, shown in Fig. 6, as well as tilting frequencies of SGs, shown in Fig. 8. We found 3 O–H stretching frequency bands in the 100–4000 cm1 range: the most intense band occurs in the 130–1330 cm1 range with a peak at around 675 cm1, a lower intensity band lies between 1450–2260 cm1 with a peak at around 1720 cm1 and the lowest intensity band is observed in the 2450–3600 cm1 range. The frequency spectra for all three SG separations are
Fig. 5 Full frequency spectra at different SG separations obtained from Fourier transform of the velocity autocorrelation function.
Fig. 6 Stretching frequency spectra of hydronium O–H bonds at different SG separations.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24103
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
Paper
PCCP
overlapping almost completely indicating that the variation in SG density exerts only a small effect on the stretching frequencies of O–H bonds. This is an important finding in this work since it shows that in this range of SG separations, the length and strength of hydrogen bonds is essentially unaltered. Buzzoni et al.38 studied IR and Raman spectra of water molecules interacting with Nafion. In order to understand the basic features of the IR spectra and assign absorption peaks, they proposed an empirical model for a solvated Zundel ion. They categorized hydrogen bonds into ‘internal’, ‘external’ and ‘solvated’ bonds based on O–O distance. For internal hydrogen bonds with O–O distance in the range of 2.45–2.57 Å, OH–O stretching frequencies are expected to be less than 2000 cm1. The external hydrogen bonds with O–O distance in the range of 2.6–2.7 Å are expected to show absorption within the range from 2500 to 3200 cm1. Finally for solvated hydrogen bonds, for which the O–O distance is larger than 2.7 Å, the frequencies are expected in the range of 3200 to 3600 cm1. We have observed the peaks roughly in similar frequency ranges. Comparing the relative intensities of the peaks we can say that the O–O distances mostly fall in the range of ‘internal’ H-bonds. To understand the effect of bond length fluctuations on stretching frequencies, we calculated the average O–O distances of different hydrogen bonds and the range of fluctuations around the average value. We divided hydrogen bonds into 3 categories: (1) donor hydrogen bonds, (2) hydrogen bonds of spectator SGs directly bonded to the transferring hydronium ion, and (3) hydrogen bonds of SGs farther away from the transferring hydronium ion. The range of average O–O distances of various H-bonds is given in Table 1. The time averaged O–O distances of different H-bonds lie in the range of 2.57–2.7 Å, which corresponds to a 2.5% variance. However, this dispersion might not be relevant because time-fluctuations of the O–O distance of any single H-bond is around 5%. All O–O distances in the system remain fairly constant upon densification of the SG array. This is achieved by a relaxation phenomenon where hydronium ions shift perpendicular to the plane spanned by SGs. Fluctuations in the height of the hydronium oxygen with respect to the plane drawn through S atoms of SGs is shown in Fig. 7. It shows that at shorter SG separation, the hydronium oxygen remains predominantly above the plane, whereas at longer SG separation, hydronium ions descend closer to the plane. We also calculated SG frequency spectra at 3 SG separations as depicted in Fig. 8. As shown in the figure, the SG separation
Fig. 7 Fluctuations in the height of O atom of a spectator hydronium ion. Values are averaged over intervals of 1000 time steps to reduce the noise.
Fig. 8 SG frequency spectra calculated at 3 SG separations.
in this small range does not alter the frequencies of the peaks. To interpret the spectra, let us consider small oscillations of SGs and hydronium ions about their equilibrium positions. As discussed in ref. 25, we consider only pair interactions in the nearest neighbor approximation. In general, in a two-component system, there are three distinct pairs: hydronium ion–hydronium ion, SG–SG, and hydronium ion–SG. The corresponding elastic constants are KHH, KSS, and KHS. These constants appear explicitly in theoretical models, see ref. 23–25. Additionally, SGs are fixed at the basal plane, which gives rise to the force returning SGs to their equilibrium position, described by an elastic constant KF. For such a system a standard approach for diatomic crystal lattices can be employed, with an additional term describing the flexibility of SGs. A schematic of the model is shown in Fig. 9.
Table 1 Influence of the hydronium ion transition on bond lengths of surrounding groups. The ranges of average O–O distances at different SG densities are given in the table. The nearest spectator groups are the SGs directly bonded to the transferring hydronium ion; the SGs further away from the transferring hydronium ion are designated as far spectators in the table
dCC (Å)
Range of average O–O distances (Å)
6.4
6.6
6.8
Donor Nearest spectators Far spectators
2.59 2.61–2.62 2.60–2.62
2.57 2.64–2.65 2.61–2.62
2.64 2.64–2.70 2.61–2.64
24104 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
Fig. 9 A schematic representation of SGs (black circles) connected to hydronium ions (red circles) by springs with spring constant KHS. The coupling constant between SGs is KSS and that between hydronium ions is KHH. Finally, KF describes the flexibility of SGs, which tends to return them to their equilibrium positions.
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
Paper
For simplicity, we will consider only a one-dimensional string of masses connected with springs. Such a 1D system captures essential features of the 2D lattice due to the hexagonal symmetry of the lattice. Because we consider only small oscillations of hydronium ions, each hydronium and each SG is subjected to a circularly symmetric potential, which is a function of the displacement from the equilibrium position. In 1D approximation, the displacement of the j-th atom is u2j = BH exp(i(2jkdCC ot)),
(4)
u2j+1 = BS exp[i((2j + 1)kdCC ot)],
(5)
where BH and BS are unknown amplitudes and k = 2p/l is a wave vector. Since the problem is 1D, we do not use vector notations. The value of l is defined by the boundary conditions in our simulations, l = 6 dCC/O3. According to Hook’s law, the force acting on the j-th atom is F2j = KHS(u2j+1 + u2j1 2u2j),
(6)
F2j+1 = KHS(u2j+2 + u2j 2u2j+1) KFu2j+1.
(7)
Inserting these equations into Newton’s second law and eliminating BH and BS, we get 1 1 KF þ þ o1;22 ¼ KHS mH mS 2mS KHS sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 1 1 KF 2KF 4 sin2 kdCC KHS þ þ mH mS 2mS KHS mH mS KHS mH mS (8) where mH and mS represent the masses of hydronium ions and SGs, respectively. The signs correspond to optical and acoustic modes. The third peak in Fig. 8 corresponds to SG–SG interactions. It can be modeled by a string of particles with the same mass connected by springs. The mathematical formalism is similar to the one shown above. The resulting frequency is 1 kdCC o32 ¼ þ KF : 4KSS sin2 (9) mS 2 The argument of the sin is different from the one in eqn (8) since the distance between SGs is twice as large as that between SGs and hydronium ions. To assign the frequencies o1,2,3 to the peaks in Fig. 8, one has to use the corresponding spectrum of hydronium ion oscillations. The frequencies o1,2 should be the same for SG and hydronium ions. The third peak, corresponding to oscillations in the hydronium subsystem only, is 1 kdCC 0 o3 2 ¼ ; (10) 4KHH sin2 mH 2 which is analogous to eqn (9). Unfortunately, it is impossible to single out the peaks in the hydronium spectrum because peaks are broad and merged together as shown in Fig. 10. However, it is clear that there are no peaks at frequencies above 700 cm1. Therefore, the peak observed at 1000 cm1 in Fig. 8 corresponds to o3. Hence, the leftmost peak is an acoustic mode, and the middle peak is the optical mode.
This journal is © the Owner Societies 2014
Fig. 10 Frequency spectrum of a hydronium ion.
The broadness of the peaks for hydronium ions, in contrast to well-distinguished peaks for SGs, can be explained based on the mechanism of energy loss suggested in ref. 30. When a hydronium ion or an SG shifts, it pulls on its nearest neighbors in transverse direction, which transfers energy to the remaining lattice. A heavy SG can easily pull a hydronium ion, whereas it is difficult for a light hydronium ion to pull a SG, which is fixed at the basal plane. Therefore, hydronium ions lose energy faster than SGs, which broadens the hydronium ion peaks. Since we can extract values of the characteristic mode frequencies from Fig. 8 and 10, eqn (8) and (9) can be solved for the elastic constants. Using o1 = 200 cm1, o2 = 625 cm1, o3 = 995 cm1, o30 = 700 cm1, mH = 19 amu, mS = 80 amu, k = 2.75 109 m1 and dCC = 6.6 1010 m, we obtain KF = 1808.68 N m1, KHS = 22.58 N m1 and KSS = 1165.02 N m1, KHH = 222.6 N m1. The comparatively small value of KHS can be explained based on the relaxation mode of hydronium ions relative to the plane spanned by SGs, discussed above. To accommodate a change of dCC, hydronium ions tend to move away from the basal plane when SGs are squeezed together and approach the basal plane on increasing the SG separation. This relaxation of hydronium ions perpendicular to the plane of SGs keeps the H-bond length relatively constant and softens the coupling between hydronium ions and SGs. The strength of coupling between hydronium ions plays a pivotal role for the collectiveness of the hydronium ion transitions. In ref. 25 and 30, we have developed a soliton theory of collective hydronium ion motion at the minimally hydrated interface of triflic acid SGs. In the fully dissociated condensed state of the interfacial system,26 the approach distinguishes two sublattices: the sublattice of sulfonate anions that fluctuate about their equilibrium positions and the sublattice of mobile hydronium ions. The two sublattices are coupled by strong hydrogen bonds. In the genuine soliton Hamiltonian,22–24 discussed in ref. 25, harmonic coupling terms represent interactions in the system of SGs and hydronium ions. To simplify the mathematical formalism, interaction terms for SG–SG coupling and SG–hydronium ion coupling are replaced by an effective potential, whereas hydronium–hydronium coupling is still treated explicitly in harmonic approximation. The resulting Hamiltonian for single file soliton motion, illustrated in Fig. 11, is X m k H¼ (11) u_ i2 þ ðuiþ1 ui Þ2 þV ðui Þ 2 2 i
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24105
View Article Online
Paper
PCCP
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
system. In the future, further theoretical studies and molecular simulations will be performed specifically to understand mechanisms of soliton creation and the energetic cost associated with it.
5. Conclusions Fig. 11 Lattice configurations of the 2D two-component lattice of H3O+ and SO3 ions. The native 2D crystal structure in (a) is shown at the critical SG density of dCC E 6.6 Å. In the introduced symbolic representation, blue triangles indicate hydronium positions and red dots represent positions of sulfonate anions. The structure for a single file hydronium ion transfer, which conserves the number of H-bonds per SG, is shown in (b).
where ui is the displacement of the ith hydronium ion, m is the mass of a hydronium ion, k is the coupling constant between neighboring hydronium ions and V(ui) is the effective potential. The Hamiltonian equations of motion of this system were solved analytically in ref. 25 for different potential functions, V(ui), in the continuum limit, assuming that the soliton width is much larger than dCC. The size of the soliton is given by sffiffiffiffiffiffiffi ka2 N0 ¼ (12) 2e It thus depends on hydronium–hydronium coupling constant, k, the hydronium relocation distance, a and the depth of the potential well, e. The effective potential V(ui) in the soliton Hamiltonian can be parameterized using free energy curves obtained from Metadynamics simulations, as reported in ref. 29. The main parameters, which define this potential and thus the static properties of solitons, are the hydronium relocation distance, a, and the depth of the potential well, e. The main parameter for the soliton dynamics is the hydronium–hydronium coupling constant, k. Using estimates for the coupling constants and scaling arguments, an explicit expression for the soliton mobility was derived in ref. 30 and a rough estimate of the mobility value was obtained. Except for the coupling constant k, values for other parameters were obtained from Metadynamics simulations. The value of the coupling constant was estimated to be B2000 N m1, based on eqn (12) and assuming that the soliton consists of 50 hydronium ions. The coupling constant KHH reported here is derived from frequency spectra, which consider only small displacements of hydronium ions, treated in harmonic approximation, whereas for soliton-like transport the coupling constant k is associated with large displacement of hydronium ions. In the limit of weak coupling, KHH would be much smaller than k and hence KHH would not provide a reasonable estimate of the soliton size. However, in the case of strong coupling, KHH could be used in eqn (12) to estimate the size of solitons. Thus, based on the KHH value obtained here, the static soliton would involve 20–25 hydronium ions. This size of soliton justifies a posteriori the continuum approximation used in the soliton theory in ref. 25 and 30. It implies a strong coupling of hydronium ions in our
24106 | Phys. Chem. Chem. Phys., 2014, 16, 24099--24107
We have performed ab initio metadynamics simulation to determine an activation barrier for the transition of hydronium ions at different SG spacing. We observed that at a SG separation of 6.4 and 6.6 Å, the energy barrier for the hydronium ion transition is high and gives rise to the formation of a hydrogen bond defect whereas at a SG separation of 6.8 Å, a spontaneous transition of a second hydronium ion remedies the hydrogen bond defect and conserves the number of hydrogen bonds in the unit cell. The activation energy in this case (0.25 eV) was found to be much lower than that obtained for the other transitions that result in the formation of hydrogen bond defects (Z0.6 eV). In spite of this change in activation energy, the stretching frequency spectra showed that in this range of SG separations, hydrogen bond strength remains essentially unaltered. Coupling constants between surface groups and hydronium ions were calculated from frequency spectra. They represent the coupling strength of the system for small fluctuations relative to the equilibrium configuration. The coupling constant between neighboring hydronium ions has been employed to estimate the size of collective, soliton-like excitations of hydronium ions. We determined that a soliton include about 20–25 hydronium ions, which would be characteristic of a strongly coupled system. We believe that the insensitivity of frequency spectra, hydrogen bond strength and coupling constants to the density of surface groups is owed to a peculiar relaxation motion of hydronium ions perpendicular to the plane of SGs.
References 1 O. Markovitch and N. Agmon, J. Phys. Chem. A, 2007, 111, 2253. 2 O. Markovitch, H. Chen, S. Izvekov, F. Parsani, G. A. Voth and N. Agmon, J. Phys. Chem. B, 2008, 112, 9456. 3 D. Marx, M. E. Tuckerman, J. Hutter and M. Parrinello, Nature, 1999, 397, 601. 4 M. E. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Phys. Chem., 1995, 99, 5749. 5 M. E. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Chem. Phys., 1995, 103, 150. 6 A. Y. Mulkidjanian, J. Heberle and D. A. Cherepanov, Biochim. Biophys. Acta, 2006, 1757, 913. 7 H. Morgan, D. M. Taylor and O. N. Oliviera Jr., Biochim. Biophys. Acta, 1991, 1062, 149. 8 B. Gabriel and J. Teissie, Proc. Natl. Acad. Sci. U. S. A., 1996, 93, 14521. 9 D. Marx, M. E. Tuckermann and M. J. Parrinello, J. Phys.: Condens. Matter, 2000, 12, A153.
This journal is © the Owner Societies 2014
View Article Online
Published on 29 September 2014. Downloaded by Christian Albrechts Universitat zu Kiel on 28/10/2014 13:11:08.
PCCP
10 J. Heberle, G. Riesle, D. Thiedemann, D. Oesterhelt and N. A. Dencher, Nature, 1994, 370, 379. 11 J. Zhang and P. R. Unwin, J. Am. Chem. Soc., 2002, 124, 2379. 12 S. Serowy, S. M. Saparov, Y. N. Antonenko, W. Kozlovsky, V. Hagen and P. Pohl, Biophys. J., 2003, 84, 1031. 13 I. Sakurai and Y. Kawamura, Biochim. Biophys. Acta, 1987, 904, 405. 14 K. Leberle, I. Kempf and G. Zundel, Biophys. J., 1989, 55, 637. 15 V. B. P. Leite, A. Cavalli and O. N. Oliviera Jr, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 6835. 16 J. Teissie, M. Prats, S. Soucaille and J. F. Tocanne, Proc. Natl. Acad. Sci. U. S. A., 1985, 82, 3217. 17 A. Polle and W. Junge, Biophys. J., 1989, 56, 27. 18 L. I. Krishtalik, Biochim. Biophys. Acta, 2000, 1458, 6. 19 C. J. Slevin and P. R. Unwin, J. Am. Chem. Soc., 2000, 122, 2507. 20 S. Fujita, A. Koiwai, M. Kawasumi and S. Inagaki, Chem. Mater., 2013, 25, 1584. 21 M. A. Ilhan and E. Spohr, J. Phys.: Condens. Matter, 2011, 23, 234104. 22 S. Pnevmatikos, Phys. Rev. Lett., 1988, 60, 1534. 23 A. V. Zolotaryuk, S. Pnevmatikos and A. V. Savin, Phys. Rev. Lett., 1991, 67, 707. 24 M. Peyrard, S. Pnevmatikos and N. Flytzanis, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 903.
This journal is © the Owner Societies 2014
Paper
25 A. Golovnev and M. Eikerling, J. Phys.: Condens. Matter, 2013, 25, 045010. 26 A. Roudgar, S. P. Narasimachary and M. Eikerling, J. Phys. Chem. B, 2006, 110, 20469. 27 A. Roudgar, S. P. Narasimachary and M. Eikerling, Chem. Phys. Lett., 2008, 457, 337. 28 S. P. Narasimachary, A. Roudgar and M. Eikerling, Electrochim. Acta, 2008, 53, 6920. 29 S. Vartak, A. Roudgar, A. Golovnev and M. Eikerling, J. Phys. Chem. B, 2013, 117, 583. 30 A. Golovnev and M. Eikerling, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 062908. 31 A. Laio and F. L. Gervasio, Rep. Prog. Phys., 2008, 71, 126601. 32 A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562. 33 http://www.cp2k.org. 34 J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167(2), 103. 35 J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105. 36 S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703. 37 A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098. 38 R. Buzzoni, S. Bordiga, G. Ricchiardi, G. Spoto and A. Zecchina, J. Phys. Chem., 1995, 99, 11937.
Phys. Chem. Chem. Phys., 2014, 16, 24099--24107 | 24107
