Article pubs.acs.org/JPCB
Ab Initio Molecular Dynamics Study of the Aqueous HOO− Ion Zhonghua Ma,*,† David Anick,*,‡ and Mark E. Tuckerman*,§ †
Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, United States Laboratory for Water and Surface Studies, Department of Chemistry, Tufts University, Medford, Massachusetts 02155, United States § Department of Chemistry, Courant Institution of Mathematical Sciences, New York University, New York, New York 10003, United States ‡
ABSTRACT: The properties of the hydroperoxide anion, HOO−, in water play a key role in many biological systems and industrial processes. However the dynamics of HOO− and its solvation shell are largely unknown. We have undertaken an ab initio molecular dynamics study of aqueous HOO− at ambient temperature in liquid water. Two solvation structures for the hydroperoxide anion account for 90% or more of the configurations in a 25 ps NVT run at 300 K: these have four hydrogen bond donors to the terminal oxygen atom of HOO− and either one or two hydrogen bond donors to its middle oxygen atom. The H of HOO− is essentially always a donor in an H-bond to a water molecule. Two structures with three donors to the terminal O and either one or two at the middle O are also important. A set of five NVE runs totaling 74 ps found considerable variability in the proportions of time spent in each type of solvation pattern. Mean lifetimes of these patterns ranged from 54 to 109 fs, after which the complexes were observed to transform into different, sometimes less favorable structures. Analysis of the electronic structure associated with different solvation patterns indicates that a traditional Lewis-type picture of hydrogen bonding at the middle oxygen and non-Lewis behavior at the terminal oxygen coexist in aqueous HOO−. The non-Lewis character of the terminal oxygen is compared to similar observations of the oxygen in the hydrated hydroxide ion from previous ab initio molecular dynamics of OH−(aq). corresponding tendencies for (HOOH)(OH−)(H2O)n−1 clusters.24 Let us designate the atoms of HOO− as follows: the H is h1, o1 is covalently bonded to h1, and o2 is the other O. We adopt here the solvation pattern notation “mAnA” to indicate that in an (HOO−)(H2O)n cluster or solvation shell, o1 is an m-fold acceptor and o2 is an n-fold acceptor, regardless of whether h1 is a donor. The term “complex” will generally be used to denote the structure corresponding to a particular pattern. Expressed using this nomenclature, the findings of Anick24 included (a) a propensity for HOO− to occupy a central position in a (HOO−)(H2O)n cluster with solvation on all sides, (b) an optimal solvation pattern of 2A4A with the second and third best being 3A4A and 1A4A, respectively, (c) because of the strength of the bonds at o2 and the unfavorability of either 3- or 5-coordination at o2, a prediction that donors to o2 would be difficult to dislodge and slow to exchange, and (d) conversely, because of the relative weakness of H-bonds at o1 and the favorability of 1A4A and 3A4A, a prediction of relatively rapid first shell exchange at o1. In addition, because optimal solvation patterns for HOOH donating to OH− and for HOO− are so different, it was reasoned that a “simple” proton transfer either way would not stick, and if a proton transferred, it would immediately transfer back. In order to endure, a proton transfer would require
1. INTRODUCTION Hydroperoxide (HOO−), the conjugate base of the weak oxoacid hydrogen peroxide (HOOH, pKa = 11.62), plays an important role in atmospheric chemistry,1 cellular biology,2 and numerous industrial processes.3 For example, the aqueous reaction of HOO− with SO2 is a step in the formation of acid rain.4,5 Immersion in a hot concentrated mixture of NH4OH + HOOH, in which HOO− is the principal oxidant, is step 1 of the “standard recipe” for cleaning semiconductors.6−8 Oxidation of lipids and proteins by HOO− and other reactive oxygen species is the principal biochemical mechanism behind aging and many degenerative diseases.9−13 Static14−19 and dynamic20 ab initio studies of aqueous HOOH have provided an understanding of how HOOH interacts with water solvent. In general, HOOH has one Hbond donor and one H-bond acceptor at each O, and the molecule as a whole is a more potent donor than acceptor. Studies that include multiple HOOH molecules21 also provide insight into the fundamental HOOH−HOOH interaction. Dynamical studies of the neutral hydroperoxy radical (HOO·) in water have found that this neutral radical tends to engage in just two H-bonds with the surrounding water: one in which H of HOO· is the donor to an H2O and one in which the terminal O of HOO· is the acceptor from an H2O.22,23 However, to our knowledge, no ab initio molecular dynamics (AIMD) studies of aqueous hydroperoxide anion have yet been carried out. In a comprehensive study of gas phase (HOO−)(H2O)n clusters, for n = 1−8 and n = 20, Anick previously noted some interesting trends and patterns and contrasted these with the © 2014 American Chemical Society
Special Issue: James L. Skinner Festschrift Received: January 23, 2014 Published: March 17, 2014 7937
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
3. RESULTS AND DISCUSSION 3.1. Hydrogen Bond Analysis. Understanding the hydrogen bond network is crucial in identifying the solvation structure of aqueous HOO−. We employ two distinct hydrogen bond (HB) definitions in our analysis. The geometric HB definition is based on an optimal geometric criterion extracted from oxygen−oxygen and oxygen−hydrogen potentials of mean force.45 This analysis yields the following HB criteria: RHd···Oa < 2.5, ROd···Oa < 3.4 Å, ∠HdOdOa < 30°, where “d” and “a” denote donor and acceptor, respectively. The same HB criteria are used for HBs donated or accepted by HOO−. The topological HB definition is the adaptation to our system of the definition given by Henchman and Irudayam:46 Od−Hd···Oa is considered to be an HB if Oa is the closest O to Hd other than Od, and RHd···Oa is the shortest H···O distance between the two molecules (which may now be H2O or HOO−). Table 1 shows
simultaneous (or prior, cf. ref 25) stabilization by substantial rearrangement of the entire first shell. Since this would be a statistically rare event, solvation actually reduces the rate of proton transfer (or equilibration) reactions, a phenomenon dubbed “solvation shell inertia”. This should be contrasted with the case of OH−(aq), which readily undergoes structural diffusion via individual proton transfer and solvation shell relaxation events.26 The aforementioned gas-phase study now sets the stage for the bulk solvation study to be presented here. Some of the motivating questions are the following: (1) What are the predominant solvation patterns for HOO− in aqueous solution, and what are their average durations/lifetimes? (2) Can AIMD verify that first-shell exchange is relatively rare at o2 but more frequent at o1? (3) What insights can be gleaned from the electronic structure of solvated HOO− that relate to its solvation behavior? We present here the combined results from six AIMD simulations of aqueous HOO−, encompassing 25 ps of NVT and 74 ps of NVE runs in all. Section 2 provides details of the ab initio molecular dynamics methodology. Section 3.1 presents a hydrogen bond analysis of the trajectories and identifies four noteworthy solvation patterns and their populations. Radial distribution functions are discussed in section 3.2, while section 3.3 examines the electronic properties of solvated HOO− that serve to explain the geometry of the major solvation complexes. Section 3.4 explores interconversions among the various solvation patterns. Conclusions and future directions are given in section 4.
Table 1. Percentage of h1, o1, and o2 Having the Specified Number of Hydrogen Bonds Obtained from the 25 ps Constant Temperature (NVT) Run Results for Geometric HB Criteria nHB 0 HB of o2 HB of o1 HB of h1
1
0 0.8 1.5 Results
2
3
0 0.3 7.8 28.8 68.0 2.2 98.5 0 0 for Topological HB Criteria
4
5
91.1 0 0
0.7 0 0
nHB
2. COMPUTATIONAL DETAILS All calculations were performed within the Car−Parrinello ab initio molecular dynamics (CPAIMD) framework.27 Electronic structure was represented within the Kohn−Sham (KS) density functional formalism using the Becke−Lee−Yang−Parr (BLYP) generalized gradient approximation for the exchangecorrelation functional and Troullier−Martins type pseudopotentials.28−30 The Grimme empirical dispersion correction was also included.31,32 The KS orbitals were expanded in a discrete variable representation (DVR) basis set.33−38 A canonical (NVT) ensemble simulation was conducted for a system of one HOO− in 31 water molecules in a periodic cubic box of length 18.635 au. A well converged DVR basis set corresponding to 75 grid points/75 basis functions along each coordinate direction was employed, which leads to fully converged energies and interatomic forces.34−39 A fictitious mass of μ = 577 au and a time step of 0.05 fs were employed in the simulation. A Nosé−Hoover chain thermostat attached to each ionic degree of freedom was used to maintain the system at 300 K.40−42 Hydrogen atoms were assigned the mass of deuterium in order to reduce the influence of nuclear quantum effects. This simulation was carried out for 25 ps, with configurations written out every 1 fs. Five microcanonical (NVE) ensemble CPAIMD simulations were performed with different initial configurations after the system was equilibrated in the NVT ensemble for 11 ps. One NVE trajectory was allowed to evolve for just over 26 ps, while the other four NVE simulations evolved for 12.125, 11.75, 12, and 12 ps. All simulations were carried out using the PINY_MD simulation package.43,44
HB of o2 HB of o1 HB of h1
0
1
2
3
4
5
0 0 0
0.2 22.1 100.0
0.3 76.5 0
4.0 1.4 0
95.0 0 0
0.5 0 0
that similar HB statistics result from the geometric and topological HB definitions for the 25 ps constant-temperature run. For the analysis of the run, we used configurations after the first 1 ps, which was taken as an equilibration phase. Using the first 5 ps as an equilibration phase does not qualitatively change the conclusions to be reported. The population of fourcoordinated o2 is 91−95%. The oxygen o1 accepts two HBs at 68−77%, and h1 donates in an HB at 99−100% during the constant-temperature simulation. In the proceeding discussion, we will refer to a “solvation complex” of HOO− as the HOO− ion together with a subset of water molecules that form either donor or acceptor HBs directly with the ion. Table 2 shows that the majority of HOO− complexes have the pattern 2A4A, followed by 1A4A, with 2A3A and 1A3A placing a distant third and fourth, respectively. The 1A4A and 2A4A patterns together comprise 90% of HOO− solvation complexes (by the geometric HB definition) or 94% (by the topological HB definition). Representative snapshots of the three most important solvation complexes are shown in Figure 1. In his analysis of the lowest energy (HOO−)(H2O)n gas phase clusters at the B3LYP and MP2 levels of theory, Anick24 similarly found that 4-fold coordination at o2 was strongly preferred. This author ranked 2A4A as having the lowest free energy using the vibrational self-consistent field method.47−49 After 2A4A, a 3A4A configuration was slightly favored over the lowest 1A4A for HOO− enclosed by the (H2O)20 dodecahedral cage. The NVT results here differ from the above findings 7938
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
hydrogen bonds comprise just 8% (geometric) or 4% (topological) of all configurations. 3.2. Radial Distribution Functions. Figures 2 and 3 show the o1-oxygen (go1O(r)), o2-oxygen (go2O(r)), o1-hydrogen
Table 2. Percentage of Oxygen Atoms Having the Specified Number of Hydrogen Bondsa Results for Geometric HB Criteria o2 2A o1 1A o1 2A o1 3A
o1 1A o1 2A o1 3A
o2 3A
o2 4A
0.8 1.9 27.6 0.2 3.9 62.0 0.1 1.8 0.3 Results for Topological HB Criteria
o2 5A 0.7 0 0
o2 2A
o2 3A
o2 4A
o2 5A
0.2 0.2 0
2.1 2.3 0.1
22.2 72.2 0.4
0 0.4 0
Column heads “o2 nA” indicate that o2 accepts “n” hydrogen bonds. Row labels “o1 nA” mean that o1 accepts “n” hydrogen bonds. a
Figure 2. o1-oxygen (go1O(r)) and o2-oxygen (go2O(r)) radial distribution functions obtained from the 25 ps NVT simulation. Here “O” is the oxygen atom of an H2O or either oxygen atom of the HOO−, i.e., it is “o2” or a water oxygen in go1O(r), and it is “o1” or a water oxygen in go2O(r).
Figure 3. o1-hydrogen (go1H(r)), o2-hydrogen (go2H(r)), and h1oxygen (gh1O(r)) radial distribution functions obtained from the 25 ps NVT simulation. Here, “h1” is the hydrogen atom of HOO−, “o1” and “o2” are as defined in Figure 1, “O” is any oxygen atom of an H2O or of the HOO−, and “H” is any hydrogen atom of an H2O or of the HOO−.
Figure 1. Snapshots from the NVT simulation. The three most common species are shown. Oxygens belonging to waters in the first solvation shell of HOO− are yellow. The atom “o1” is the oxygen atom that is covalently bonded to the hydrogen atom of HOO−, and “o2” is the terminal oxygen atom of HOO− (these definitions are also given in the text). The atoms o1 and o2 of HOO− are highlighted as blue. Other oxygen atoms are red, and hydrogen atoms are white. Green dashed lines indicate hydrogen bonds (HBs). (A) o1 accepts two HBs, o2 accepts four HBs (2A4A). (B) o1 accepts one HB, o2 accepts four HBs (1A4A). (C) o1 accepts two HBs, o2 accepts three HBs (2A3A).
(go1H(r)), o2-hydrogen (go2H(r)), and h1-oxygen (gh1O(r)) radial distribution functions (RDFs) obtained from the NVT simulation. Each RDF was calculated with a bin size of 0.03 Å. The sharp peak near 1.5 Å obviously represents the o1−o2 covalent bond and is identical for go1O(r) and go2O(r). Coordination numbers, obtained by integrating g(r) from r = 2.0 Å to the subsequent minimum, and the position and height of the second maxima for go1O(r) and go2O(r), are summarized in Table 3. The average coordination number for o1 is 2.6, whereas for o2 it is 4.2, both of which are consistent with the results of the hydrogen bonding analysis in section 3.1. As Figure 2 shows, the average distance for o2 to the oxygens in its first solvation shell is shorter than that for o1, but o1 has a much closer second solvation shell than o2. This may be due to there being one or two more water molecules in the first
mainly in that the 1A4A configuration has a significant population while 3A4A is rare. This likely indicates that the dodecahedral cage is “artificial” in that it allows a simultaneous close approach by eight water molecules whereas this occurs only rarely in a finite-temperature bulk aqueous environment. In the gas phase, 3A4A has the advantage of providing two additional HBs (compared with 1A4A) by turning the protons inward and using the HOO− as their acceptor. However, in the bulk solution this is not necessary, as there are other water layers to receive outward-pointing hydrogens from the first and second shells. The fact that the end oxygen o2 of HOO− accepts primarily four HBs is also similar to that of the oxygen of OH−(aq), which primarily exists in a 4-fold coordinated state.25,26,39,50 Structures in which o2 accepts just three
Table 3. Structural Data of the Second Peaks of go1O(r) and go2O(r)
o1 o2 7939
av coord
gmax
rmax (Å)
gmin
rmin (Å)
2.6 4.2
2.25 3.42
2.79 2.76
1.04 0.32
3.09 3.25
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
⎛ |∇n(r)|2 ⎞ ⎟ χ (r) = ⎜⎜∑ |∇ϕi(r)|2 − 4n(r) ⎟⎠ ⎝ i
solvation shell of o2 (typically 4A or 3A for o2 vs 2A or 1A for o1). Keeping four waters in the first solvation shell of o2 hampers the approach of waters from the second solvation shell of o2. Similar features can be found in go1H(r) and go2H(r). Between 1.5 Å and 2.2 Å go2H(r) has two peaks that are due to h1 and hydrogens that donate to o2. The average distance between o2 and h1 is 1.98 Å . Similarly, the peak of gh1O(r) near 1.7 Å is due to the oxygen atom of the H2O that is H-bonded to o1 via h1 and the peak near 2.0 Å is for o2. The RDFs provide a measure of the average structure of the liquid. In order to capture the detailed arrangement of HOO−, Figure 4 shows the distributions of the o1−o2 distance (ro1o2)
⎞ ⎛3 2 2/3 5/3 ⎜ (6π ) n (r)⎟ ⎠ ⎝5 (2)
where n(r) is the electron density and {ϕi(r)} are the KS orbitals. The ELF is scaled so that f(r) ∈ [0, 1], with perfect localization at r corresponding to f(r) = 1. Figure 5 shows the ELF for gas phase HOO− and for three representative snapshots from the NVT simulation. Isosurfaces
Figure 4. Distribution of o1−o2 distance (red), o1−h1 distance (blue), and angle ∠h1o1o2 (purple).
and the o1−h1 distance (ro1h1) with a bin size of 0.013 Å and the angle ∠h1o1o2 with a bin size 1.75° obtained from the NVT simulation. In Table 4, calculated geometry data from high-level optimizations of HOOH and of solvated HOO· (neutral HOO radical) and HOO− are summarized, along with the methods used.24,51,52 The ro1o2 value found for HOO− from an average of the 300 K simulation is 1.50 ± 0.05 Å, which exceeds the 0 K o1−o2 distances for solvated HOO−, isolated HOOH, and solvated HOO· by 0.02, 0.05, and 0.17 Å, respectively. The average ro1h1 value at 300 K is found to be 0.98 ± 0.04 Å, which is similar to the ro1h1 value found for solvated HOO− and HOO· clusters at 0 K. The average angle ∠h1o1o2 is found to be 101° ± 6°, which varies in a larger range relative to other reported results.24,51,52 3.3. Electron Localization. In order to understand the solvation behavior of HOO−, the spatial distribution of electronic charge was investigated. We computed the electron localization function (ELF)53,54 of both gas phase HOO− and its hydration complexes. The ELF f(r) is defined in terms of densities and takes the form f (r) = 1/(1 + χ 2 (r))
Figure 5. Isosurfaces of the electron localization function (ELF), f(r) = 0.954: (A) gas-phase HOO−; (B) 2A4A; (C) 1A4A; (D) 2A3A.
at 0.954 are plotted, and for clarity, only waters in the first solvation shell are included. The ELF of the gas phase HOO− has a ring attractor around o2 and two overlapping lone pairs situated at o1, and this same picture occurs when solvating water molecules are added. Thus, the HOO− complex can be regarded as containing two types of oxygens: o1 corresponds very closely to what would be expected from a traditional Lewis-type picture of the anion, while the ring attractor resulting from delocalization of the lone pairs around o2 indicates a significant departure from the Lewis picture. Note that the ELF of o2 is very similar to that of the oxygen of OH−(aq).50,55,56 The ELF provides an explanation for the geometry of the HOO− complexes and, in particular, for their strong preference for the 2A4A arrangement. The traditional Lewis-type picture of hydrogen bonding suggests one covalent O−H bond, two lone pairs around the middle oxygen atom, one covalent O−O bond, and three lone pairs around the end oxygen atom. On the
(1)
Table 4. Geometry Data for a Free HOOH Molecule, HOO−(H2O)6, HOO·(H2O)2, and from the NVT Simulation species
solvation
method
ro1o2 (Å)
ro1h1 (Å)
HOOH HOOH HOO·(H2O)2 HOO−(H2O)6 HOO− + 31H2O
free free 0A1A 1A4A varies
CCSD(T)/cc-pVQZ exp QCISD MP2/aug-cc-pVTZ CPMD
1.4525 1.452 1.331 1.483 1.50 ± 0.05
0.9627 0.965 0.988 0.984 0.98 ± 0.04
7940
∠h1o1o2 (deg) 99.91 100 100.1 101 ± 6
ref 51 51 52 24 this work
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
basis of this picture, one would predict that HOO− accepts three HBs via its oxygen lone pairs at the end oxygen atom o2 and two HBs at the middle oxygen atom o1 with the possibility of a donor HB the h1 site, i.e., 2A3A. For o1 and h1, the traditional picture agrees with the ELF calculation. However, the ELF exhibits non-Lewis behavior at o2. Instead of three approximately distinct lone pairs, there is a non-Lewis picture of strong delocalization of the lone pairs, leading to the formation of a uniform ring of enhanced electron pair probability. The ring preferentially accepts four HBs in an approximately planar configuration, a configuration that is also found for OH−(aq).26,39,50,57,58 As mentioned above, 90−95% of the configurations have four HBs at o2, corresponding to clear non-Lewis behavior. The ring persists even for the 2A3A cases and in the gas phase, indicating that it is intrinsic to the HOO− electronic structure rather than an effect of the solvation shell. 3.4. Transitions among HOO− Complexes. A “transition” between two HOO− solvation complexes (i.e., HOO− together with its first shell) is said to occur when the local Hbonding pattern (i.e., mAnA) changes. The “duration” of a complex is the length of time a pattern persists until a transition occurs. An “exchange” comprises two or more transitions that start and end with the same pattern, in the course of which at least one first-shell water is replaced by a different one. We employed the geometric criterion to determine transitions and exchanges. Transitions and exchanges both provide insight into solvation dynamics. Typically, there will be many more transitions than exchanges. For example, suppose a first-shell water whose RH···O is near the H-bonding cutoff of 2.5 Å moves slightly further away so that RH···O briefly exceeds 2.5 Å and then returns. Two transitions have occurred but no exchanges. The majority of transitions in our HOO− system were of this “rattling” nature. It could be argued that determination of the dwell time of waters in the first solvation shell should ignore small-amplitude rattles and would be better measured by exchange events. On the other hand, exchanges are not a perfect measure of dwell time either, as they can also include “rattling” when there are two water protons hovering near o1 or near o2 and taking turns venturing briefly across the 2.5 Å threshold. In addition, focusing solely on exchanges requires making an arbitrary choice of one solvation pattern as the “home” pattern from which the system departs and returns. We will discuss the findings for exchanges briefly and then focus primarily on transitions. In the constant temperature (NVT) simulation, there was one first solvation shell exchange event at o1 and no first solvation shell exchange events at o2, possibly suggesting that exchange at o1 may be more probable than at o2, consistent with an extrapolation from the earlier gas phase findings.24 However, such a small number of events is not statistically significant. By contrast, in the constant energy (NVE) runs, of which 69 ps in all were used in the analysis, many exchanges occurred at both the o1 and o2 sites. To detect exchanges at o2, all configurations with exactly four waters donating to o2 were arranged chronologically together with the labels of the donating waters. Whenever the set of four donors was seen to lose one member and gain a different member, an exchange was deemed to have occurred. As noted previously, fourcoordination at o2 occurs in 90% of the configurations. Because o1-1A and o1-2A are both seen frequently, the analysis of exchanges at o1 was carried out for both one- and two-
coordinated o1. Results are listed in Table 5. Exchanges occurred 3−4 times more frequently at o1 than at o2, in keeping with Anick’s prediction.24 Table 5. Picoseconds Analyzed and Number of Exchange Events in the Five NVE Simulations NVE NVE NVE NVE NVE total
1 2 3 4 5
ps
o1-1A
o1-2A
o2-4A
25.130 11.125 10.750 11.000 11.000 69.005
40 9 17 14 18 98
61 15 22 13 13 124
14 4 6 2 5 31
Table 6. Occurrence Percentages for the Top Five Solvation Patterns in the NVE Simulations NVE 1 NVE 2 NVE 3 NVE 4 NVE 5 av std dev
1A3A
1A4A
2A3A
2A4A
0A4A
13.1 0.2 14.8 1.2 14.5 8.8 7.4
46.2 14.3 20.4 8.3 20.0 21.8 14.5
15.1 6.4 6.9 37.3 18.0 16.7 12.6
7.6 71.5 50.2 43.7 38.6 42.3 23.1
14.1 0.1 1.3 0.4 5.3 4.2 5.9
In order to analyze transitions among HOO− complexes at 300 K, the populations of the five most frequent solvation patterns seen in the combined NVE simulations are listed in Table 6, along with their occurrence frequencies for each NVE run. Table 6 shows that the proportions of “2A4A”, “1A4A”, “2A3A”, “1A3A”, and “0A4A” patterns varied widely, depending on the initial configuration. Nevertheless, their average occurrence rates over all five NVE simulations are consistent with what is observed in the NVT simulation. The 1A4A and 2A4A patterns are the predominant structures in the NVE ensemble just as they are in the NVT ensemble. The average populations for the 1A3A and 2A3A patterns indicate that they also play an important role. Overall, the 0A4A pattern is minor. All occurrence rates exhibit large fluctuations from one run to the next because of the relatively short run lengths: In each run, certain solvation patterns have very low occurrence frequencies in the sense that the number of independent occurrences of these complexes can be low because of the fact that either they have long durations or they have low occurrence probabilities. Consequently, statistics on the average durations of the solvation complexes is somewhat poor, and the average durations from these runs must therefore be taken as very approximate. Figure 6 shows the number of HBs accepted and donated by HOO− as a function of time for the run denoted NVE1. We focused on transitions among the top four solvation patterns for a more detailed analysis. Transitions occur via HB forming and breaking events between HOO− and solvating waters. Table 7 shows the total number of transitions among the top four patterns from the full set of NVE configurations. In Table 7, the left column indicates the starting pattern, and the top row designates the ending pattern. Interestingly, we find that the major transition events can be characterized as a “three-level” structure. One terminal level is the solvation pattern 1A3A. One terminal level is 2A4A. The other two 1A4A and 2A3A are in the middle. The transition only occurs 7941
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
The correlation functions CcmAnA(t) for the five NVE trajectories are shown in Figure 7. The five mean durations and their averages and standard deviations are listed in Table 8 (compare Figure 7 to Figure 6 of ref 62, which shows the difference in the decay of Cc(t) for three- and four-coordinated hydroxide ions in water). Considering all computed mean durations, we conclude that the interconversion of solvation patterns occurs on time scales roughly in the range of 54−109 fs. The 2A4A is the most stable solvation pattern and has the longest mean duration of 109 fs. The 1A4A and 2A3A are the two next most stable solvation patterns and have mean durations of 79 and 84 fs, respectively. Considering their large standard deviations, these two species could be considered as being equally long-lived. The shortest average mean duration, 54 fs, belongs to 1A3A. To some extent, mean durations reflect the gas-phase stability of a complex with longer durations corresponding to lower free energies. Using optimized geometries for (HOO−)(H2O)20 clusters at the B3LYP and MP2 levels of theory and with correction for 298 K,24 the free energies for 2A4A, 1A4A, and 2A3A were found to be −32.97, −30.09, and −29.16 kcal/ mol, respectively (see ref 24 for details of the analysis). The 1A3A cluster would clearly have a higher free energy and was not explicitly computed. We find that the gas-phase free energy ordering 2A4A < 1A4A ≲ 2A3A < 1A3A corresponds to the ordering by mean durations, τ2A4A > τ1A4A ≈ τ2A3A > τ1A3A. Moreover, this ordering supports our observation that transitions between the solvation complexes proceed from the 1A3A pattern to 2A4A via the neighboring intermediate structures.
Figure 6. Number of HBs as a function of time for 24 ps of the NVE1 simulation. Geometric and topological criteria for defining H-bonds were used in the upper and lower panels, respectively.
Table 7. Total Number of Transition Events among the Top Four Complexes 2A4A 1A4A 2A3A 1A3A
2A4A
1A4A
2A3A
1A3A
0 94 72 0
98 0 86 34
66 74 0 61
0 38 56 0
4. CONCLUSION In this paper, we reported the structural and dynamical properties of a hydroperoxide anion in water at ambient temperature obtained from extensive AIMD simulations. Our analysis highlights the importance of four HOO− hydration species in water. The terminal oxygen atom of the hydroperoxide anion accepts four hydrogen bonds with water molecules 91% (respectively 95%) of the time, based on a geometric (respectively a topological) criterion for hydrogen bonding. Thus, the HOO− contains one oxygen that adheres closely to the traditional Lewis model while one departs significantly from this picture. Roughly 68% (respectively 77%) of the configurations are 2-fold around the middle oxygen atom of the hydroperoxide anion, while single-fold patterns represent about 29% (respectively 22%) of the configuration. The hydroperoxide anion almost always donates one hydrogen bond. The hydration shell structure of HOO− is characterized by radial distribution functions. We have demonstrated the geometrical parameters: the equilibrium bond length of rOO of HOO− is about 1.5 ± 0.05 Å; the average rOH distance of HOO− is about 0.98 ± 0.04 Å; the bend angle of HOO− is about 101° ± 6°. Analysis of HOO− based on the electron localization function further confirms that a Lewis and a non-Lewis structure exist simultaneously in HOO−. The hypercoordination of a first shell water at site of the terminal oxygen atom is due to the non-Lewis structure, namely, a strong delocalization of the lone pairs and formation of a uniform ring of enhanced electron pairing probability. The constant energy simulation results indicate that transitions between four major solvation complexes occur via specific events involving the 1A4A and
between two neighboring levels or within the level. There are no transition events between the two terminal levels. As we will argue shortly, this multilevel structure accords well with the free energy ordering of gas-phase solvation complexes reflecting these patterns. In order to determine mean duration times for the top four solvation patterns, we employ a strategy based on population autocorrelation functions.59−62 In the present adaptation of the method, we construct, from an NVE trajectory, two population functions, h(t) and H(t), which depend implicitly on a set of initial conditions at t = 0 and a target or reference pattern mAnA. The population function h(t) equals 1 if HOO− has the solvation pattern “mAnA” at time t, and it equals 0 otherwise. The function H(t) is 1 if the solvation pattern remains intact for a total time t and is 0 otherwise. For each solvation pattern “mAnA”, we define the “continuous transfer correlation function” CmAnA (t) and mean duration τmAnA by c Ccm An A(t ) =
⟨h(0) H(t )⟩ ⟨h⟩
(3)
and τ m An A =
∫0
∞
Ccm An A(t ) dt
(4)
where the averages are taken over all initial conditions. 7942
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
Figure 7. Continuous transfer correlation functions for the five NVE simulations. The geometric H-bonding criteria were used to identify solvation patterns.
Table 8. Mean Durations τmAnA for the Top Four Solvation Patterns from the Five NVE Runsa
complexes predicted in this study could be probed by techniques such as neutron scattering in much the same way as has been done for OH−(aq),63−66 and we hope that the present study might inspire future experiments in this direction.
τ (fs) NVE 1 NVE 2 NVE 3 NVE 4 NVE 5 av std dev
2A4A
1A4A
2A3A
1A3A
59 177 97 92 121 109 44
128 57 65 37 110 79 38
66 26 32 224 74 84 80
106 8 76 12 68 54 43
■
AUTHOR INFORMATION
Corresponding Authors
*Z.M.: e-mail,
[email protected]. *D.A.: e-mail,
[email protected]. *M.E.T.: e-mail,
[email protected]. Notes
The authors declare no competing financial interest.
Averages over the five runs and corresponding standard deviations are also provided. a
■
ACKNOWLEDGMENTS M.E.T. and Z.M. acknowledge support from the National Science Foundation, Grant CHE-1012545. M.E.T. acknowledges additional support from Grant CHE-1301314.
2A3A patterns as intermediates. We found that interconversion of HOO− solvation complexes occurs on time scales in the range of 54−109 fs. The most stable solvation structure is that the terminal oxygen atom of the hydroperoxide anion accepts four hydrogen bonds and its middle oxygen atom accepts two hydrogen bonds. The large difference in exchange times between o1 and o2 indicates that the latter is likely the rate-limiting process in the overall diffusion of HOO− in water. Indeed, the relative stability of the 4-fold coordination around o2 suggests that the
■
REFERENCES
(1) Vione, D.; Maurino, V.; Ninero, C.; Pelizzetti, E. The Atmospheric Chemistry of Hydrogen Peroxide: A Review. Ann. Chim. 2003, 93, 477−488. (2) Stone, J. R.; Yang, S. Hydrogen Peroxide: A Signaling Messenger. Antioxid. Redox Signaling 2006, 8, 243−270.
7943
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
Complexes: Quantum Chemical Study. J. Chem. Phys. 2006, 124, 214309. (22) Iyengara, S. S. Dynamical Effects on Vibrational and Electronic Spectra of Hydroperoxyl Radical Water Clusters. J. Chem. Phys. 2005, 123, 084310. (23) Chalmet, S.; Ruiz-Lopez, M. F. The Structures of Ozone and HOx Radicals in Aqueous Solution from Combined Quantum/ Classical Molecular Dynamics Simulations. J. Chem. Phys. 2006, 124, 194502. (24) Anick, D. J. Comparison of Hydrated Hydroperoxide Anion (HOO − )(H 2 O) n Clusters with Alkaline Hydrogen Peroxide (HOOH)(OH−)(H2O)n−1 Clusters, n = 1−8, 20: An ab Initio Study. J. Phys. Chem. A 2011, 115, 6327−6338. (25) Tuckerman, M. E.; Chandra, A.; Marx, D. Structure and Dynamics of OH−(aq). Acc. Chem. Res. 2006, 39, 151−158. (26) Marx, D.; Chandra, A.; Tuckerman, M. E. Aqueous Basic Solutions: Hydroxide Solvation, Structural Diffusion, and Comparison to the Hydrated Proton. Chem. Rev. 2010, 110, 2174−2216. (27) Car, R.; Parrinello, M. Unified Approach for MolecularDynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (28) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic-Behavior. Phys. Rev. A 1988, 38, 3098−3100. (29) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle− Salvetti Correlation-Energy Formula into a Functional of the ElectronDensity. Phys. Rev. B 1988, 37, 785−789. (30) Troullier, N.; Martins, J. Efficient Pseudopotentials for PlaneWave Calculations. Phys. Rev. B 1991, 43, 1993−2006. (31) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (32) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (33) Liu, Y.; Yarne, D. A.; Tuckerman, M. E. Ab Initio Molecular Dynamics Calculations with Simple, Localized, Orthonormal RealSpace Basis Sets. Phys. Rev. B 2003, 68, 125110. (34) Lee, H.-S.; Tuckerman, M. E. Ab Initio Molecular Dynamics Simulation with Discrete Variable Representation Basis Sets: Techniques and Application to Liquid Water. J. Phys. Chem. A 2006, 110, 5549−5560. (35) Lee, H.-S.; Tuckerman, M. E. Structure of Liquid Water at Ambient Temperature from ab Initio Molecular Dynamics Performed in the Complete Basis Set Limit. J. Chem. Phys. 2006, 125, 154507. (36) Lee, H.-S.; Tuckerman, M. E. Dynamical Properties of Liquid Water from ab Initio Molecular Dynamics in the Complete Basis Set Limit. J. Chem. Phys. 2007, 126, 164501. (37) Ma, Z.; Tuckerman, M. Constant Pressure ab Initio Molecular Dynamics with Discrete Variable Representation Basis Sets. J. Chem. Phys. 2010, 133, 184110. (38) Ma, Z.; Zhang, Y.; Tuckerman, M. E. Ab Initio Molecular Dynamics Study of Water at Constant Pressure Using Converged Basis Sets and Empirical Dispersion Corrections. J. Chem. Phys. 2012, 137, 044506. (39) Ma, Z.; Tuckerman, M. E. On the Connection between Proton Transport, Structural Diffusion, and Reorientation of the Hydrated Hydroxide Ion as a Function of Temperature. Chem. Phys. Lett. 2011, 511, 177−182. (40) Martyna, G. J.; Klein, M. L.; Tuckerman, M. E. Nosé Hoover Chains: The Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 1992, 97, 2635−2643. (41) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. Efficient Molecular Dynamics and Hybrid Monte Carlo Algorithms for Path Integrals. J. Chem. Phys. 1993, 99, 2796−2808. (42) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117−1157.
(3) Jones, C. W. In Applications of Hydrogen Peroxide and Derivatives; Clark, J. H., Braithwaite, M. J., Eds.; RSC Clean Technology Monographs; The Royal Society of Chemistry: Cambridge, U.K., 1999. (4) Vincent, M. A.; Palmer, I. J.; Hillier, I. H.; Akhmatskaya, E. Exploration of the Mechanism of the Oxidation of Sulfur Dioxide and Bisulfite by Hydrogen Peroxide in Water Clusters Using ab Initio Methods. J. Am. Chem. Soc. 1998, 120, 3431−3439. (5) Smith, A.; Vincent, M. A.; Hillier, I. H. Mechanism of Acid Dissociation in Water Clusters: Electronic Structure Studies of (H2O)nHX (n = 4, 7; X = OH, F, HS, HSO3, OOSO2H, OOH· SO2). J. Phys. Chem. A 1999, 103, 1132−1139. (6) Yamamoto, K.; Nakamura, A.; Hase, U. Control of Cleaning Performance of an Ammonia and Hydrogen Peroxide Mixture (APM) on the Basis of a Kinetic Reaction Model. IEEE Trans. Semicond. Manuf. 1999, 12, 288−294. (7) Wang, H.; Song, Z.; Liu, W.; Kong, H. Effect of Hydrogen Peroxide Concentration on Surface Micro-roughness of Silicon Wafer after Final Polishing. Microelectron. Eng. 2011, 88, 1010−1015. (8) Siddiqui, S.; Keswani, M.; Brooks, B.; Fuerst, A.; Raghavan, S. A Study of Hydrogen Peroxide Decomposition in Ammonia-Peroxide Mixtures. Microelectron Eng. 2011, 102, 68−73. (9) Fuchs, H. J.; Borders, C. L. J. Affinity Inactivation of Bovine Cu,Zn Superoxide Dismutase by Hydroperoxide Anion HO2−. Biochem. Biophys. Res. Commun. 1983, 116, 1107−1113. (10) Tuan, L. Q.; Umakoshi, H. Toshinori, Liposome Membrane Can Act Like Molecular and Metal Chaperones for Oxidized and Fragmented Superoxide Dismutase. Enzyme Microb. Technol. 2009, 44, 101−106. (11) Babizhayev, M. A.; Vishnyakova, K. S.; Yegorov, Y. E. TelomereDependent Senescent Phenotype of Lens Epithelial Cells as a Biological Marker of Aging and Cataractogenesis: The Role of Oxidative Stress Intensity and Specific Mechanism of Phospholipid Hydroperoxide Toxicity in Lens and Aqueous. Fundam. Clin. Pharmacol. 2011, 25, 139−162. (12) Qiao, L.; Lu, Y.; Liu, B.; Girault, H. H. Copper-Catalyzed Tyrosine Nitration. J. Am. Chem. Soc. 2011, 133, 19823−19831. (13) Aung-Htut, M. T.; Ayer, A.; Breitenbach, M.; Dawes, I. W. Oxidative Stresses and Ageing. In Ageing Research in Yeast; Breitenbach, M., Jazwinski, S., Laun, P., Eds.; Subcellular Biochemistry, Vol. 57; Springer Science+Business Media B.V.: Dordrecht, The Netherlands, 2012; pp 13−54. (14) Dobado, J. A.; Molina, J. M. Ab Initio Molecular Orbital Study of the Hydrogen Peroxide−Water Complex (HOOH···H2O). J. Phys. Chem. 1994, 98, 1819−1825. (15) Gonzalez, L.; Mo, O.; Yanez, M. High-Level ab Initio versus DFT Calculations on (H2O2)2 and H2O2−H2O Complexes as Prototypes of Multiple Hydrogen Bond Systems. J. Comput. Chem. 1997, 18, 1124−1135. (16) Ju, X.-H.; Xiao, J.-J.; Xiao, H.-M. Theoretical Study on Intermolecular Interactions and Thermodynamic Properties of Water−Hydrogen Peroxide Clusters. J. Mol. Struct.: THEOCHEM 2003, 626, 231−238. (17) Kulkarni, A. D.; Pathak, R. K.; Bartolotti, L. J. Structures, Energetics, and Vibrational Spectra of H2O2···(H2O)n, n = 1−6 Clusters: Ab Initio Quantum Chemical Investigations. J. Phys. Chem. A 2005, 109, 4583−4590. (18) Du, D.; Fu, A.; Zhou, Z. Theoretical Study of the Rotation Barrier of Hydrogen Peroxide in Hydrogen Bonded Structure of HOOH−H2O Complexes in Gas and Solution Phase. J. Mol. Struct.: THEOCHEM 2005, 717, 127−132. (19) Sennikov, P. G.; Ignatov, S. K.; Schrems, O. Complexes and Clusters of Water Relevant to Atmospheric Chemistry: H2O Complexes with Oxidants. ChemPhysChem 2005, 6, 392−412. (20) Martins-Costa, M. T. C.; Ruiz-Lopez, M. F. Molecular Dynamics of Hydrogen Peroxide in Liquid Water Using a Combined Quantum/Classical Force Field. Chem. Phys. 2007, 332, 341−347. (21) Kulkarni, A. D.; Pathak, R. K.; Bartolotti, L. J. Effect of Additional Hydrogen Peroxide to H2O2···(H2O)n, n = 1 and 2 7944
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
The Journal of Physical Chemistry B
Article
(43) Tuckerman, M.; Yarne, D.; Samuelson, S.; Hughes, A.; Martyna, G. Exploiting Multiple Levels of Parallelism in Molecular Dynamics Based Calculations via Modern Techniques and Software Paradigms on Distributed Memory Computers. Comput. Phys. Commun. 2000, 128, 333−376. (44) PINY MD is available as an open-source code under the Common Public License and capable of performing a wide variety of molecular dynamics, AIMD, and path-integral simulations under different ensembles. It can be downloaded from the homepage of Prof. Tuckerman, New York University: http://homepages.nyu.edu/ ∼mt33/PINY_MD/PINY.html. (45) Kumar, R.; Schmidt, J. R.; Skinner, J. L. Hydrogen Bonding Definitions and Dynamics in Liquid Water. J. Chem. Phys. 2007, 126, 204107. (46) Henchman, R. H.; Irudayam, S. J. Topological Hydrogen-Bond Definition To Characterize the Structure and Dynamics of Liquid Water. J. Phys. Chem. B 2010, 114, 16792−16810. (47) Bowman, J. M. Self-Consistent Field Energies and Wave Functions for Coupled Oscillators. J. Chem. Phys. 1978, 68, 608−610. (48) Bowman, J. M. The Self-Consistent-Field Approach to Polyatomic Vibrations. Acc. Chem. Res. 1986, 19, 202−208. (49) Carter, S.; Culik, S.; Bowman, J. M. Vibrational Self-Consistent Field Method for Many-Mode Systems: A New Approach and Application to the Vibrations of CO Adsorbed on Cu(100). J. Chem. Phys. 1997, 107, 104508−10469. (50) Tuckerman, M. E.; Marx, D.; Parrinello, M. The Nature and Transport Mechanism of Hydrated Hydroxide Ions in Aqueous Solution. Nature 2002, 417, 925−929. (51) Koput, J. An ab Initio Study on the Equilibrium Structure and Torsional Potential Energy Function of Hydrogen Peroxide. Chem. Phys. Lett. 1995, 236, 516−520. (52) Tachikawa, H.; Abe, S. Direct ab Initio MD Study on the Interaction of Hydroperoxy Radical (HOO) with Water Molecules. Phys. Chem. Chem. Phys. 2010, 12, 3904−3909. (53) Becke, A. D.; Edgecombe, K. E. A Simple Measure of Electron Localization in Atomic and Molecular Systems. J. Chem. Phys. 1990, 92, 5397−5403. (54) Savin, A.; Jepsen, O.; Flad, J.; Andersen, O. K.; Preuss, H.; von Schnering, H. G. Electron Localization in Solid-State Structures of the Elements: The Diamond Structure. Angew. Chem., Int. Ed. Engl. 1992, 31, 187−188. (55) Liu, Y.; Tuckerman, M. E. Protonic Defects in Hydrogen Bonded Liquids: Structure and Dynamics in Ammonia and Comparison with Water. J. Phys. Chem. B 2001, 105, 6598−6610. (56) Mundy, C. J.; Kuo, I.-F. W.; Tuckerman, M. E.; Lee, H.-S.; Tobias, D. J. Hydroxide Anion at the Air−Water Interface. Chem. Phys. Lett. 2009, 481, 2−8. (57) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Molecular Dynamics Simulation of the Solvation and Transport of Hydronium and Hydroxyl Ions in Water. J. Chem. Phys. 1995, 103, 150−161. (58) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Molecular Dynamics Simulation of the Solvation and Transport of H3O+ and OH− Ions in Water. J. Phys. Chem. 1995, 99, 5749−5752. (59) Luzar, A.; Chandler, D. Hydrogen-Bond Kinetics in Liquid Water. Nature 1996, 379, 55−57. (60) Luzar, A. Resolving the Hydrogen Bond Dynamics Conundrum. J. Chem. Phys. 2000, 113, 10663−10675. (61) Chandra, A.; Tuckerman, M. E.; Marx, D. Connecting Solvation Shell Structure to Proton Transport Kinetics in Hydrogen-Bonded Networks via Population Correlation Functions. Phys. Rev. Lett. 2007, 99, 145901. (62) Tuckerman, M. E.; Chandra, A.; Marx, D. A Statistical Mechanical Theory of Proton Transport Kinetics in HydrogenBonded Networks Based on Population Correlation Functions with Applications to Acids and Bases. J. Chem. Phys. 2010, 133, 124108. (63) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M. A.; Soper, A. K. Solvation of Hydroxyl Ions in Water. J. Chem. Phys. 2003, 119, 5001− 5004.
(64) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M. A.; Soper, A. K. Ions in Water: The Microscopic Structure of Concentrated NaOH Solutions. J. Chem. Phys. 2004, 120, 10154−10162. (65) Imberti, S.; Botti, A.; Bruni, F.; Cappa, G.; Ricci, M.; Soper, A. K. Ions in Water: The Microscopic Structure of Concentrated Hydroxide Solutions. J. Chem. Phys. 2005, 122, 194509. (66) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M. A.; Soper, A. K. Solvation Shell of OH− Ions in Water. J. Mol. Liq. 2005, 117, 81−84.
7945
dx.doi.org/10.1021/jp5008335 | J. Phys. Chem. B 2014, 118, 7937−7945
