Nuclear quantum effect and H/D isotope effect on Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions

Abstract

Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions are fundamental and important ones in atmospheric chemistry. In this study, we focused on the nuclear quantum effect (NQE) of the hydrogen nucleus on these reactions with the aid of the multicomponent quantum mechanics (MC_QM) method, which can directly take account of NQE of light nuclei. Our study reveals that the NQE of the hydrogen nucleus lowers the activation barriers of the reactions and enhances the catalytic effects of second and third water molecules. In particular, we find that (i) the NQE of the proton removes the activation barrier of the reverse reaction of HCl + OH· → Cl· + H₂O, and (ii) the catalytic effect of the third water molecule appears in only our MC_QM calculation. We also analyze the H/D isotope effects on these reactions by using the MC_QM method.

---

1. Introduction

Water is indispensable to living organisms and also has important roles in various chemical reactions as a solvent. For instance, it is known that solvent water molecules lower the activation energy of proton transfer reactions by acting as a proton shuttle in the relay mechanism.¹ Thus, reactions involving water molecules have attracted a lot of attention in a wide variety of branches of chemistry, such as atmospheric chemistry, catalytic chemistry, surface chemistry, biochemistry, and so on.^2–4 In particular, reactions with halogen atoms have attracted special interest. For example, the reactions with Cl·, which is known as the major cause of the serious environmental problem of “ozone holes” in Antarctica,⁵ have been earnestly studied by several groups.^5–7 Recently, Li and coworkers reported the potential energy profiles for Cl· + (H₂O)_n (n = 1 and 2) reactions using the “gold-standard” CCSD(T)/cc-pVQZ method. As Li and coworkers mentioned,⁸ investigations of the reactions with a few water molecules are important as an initial step to understand chemical reactions in water at the molecular level. Therefore, it is important to analyze such small basic chemical reactions with a high-accuracy method, such as the CCSD(T)/cc-pVQZ method. They reported that the second water molecule acted as a catalyst, that is, the second water molecule lowered the activation energy for the hydrogen transfer reaction in the Cl· + H₂O reaction.

On the other hand, the importance of the nuclear quantum effect (NQE) of a light nucleus, such as a proton, has been recognized in various fields in tandem with the recent advances of experimental techniques and computational methods. The NQE of a proton provides significant effects in several reactions and systems, in which the hydrogen atom takes a central role, such as hydrogen transfer reactions, hydrogen absorption reactions, hydrogen-bonded systems, and so on. In addition, NQE is one of the main contribution factors to form low-barrier hydrogen-bonded (LBHB) systems,^9–14 which have been considered to be important for functional expression of some proteins. Also, it is known that deuterium substitution often induces significant hydrogen/deuterium (H/D) isotope effect, since a deuterium is a twice-heavier isotope of hydrogen. For example, the phase transition temperature of hydrogen-bonded ferroelectric material drastically (>100 K) changes by replacing a hydrogen-bonded hydrogen atom with deuterium,¹⁵ and the reaction rate constant of 1,5-sigmatropic hydrogen transfer reaction in 1,3-pentadiene is more than 10 times faster than the transfer in deuterated system.² In addition, deuterium substitution of solvent molecules (solvent H/D isotope effect) is also important for molecular properties. For example, the critical temperature of lens protein of γB-crystallin in D₂O is 16 K higher than that in H₂O.¹⁶ Since Cl· + (H₂O)_n reactions are important as an initial step to understand chemical reactions in solution, as mentioned above, we believe that the analyses of the Cl· + (D₂O)_n reactions are also important toward understanding the chemical reactions involving heavy water.

However, the conventional quantum mechanical (QM) calculations usually can not take account of nuclear motion based on the Born–Oppenherimer (BO) approximation.¹⁷ In the framework of BO approximation, only electronic Schrödinger equation is solved under the field of clamped nuclear charges. Therefore, it is difficult for the conventional QM methods to directly reflect NQE on electronic structure. On the other hand, we have recently proposed multicomponent QM (MC_QM) methods,^18–21 which can directly reflect NQE of light nuclei on electronic structures. We have already successfully analyzed H/D isotope effect on geometries (H/D geometrical isotope effect: GIE) in various hydrogen-bonded systems and H/D kinetic isotope effect (KIE) in several hydrogen transfer reactions using our MC_QM methods.^22,23 In addition, Ishimoto and Koyama successfully revealed the dynamic behavior of nuclear wavefunction of hydrogen nucleus in H₅O₂⁺ and its isotopomers by performing molecular dynamics simulation on the MC_QM effective potential energy hypersurface.²⁴ Quite recently, we have proposed MC_QM-climbing image-nudged elastic band (CI-NEB) method¹⁸ by combining our MC_QM method with CI-NEB method,²⁵ and have successfully analyzed H/D isotope effect in proton/hydrogen transfer reactions^18,22,23 and H/D isotope effect in isomerization reactions²³ including NQE of light nuclei with the aid of MC_QM-CI-NEB method.

Although CCSD(T) calculation can predict electronic energy and molecular geometries quite accurately, NQEs of light nuclei are not taken into account within the framework of BO approximation. For the systems and the reactions, in which NQEs of light nuclei are remarkable, we have to not only improve the accuracy of electronic structure calculations but also adequately include NQE of light nuclei beyond the BO approximation. Thus, the main subject of the present study is to reveal NQE of hydrogen nuclei on Cl· + (H₂O)_n (n = 1–3) reactions using MC_QM-CI-NEB method. First, we investigate the performance of several major exchange–correlation density functionals to obtain the reliable optimized geometries of the stationary points structures in Cl· + (H₂O)_n (n = 1–3) reactions. Because one of the efficient ways to obtain reliable energies is calculating CCSD(T) energies at the geometries optimized by DFT method with the appropriate exchange–correlation functional. Although it is known that DFT calculations can provide reasonable energy and molecular geometries, the accuracy of DFT calculations strongly depends on the selection of exchange–correlation functional. We, thus, try to determine the best exchange–correlation functional for Cl· + (H₂O)_n (n = 1–3) reactions by calculating CCSD(T) energy at the geometry optimized by DFT (denoted as CCSD(T)//DFT), and compare it with the energy obtained by pure CCSD(T) calculation. Then the NQE on Cl· + (H₂O)_n and Cl· + (D₂O)_n (n = 1–3) reactions have been analyzed using MC_CCSD(T) and MC_DFT calculations.

In the next section, we briefly introduce our MC_QM method and CI-NEB method. Computational details to analyze Cl· + (H₂O)_n (n = 1–3) reactions are given in Sec. 3. The computational details for MC_QM and MC_QM-CI-NEB calculations are also given in Sec. 3. Performance of the several density functionals to reproduce CCSD(T) stationary point geometries is discussed in Sec. 4-1. The NQE and H/D isotope effect on the reactions are discussed in Sec. 4-2. Finally, some concluding remarks are given in Sec. 5.

2. Theory

2.1. Multicomponent quantum mechanics (MC_QM) method

In this section, we would like to introduce our MC_QM method briefly (see ref. 19–21 and references therein for more detailed information). In our MC_QM approach, the total Hamiltonian for system containing N_e-electrons, N_p-quantum nuclei, and M-classical nuclei, is expressed aswhere the indices of i and j refer to the electrons, p and q to the quantum nuclei, and A and B to the classical nuclei. In addition, Z_A represents nuclear charge of Ath nucleus, and M_p is the mass of pth quantum nucleus. The first to third terms are the conventional terms for electrons, and the fourth to sixth terms are those for quantum nuclei. The seventh term represents the coulombic interaction between electron and quantum nucleus, and the last term represents the classical nuclear repulsion.

The total wavefunction Ψ₀ is represented by a simple product of electronic and nuclear wavefunction in the Hartree–Fock level of MC_QM method,

Ψ₀ = Φ^e₀Φ^p₀, (2)

where Φ^e₀ and Φ^p₀ are the electronic and the nuclear wavefunctions, respectively. The effective one-particle (one-electron and one-nucleus) operators arewhere ĥ_e and ĥ_p are one-particle operators for electron and quantum nucleus, Ĵ_e and are the electronic Coulomb and exchange operators, and Ĵ_p is the Coulomb operator for quantum nuclei. The nuclear exchange term is ignored in this study due to its small contribution.²⁶ To solve the electronic Fock equation, we use linear combination of Gaussian-type functions (LCGTFs) for expanding electronic molecular orbitals (MOs),where C^e_μi is the LCGTF coefficients for an electron, and χ is the Gaussian-type function (GTF). Although quantum nuclear MOs should be expanded using a suitable number of GTFs as well as the electronic MOs, we have already demonstrated that nuclear quantum effect is adequately taken into account by using only single s-type GTF as nuclear basis function.^18,20,22,23

It should be noted here that new types of many-body effects, such as electron–quantum nucleus and quantum nucleus–quantum nucleus correlation effects arise in MC_QM scheme. To improve accuracy, we have to evaluate these correlation effects as well as the conventional electron–electron one. However, it is worth mentioning that the contributions from electron–quantum nucleus and quantum nucleus–quantum nucleus correlations are much smaller than that from electron–electron one for typical molecular systems, and we have shown that H/D isotope effects in various hydrogen-bonded systems and hydrogen transfer reactions can be analyzed evaluating only electron–electron correlation effect in MC_QM calculations.^18–23 In particular, we have already demonstrated that H/D isotope effects in similar reactions of F· + (H₂O)_n → HF + OH·(H₂O)_n−1 (n = 1–3) can be adequately analyzed even if only electron–electron correlation effect is evaluated using established density functionals for electrons and CCSD(T) method.²² Therefore, we also evaluate only electron–electron correlation effect in the present MC_QM calculations.

2.2. Climbing image-nudged elastic band (CI-NEB) method

To find out the transition state (first order saddle point) structure on MC_QM effective potential energy surface (PES), we have carried out MC_QM-CI-NEB calculations. Here, we would like to briefly introduce the algorithm of the CI-NEB method.^25,27,28

In the NEB calculations, the total force acting on an image is calculated as the sum of the spring force along the tangent and the true force perpendicular to the tangent

F_i = F^s_i|_|| − ∇V(R_i)|_⊥. (6)

The first term of the right-hand side of eqn (6) corresponds to the spring force, which is given as

F^s_i|_|| = k(|R_i+1 − R_i| − |R_i − R_i−1|)[small tau, Greek, circumflex]_i, (7)

where k is the spring constant. We adopted the variable spring constants²⁷ to efficiently improve the resolution of the vicinity of the transition statewhere Δk = k_max − k_min, V_max is the highest energy for the whole elastic band, V_i is the higher energy of the two images connected by spring i. According to Henkelman's original paper,²⁷ we chose the energy of the higher energy endpoint of the MEP as V_ref.

The second term of the right-hand side of eqn (6) is given as

∇V(R_i)|_⊥ = ∇V(R_i) − ∇V(R_i)[small tau, Greek, circumflex]_i[small tau, Greek, circumflex]_i. (9)

[small tau, Greek, circumflex] in eqn (7) and (9) is the unit vector along the reaction path. We used Henkelman's revised tangent²⁸ in this study. The revised tangents are defined aswhere R_i and V_i are the position vector and energy of i-th image, respectively. For the image at an energy minimum (V_i+1 > V_i < V_i−1) or at an energy maximum (V_i+1 < V_i > V_i−1), the tangent becomeswhereand

For the image with the highest energy, the total force is evaluated according to the following expression instead of eqn (6)

F_imax = ∇V(R_imax) − 2∇V(R_imax)[small tau, Greek, circumflex]_imax[small tau, Greek, circumflex]_imax. (14)

It should be noted here that we can obtain both minimum energy path and transition state structure on the MC_QM effective PES by minimizing (CI-)NEB forces in eqn (6) and (14) using MC_QM method. We have already demonstrated that MC_QM-CI-NEB approach adequately provides effective transition state structure on the MC_QM PES for several reactions.^18,22,23

3. Computational detail

To investigate the performance of density functionals for geometry optimization, all stationary points of Cl· + (H₂O)_n → HCl + OH(H₂O)_n−1 (n = 1–3) reactions were optimized using five major density functionals (B3LYP, CAM-B3LYP, M06, ωB97XD, and MPW1K) with correlation-consistent cc-pVTZ electronic basis set. Then, we calculated CCSD(T)/cc-pVQZ energy at the geometry optimized by DFT (CCSD(T)/cc-pVQZ//DFT/cc-pVTZ). It is known that diffuse functions are often important to describe hydrogen-bonded systems. Thus, first, we investigated the importance of diffuse functions for calculating relative potential energy profile for Cl· + H₂O reaction. The calculated relative energies and the optimized geometrical parameters obtained with cc-pVXZ and aug-cc-pVXZ (X = T or Q) electronic basis sets are shown in Fig. S1 in ESI.† The activation energy obtained by CCSD(T)/aug-cc-pVQZ//ωB97XD/aug-cc-pVTZ calculations (20.7 kcal mol⁻¹) is almost equal to that by CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculations (20.6 kcal mol⁻¹). Therefore, the effect of diffuse functions is also negligible in Cl· + H₂O reactions, as in the case of F + H₂O ones.²²

Next, to analyze HCl + OH·(H₂O)_n−1 → Cl· + (H₂O)_n and DCl + OD·(D₂O)_n−1 → Cl· + (D₂O)_n (n = 1–3) reactions including NQE of proton and deuteron, we performed MC_CCSD(T)/cc-pVQZ//MC_DFT/cc-pVTZ calculations. The effective transition state structure on MC_QM effective potential energy hypersurface is obtained by MC_DFT-climbing image-nudged elastic band (CI-NEB)¹⁸ calculations. In MC_QM calculations, all protons and deuterons were treated as quantum wavefunction. We adopted single s-type Gaussian-type functions (GTFs), exp{−α(r − R)²}, as nuclear wavefunction. The α value in GTF is the orbital exponent value and determines the spatial distribution of nuclear wavefunction. We used average α values (a^H_ave = 24.1825 and a^D_ave = 35.6214) for geometry optimizations.^29–31 To refine the spatial distribution of nuclear wavefunctions and the total energy, we optimized α (α_opt) values by MC_ωB97XD/cc-pVQZ//MC_ωB97XD/cc-pVTZ calculations. Then we calculated the more reliable MC_CCSD(T)/cc-pVQZ//MC_ωB97XD/cc-pVTZ energies with these α_opt values. We have confirmed that α_opt values obtained by MC_ωB97XD/cc-pVQZ and MC_CCSD(T)/cc-pVQZ calculations are similar to each other and MC_CCSD(T)/cc-pVQZ energies obtained using MC_ωB97XD/cc-pVQZ α_opt values and using MC_CCSD(T)/cc-pVQZ α_opt values are also quite similar to each other (see, Fig. S2 in ESI†). All calculations were performed with the modified version of GAUSSIAN 09 program package.³²

4. Results and discussion

4.1. DFT calculations for Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions

To determine the best density functional for the reactions of Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3), we optimized the stationary points structures of Cl· + H₂O → HCl + OH· reaction using five major density functionals (B3LYP, CAM-B3LYP, M06, ωB97XD, and MPW1K) and compared the optimized interatomic distances in DFT-optimized geometries with those in CCSD(T)/cc-pVQZ-optimized geometry.

Table 1 lists the optimized geometrical parameters in the stationary points structures of Cl· + H₂O → HCl + OH· reaction obtained by DFT/cc-pVTZ calculations. Guo's CCSD(T)/cc-pVQZ values⁶ are also listed in Table 1 for comparison. We can find that all five density functionals can reproduce the CCSD(T)/cc-pVQZ interatomic distance well, except for transition state. The O⋯H1 distances in transition state are 1.515 Å and 1.413 Å in B3LYP and M06 geometries, respectively. Both of them were clearly longer than that in CCSD(T) geometry (1.324 Å). Also, Cl⋯H1 distance of 1.345 Å in B3LYP geometry is slightly shorter than that in CCSD(T) one (1.391 Å).

Table 1 Optimized geometrical parameters (Å) of the stationary point structures of Cl· + H₂O → HCl + OH· reaction obtained by DFT/cc-pVTZ methods

Method	Reactant	Entrance complex			Transition state			Exit complex			Products
	R(O–H1)	R(O–H1)	R(O⋯Cl)	R(Cl⋯H1)	R(O–H2)	R(O⋯H1)	R(Cl⋯H1)	R(O–H2)	R(O⋯H1)	R(Cl⋯H1)	R(O–H2)	R(Cl–H1)
a Ref. 6.
B3LYP	0.961	0.965	2.508	2.653	0.975	1.515	1.325	0.976	1.979	1.297	0.974	1.283
CAM-B3LYP	0.960	0.963	2.482	2.681	0.973	1.396	1.368	0.975	1.959	1.293	0.973	1.280
M06	0.957	0.962	2.524	2.678	0.971	1.413	1.371	0.972	1.988	1.297	0.970	1.283
ωB97XD	0.957	0.960	2.501	2.698	0.970	1.386	1.374	0.971	1.977	1.294	0.969	1.280
MPW1K	0.950	0.953	2.498	2.708	0.963	1.303	1.393	0.964	1.979	1.285	0.962	1.273
CCSD(T)/cc-pVQZ^a	0.958	0.960	2.601	2.845	0.971	1.324	1.391	0.971	2.045	1.285	0.970	1.277

---

Next, we evaluated CCSD(T)/cc-pVQZ energy at the geometry optimized with DFT/cc-pVTZ (CCSD(T)/cc-pVQZ//DFT/cc-pVTZ) to obtain high-accuracy energy. The relative energy diagram for the reaction is depicted in Fig. 1. We calculated the CCSD(T)/cc-pVQZ energy at only the CAM-B3LYP, ωB97XD, and MPW1K/cc-pVTZ optimized geometries, since B3LYP and M06 calculations cannot reproduce transition state structure well, as mentioned above. As shown in Fig. 1, all of the CCSD(T)/cc-pVQZ//DFT/cc-pVTZ calculations reproduced CCSD(T)/cc-pVQZ potential energy diagram for Cl· + H₂O → HCl + OH· reaction well. Difference between the pure CCSD(T) and the CCSD(T)//DFT relative energies are less than 0.3 kcal mol⁻¹ for all stationary point structures.

Fig. 1 The potential energy profile (kcal mol⁻¹) and stationary point structures for Cl· + H₂O → HCl + OH· reaction obtained by CCSD(T)/cc-pVQZ//DFT/cc-pVTZ calculations. Relative energies are calculated as the energy difference between the energy of each stationary point structure and sum of energies of Cl radical and H₂O (Reactants).

Next, we would like to focus on Cl· + (H₂O)₂ reaction. Fig. 2 shows the optimized structures and Table 2 and 3 list the optimized geometrical parameters of the entrance complex and transition state for the reaction, respectively. The optimized geometrical parameters of other stationary point structures are shown in Tables S1–S3 in ESI.† For Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction, intermolecular O2⋯H2 distance of entrance complex with B3LYP and CAM-B3LYP calculation is 0.066 Å and 0.081 Å shorter than that with CCSD(T)/cc-pVQZ, respectively. In addition, intermolecular Cl⋯H3 distances with B3LYP, CAM-B3LYP, and M06 calculations are 0.057–0.100 Å shorter than that with CCSD(T)/cc-pVQZ. Thus, these functionals are found to overestimate the interaction energy between Cl(H₂O) and the second water molecule. On the other hand, the covalent O1–H2 bond lengths in transition state with B3LYP and M06 calculations are 0.016 Å and 0.014 Å longer than that with CCSD(T)/cc-pVQZ, respectively. Meanwhile, the O1–H2 bond length in MPW1K transition state geometry is 0.987 Å, which is almost the same with CCSD(T)/cc-pVTZ value (0.988 Å). However, this is only the accidental coincidence, since MPW1K 0.008 Å underestimate the covalent O–H bond length in water monomer molecule than CCSD(T) (see, reactant in Table 1). Also, B3LYP, M06, and MPW1K calculations gave slightly (about 0.06 Å) shorter intermolecular O2⋯H2 distance compared to CCSD(T). Therefore, we conclude that ωB97XD functional is the best choice to reproduce the CCSD(T)/cc-pVQZ stationary point structures for Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction.

Fig. 2 The potential energy profile (kcal mol⁻¹) and stationary point structures for Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction obtained by CCSD(T)/cc-pVQZ//DFT/cc-pVTZ and CCSD(T)/cc-pVQZ⁷ calculations. Relative energies are calculated as the energy difference between the energy of each stationary point structure and sum of energies of Cl radical and (H₂O)₂ (reactants).

Table 2 Optimized geometrical parameter (Å) of entrance complex for Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction obtained by DFT/cc-pVTZ methods

Method	R(O1–H1)	R(O1–H2)	R(O2–H3)	R(O2–H4)	R(O1⋯Cl)	R(O2⋯H2)	R(Cl⋯H1)	R(Cl⋯H2)	R(Cl⋯H3)
a Ref. 7.
B3LYP	0.965	0.987	0.968	0.962	2.453	1.765	2.563	2.607	2.654
CAM-B3LYP	0.963	0.983	0.966	0.961	2.410	1.750	2.655	2.590	2.554
M06	0.962	0.980	0.965	0.959	2.476	1.806	2.676	2.622	2.511
ωB97XD	0.960	0.979	0.963	0.958	2.431	1.784	2.659	2.599	2.613
MPW1K	0.953	0.972	0.956	0.951	2.399	1.769	2.650	2.584	2.614
CCSD(T)/cc-pVQZ^a	0.961	0.974	0.963	0.959	2.452	1.831	2.725	2.637	2.611

---

Table 3 Optimized geometrical parameter (Å) of transition state for Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction obtained by DFT/cc-pVTZ methods

Method	R(O1–H2)	R(O2–H3)	R(O2–H4)	R(O1⋯H1)	R(O2⋯H2)	R(Cl⋯H1)	R(Cl⋯H3)
a Ref. 7.
B3LYP	1.004	0.969	0.962	1.268	1.710	1.490	2.599
CAM-B3LYP	0.996	0.966	0.961	1.344	1.731	1.406	2.677
M06	1.002	0.969	0.958	1.240	1.715	1.507	2.452
ωB97XD	0.992	0.962	0.958	1.343	1.764	1.407	2.740
MPW1K	0.987	0.956	0.951	1.262	1.717	1.435	2.739
CCSD(T)/cc-pVQZ^a	0.988	0.963	0.959	1.293	1.777	1.420	2.684

---

The calculated potential energy diagram for Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction is shown in Fig. 2. Based on the above results, we calculated CCSD(T)/cc-pVQZ energies at the ωB97XD/cc-pVTZ-optimized geometries. Meanwhile, Li and coworkers claimed that MPW1K is the best choice to analyze a similar reaction of F· + H₂O → HF + OH·.³³ Thus, we also calculated CCSD(T)/cc-pVQZ energies at the geometries optimized by MPW1K for comparison. Both of CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ and CCSD(T)/cc-pVQZ//MPW1K calculations reproduced CCSD(T)/cc-pVQZ potential energy diagram⁷ well. It should be noted here that in F· + (H₂O)₂ → HCl + OH·(H₂O) reaction, although the CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculation reproduced the CCSD(T)/cc-pVQZ potential energy diagram well, the CCSD(T)/cc-pVQZ//MPW1K/cc-pVTZ calculation failed to reproduce CCSD(T)/cc-pVQZ potential energy diagram even qualitatively.²² Therefore, the performance of exchange–correlation density functional is different even for similar reactions of X· + (H₂O)_n (X = F and Cl).

The relative energy of transition state lies about 16 kcal mol⁻¹ above the sum of energies of separate reactants, and is about 4 kcal mol⁻¹ lower than that of the reaction with water monomer (Fig. 1). Thus, the second water molecule lowers the barrier of Cl· + H₂O → HCl + OH· reaction, and acts as a catalyst as is the case of F· + (H₂O)₂ reaction.²² Since the catalytic effect of the second water molecule observed in CCSD(T)/cc-pVQZ calculation is adequately reproduced by CCSD(T)/cc-pVQZ//DFT/cc-pVTZ calculation, evaluating CCSD(T)/cc-pVQZ energy at the geometry optimized DFT/cc-pVTZ is a good way to analyze Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 reactions efficiently. In addition, the relative energy of transition state is 1.2 kcal mol⁻¹ lower than products in CCSD(T)/cc-pVQZ results. In other words, the reverse reaction, HCl + OH·(H₂O) → Cl· + (H₂O)₂, has negative activation barrier. Such negative activation barrier of the reverse reaction is also adequately predicted by CCSD(T)/cc-pVQZ//DFT/cc-pVTZ calculations.

The potential energy diagram and optimized structures for Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ reaction are shown in Fig. 3, and the optimized geometrical parameters of each stationary point structure are listed in Tables S4–S8 in ESI.† As far as we know, CCSD(T)/cc-pVQZ study on Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ reaction has not been reported so far. Since CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculations gave the reliable potential energy diagrams for the reaction of Cl· with H₂O and (H₂O)₂, as mentioned above, we would like to analyze the reaction of Cl· with (H₂O)₃ in the same manner.

Fig. 3 The potential energy profile (kcal mol⁻¹) and stationary point structures for Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ reaction obtained by CCSD(T)/cc-pVQZ//DFT/cc-pVTZ calculations. Relative energies are calculated as the energy difference between the energy of each stationary point structure and sum of energies of Cl radical and (H₂O)₃ (reactants).

In Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ reaction, the relative energy of entrance complex and products are slightly higher than those in Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction, whereas that of transition state is unchanged by existence of the third water molecule. Therefore, the third water molecule hardly act as a catalyst in Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ reaction. It should be mentioned that the relative energy of transition state in F· + (H₂O)₃ → HF + OH·(H₂O)₂ reaction is 1.5 kcal mol⁻¹ lower than that in the reaction of F· + (H₂O)₂. Thus, unlike Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ reaction, the catalytic effect of the third water molecule is observed in F· + (H₂O)₃ → HF + OH·(H₂O)₂ reaction.²²

Let us briefly summarize the results of the conventional DFT and CCSD(T)//DFT calculations of Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions. First of all, we have investigated the performance of several density functionals for geometry optimization calculations of the reactions of Cl· with (H₂O)_n (n = 1 and 2). The B3LYP and M06 calculations tended to overestimate the interatomic O1⋯H1 distance in transition state of Cl· + H₂O → HCl + OH· reaction, and the MPW1K calculation slightly underestimated covalent O–H bond length in H₂O monomer. Meanwhile, CAM-B3LYP/cc-pVTZ calculations gave shorter R(O2⋯H2) and R(Cl⋯H3) interatomic distances compared to CCSD(T)/cc-pVQZ ones in entrance complex of Cl· + (H₂O)₂ → HCl + OH·(H₂O) reaction. The CAM-B3LYP calculations, thus, overestimated the interaction between Cl(H₂O) and the second water molecule. Therefore, our results indicate that the ωB97XD/cc-pVTZ method is the best choice to reproduce the geometries of CCSD(T)/cc-pVQZ stationary point structure in Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions. The potential energy diagrams obtained by CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculations are indeed close to those obtained by CCSD(T)/cc-pVQZ calculations.^6,7 Therefore, we have clearly demonstrated that the reliable high-accuracy energies can be obtained by evaluating CCSD(T)/cc-pVQZ energy at the geometries optimized by ωB97XD/cc-pVTZ instead of calculating time-consuming CCSD(T)/cc-pVQZ geometry optimization calculation. Moreover, we have found that the third water molecule hardly catalyze the reaction unlike the second water molecule in Cl· + (H₂O)₂ reaction and the third water molecule in F· + (H₂O)₃ reaction.²²

4.2. Nuclear quantum effect and H/D isotope effect on Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions

Next, we performed MC_CCSD(T)/cc-pVQZ//MC_ωB97XD/cc-pVTZ calculations to analyze Cl· + H₂O → HCl + OH· reaction including NQE of proton. The optimized structures and the potential energy diagram are shown in Fig. 4 (quantum H). We can find that quantum mechanical treatment of proton affects geometrical parameters. For example, the covalent O–H2 and Cl–H1 bond lengths in exit complex of quantum H are 0.995 Å and 1.323 Å, respectively, whereas those are 0.971 Å and 1.294 Å in the conventional ωB97XD result (denoted as classical H). Fig. 4 also shows the results of the reaction with D₂O (quantum D). In quantum D results, the covalent O–D2 and Cl–D1 bond lengths in exit complex are 0.988 Å and 1.314 Å, respectively, which are shorter than the counterpart of quantum H. These changes of covalent bond lengths are due to the direct inclusion of the anharmonicity of potential energy curve along the covalent bond direction by MC_QM method. On the other hand, the intermolecular O1⋯D1 distance (1.931 Å) in quantum D is slightly longer than the intermolecular O1⋯H1 distance (1.910 Å) in quantum H. We have already revealed that the elongation of intermolecular distances can be understood by the difference between localized nature of protonic and deuteronic wavefunctions.²⁰ Since deuteron is twice heavier than proton, deuteronic wavefunction is more localized compared to protonic one. Consequently, deuteronic wavefunction more strongly attracts surrounding electrons, and deuterium becomes less positive and acts as weaker hydrogen nucleus-donor than hydrogen. In fact, the optimized orbital exponent α values for H1 and H2 of exit complex are 20.6 and 22.8, and those for D1 and D2 are 30.9 and 34.1, respectively. (see, Fig. S3 in ESI†). Since larger α value represents more localized nuclear wavefunction, deuteronic wavefunctions are indeed more localized than protonic ones in our MC_QM calculations. In addition, the NBO charges on proton and deuteron are 0.154 and 0.150 in quantum H and quantum D, respectively. Deuteron is indeed less positive than proton.

Fig. 4 The potential energy profile (kcal mol⁻¹) and the optimized geometrical parameters (Å) for Cl· + H₂O → HCl + OH· and Cl· + D₂O → DCl + OD· reaction obtained by MC_CCSD(T)/cc-pVQZ//MC_ωB97XD/cc-pVTZ and CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculations. Relative energies are calculated as the energy difference between the energy of each stationary point structure and sum of energies of Cl radical and H₂O or D₂O (reactants).

Next, we would like to focus on NQE on the potential energy profiles of Cl· + H₂O → HCl + OH· reaction (Fig. 4). The relative energy of quantum H is always lower than that of quantum D and classical H. In other words, NQE of proton lowered the relative energies of all stationary point structures. In particular, larger stabilization by NQE is observed in transition state and exit complex due to the hydrogen-bonded interaction. Larger stabilization appears in quantum H rather than in quantum D, because NQE of heavier nucleus (deuteron) is smaller than that of lighter one (proton). In particular, the relative energy of transition state of quantum H is 0.4 kcal mol⁻¹ lower than products. The reverse reaction HCl + OH· → Cl· + H₂O, thus, has a negative activation barrier when considering NQE of proton, even only one water molecule participate in the reaction. We should note here that harmonic zero-point vibration energy correction (also shown in Fig. 4) lower not only transition state but also products. The harmonic ZPVE correction, thus, did not remove the activation barrier of the reverse reaction.

To analyze the stabilization by NQE in detail, we additionally performed the MC_CCSD(T)/cc-pVQZ energy calculation at the geometry optimized by ωB97XD/cc-pVTZ. (MC_CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ) for Cl· + H₂O → HCl + OH· reaction. In such calculation, the geometrical relaxation induced by NQE is not taken into account, whereas the nuclear quantum fluctuation of proton is represented in MC_CCSD(T) energy calculation. The calculated energy diagram is shown in Fig. S4 in ESI.† Although the relative energies obtained in MC_CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculations are slightly (0.5–0.7 kcal mol⁻¹) higher than those obtained in MC_CCSD(T)/cc-pVQZ//MC_ωB97XD/cc-pVTZ calculations, the potential energy profiles obtained in MC_CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ and MC_CCSD(T)/cc-pVQZ//MC_ωB97XD/cc-pVTZ calculations are almost parallel to each other. In addition, the activation barrier for the reverse reaction is also negative in MC_CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ result. These results show that the direct quantum mechanical treatment of light nuclei by appropriate method, such as our MC_QM, is indispensable to analyze potential energy profiles for Cl· + H₂O → HCl + OH· reaction adequately.

The optimized structure and potential energy profiles for Cl· + (H₂O)₂ → HCl + OH·(H₂O) and Cl· + (D₂O)₂ → DCl + OD·(D₂O) reactions are shown in Fig. 5. The covalent bond lengths and the intermolecular distances are the shortest and the longest in quantum H, respectively, as well as the reaction with one H₂O molecule. Meanwhile, NQE of proton and deuteron seemingly did not affect the intermolecular O1⋯H1 and O1^…D1 distances in transition state, since the intermolecular O1⋯H1(D1) distances in transition state in quantum H, quantum D, and classical H are very similar to each other. Since H1 atom is located in the almost midpoint between O1 and Cl, H1 atom is loosely bounded to both O1 and Cl in transition state. Hence, the O1⋯H1 distances do not show the typical H/D geometrical isotope effect in transition state. We would like to advocate that NQE of proton and deuteron also provide important effect on transition state because other parameters, such as O1–H2, O2–H3, O2–H4, O2⋯H2, and Cl⋯H3, are affected by NQE, and the spatial distributions of protonic and deuteronic wavefunctions are different from each other (see, optimized exponent α_opt value in Cl· + (H₂O)₂ → HCl + OH·(H₂O) and Cl· + (D₂O)₂ → DCl + OD·(D₂O) reactions in Fig. S5 in ESI†).

Fig. 5 The potential energy profile (kcal mol⁻¹) and the optimized geometrical parameters (Å) for Cl· + (H₂O)₂ → HCl + OH·(H₂O) and Cl· + (D₂O)₂ → DCl + OD·(D₂O) reaction obtained by MC_CCSD(T)/cc-pVQZ//MC_ωB97XD//cc-pVTZ and CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculations. Relative energies are calculated as the energy difference between the energy of each stationary point structure and sum of energies of Cl radical and (H₂O)₂ or (D₂O)₂ (reactants).

Next, we would like to focus on the potential energy profile of Cl· + (H₂O)₂ → HCl + OH·(H₂O) and Cl· + (D₂O)₂ → DCl + OD·(D₂O) reactions. As in the case of classical H, the second water molecule lowered the relative energies of all stationary points compared to the reaction with one H₂O (or D₂O) molecule. NQEs of proton and deuteron stabilize all the stationary point structures in quantum H and quantum D. The relatively large stabilization effects appear in transition state and exit complex due to the number of hydrogen-bonded interactions, as well as the reactions with one water molecule. Consequently, NQE of proton and deuteron lower the activation barrier of Cl· + (H₂O)₂ → HCl + OH·(H₂O) and Cl· + (D₂O)₂ → DCl + OD·(D₂O) reactions.

The optimized structures and potential energy profiles for Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ and Cl· + (D₂O)₃ → DCl + OD·(D₂O)₂ are shown in Fig. 6 (see, Fig. S6 in ESI† for the optimized exponent α_opt values). NQE also impacts the geometrical parameters of the stationary point structures, like the reactions with one or two H₂O (or D₂O) molecule(s). The covalent bond length and interatomic distances were elongated and shortened by including NQE, respectively.

Fig. 6 The potential energy profile (kcal mol⁻¹) and the optimized geometrical parameters (Å) for Cl· + (H₂O)₃ → HCl + OH·(H₂O)₂ and Cl· + (D₂O)₃ → DCl + OD·(D₂O)₂ reaction obtained by MC_CCSD(T)/cc-pVQZ//MC_ωB97XD/cc-pVTZ and CCSD(T)/cc-pVQZ//ωB97XD/cc-pVTZ calculations. Relative energies are calculated as the energy difference between the energy of each stationary point structure and sum of energies of Cl radical and (H₂O)₃ or (D₂O)₃ (reactants).

Due to the difference of the number of hydrogen-bonded interactions in the structure, the greater stabilizations are again found in transition state and exit complex rather than in entrance complex. As mentioned in Sec. 4-1, the relative energies of transition state in Cl· + (H₂O)₂ and Cl· + (H₂O)₃ reactions of classical H are the same (16.3 kcal mol⁻¹). The catalytic effect of the third water molecule is, thus, not observed in classical H calculation. However, the relative energies of transition state in quantum H and quantum D are 1.7 kcal mol⁻¹ and 1.1 kcal mol⁻¹ lower than those in the reactions with two water molecules, respectively. Therefore, we can find that the clear catalytic effect of the third water molecule appears when NQEs of proton and deuteron are adequately taken into account. The potential energy diagrams of Cl· + (H₂O)₂ and Cl· + (H₂O)₃ reactions corrected for ZPVE are shown in Fig. S7 and S8 in ESI† for comparison. The ZPVE-corrected relative energy of transition state in Cl· + (H₂O)₃ reaction is lower than that in Cl· + (H₂O)₂ one. However, ZPVE correction destabilizes the relative energy of entrance complex. This tendency is opposite to the MC_QM results. Thus, we claim to that the geometrical relaxation effect is also important to adequately analyze the catalytic effect of water molecule(s) on the reaction.

Let us summarize NQE on the reactions of Cl· + (H₂O)_n (n = 1–3) briefly. We have analyzed the influence of NQE of proton and deuteron on the stationary point structures and potential energy profiles of Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 and Cl· + (D₂O)_n → DCl + OD·(D₂O)_n−1 (n = 1–3) reactions using MC_QM method. We have clearly demonstrated that NQEs of hydrogen nuclei lower the relative energies of all stationary point structures in Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) and its deuterated reactions. Furthermore, the activation barriers for these reactions (relative energy of transition state) are also lowered by NQE. In particular, the relative energy of transition state of Cl· + H₂O → HCl + OH· reaction lies below the separated products molecules in quantum H. Thus, the reverse reaction of Cl· + H₂O → HCl + OH· has negative activation barrier when considering NQE of proton. In addition, our calculation suggests that the catalytic effect of the third water molecule only appears when NQEs of proton and deuteron are adequately taken into account.

5. Conclusion

In this study, we analyzed Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 and Cl· + (D₂O)_n → DCl + OD·(D₂O)_n−1 (n = 1–3) reactions using MC_QM method, which can directly reflect NQEs of light nuclei on electronic states. First, we investigated the performance of five major density functionals to efficiently reproduce the CCSD(T)/cc-pVQZ structure. As a result, ωB97XD/cc-pVTZ calculation gave the best structure among these five functionals. Next, we investigated the NQE of proton and deuteron on the optimized structures and energy diagram of the reactions using MC_QM method, and find that NQE of hydrogen nuclei lower the relative energies of the stationary point structures. In particular, NQE significantly lowers the activation barrier of Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions, and consequently the relative energy of transition state in Cl· + H₂O → HCl + OH· is lower than that of separated products when NQE of proton is directly considered. Furthermore, although the conventional DFT calculations gave the same activation energies for the reactions with (H₂O)₂ and (H₂O)₃, the activation energy in the reaction with (H₂O)₃ is lower than that with (H₂O)₂ in the results obtained with MC_QM method. In other words, the catalytic effect of the third water molecule is observed in only MC_QM calculations. Therefore, we clearly demonstrate the importance of direct inclusion of protonic and deuteronic quantum natures to analyze Cl· + (H₂O)_n → HCl + OH·(H₂O)_n−1 (n = 1–3) reactions, as well as the reactions of F· + (H₂O)_n.²²

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work was supported by JSPS KAKENHI Grants No. 16K17851 (to T. U.) and 15KT0067 (to M. T. and T. U.). Portion of these computations was performed at Research Center for Computational Science (RCCS), Okazaki.

References

1. A. Petrone, P. Cimino, G. Donati, H. P. Hratchain, M. J. Frisch and N. Rega, J. Chem. Theory Comput., 2016, 12, 4925.
2. W. R. Roth and J. König, Justus Liebigs Ann. Chem., 1966, 113, 24.
3. W. Yin, S. Yan, M. Qin, Z. Li and Y. Bu, Int. J. Quantum Chem., 2006, 106, 1528–1543.
4. E. Brunk, J. S. Arey and U. Rothlisberger, J. Am. Chem. Soc., 2012, 134, 8608–8616.
5. R. J. Buzerk, J. R. Barker and J. S. Francisco, J. Phys. Chem. A, 2012, 116, 4712.
6. Y. Guo, M. Zhang, Y. Xie and H. F. Schaefer III, J. Chem. Phys., 2013, 139, 041101.
7. G. Li, H. Wang, Q. Li, Y. Xie and H. F. Schaefer III, J. Chem. Phys., 2014, 140, 124316.
8. G. Li, Y. Xie and H. F. Schaefer III, Chem. Phys. Lett., 2016, 648, 1.
9. K. Saito and H. Ishikita, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 167.
10. K. Saito and H. Ishikita, Biochemistry, 2012, 51, 1171.
11. Y. Kanematsu and M. Tachikawa, J. Chem. Phys., 2014, 141, 185101.
12. M. Nadal-Ferret, R. Gelabert, M. Moreno and J. M. Lluch, J. Am. Chem. Soc., 2014, 136, 3542.
13. D. A. Nichols, J. C. Hargis, R. Sanishvili, P. Jaishankar, K. Defrees, E. W. Smith, K. K. Wang, F. Prati, A. R. Renslo, H. L. Woodcock and Y. Chen, J. Am. Chem. Soc., 2015, 137, 8086.
14. T. Graen, L. Inhester, M. Clemens, H. Grubmüller and G. Groenhof, J. Am. Chem. Soc., 2016, 138, 16620.
15. N. Dalal, A. Klymachyov and A. Bussmann-Holder, Phys. Rev. Lett., 1998, 81, 5924.
16. S. Bucciarelli, N. Mahmodi, L. C. Dujat, M. Jehannin, C. Jud and A. Stradner, J. Phys. Chem. Lett., 2016, 7, 1610.
17. M. Born and J. R. Oppenheimer, Ann. Phys., 1927, 84, 457.
18. T. Udagawa, K. Suzuki and M. Tachikawa, ChemPhysChem, 2015, 16, 3156.
19. M. Tachikawa, K. Mori, H. Nakai and K. Iguchi, Chem. Phys. Lett., 1998, 290, 437.
20. T. Udagawa and M. Tachikawa, J. Chem. Phys., 2006, 125, 244105.
21. T. Udagawa, T. Tsuneda and M. Tachikawa, Phys. Rev. A, 2014, 89, 052519.
22. T. Udadgawa and M. Tachikawa, J. Chem. Phys., 2016, 145, 164310.
23. T. Udagawa, K. Sugiura, K. Suzuki and M. Tachikawa, RSC Adv., 2017, 7, 9328.
24. T. Ishimoto and M. Koyama, Chem. Phys., 2012, 392, 166.
25. D. Sheppard, R. Terrell and G. Henkelman, J. Chem. Phys., 2008, 128, 134106.
26. H. Nakai, Int. J. Quantum Chem., 2002, 86, 411.
27. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901.
28. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978.
29. M. Hashimoto, T. Ishimoto, M. Tachikawa and T. Udagawa, Int. J. Quantum Chem., 2016, 116, 961.
30. T. Ishimoto, M. Tachikawa and U. Nagashima, Int. J. Quantum Chem., 2006, 106, 1465.
31. T. Ishimoto, M. Tachikawa and U. Nagashima, J. Chem. Phys., 2006, 125, 144103.
32. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, N. J. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford CT, 2010.
33. G. Li, L. Zhou, Q. Li, Y. Xie and H. F. Schaefer III, Phys. Chem. Chem. Phys., 2012, 14, 10891.

---

Footnote

† Electronic supplementary information (ESI) available: See supplementary material for some test calculations, the optimized exponent α_opt values obtained in MC_QM calculations, and the optimized geometrical parameters of all stationary point structures obtained by DFT/cc-pVTZ methods. See DOI: 10.1039/c8ra02679c

---