Quantum Proton Transfer in Hydrated Sulfuric Acid Clusters: A

ARTICLE pubs.acs.org/JPCA

Quantum Proton Transfer in Hydrated Sulfuric Acid Clusters: A Perspective from Semiempirical Path Integral Simulations Shuichi Sugawara, Takehiro Yoshikawa, and Toshiyuki Takayanagi* Department of Chemistry, Saitama University, 255 Shimo-Okubo, Sakura-ku, Saitama City, Saitama 338-8570, Japan

Motoyuki Shiga Center for Computational Science and E-systems, Japan Atomic Energy Agency, Higashi-Ueno 6-9-3, Taito-ku, Tokyo 110-0015, Japan

Masanori Tachikawa Quantum Chemistry Division, Graduate School of Nanobioscience, Yokohama-City University, Seto 22-2, Kanazawa-ku, Yokohama 236-0027, Japan ABSTRACT: We have carried out path-integral molecular dynamics simulations for hydrated sulfuric acid clusters to understand acid-dissociation and hydrogen-bonded structural rearrangement processes in these clusters from a quantum mechanical viewpoint. The simulations were performed using the PM6 semiempirical electronic structure level whose parameters were modified on the basis of the specific reaction parameters strategy so that relative energies of optimized structures, as well as water binding energies reproduce ab initio and density-functional theory calculations. We have found that the acid dissociation processes, first and second deprotonation, effectively occur in a hydrated cluster with a specific cluster size. The mechanisms of the proton-transfer processes were analyzed in detail and it was found that the distance between O in sulfuric acid and O in the proton-accepting water is playing an important role. We also found that the water coordination number of the poton-accepting water is important in the proton-transfer processes.

1. INTRODUCTION It is well-known that hydrated sulfuric acid clusters are key components in atmospheric chemistry.1
6 In fact, sulfuric acid is playing a significant role in acid rain formation, aerosol formation affecting climate change, and ozone depletion in the polar stratosphere. Therefore, the hydrated sulfuric acid system has been extensively investigated from both experimental and theoretical viewpoints. Electronic structure calculation is a useful approach to understand detailed structures of hydrated acid clusters at the molecular level.7
20 In particular, a most fundamental question addressed by previous theoretical studies is: how many water molecules are necessary to cause acid dissociation. Using densityfunctional theory (DFT), Arstila et al.7 have shown that the ionic dissociation does not occur in small H2SO4 3 H2O and H2SO4 3 (H2O)2 clusters while the relative energy between the ionic HSO4
3 H3O+ 3 (H2O)2 structure and the neutral (nonionic) H2SO4 3 (H2O)3 structure is small in the trihydrate sulfuric acid cluster. Bandy and Ianni8 studied the structures and hydration free energies for H2SO4 3 (H2O)n (n = 1
7) clusters at the B3LYP DFT level. Re et al.9 have performed more systematic B3LYP-level calculations on the structures and energetics of H2SO4 3 (H2O)n (n = 1
5) clusters. They have found that this r 2011 American Chemical Society

system has many hydrogen-bonded configurations including both neutral and ionic structures within a several kcal/mol energy range. Moreover, Ding and Laasonen15 have carried out electronic structure calculations for larger H2SO4 3 (H2O)n (n = 6
9) clusters to understand the complete deprotonation processes to form the SO42
moiety. They have concluded that sulfuric acid can be fully deprotonated in clusters containing eight or more water molecules. More recent theoretical studies have focused on dynamical aspects of the hydrated sulfuric acid system using first-principles molecular dynamics approach.21
26 Choe et al.22 have applied Car
Parrinello molecular dynamics (CPMD) simulations to aqueous solutions of sulfuric acid in the condensed phase and discussed the proton diffusion mechanisms. Anderson et al.24 have recently reported the CP2K study for (H2SO4)m 3 B 3 (H2O)6 clusters, where B is a base molecule, in order to understand the effect of base molecules on the deprotonation processes of H2SO4. Hammerich et al.25 have performed ab initio molecular dynamics simulations for a solution of a single H2SO4 Received: March 14, 2011 Revised: September 9, 2011 Published: September 12, 2011 11486

dx.doi.org/10.1021/jp202380h | J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A molecule in a periodic box with 63 water molecules. Interestingly, they have observed the fully deprotonated SO42
moiety in their calculated trajectories. This result is in high contrast with the CPMD result of Choe et al., where no such species were found.22 Our research group is currently interested in nuclear quantum effects on the proton transfer, as well as the structural rearrangement processes, in hydrated acid clusters. Path-integral molecular dynamics (PIMD) technique is a useful tool to simulate quantum molecular systems in thermal equilibria.27 First-principle electronic structure methods, such as ab initio molecular orbital and DFT approaches, can be combined with PIMD simulations where the Born
Oppenheimer potential is computed on the fly;28,29 however, actual applications are currently limited to small chemical systems and relatively short simulation runs because of huge computational demands. As an alternative, we have been recently employing semiempirical molecular orbital methods in the combination with PIMD.30
34 It should be mentioned that a typical CPU time for obtaining the potential energy and its gradients at a given molecular geometry with the semiempirical method is much faster than that with the standard DFT method by a factor of 102
103. In our previous study,30 we have performed PIMD simulations for the H2SO4 3 (H2O)n (n = 1
6) clusters using a semiempirical model called PM6,35 for which the description of hydrogenbonded interaction is significantly improved over the widely used AM1 and PM3 methods. We were able to find important quantummechanical natures in proton transfer processes, ionic dissociation, and hydrogen-bonded network rearrangements in the H2SO4 3 (H2O)n clusters. However, we also noticed that the PM6 potential energy surface needs improvement to study the hydrated acid cluster system in more detail. For example, from PM6 geometry optimization, we found some structures of small H2SO4 3 (H2O)n (n = 1
3) clusters having bifurcated hydrogenbonds (i.e., the structure in which two protons of water are weakly bound to the O-site of H2SO4), which is not optimized with ab initio calculations. In addition, we found that water binding energies obtained at the PM6 level are underestimated in light of previous DFT calculations. Therefore, the potential energy landscape should be more accurate to describe thermal and quantum fluctuations of hydrogen bonded network in H2SO4 3 (H2O)n clusters. Fortunately, since semiempirical methods contain various parameters that can be tuned to improve the accuracy of the potential energy surface for a specific molecular system. This computational strategy was originally developed by Truhlar and his co-workers and is called the SRP (specific reaction parameters) method.36 In this work we have redetermined PM6 parameters for the H2SO4 3 (H2O)n cluster system referring to the energy values obtained from ab initio and DFT calculations at various local minimum geometries. Using the newly developed PM6/SRP potential energy surfaces, PIMD calculations have been performed.

2. PM6/SRP PARAMETRIZATION In this work, the semiempirical PM6 Hamiltonian35 has been used as a starting point for the SRP adjustment, as mentioned above. In an early stage of the present investigation, we have optimized the PM6 parameters using the standard minimization scheme,36 where semiempirical parameters are optimized by minimizing the target functional value representing a measure of the difference between ab initio electronic structure properties (energy, geometric parameters, or vibrational frequencies), and

ARTICLE

Table 1. Original PM6 Parameters and Specific Reaction Parameters (SRP) Determined in This Work for the H2SO4(H2O)n Cluster System atom H

O

S

parameter

PM6

PM6/SRP1

PM6/SRP2
11.3177

USS [eV]


11.247


11.2938

ζS [a0
1]

1.268641

1.294274

1.322400

βS [eV]


8.35298


8.31245


8.28630

GSS [eV]

14.44869

14.39617

14.53310

a1 [none] b1 [Å
2]

0.024184 3.055953


0.03042 3.035113


0.03000 3.05320

c1 [Å]

1.786011

1.781222

1.797200

USS [eV]


91.6788


91.6116


93.2687
70.7399

UPP [eV]


70.4609


70.4604

ζS [a0
1]

5.421751

5.346753

5.361100

ζP [a0
1]

2.27096

2.205743

2.232800

βS [eV]


65.6351


65.5645


64.0678

βP [eV] GSS [eV]


21.6226 11.30404


21.6739 11.24599


21.9825 11.20350 13.87990

GPP [eV]

13.61821

13.57407

GSP [eV]

15.80742

15.73668

15.96090

GP2 [eV]

10.33277

10.26957

10.20460

HSP [eV]

5.010801

5.107848

5.20150

a1 [none]


0.01777

0.005396

0.006000

b1 [Å
2]

3.05831

2.985559

2.984500

c1 [Å] USS [eV]

1.896435
47.5307

1.852654
47.5952

1.839100
48.3073

UPP [eV]


39.191


39.2429


39.9908

UDD [eV]


46.3069


46.3376


46.8495

ζS [a0
1]

2.192844

2.041152

2.051500

ζP [a0
1]

1.841078

1.613646

1.610200

ζD [a0
1]

3.109401

2.851224

2.817300

βS [eV]


13.8274


13.7234


13.8806

βP [eV] βD [eV]


7.66461
9.98617


7.59136
10.2853


7.58930
10.1559

GSS [eV]

9.17035

9.35623

9.39650

GPP [eV]

8.165473

8.388687

8.459700

GSP [eV]

5.944296

6.037333

6.046200

GP2 [eV]

7.301878

7.085321

7.191200

HSP [eV]

5.005404

5.105629

5.143100

the corresponding values obtained with the target semiempirical level. We have applied this technique using various reference ab initio and DFT data for the small H2SO4 3 H2O and H2SO4 3 (H2O)2 clusters. However, we have found that such parametrization does not work for larger hydrated clusters. More specifically, we have noticed that the parameters optimized for these small clusters do not reproduce the important quantities such as water binding energies for larger hydrated clusters. The above-mentioned failure is presumably because of the difficulty for obtaining good initial parameters for optimization. Therefore, we have employed the following screening procedure. First, we have randomly generated 105 sets of parameters within (30% range of the original PM6 parameters. For each set of parameters we have done geometry optimizations for small H2SO4 3 H2O and H2SO4 3 (H2O)2 clusters. Then, we have chosen 200 sets of parameters which reproduce relative energies and water binding energy within 2 kcal/mol, hydrogen-bond 11487

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A

ARTICLE

Figure 2. Water binding energies of the most stable structures as a function of the number of water molecule in the H2SO4 3 (H2O)n cluster. B3LYP, BLYP, and PW91 results are taken from refs 9 and 15, respectively.

Table 2. Comparison of Water Binding Energies (in kcal/ mol) for the H2SO4 3 (H2O)n (n = 1
9) Clusters

Figure 1. Comparison of optimized H2SO4 3 H2O and H2SO4 3 (H2O)2 clusters obtained at B3LYP, PM6/SRP1, PM6/SRP2, and PM6 levels of theory. The B3LYP results are taken from ref 9. Some selected bond distances (in Å) and angles (in deg.) are shown. Relative energies are also shown in the unit of kcal/mol.

distances within 0.1 Å, and hydrogen-bonding angles within 10 degrees (comparing to the B3LYP results,37 see Figure 1). For each set of the chosen parameters, we have done extensive geometry optimization for larger H2SO4 3 (H2O)n (n = 3
12) clusters to obtain water binding energy as a function of the cluster size n. Here, we have performed a comprehensive geometry optimization for H2SO4 3 (H2O)n starting from various randomly generated geometries (1000 initial structures for each cluster size) with atom
atom distance, 2.5 Å e ROO e 3.5 Å for intermolecular species, and with atom
atom distances, 1.4 Å e RSO e 1.7 Å and 0.7 Å e ROH e 1.5 Å, for intramolecular species. In addition, we have calculated the potential energy profile of the proton-transfer process for the H2SO4 3 (H2O)2 cluster. From the detailed comparison, we have finally chosen the parameters which reproduce the water binding energies obtained at the MP2, B3LYP, and PW91/DNP levels and reasonably reproduce the proton-transfer barrier for the H2SO4 3 (H2O)2 cluster. It should be mentioned that Nadykto et al.19
21 have found good agreement between the water-binding free energies computed at the PW91 level and experimental data for larger hydrated clusters. The parameter set, which has thus been obtained with this screening procedure, is presented in Table 1 along with the original PM6 parameters. The obtained parameters are denoted as PM6/SRP1 hereafter. Using the obtained PM6/SRP1 parameter set as initial guess, we have further performed parameter optimization on the basis of the B3LYP/D95++(d, p) results (including geometries and relative energies for local minima, proton-transfer barrier profile for the dimer, water binding energies) for the small H2SO4 3

n

PM6/SRP1

PM6/SRP2

PM6

MP2a

B3LYPb

1 2

12.9 25.7

12.5 25.0

10.8 22.1

12.7 25.2

12.8 25.6

3

38.6

38.5

31.9

38.7

38.5

4

55.6

54.3

40.7

53.1

5

69.8

68.6

51.3

65.1

6

82.9

80.7

61.9

PW91c

BLYPc

69.4

55.9

83.3

67.1

7

95.7

93.6

70.6

96.7

78.3

8

110.0

105.4

81.4

108.2

87.3

9

121.5

119.2

87.7

121.8

97.0

a

Values were obtained at the MP2/aug-cc-pVDZ level. b Taken from ref 9. c Taken from ref 15.

H2On (n = 1, 2, and 3) clusters. The optimized parameters, denoted as PM6/SRP2, are also summarized in Table 1. Notice that the PM6/SRP2 parameters differ by only small amounts from the PM6/SRP1 ones. It is generally difficult to obtain fully optimized parameters giving a global minimum due to a very high dimensional optimization landscape. Although our PM6/SRP2 parameter set may correspond to a local minimum, it is encouraging that the obtained parameters reproduce DFT results at a fairly good level. First, we compare low-lying optimized structures of H2SO4 3 H2O obtained by the B3LYP, PM6/SRP1, PM6/SRP2, and original PM6 methods in Figure 1. The B3LYP result shows that the most stable structure has a cyclic hydrogen-bonded form. The PM6/SRP1 and PM6/SRP2 methods give structures very similar to B3LYP as the most stable one, while the original PM6 method predicts a structure largely distorted compared with that of B3LYP. The B3LYP, PM6/SRP1, and PM6/SRP2 methods predict the second-lowest energy structures, where two hydrogen atoms of sulfuric acid contribute to hydrogen-bonding with a single water molecule. A similar structure is obtained by the PM6 method as well; but the relative stability of the two lowest-energy structures were found to be opposite to that of the B3LYP method. Figure 1 also compares geometric structures of the H2SO4 3 (H2O)2 complex optimized with the B3LYP, PM6/ 11488

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A

ARTICLE

Table 3. Energy difference (in kcal/mol) between the most stable neutral and partially deprotonated cluster structures for the H2SO4 3 (H2O)n (n = 3
9) clusters n

PM6/SRP1

PM6/SRP2

PM6

B3LYP/D95(d, p)b

B3LYP/D95++(d, p)b

BLYP/DNPc

B3LYP/6-311++G(2d, 2p)d

ΔE = E(HSO4
)
E(H2SO4)
1.3

3

0.05 (
0.89)a

1.37 (
0.60)

0.6 (0.4)a

0.3 (
0.9)a

4

2.34 (1.07)

4.78 (1.48)

0.1

1.6 (1.6)

0.9 (0.6)

5

5.20 (3.51)

4.33 (4.50)

1.8

4.5 (3.7)

2.7 (2.4)

6 7

6.45 (4.87) 5.35 (4.51)

5.85 (5.20) 10.27 (7.72)

3.5 3.2

8

10.70 (8.99)

4.72 (3.99)

3.8

9

8.90 (7.11)

13.18 (9.61)

4.2


1.4 (
0.2)a 0.2 (0.4) 0.5 (0.4) 4.9

3.2 (1.2)

a

Numbers in parentheses indicate relative energies corrected with zero-point vibrational energies within harmonic approximation. b Taken from ref 9. c Taken from ref 15. d Taken from ref 8.

SRP1, PM6/SRP2 and PM6 methods. It is seen that the structures optimized with the PM6/SRP1 and PM6/SRP2 methods reasonably agree with those optimized at the B3LYP level. In the case, where a water dimer forms a hydrogen-bonding cyclic structure with sulfuric acid, the original PM6 method gives a bifurcated hydrogen-bonded structure. This is in contrast to a single hydrogen-bonded structure obtained by the B3LYP and PM6/SRP methods. Figure 2 and Table 2 compare water binding energies for the lowest-energy structures as a function of the cluster size n, that is, the number of water molecule in the cluster. It is seen that the present PM6/SRP1 and PM6/SRP2 results fairly agree with the previously reported PW91/DNP15 and the B3LYP/D95++(d, p) results.9 Notice that the original PM6 method, as well as the BLYP/DNP method,15 tend to underestimate water binding energies systematically as the cluster size increases. Table 3 summarizes the relative energy difference between the lowest undissociated (H2SO4 or neutral) and dissociated (HSO4
) structures obtained for each cluster size. The energies of the undissociated and dissociated structures are found to become comparable for n = 3
4 . This is also true after the harmonic zero-point energy correction is added to the energy of each structure. For larger clusters with n > 4, the dissociated structures is lower in potential energy than the undissociated ones both by the PM6 and PM6/SRP methods. This trend is qualitatively consistent with previous DFT calculations. Table 4 summarizes the energy difference between the lowest singly (HSO4
) and doubly (SO42
) deprotonated structures optimized for each cluster size. From the PM6/SRP1 and PM6/SRP2 results presented in Table 4, we can see that the energies of the singly and doubly deprotonated structures become energetically comparable at the n = 9
10 cluster size. Again, this trend is qualitatively similar to previous DFT results.15 The proton-transfer barrier is a key quantity that is related to the extent of structural rearrangement via proton exchange processes. However, as far as we are aware, proton-transfer barriers in the hydrated sulfuric acid clusters have not been fully addressed in previous electronic structure studies from the quantum chemical viewpoint. In Figure 3, potential energy profiles obtained from several DFT, PM6/SPR1 and PM6/ SRP2 calculations are plotted along the triple proton-transfer reaction coordinate of H2SO4 3 (H2O)2 cluster. The DFT results presented in Figure 3 were obtained from the intrinsic reaction coordinate (IRC) calculations implemented in the Gaussian03 program package.38 The present PM6/SRP1 and PM6/SRP2

methods are found to give barriers of 4.91 and 4.95 kcal/mol, respectively. These values are slightly smaller than the B3LYP/ D95++(d, p) result but is close to the B3LYP/D95(d, p) result. It is seen that the PW91 method39 gives relatively smaller barrier heights for both the D95(d, p) and D95++(d, p) basis sets.40
42 To understand the reliability of the DFT values, we have also performed ab initio MP2-level calculations. The initial and transition-state geometries were determined at the MP2/augcc-pVDZ level, although we did not perform IRC analyses. The MP2 calculations give barriers of 8.29, 6.92, 6.72, and 5.68 kcal/ mol with aug-cc-pVDZ,43 aug-cc-pVTZ,44 6-31++G(3df, 3pd), and 6-31++G(3df, 3p2d) basis sets,45 respectively. Thus, it is found that the MP2 calculation with large basis sets yields a barrier height comparable to the B3LYP values. In the case of the DFT results, the barrier height somewhat depends on the choice of DFT functional as well as the quality of the basis set used. Nevertheless, it is seen that the PM6/SRP1 and PM6/SRP2 results are in the range of the results obtained with the commonly used DFT-B3LYP methods. Since we found that the present PM6/SRP1 and PM6/SRP2 methods show reasonably good agreement with the DFT results for small clusters, it may be interesting to report some characteristic features of the low-lying optimized structures for larger clusters. Figure 4 presents typical geometric structures optimized with the PM6/SRP1 method for the H2SO4 3 (H2O)10 cluster. We found that the fully deprotonated SO42
moiety is seen only when all the oxygen atoms in sulfuric acid are hydrogen-bonded by water molecules. Notice that the counter H3O+ cations are located just next to the SO42
moiety, indicating that the corresponding structure can be categorized into a so-called contaction-pair type. Figure 4 also shows another optimized structure but slightly higher in relative energy. As this structure includes the HSO4
and H3O+ moieties bridged by the H2O molecules in between, it can be called solvent-separated-ion-pair type.

3. METHOD OF PIMD SIMULATIONS In the PIMD method, the nuclei are quantized following Feynman’s path-integral formalism,46 where the quantum mechanical character of nuclei is described by cyclic bead chains consisting of classical quasiparticles. The PIMD accounts not only for thermal fluctuations and bond rearrangements, but also for nuclear tunneling and quantized vibrational effects. Thermal equilibrium averages are computed from time averages over fictitious molecular dynamics trajectories generated from the 11489

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A

ARTICLE

Table 4. Energy Difference (in kcal/mol) between the Most Stable Singly Deprotonated and Doubly Deprotonated Cluster Structures for the H2SO4 3 (H2O)n (n = 6
12) Clusters n

PM6/SRP1a

PM6/SRP2a

BLYP/DNPb,c

PW91/DNPb,c

RI-MP2/TZVPb,c

ΔE = E(SO42
)
E(HSO4
)
3.58 (
4.96)


8.1 (
10.2)


4.6 (
6.7)

7 8


2.42 (
3.94)
3.48 (
4.20)


3.05 (
4.59)
0.07 (
0.88)


5.0 (
5.8)
2.4 (0.1)


1.1 (
3.9) 0.0 (2.9)


2.4 (0.3)

9


0.75 (
1.66)


1.17 (
0.73)


1.2 (
0.4)

1.6 (3.4)


1.5 (0.2)

10

0.07 (
0.65)

1.16 (0.57)

11

0.76 (
1.64)

2.73 (2.21)

12

2.16 (1.24)

1.42 (0.86)

6

a

Numbers in parentheses for the PM6/SRP1 and PM6/SRP2 results indicate relative energies corrected with zero-point vibrational energies within harmonic approximation. b Taken from ref 15. c Numbers in parentheses for the B3LYP, PW91, and RI-MP2 results indicate relative free energies at 298 K.

Figure 3. Potential energy profiles for the proton-transfer process in the H2SO4 3 (H2O)2 cluster as a function of the reaction coordinate.

equation-of-motions described in the normal mode coordinate representation. The PIMD code is developed based on velocity
Verlet time integrator47,48 with multiple time step algorithm49 using normal mode coordinates. The system temperature is controlled with a massive N
ose
Hoover chain technique.50,51 Note that the only important assumptions made here are Born
Oppenheimer approximation and neglect of the exchange between hydrogen nuclei, which should work well in the present system. Our PIMD computational code was linked to the semiempirical MOPAC2007 code52 so that each input data of a molecular configuration returns the output data of the corresponding potential energy and atomic forces at the semiempirical PM6/SRP level. The time increment was set to Δt = 5 atomic unit (0.12 fs). The PIMD simulations were carried out with 48 imaginary-time slices (number of beads) at T = 300 K. After equilibration, the PIMD simulations have been run for 120
200 ps for the respective clusters depending on their statistical convergence. The PIMD simulations were performed for the

Figure 4. Typical structures optimized with the PM6/SRP1 method for the H2SO4 3 (H2O)10 cluster. The arrows indicate the H3O+ hydronium ions in the cluster. The relative energies with respect to the most stable structure are also shown in the unit of kcal/mol.

H2SO4 3 (H2O)n (n = 4, 6, 8, 10, and 12) clusters. At this temperature, various structural rearrangements of hydrogenbonded network are found to occur, while the water molecules rarely evaporate from the cluster. The details of our computational procedure are also described in refs 28
30.

4. RESULTS AND DISCUSSION Figure 5 displays the atom
atom distribution function for the Os
H pair obtained from the PIMD simulations at T = 300 K, where Os is one of the four oxygen atoms in H2SO4 (HSO4
or SO42
), while H indicates one of all hydrogen atoms in the H2SO4 3 (H2O)n system. Since we found that both the PM6/ SRP1 and PM6/SRP2 parameter sets give very similar results, only the PM6/SRP1 result is presented in Figure 5. The first peaks around 1 Å indicate the distributions of the chemically bound Os
H covalent bond in H2SO4 or HSO4
in the H2SO4 3 (H2O)n 11490

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A

ARTICLE

Figure 5. Atom-atom distribution functions (solid lines) of the Os
H pair for the H2SO4 3 (H2O)n (n = 4, 6, 8, 10, and 12) clusters obtained from the PIMD simulations at T = 300 K on the PM6/SRP1 potential energy surfaces. Dashed lines correspond to integrations of the distribution functions.

Table 5. Fraction of Each Species Determined from PIMD Simulations at T = 300 K for the H2SO4 3 (H2O)n (n = 4, 6, 8, 10, and 12) Clusters on the Semiempirical PM6/SRP1 Potential Energy Surfaces n

H2SO4

HSO4


SO42


4

0.62

0.38

0.00

6

0.01

0.99

0.00

8

0.00

1.00

0.00

10

0.00

0.92

0.08

12

0.00

0.91

0.09

clusters. The second broad peaks around 1.8
1.9 Å are the contribution from hydrogen-bonds between Os in sulfuric acid and H in water molecule. It is interesting to notice that the intensity of the first peak decreases as the cluster size is larger. This trend clearly indicates that the fractions of the HSO4
and SO42
protondissociated moieties increase with the number of surrounding water molecules. In addition, density minima can be seen around ROs
H = 1.2
1.4 Å for all the clusters studied here. To understand this region more quantitatively, we also plot integrations of the distribution functions in Figure 5. This integration value gives information on the coordination number at given Os
H distance. A flat behavior in the integration curves is clearly seen around 1.3 Å. Thus, we take 1.3 Å as the critical distance to define the species to be either ionic or nonionic. When the Os
H bond length is longer than 1.3 Å, the corresponding distance is assumed to be dissociated. More specifically, if one of the four Os
H distances is shorter than 1.3 Å, while other Os
H distances are longer than 1.3 Å, the ionic HSO4


Figure 6. Two-dimensional contour plots of the probability distribution obtained from the PIMD simulation for the H2SO4 3 (H2O)4 cluster on the PM6/SRP1 and PM6/SRP2 potential energy surfaces as a function of (a) ΔR and R(Os
Ow), (b) ΔR and θ, and (c) ΔR and Cw. The definitions of these coordinates are also shown.

species is assumed to exist in the cluster. Similarly, if two Os
H bonds in H2SO4 are simultaneously longer than 1.3 Å, the doubly ionized SO42
species is assumed to exist in the cluster. The fractions thus determined are summarized in Table 5. In the smallest cluster (n = 4) studied here, both nonionic H2SO4 and ionic HSO4
species coexist in the cluster. In the n = 6 and 8 clusters, the fraction of the H2SO4 neutral species is quite small indicating that H2SO4 is nearly dissociated into HSO4
. The PIMD calculation gives the fraction of the fully deprotonated SO42
species to be 8 and 9% for the n = 10 and 12 clusters, respectively. From the proton distributions presented in Figure 5, it is interesting to notice that the value of the distribution around ROs
H ∼ 1.3 Å is significantly dependent on the cluster size. In fact, the corresponding density is nearly zero for the n = 6 and 8 clusters while the densities show finite (nonzero) values for the n = 4, 10, and 12 clusters. This result clearly suggests that the proton-transfer processes effectively occur in these clusters. More specifically, in the H2SO4 3 (H2O)4 cluster, the protontransfer process causing the H2SO4 T HSO4
exchange with the nearest water molecule should occur. On the other hand, in the n = 10 and 12 clusters, it is suggested that the HSO4
T SO42
proton-exchange process with the nearest water molecule occurs. To understand these proton-transfer processes more quantitatively, 11491

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A

ARTICLE

we have analyzed these processes with two-dimensional contours plots of nuclear density distributions as a function of appropriate coordinates. Figure 6a shows the two-dimensional contour plots of the density distribution of H2SO4 3 (H2O)4 as a function of ΔR = rs
rw and R(Os
Ow), where rs and rw are distances between H and Os and between H and Ow, respectively (defined also in the figure). Here, Ow is the oxygen atom of the proton-accepting water having the shortest distance between Os and the oxygen atom in water. ΔR corresponds to the proton-transfer reaction coordinate; The negative region of ΔR corresponds to the nonionic H2SO4 species, while the positive region to the ionic HSO4
species. From these plots, we can see that the protontransfer process is strongly correlated with the Os
Ow distance, as frequently pointed out in a number of previous theoretical studies; the transfer process is promoted when the Os
Ow distance becomes shorter. This result is also in contrast with our previous PIMD study using the original PM6 parameters, where we did not find such a strong correlation between the proton-transfer coordinate and the Os
Ow distance.30 Presumably, the proton-transfer transition-state geometries at the original PM6 level are not reasonable since the bifurcated hydrogenbond is stabilized too much. Figure 6b displays another kind of contour plot as a function of ΔR and angle θ (= — Os
H
Ow). It is interesting to note that proton-transfer is promoted when θ is large although the corresponding angle is somewhat smaller than 180 degrees. Figure 6c displays the contour plot as a function of ΔR and the collective coordinate Cw, where Cw is a dimensionless coordination number of water with respect to the proton-accepting water (which is directly attached to H2SO4). To obtain a smooth distribution with respect to this coordinate, Cw is defined as Cw ¼

∑i

1 þ tanh½αðr0
rOw
Oi 0 Þ w

2

ð1Þ

where the cutoff distance r0 was set to 3.2 Å and damping parameter α to 1.0 Å
1 in this work. Ow is the oxygen atom of the proton-accepting water and Owi is the oxygen atom of all the other water molecules. The functional form of the coordination number Cw is different from that used in Reference;53 however, its physical implication is essentially the same. From the result in Figure 6c, it is seen that the proton-transfer process is strongly correlated with Cw; the ionic dissociation to form HSO4
effectively occurs when the water coordination number Cw is large. Another interesting point seen in Figure 6c is that the distribution for ΔR < 0 has two broad peaks at Cw ≈ 1.4 and 0.8. Our detailed analysis shows that these two peaks correspond to different hydrogen-bonded networks in the water cluster moiety. More specifically, the former distribution (Cw ≈ 1.4) qualitatively corresponds to the branched hydrogen-bonded water cluster structures while the latter distribution (Cw ≈ 0.8) to the linear hydrogen-bonded water cluster structures. To understand nuclear quantum effects in proton-transfer processes for the H2SO4 3 (H2O)4 cluster, we have further carried out classical mechanics simulations as well as the PIMD simulation but with a smaller bead number of 8 (a nonconverged bead number). Figure 7 shows resulting contour plots similar to Figure 6a. In the classical mechanics simulation result presented in Figure 7a, only the ionic HSO4
species is seen. The fraction of the HSO4
species was estimated to be 100% in this simulation (see Table 5). This indicates that the ionic structures are energetically

Figure 7. Two-dimensional contour plots of the probability distribution for the H2SO4 3 (H2O)4 cluster on the PM6/SRP1 potential energy surface as a function of ΔR and R(Os
Ow): (a) classical result and (b) PIMD result with 8 beads.

stable than the neutral structures even after temperature effects were taken into account. This trend is qualitatively understood from the potential energy difference between ionic and neutral structures for the H2SO4 3 (H2O)4 cluster presented in Table 3, where the deprotonated ionic structure is lower in potential energy than the neutral structure. It is interesting to notice that the energy difference between the neutral and ionic structures becomes smaller after zero-point vibrational energy correction. Thus, vibrational quantization significantly affects the thermodynamic stability between the ionic and neutral structures in the cluster. The PIMD simulation with 8 beads show an intermediate result between Figure 6a and Figure 7a, as is presented in Figure 7b. The fraction of the HSO4
ionic species was obtained to be 50%, which is larger than the converged result of 38% presented in Table 5. These computational results thus presented in Figure 7 clearly indicate that nuclear quantum effects should be taken into account for understanding of free-energy profiles for proton-transfer processes at a quantitative level. Figure 8 displays contour plots similar to Figure 6 but for the larger H2SO4 3 (H2O)10 and H2SO4 3 (H2O)12 clusters on the PM6/SRP1 potential energy surfaces. As mentioned previously, we can understand the detailed proton-exchange mechanism, HSO4
T SO42
, from the results presented in Figure 8. Interestingly, a general trend is quite similar to the first proton-dissociation mechanism seen in the H2SO4 3 (H2O)4 cluster case. Again, both the Os
Ow distance and — Os
H
Ow angle are coupled with the proton-transfer reaction coordinate. In addition, the coordination number around the proton-accepting water is also seen to play an important role in this second deprotonation process. It is interesting to notice that the distribution around Cw ≈ 0.7
0.8 (linear hydrogen-bonded water cluster structures) is much larger than that around Cw 1.4
1.5 (branched hydrogenbonded water cluster structures). 11492

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A

ARTICLE

Voth59 have found that nuclear quantum effects significantly lower the proton-transfer barrier in their simulations for an excess proton in bulk phase water. This also indicates that quantization of the nuclei leads to significant differences in free-energy surfaces obtained from classical mechanics simulations. In our previous path-integral simulations for HCl(H2O)n and H2SO4(H2O)n clusters,30,31 we have found that quantum simulations show liquid-like cluster structures while classical simulations show solid-like structures at the same temperature presumably due to lower quantum proton-transfer barriers. Interestingly, Schwelger et al.60 have stated that including proton quantum effects is approximately equivalent to a temperature increase of 50 and 100 K.

Figure 8. Two-dimensional contour plots of the probability distribution obtained from the PIMD simulation for the H2SO4 3 (H2O)10 and H2SO4 3 (H2O)12 clusters on the PM6/SRP1 potential energy surfaces. Plot conventions are the same as Figure 6. Left and right panels show the result for the H2SO4 3 (H2O)10 and H2SO4 3 (H2O)12 clusters, respectively.

From the results presented in Figures 6 and 8, the acid dissociation processes of a sulfuric acid H2SO4 molecule, including first and second deprotonation, is strongly correlated with the water coordination number around the proton-accepting water molecule. It should be mentioned that the importance of the water coordination number has been pointed out in other acid systems.54
58 For example, Hynes and his co-workers54
58 have carried out molecular dynamics simulations of nitric acid dissociation at an aqueous interface and discussed the acid dissociation mechanisms of nitric acid from a microscopic viewpoint in some detail. They have found that the coordination number for the proton-accepting water plays a very important role in acid dissociation although they have mentioned that their simulations did not suffice to provide a complete statistical picture due to huge computational demands. Finally, we would like to give some comments on quantum effects for proton-transfer processes. Extensive theoretical studies demonstrating the importance of quantum effects in proton-transfer processes have already been available. For example, Schmitt and

5. CONCLUDING REMARKS In this work, we have studied the acid dissociation of a sulfuric acid H2SO4 molecule in hydrated clusters using the quantum PIMD simulations combined with an on-the-fly semiempirical direct dynamics technique. Starting an original PM6 semiempirical level, we have modified its parameters so as that the optimized structures, their relative energies, water binding energies, and proton-transfer barriers give reasonable agreement with previously reported DFT and ab initio electronic structure results. Because the semiempirical PM6/SRP electronic structure calculation is computationally efficient, the “semiempirical PIMD” calculation allows one to collect enough statistics to investigate largeamplitude and collective motions of complex quantum many-body systems. This helps the analyses of structural rearrangements and proton-transfer processes of molecular clusters in different sizes, which would require a huge computational cost if an ab initio PIMD simulation were to carry out instead. From our PIMD simulations, we have been able to see several important features in the acid dissociations as well as the structural rearrangements in hydrated sulfuric acid clusters. We have found that the proton-exchange process with water, H2SO4 T HSO4
or HSO4
T SO42
, mainly occur in a specific cluster sizes; n = 4 for the former, 10 and 12 for the latter. This is quite reasonable since the relative free energy difference between these nonionic and ionic species should be dependent on the number of water molecules. It has been found that the distance between the sulfuric acid oxygen atom and the oxygen atom in the proton-accepting water molecule is playing an important role in these proton-transfer processes. In addition, the coordination number of water around the proton-accepting water is another important factor for controlling the proton-exchange process. Currently, we are planning to apply the present method to very large clusters as well as bulk solutions to explore the mechanism of acid dissociation at aqueous surface. This is relevant to aerosol uptake processes that is of fundamental importance in the field of atmospheric chemistry. ’ AUTHOR INFORMATION Corresponding Author

*Fax: +81-48-858-3700. E-mail: [email protected].

’ ACKNOWLEDGMENT We would like to thank Grant-in-Aid for Scientific Research and for the priority area by Ministry of Education, Culture, Sports, Science and Technology, Japan. 11493

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494

The Journal of Physical Chemistry A

’ REFERENCES (1) Calhoun, J. A.; Charlson, R. J.; Bates, T. S. Geophys. Res. Lett. 1991, 18, 1877–1880. (2) Charlson, R. J.; Schwartz, S. E.; Hales, J. M.; Cess, R. S.; Coakley, J. A.; Hansen, J. E.; Hofmann, D. J. Science 1992, 255, 423–430. (3) Laaksonen, A.; Talanquer, V.; Oxtoby, D. W. Annu. Rev. Phys. Chem. 1995, 46, 489–524. (4) Bianco, R.; Hynes, J. T. Acc. Chem. Res. 2006, 39, 159–165. (5) Vaida, V.; Kjaergaard, H. G.; Hintze, P. E.; Donaldson, D. J. Science 2003, 299, 1566–1568. (6) Donaldson, D. J.; Tuck, A. F.; Vaida, V. Chem. Rev. 2003, 103, 4717–4730. (7) Arstila, H.; Laasonen, K.; Laaksonen, A. J. Chem. Phys. 1998, 108, 1031–1039. (8) Bandy, A. R.; Ianni, J. C. J. Phys. Chem. 1998, 102, 6533–6539. (9) Re, S.; Osamura, Y.; Morokuma, K. J. Phys. Chem. A 1999, 103, 3535–3547. (10) Loerting, T.; Liedl, K. R. J. Phys. Chem. A 2001, 105, 5137– 5145. (11) Ding, C.-G.; Taskila, T.; Laasonen, K.; Laaksonen, A. Chem. Phys. 2003, 287, 7–19. (12) Standard, J. M.; Buckner, I. S.; Pulsifer, D. H. J. Mol. Struct. (THEOCHEM) 2004, 673, 1–16. (13) Havey, D. K.; Feierabend, K. J.; Vaida, V. J. Mol. Struct. (THEOCHEM) 2004, 680, 243–247. (14) Bianco, R.; Hynes, J. T. Theo. Chem. Acc. 2004, 111, 182–187. (15) Ding, C.-G.; Laasonen, K. Chem. Phys. Lett. 2004, 390, 307–313. (16) Arrouvel, C.; Viossat, V.; Minot, C. J. Mol. Struct. (THEOCHEM) 2005, 718, 71–76. (17) Bianco, R.; Wang, S.; Hynes, J. T. J. Phys. Chem. B 2005, 109, 21313–21321. (18) Ding, C.-G.; Laasonen, K.; Laaksonen, A. J. Phys. Chem. A 2003, 107, 8648–8658. (19) Natsheh, A. A.; Nadykto, A. B.; Mikkelsen, K. V.; Yu, F.; Ruuskanen, J. J. Phys. Chem. A 2004, 108, 8914–8929. (20) Nadykto, A. B.; Yu, F. Chem. Phys. Lett. 2007, 435, 14–18. (21) Nadykto, A. B.; Yu, F.; Herb., J. Chem. Phys. 2009, 360, 67–73. (22) Choe, Y.-K.; Tsuchida, E.; Ikeshoji, T. J. Chem. Phys. 2007, 126, 154510–154518. (23) Vchirawongkwin, V.; Rode, B. M. Chem. Phys. Lett. 2007, 443, 152–157. (24) Anderson, K. E.; Siepmann, J. I.; McMurry, P. H.; VandeVondele, J. J. Am. Chem. Soc. 2008, 130, 14144–14147. (25) Hammerich, A. D.; Buch, V.; Mohamed, F. Chem. Phys. Lett. 2008, 460, 423–431. (26) Vchirawongkwin, V.; Kritayakornupong, C.; Rode, B. M. J. Phys. Chem. B 2010, 114, 11561–11569. (27) Berne, B. J.; Thirumalai, D. Annu. Rev. Phys. Chem. 1986, 37, 401–424. (28) Shiga, M.; Tachikawa, M.; Miura, S. Chem. Phys. Lett. 2000, 332, 396–402. (29) Shiga, M.; Tachikawa, M.; Miura, S. J. Chem. Phys. 2001, 115, 9149–9159. (30) Kakizaki, A.; Motegi, H.; Yoshikawa, T.; Takayanagi, T.; Shiga, M.; Tachikawa, M. J. Mol. Struct. (THEOCHEM) 2009, 901, 1–8. (31) Takayanagi, T.; Takahashi, K.; Kakizaki, A.; Shiga, M.; Tachikawa, M. Chem. Phys. 2009, 358, 196–202. (32) Yoshikawa, T.; Motegi, H.; Kakizaki, A.; Takayanagi, T.; Shiga, M.; Tachikawa, M. Chem. Phys. 2009, 365, 60–68. (33) Yoshikawa, T.; Sugawara, S.; Takayanagi, T.; Shiga, M.; Tachikawa, M. Chem. Phys. Lett. 2010, 496, 14–19. (34) Sugawara, S.; Yoshikawa, T.; Takayanagi, T.; Tachikawa, M. Chem. Phys. Lett. 2011, 501, 238–244. (35) Stewart, J. J. P. J. Mol. Model 2007, 13, 1173–1213. (36) Rossi, I.; Truhlar, D. G. Chem. Phys. Lett. 1995, 233, 231–236. (37) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652.

ARTICLE

(38) Frisch, M. J. et al. Gaussian 03, revision D.02; Gaussian Inc.: Wallingford, CT, 2004. (39) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671–6687. (40) Dunning Jr, T. H.; Hay, P. J. Modern Theoretical Chemistry, 3rd; Schaefer III, H. F., ed.; Plenum: New York, 1976, 1
28. (41) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984, 80, 3265–3269. (42) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R. J. Comput. Chem. 1983, 4, 294–301. (43) Dunning, T. H., Jr J. Chem. Phys. 1989, 90, 1007–1023. (44) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796–6806. (45) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257–2261. (46) Feyman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (47) Frenkel, D.; Smit, B.; Understanding Molecular Simulation: From Algorithms to Applications, Academic Press: San Diego, CA, 2002. (48) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press; Oxford, U.K., 2005. (49) Streett, W. B.; Tildesley, D. J.; Saville, G. Mol. Phys. 1978, 35, 639–648. (50) Hoover, W. G. Phys. Rev. A 1985, 31, 1695–1697. (51) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635–2643. (52) MOPAC2007; Stewart Computational Chemistry: Colorado Springs, CO, 2007; http://OpenMOPAC.net. € Masia, M.; Kaczmarek, A.; (53) Gutberlet, A.; Schwaab, G.; Birer, O.; Forbert, H.; Havenith, M.; Marx, D. Science 2009, 324, 1545–1548. (54) Ando, K.; Hynes, J. T. J. Phys. Chem. B 1997, 101, 10464–10478. (55) Bianco, R.; Wang, S.; Hynes, J. T. J. Phys. Chem. A 2007, 111, 11033–11042. (56) Bianco, R.; Wang, S.; Hynes, J. T. J. Phys. Chem. A 2008, 112, 9467–9476. (57) Wang, S.; Bianco, R.; Hynes, J. T. J. Phys. Chem. A 2009, 113, 1295–1307. (58) Wang, S.; Bianco, R.; Hynes, J. T. Phys. Chem. Chem. Phys. 2010, 12, 8241–8249. (59) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361–9381. (60) Schwegler, E.; Grossman, J. C.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 121, 5400–5409.

11494

dx.doi.org/10.1021/jp202380h |J. Phys. Chem. A 2011, 115, 11486–11494