Mechanism of ionic dissociation of HCl in the smallest water clusters

Abstract

The dissociation of strong acids into water is a fundamental process in chemistry and biology. Determining the minimum number of water molecules that can result in an ionic dissociation of hydrochloric acid (HCl → H⁺ + Cl⁻) remains a challenging subject. In this study, the reactions of H₂O with HCl(H₂O)_n−1 (HCl–H₂O cluster), i.e., HCl(H₂O)_n−1 + H₂O (n = 3–7), were investigated by using the direct ab initio molecular dynamics (AIMD) method. Direct AIMD calculations were performed to set the collision energy of H₂O to zero for all trajectories. For n = 3, no reaction occurred. In contrast, HCl dissociated to H⁺ + Cl⁻ at n = 4, forming a contact ion pair (cIP) and solvent-separated ion pair (ssIP) as products. The reactions were expressed as HCl(H₂O)₃ + H₂O → H₃O⁺(H₂O)₂Cl⁻ (ssIP), and HCl(H₂O)₃ + H₂O → H₃O⁺(Cl⁻)(H₂O)₂ (cIP). The ion pair (IP) products were dependent on the collision site of H₂O relative to HCl(H₂O)₃. For n = 5–7, both IPs were formed through the reaction between H₂O and HCl(H₂O)_n−1 (n = 5–7). The reaction between HCl and (H₂O)₄ (HCl + (H₂O)₄ → HCl(H₂O)₄) was non-reactive in IP formation. The reaction mechanism was discussed based on the theoretical results.

---

1. Introduction

Dissociation reactions of chemicals into solvents are a fundamental process in chemistry and biology. Solubility phenomena can be of two types: (1) simple solvation of the solute by the solvent and (2) dissociation of solute molecules into the solvent, followed by solvation of dissociated positive and negative ions. The ionic dissociation of a strong acid, HX, into water (HX → H⁺ + X⁻ (X = F and Cl)) is an example of these phenomena; these ions are solvated by the solvent.

Microhydrated hydrogen chloride, HCl(H₂O)_n, is an experimentally^1–11 and theoretically^12–20 well-studied cluster. This cluster is a model system in which a strong acid dissolves in water. The number of water molecules required for an HCl molecule to dissociate into the ionic state of H⁺ + Cl⁻ is of particular interest.^21–27

Re et al.²⁸ theoretically investigated the dissociation of acids into water clusters. They evaluated the structures of HCl(H₂O)_n (n = 1–5) using the density functional theory method, revealing that the HCl(H₂O)_n clusters have two structural types: non-proton transferred (non-PT) and proton transferred (PT) structures. Scheme 1 displays these structures. In the non-PT structure, HCl remains a neutral molecule (HCl), whereas in the PT structure, HCl dissociates into H⁺ + Cl⁻. Calculations show that the most stable structures for n = 1–3 are non-PT structures. However, at n = 4, the non-PT and PT structures exhibit similar energies, whereas at n = 5, the stability of the energy of the PT structure is greater than that of the non-PT structure. Therefore, the ionic dissociation of HCl occurs at n = 4 to n = 5. The non-PT and PT structures are schematically expressed as HCl(H₂O)_n and H₃O⁺(H₂O)_n−1Cl⁻, respectively.

Scheme 1 Structural types of the non-proton transferred (non-PT) and proton transferred (PT) structures of HCl(H₂O)₃.

Odde et al.²⁹ performed systematic calculations for HX(H₂O)_n clusters (X = F, Cl, Br, and I) up to n = 6 using the B3LYP and second-order-perturbation Møller–Plesset (MP2) levels of theory. The non-PT and PT structures exhibited different vibrational spectra, suggesting that the analysis of the IR spectrum can be used to monitor the ion dissociation reaction, HCl(H₂O)_n → H₃O⁺(H₂O)_n−1Cl⁻.

To experimentally understand the dissociation of HCl into water clusters, Gutberlet et al. observed microhydrated HCl in helium droplets at below 1 K.³⁰ High-resolution mass-selective infrared laser spectroscopy revealed that the aggregation of HCl with water molecules, HCl(H₂O)_n, results in the formation of hydronium ions (H₃O⁺). Computer simulations were conducted to determine the smallest cluster size to be n = 4.^31,32

In 2019, Mani et al. experimentally examined two reaction pathways of HCl⁻ and H₂O reactions using mass-selective infrared spectroscopy.³³ The reaction pathways are expressed as follows:

HCl(H₂O)_m + H₂O → H₃O⁺(H₂O)_mCl⁻, (1)

HCl + (H₂O)_m → HCl(H₂O)_m. (2)

The stepwise addition of H₂O to HCl(H₂O)_m (reaction (1)) led to an ion pair (IP) state (m = 3) composed of H₃O⁺(H₂O)₃Cl⁻, whereas the addition of H₂O clusters to bare HCl did not result in the dissociation reaction (reaction (2)).

The ionic dissociation of the radical cation, HCl⁺, in water clusters has been previously investigated. It involves the ionization states of HCl(H₂O)_n. Using static ab initio calculations, Calatayud et al. investigated the ionization of HCl and HF molecules trapped in an ice lattice in which one water molecule was replaced by an HCl or HF molecule.³⁴ The potential energy curve plotted as a function of the H–Cl distance of HCl⁺ showed that HCl⁺ dissociated to H⁺ + Cl without an activation barrier in the ice lattice.³⁵ The dynamics calculations suggested that very fast proton transfer occurs in the HCl⁺ system.

In the present study, the reaction of HCl(H₂O)_n−1 with H₂O (n = 3–7) was investigated using the direct ab initio molecular dynamics (AIMD) method to elucidate the reaction mechanism of the dissociation dynamics of HCl in water clusters.^36,37 In particular, the minimum number of water molecules required to dissociate HCl was determined.

In previous studies, the reaction mechanisms between NO⁺ and water clusters in the atmosphere were investigated. The bimolecular reaction, NO⁺ + (H₂O)_n (n = 1–7), was examined by using the direct AIMD method.³⁸ The formation mechanism of nitrous acid in the atmosphere was proposed. In addition, the direct AIMD method was applied to bimolecular and intracluster S_N2 reactions.³⁹ The effects of water clusters on the S_N2 reaction mechanism were discussed.

2. Computational methods

2.1. Ab initio calculations

The geometries of HCl(H₂O)_n (n = 2–6) were optimized using the MP2 method with a 6-311++G(d,p) basis set. The CAM-B3LYP functional was also used for geometry optimization.⁴⁰ The coupled-cluster single, double, and perturbative triple excitation (CCSD(T)) method⁴¹ was used for single-point energy calculations of the MP2 geometries. Atomic and molecular charges were determined by applying natural population analysis. All static ab initio calculations were performed using the standard Gaussian 09 program package.⁴²

2.2. Direct AIMD calculations

For the direct AIMD calculations, the H₂O and HCl–water clusters, HCl(H₂O)_n, were separately optimized at the CAM-B3LYP/6-311G(d,p) level. Subsequently, H₂O was positioned around HCl(H₂O)_n at a distance of 4.0–6.0 Å; the trajectory began from this position. The collision energy was zero (E_coll = 0) at time zero. The excess energy, momentum vector, and rotational temperature of the reaction system were assumed to be zero (at 0 fs). The equations of motion for N atoms in the reaction system are given bywhere j = 1–3 N, H is the classical Hamiltonian, Q_j is the Cartesian coordinate of the j-th mode, and P_j is the conjugated momentum; these equations were numerically solved (NVE ensemble). The velocity Verlet algorithm, with a time step of 0.10–0.25 fs, was used to solve the equations of motion for the system. The maximum simulation time was 2.0 ps. The total energy drift in all trajectory calculations was <0.01 kcal mol⁻¹. The numbers of trajectory runs were 160 (n = 4), 50 (n = 5), 20 (n = 6), and 10 (n = 7), respectively. The effects of the zero-point energy (ZPE) were included by employing the classical vibrational sampling method (microcanonical ensemble).^43–46 The effects of the functional on the reaction mechanism were investigated using the wB97XD functional and compared with those of the CAM-B3LYP functional. Direct AIMD calculations were performed using custom-made AIMD codes combined with Gaussian09 (our own AIMD codes written by C language).^36,37

It should be noted that the present reaction system has no activation barrier. Therefore, the tunneling effects may be secondary matter.

3. Results

3.1. Structure of HCl(H₂O)_n (n = 2–4)

The optimized structures of HCl(H₂O)_n (n = 3–6) obtained at the MP2/6-311++G(d,p) level of theory are illustrated in Fig. 1; additionally, the optimized geometrical parameters obtained at several levels of theory are provided in Table 1 and Table S1 (ESI†). For n = 3, the intermolecular distance between HCl and W1 was r(Cl–O1) = 2.988 Å for the Cl–O distance. The distance between HCl and W3 was r(Cl–O3) = 3.308 Å, whereas the water–water distances were r(O1–O2) = 2.727 Å and r(2–O3) = 2.770 Å. The distances determined by using the CAM-B3LYP functional were r(Cl–O1) = 2.915 Å and r(Cl–O3) = 3.258 Å, which corroborated those of the MP2 calculations. Similar structures were obtained for the larger systems (n = 4–6).

Fig. 1 Optimized structures of microhydrated hydrogen chloride, HCl(H₂O)_n (n = 3–6), evaluated at the MP2/6-311++G(d,p) level. The geometrical parameters calculated at all levels of theory are listed in Table S1 (ESI†).

Table 1 Optimized geometrical parameters of HCl(H₂O)_n (n = 3), calculated at the MP2/6-311++G(d,p), CAM-B3LYP/6-311++G(d,p), CAM-B3LYP/6-311G(d,p), and CAM-B3LYP/6-31G(d) levels of theory

n = 3cyc
	MP2/6-311++G(d,p)	CAM-B3LYP/6-311++G(d,p)	CAM-B3LYP/6-311G(d,p)
r 1[Cl-O(W1)]	2.988	2.915	2.808
r 2[O(W1)-O(W2)]	2.727	2.671	2.603
r 3[O(W2)-O(W3)]	2.77	2.721	2.681
r 4[O(W3)-Cl]	3.308	3.258	3.211

---

In n = 5 and 6, HCl was located on the surface of water clusters (rather than inside). This is due to the fact that the interaction energy between H₂O–H₂O in the water cluster is slightly larger than that of HCl–H₂O.

The binding energies of HCl with the (H₂O)_n cluster were calculated at the MP2/6-311++G(d,p) level. For n = 4, 5, and 6, the binding energies corrected by the basis set superposition errors (BSSEs) were calculated as 0.3, 4.1, and 4.2 kcal mol⁻¹, respectively, indicating that the binding energies of HCl to (H₂O)_n are slightly smaller than those of H₂O–H₂O (4.6 kcal mol⁻¹).

3.2. Reaction dynamics of H₂O + HCl(H₂O)_n−1 (n = 3 + 1): solvent-separated ion pair (ssIP) channel

Using these optimized geometry of n = 3, direct AIMD calculations were performed for H₂O + HCl(H₂O)₃ collision reactions. Fig. 2 shows the collision regions (I–IV) and types of H₂O around the HCl–H₂O cluster. The collision energy of H₂O(W0) was set to zero (E_coll = 0) at time zero for all trajectories. The initial position of the colliding H₂O was distributed in 4–6 Å from the nearest molecule of the cluster, and the height of H₂O from the molecular place was 0–6 Å. The colliding H₂O was started from mostly out-of-plane. Four collision regions were examined by performing direct AIMD calculations. In region I, H₂O collided directly with HCl without collision energy (zero collision energy). In regions II, III, and IV, H₂O(W0) attacked W1, W2, and W3, respectively. Two types of collisions—the OH and oxygen collisions—were mainly examined in the calculations In the case of the OH collision, the O–H group of H₂O(W0) was almost oriented toward the HCl–H₂O cluster at time zero. In contrast, the oxygen atom of H₂O(W0) was almost oriented toward the cluster at time zero during the oxygen collision. A total of 40 trajectories were run from each region (regions I – IV); the total number of trajectories was 160 in the reaction n = 3 + 1 system. The initial orientation of the colliding water molecules was randomly generated around the OH- and oxygen-oriented forms.

Fig. 2 Collision regions of H₂O around HCl(H₂O)₃ and collision types of H₂O with HCl(H₂O)₃ at initial separation (time zero). Collision regions are classified into four regions in the present calculations. Region I: H₂O collides directly with HCl; regions II, III, and IV: H₂O collides with W1, W2, and W3, respectively. The initial orientation of H₂O relative to HCl(H₂O)₃ in time zero (OH) and oxygen orientation are given.

Snapshots of the reactive sample trajectory for HCl(H₂O)₃ + H₂O (n = 3 + 1) to form ssIPs are displayed in Fig. 3. Notation (n = 3 + 1) denotes that one H₂O molecule collides with an HCl(H₂O)_n cluster composed of three water molecules (n = 3). Direct AIMD calculations were performed at the CAM-B3LYP/6-311G(d,p) level of theory, and the optimized structures of HCl(H₂O)₃ and H₂O were used as the geometries at time zero. In this trajectory, the H₂O(W0) was located in region III and the distance between H₂O(W0) and W2 was R = 4.112 Å at time zero (0 fs). The collision energy of W0 was set to zero (0 kcal mol⁻¹) at time zero (0 fs).

Fig. 3 Sample trajectory for the solvent-separated ion pair (ssIP) channel: (A) snapshots of the HCl(H₂O)_n + H₂O (n = 3 + 1) reaction system given as a function of time and (B) time evolutions of potential energy. H₂O(W0) approaches HCl from region III.

W0 gradually approached W2 after the reaction began; the distance between W0 and W2 was R = 2.380 Å at 260 fs. The distances of r₁ and r₂ were 1.465 and 1.381 Å at 260 fs, respectively. The bond length of H–Cl was r₁ = 1.381 Å, which was marginally longer than the H–Cl equilibrium distance (1.283 Å). The collision of W0 with W2 occurred at 292 fs with R = 1.674 Å. A Zundel-type complex was formed between W1 and W2. Simultaneously, the proton of HCl was transferred to W1. The position of the proton was r₁ = 1.491 Å and r₂ = 1.289 Å, indicating that the proton was closer to W1 than it was to Cl⁻. The proton of W1 approached that of W2. W1 relied on protons from HCl to W2. At 355 fs, the proton was fully transferred to W2, which produced H₃O⁺(W2). Eigen-type complex formed between W1 and W2. HCl was ionized to the Cl⁻ ion. Concurrently, the W1 proton was transferred to W2; thus, ssIP was formed at 708 fs. After H₂O(W0) transferred (708–860 fs), ssIP was solvated by W0, and stable ssIP was formed as the product at 860 fs.

The potential energy of the reaction system is plotted as a function of time in Fig. 3B. The initial zero energy corresponded to the total energy of the HCl(H₂O)₃ + H₂O reaction system. After the reaction began, W0 gradually approached W2 and the energy decreased from 0 to −5.0 kcal mol⁻¹ (292 fs), which corresponded to the solvation energy of HCl(H₂O)₃ by H₂O. An incomplete structure of ssIP was formed at 355 fs. After the transfer of W0 to ssIP, a fully stable ssIP was formed at 860 fs owing to the complete solvation of H₂O to ssIP.

3.3. Reaction dynamics of H₂O + HCl(H₂O)_n−1 (n = 3 + 1): contact ion pair (cIP) channel

Snapshots of the sample trajectory leading to the cIP channel in HCl(H₂O)₃ + H₂O (n = 3 + 1) are presented in Fig. 4. The H₂O(W0) was located in region I, and the distance of H₂O(W0) from W1 was R = 4.200 Å at time zero (0 fs). The collision energy, W0, was set to zero (0 kcal mol⁻¹) at time zero (0 fs). This reaction occurred in region I.

Fig. 4 Sample trajectory for the contact ion pair (cIP) channel: (A) snapshots of the HCl(H₂O)_n + H₂O (n = 3 + 1) reaction system given as a function of time and (B) time evolutions of potential energy. H₂O(W0) approaches HCl from region I.

After the reaction started, W0 gradually approached HCl; the distance between W0 and W1 was R = 2.368 Å at 302 fs. The distances of r₁ and r₂ were 1.670 and 1.110 Å at 302 fs, respectively. Zundel complex was formed. The proton of HCl was fully transferred to W1, and W1 was converted to H₃O⁺(W1). Concurrently, the collision between W0 and HCl occurred, resulting in the formation of cIP as Cl⁻–H₃O⁺(W1). At 860 fs, the stability of cIP increased owing to structural deformation. W0 moved from Cl⁻ to H₃O⁺, which increased the stability of the entire structure. Stable cIP and Eigen complex (W1) were formed at 1260 fs.

The potential energy of the reaction system is plotted as a function of time in Fig. 4B. The initial zero energy corresponded to the total energy of the HCl(H₂O)₃ + H₂O reaction system. After the commencement of the reaction, W0 gradually approached HCl and the energy decreased from 0 to −5.0 kcal mol⁻¹ (302 fs), which is the solvation energy (hydrogen bonding energy) of HCl(H₂O)₃ by H₂O. The formation of cIP at 750–860 fs was temporary. After the transfer of W0 to H₃O⁺(W1), a more stable cIP was formed at 1260 fs because of its complete solvation.

3.4. Summary of reaction dynamics of H₂O + HCl(H₂O)_n−1 (n = 3 + 1)

Similar direct AIMD calculations were performed for regions I–IV. A total of 160 trajectories were run; the main products of these regions are summarized in Table 2. The main product in regions I and II was cIP, whereas that of region III was ssIP. In region IV, almost all trajectories were non-reactive. The product states were strongly dependent on the collision sites of H₂O around HCl(H₂O)₃. The reason for this is discussed in the following section.

Table 2 Summary of the main product channels in the bi-molecular reaction, HCl(H₂O)_n + H₂O (n = 3 + 1). CAM-B3LYP and wB87XD indicate that direct AIMD calculations were performed at the CAM-B3LYP/6-311G(d,p) and wB87XD/6-311G(d,p) levels of theory, respectively

Region	CAM-B3LYP	wB87XD
I	cIP	cIP
II	cIP	cIP
III	ssIP	ssIP
IV	Non-reactive	Non-reactive, cIP

---

The sample trajectories leading to ssIP (reaction in region III) are depicted in Fig. S1–S3 (ESI†). Similar features were obtained for all trajectories. The trajectories in region IV were non-reactive (a sample trajectory is shown in Fig. S4, ESI†); the trajectories in region II mainly led to cIP (a sample trajectory is shown in Fig. S5, ESI†).

Two types of channels were obtained as nonreactive trajectories: (1) HCl molecules dissociated to the ions (H⁺ + Cl⁻) and then moved back to the HCl molecule, and (2) HCl molecules never dissociated and remained as HCl.

3.5. Reaction dynamics of H₂O + HCl(H₂O)_n−1 (n = 4 + 1)

The reaction dynamics for ssIP formation from the (4 + 1) reaction are shown in Fig. 5. Fig. 5A and B show snapshots and the potential energy of the system, respectively. The ZPE level was the total energy of (H₂O + HCl(H₂O)₄). The distance between water molecules W0 and W2 was R = 4.794 Å. After the reaction began, W0 gradually approached W2; these molecules collided at 290 fs. The energy was reduced to −6.0 kcal mol⁻¹. At 515 fs, a book-form structure composed of HCl(H₂O)₄ was formed, in which all the molecules were bound by to each other through hydrogen bonds. The energy decreased to −13.0 kcal mol⁻¹. Subsequently, the proton of HCl was transferred to W1, and the H₃O⁺(Cl⁻) complex was tentatively formed at 760 fs (the energy was −12.0 kcal mol⁻¹). The proton of H₃O⁺ was further transferred to the next water molecule (W2), and W2 was converted to H₃O⁺(W2) at 927 fs. Simultaneously, the Eigen complex was formed. Concurrently, ssIP was formed in H₃O⁺(W2)–Cl⁻ (the energy was −13.0 kcal mol⁻¹). The stability of ssIP increased owing to structural deformation at 1042 fs; the energy was stabilized to −19.5 kcal mol⁻¹. The reaction and ssIP formation were completed at 1042 fs.

Fig. 5 Sample trajectory for the ssIP channel: (A) snapshots of the HCl(H₂O)_n + H₂O (n = 4 + 1) reaction system given as a function of time and (B) time evolutions of potential energy.

Similar direct AIMD calculations for H₂O + HCl(H₂O)_n−1 (n = 4 + 1) were performed for 50 trajectories. The products were ssIP and cIP in the (n = 4 + 1) dynamics calculations.

Similar direct AIMD calculations for H₂O + HCl(H₂O)_n−1 (n = 4 + 1), (n = 6 + 1), and (n = 6 + 1) were performed for 50 trajectories. The products were ssIP and cIP in the (n = 4 + 1) dynamic calculations. Some sample trajectories are provided in the ESI.†

Similar calculations were performed for a larger system (H₂O + HCl(H₂O)_n−1 (n = 5 + 1)); both ssIP and cIP were produced. The sample trajectories are shown in Fig. S6 (ESI†).

3.6. Structures of IPs

Direct AIMD calculations showed that two types of IP were formed by following the reaction: HCl(H₂O)_n−1 + H₂O (n = 4–6) → ssIP or cIP. This section describes the accurate structures of the IPs. The optimized structures of the IPs for n = 4 and 5 evaluated at the MP2/6-311++G(d,p) level are shown in Fig. 6. IP structures existed for n = 4 and 5: cIP and ssIP. In cIP, Cl⁻ and H₃O⁺ ions were in direct contact each other; the proton of H₃O⁺ interacts directly with the Cl⁻ ions. The distance between Cl⁻ and the oxygen atom of H₃O⁺ was 2.773 Å in cIP (n = 4). In contrast, solvent molecules separated Cl⁻ and H₃O⁺ in ssIP. The distances between Cl⁻ and H₃O⁺ were 3.497 Å (n = 4) and 3.675 Å (n = 5) for Cl–O distances. The dipole moment of H₃O⁺ oriented toward the Cl⁻ ion. The CAM-B3LYP/6-311++G(d,p) level exhibited similar IP structures (Table S2, ESI†).

Fig. 6 Optimized structures of the ion pair (IP) states of H₃O⁺(H₂O)_nCl⁻ (n = 4 and 5) calculated at the MP2/6-311++G(d,p) level. The geometrical parameters calculated at all levels of theory are presented in Table S2 (ESI†).

The relative energies of cIP and ssIP are listed in Table 3. For n = 4, the energetic stability of ssIP was marginally higher than that of cIP. In contrast, cIP was the more stable form at n = 5. This relationship was not affected by higher electron correlation, coupled-cluster single and double (CCSD), and CCSD(T). This occurred because of the shorter interatomic distance between Cl⁻ and H₃O⁺ in ssIP (n = 5). However, the energy difference between ssIP and cIP was significantly small. Hence, both IPs were formed after the reaction of H₂O with the HCl(H₂O)_n clusters.

Table 3 Relative energies between cIP and ssIP (in kcal mol⁻¹)

	cIP	ssIP
MP2/6-311++G(d,p)//MP2/6-311++G(d,p)
n =4	0.0	−0.62
n =5	0.0	0.82
CCSD/6-311++G(d,p)//MP2/6-311++G(d,p)
n =4	0.0	−0.11
n =5	0.0	1.16
CCSD(T)/6-311++G(d,p)//MP2/6-311++G(d,p)
n =4	0.0	−0.41
n =5	0.0	1.01

---

3.7. Effects of functional and of ZPE on the reaction mechanism

Above dynamics calculations were performed using the CAM-B3LYP functional. In this section, we will consider the effect of the functionals used on the reaction mechanism. To assess the effect of the functional on the reaction mechanism, direct AIMD calculations were performed using the wB97XD functional. A total of 60 trajectories were run. The results of the sample trajectories are shown in Fig. S7 (ESI†), and the main products are summarized in Table 2. The wB97XD functional yielded results similar to those obtained by using the CAM-B3LYP functional. Therefore, there is no effect of the functionals used on the reaction mechanism.

To elucidate the effects of ZPE on the reaction mechanism, direct AIMD calculations including ZPE were performed for the n = 3 + 1 reaction system. All calculations were performed at the CAM-B3LYP/6-31G(d) level of theory. First, the geometries of HCl(H₂O)_n−1 and H₂O were separately optimized. Subsequently, vibrational energies were applied to the reaction system. A microcanonical ensemble of the classical harmonic vibrations sampling method was applied to the reaction system. All atoms in the initial state of the reaction system had momenta that corresponded to their vibrational energies. The trajectories of the samples are shown in Fig. S8 (ESI†). At 503 fs, cIP was formed after proton transfer from HCl to W1. Thereafter, further proton transfer from H₃O⁺(W1) to W2 occurred, resulting in the formation of ssIP at 541 fs. Furthermore, proton transfer from H₃O⁺(W3) to Cl⁻ resulted in the formation of HCl, indicating that the proton was returned to the hydrogen bond network. The trajectory moved back and forth between cIP, ssIP, and HCl. This occurred because the ZPE of the reaction system was higher than the reaction energy (5 kcal mol⁻¹). These results indicated that ZPE enhanced reaction efficiency but caused the instability of the IPs.

4. Discussion

4.1. Roles of water molecules in dissociation

Fig. 7 shows the reaction model for HCl dissociation in a water cluster. The roles of each water molecule are discussed in this section. First, the collision in region III was considered. After the collision of water molecules (W0) with W2 in HCl(H₂O)₃, W0 provided energy to HCl(H₂O)₃. This energy corresponded to the hydrogen-bonding energy of H₂O to HCl(H₂O)₃ (∼5 kcal mol⁻¹). Furthermore, W0 increased the proton affinity of W2 via hydrogen bonding. HCl was dissociated to H⁺ + Cl⁻ after acceptance of the collision energy from W0. W1 received a proton from HCl (first proton transfer, PT1) and then relayed the proton to W2 (second PT: PT2); PT1 and PT2 occurred simultaneously. The proton affinity of W2 increased owing to the hydrogen bond with W0; additionally, W2 was converted to H₃O⁺ after the proton relay from HCl via W1. W3 stabilized the Cl⁻ ion caused by the dissociation of HCl. These results reveal that the four water molecules are effective in dissociating HCl from the water clusters and are required to dissociate HCl.

Fig. 7 Reaction model and roles of H₂O in the formation of IPs from the reaction for HCl(H₂O)₃ + H₂O. PT denotes the proton transfer process.

The next case involved the reaction in region II. When W0 collided with and attached to W1 (region II), the proton affinity of W1 increased because of the formation of a new hydrogen bond from W0. After the collision, the proton of HCl was transferred to W1 (PT1). The proton remained in W1, and PT2 did not occur because the proton affinity of W1 was higher than that of W2. Hence, cIP was formed between Cl⁻ and H₃O⁺(W1). For the reactions in region I, namely the collision of W0 with HCl (region I), the proton of W0 stabilized Cl⁻, resulting in the occurrence of only PT1. Therefore, cIP was formed by the collision of W0 with HCl.

4.2. Remarks

Several approximations have been introduced to calculate the reaction dynamics. First, the medium surrounding the reaction system (rare gas) was not considered. In the experiments, the reaction system was surrounded by a rare gas (helium). It is generally known that in the reaction of molecules in a helium droplet, the heat of reaction of product (PD) is transferred to the translational kinetic energy of the helium. Therefore, PD can survive without decomposition. Since PDs in this study are ion pairs, the efficiency of survival of PDs is expected to be more efficient, if rare gas is included in direct AIMD calculations. We plan to perform these calculations in the near future. Ab initio MD calculations using quantum mechanical (QM), molecular mechanics (MM), and ONIOM (our own n-layered integrated molecular orbital and molecular mechanics) methods^47–52 are effective for calculating the medium effects.

Next, the number of trajectories was limited in the present calculation, because the direct AIMD method requires a long CPU time. Therefore, we discuss only the qualitative nature of this study. A discussion of the branching ratios of ssIP and cIP was not included. By performing further calculations, detailed ratios and mechanisms can be discussed.

Finally, for larger clusters, further stabilized structures may exist. To obtain more detailed and stable structures, it is necessary to use a genetic algorithm^53–57 to produce all the stable structures. It should be noted that the reaction model was constructed based on the above assumptions. Despite the several approximations introduced herein, the results enable us to obtain valuable information on the dissociation mechanism of HCl in smallest water cluster.

5. Conclusion

Dissociation of strong acids in water is a basic process in chemistry. Therefore, determining the minimum number of water molecules required for the ionic dissociation of hydrochloric acid (HCl → H⁺ + Cl⁻) is a challenging and fundamental subject. In this study, the direct AIMD method was applied to a bimolecular reaction (HCl(H₂O)_n−1 + H₂O) to solve this problem. The calculations indicated that HCl can dissolve into four water molecules when H₂O collides with the HCl(H₂O)₃ cluster. The reaction is expressed as follows:

HCl(H₂O)₃ + H₂O(W0) → H₃O⁺(H₂O)₃(Cl⁻) (ssIP channel)

→ H₃O⁺(Cl⁻)(H₂O)₃ (cIP channel),

where H₂O(W0) represents an additional water molecule. These reactions proceeded easily without the collision energy of H₂O(W0). Both ssIP and cIP reaction channels were determined to be reactive collisions. These channels were strongly dependent on the collision region of H₂O(W0) relative to HCl(H₂O)₃. If H₂O collided with HCl and H₂O(W1), cIP was mainly formed. In contrast, when H₂O collided with H₂O(W2) at the opposite position of HCl, ssIP was formed. Thus, the reaction sites (regions I–IV) determined the product states.

The calculations also showed that the reaction of HCl with water clusters results in only cluster formation, which is expressed as

HCl + (H₂O)_n → HCl(H₂O)_n (n = 4 and 5).

HCl was not dissociated. This occurred because the water molecules in the reaction system did not function efficiently. In contrast, each water molecule in the reaction system of HCl(H₂O)₃ + H₂O functioned efficiently during the dissociation of HCl. Therefore, HCl can be dissociated via the following reaction: HCl(H₂O)₃ + H₂O(W0). These results corroborate previous experimental and theoretical calculations.^30–33

Conflicts of interest

The author declares no competing interests.

Acknowledgements

The author acknowledges the partial support from JSPS KAKENHI (Grant Numbers 21K04973 and 21H05415). This paper was written during my hospitalization at the Matsuda Orthopedics Memorial Hospital in Sapporo. I would like to thank the doctors and medical staff.

References

1. A. K. Samanta, Y. Wang, J. S. Mancini, J. M. Bowman and H. Reisler, Chem. Rev., 2016, 116, 4913–4936.
2. B. E. Rocher-Casterline, A. K. Mollner, L. C. Ch’ng and H. Reisler, J. Phys. Chem. A, 2011, 115, 6903–6909.
3. A. K. Samanta, G. Czako-, Y. Wang, J. S. Mancini, J. M. Bowman and H. Reisler, Acc. Chem. Res., 2014, 47, 2700–2709.
4. J. S. Mancini and J. M. Bowman, Phys. Chem. Chem. Phys., 2015, 17, 6222–6226.
5. M. Farnik, M. Weimann and M. A. Suhm, J. Chem. Phys., 2003, 118, 10120–10136.
6. D. Skvortsov, S. J. Lee, M. Y. Choi and A. F. Vilesov, J. Phys. Chem. A, 2009, 113, 7360–7365.
7. M. Letzner, S. Gruen, D. Habig, K. Hanke, T. Endres, P. Nieto, G. Schwaab, L. Walewski, M. Wollenhaupt and H. Forbert, et al. , J. Chem. Phys., 2013, 139, 154304.
8. D. Skvortsov, S. J. Lee, M. Y. Choi and A. F. Vilesov, J. Phys. Chem. A, 2009, 113, 7360–7365.
9. M. Ortlieb, O. Birer, M. Letzner, G. W. Schwaab and M. Havenith, J. Phys. Chem. A, 2007, 111, 12192–12199.
10. A. J. Huneycutt, R. J. Stickland, F. Hellberg and R. J. Saykally, J. Chem. Phys., 2003, 118, 1221–1229.
11. J. Zischang, D. Skyortsov, M. Y. Choi, R. A. Mata, M. A. Suhm and A. F. Vilesov, J. Phys. Chem. A, 2015, 119, 2636–2643.
12. K. Zuraski, D. Kwasniewski, A. K. Samanta and H. Reisler, J. Phys. Chem. Lett., 2016, 7, 4243–4247.
13. W. Lin and F. Paesani, J. Phys. Chem. A, 2013, 117, 7131–7141.
14. R. L. DeKock, B. M. Brandsen and J. R. Strikwerda, Theor. Chem. Acc., 2011, 130, 871–881.
15. F. M. Balci and N. Uras-Aytemiz, J. Phys. Chem. A, 2011, 115, 5943–5954.
16. T. Takayanagi, K. Takahashi, A. Kakizaki, M. Shiga and M. Tachikawa, Chem. Phys., 2009, 358, 196–202.
17. U. F. T. Ndongmouo, M.-S. Lee, R. Rousseau, F. Baletto and S. Scandolo, J. Phys. Chem. A, 2007, 111(49), 12810–12815.
18. M. E. Alikhani and B. Silvi, Phys. Chem. Chem. Phys., 2003, 5, 2494–2498.
19. M. Masia, H. Forbert and D. Marx, J. Phys. Chem. A, 2007, 111, 12181–12191.
20. M. E. Alikhani and B. Silvi, Phys. Chem. Chem. Phys., 2003, 5, 2494–2498.
21. O. I. Arillo-Flores, M. F. Ruiz-Lopez and M. I. Bernal-Uruchurtu, Theor. Chem. Acc., 2007, 118, 425–435.
22. V. Buch, A. Dubrovskiy, F. Mohamed, M. Parrinello, J. Sadlej, A. Hammerich, A. D. Hammerich and J. P. Devlin, J. Phys. Chem. A, 2008, 112, 2144–2161.
23. M. Masia, H. Forbert and D. Marx, J. Phys. Chem. A, 2007, 111, 12181–12191.
24. S. Sugawara, T. Yoshikawa, T. Takayanagi and M. Tachikawa, Chem. Phys. Lett., 2011, 501, 238–244.
25. H. Forbert, M. Masia, A. Kaczmarek-Kedziera, N. N. Nair and D. Marx, J. Am. Chem. Soc., 2011, 133, 4062–4072.
26. J. S. Mancini and J. M. Bowman, Phys. Chem. Chem. Phys., 2015, 17, 6222–6226.
27. S. Sugawara, T. Yoshikawa, T. Takayanagi and M. Tachikawa, Chem. Phys. Lett., 2011, 501, 238–244.
28. S. Re, Y. Osamura, Y. Suzuki and H. F. Schaefer III, J. Chem. Phys., 1998, 109, 973–977.
29. S. Odde, B. J. Mhin, S. Lee, H. M. Lee and K. S. Kim, J. Chem. Phys., 2004, 120, 9524–9535.
30. A. Gutberlet, G. Schwaab, O. Birer, M. Masia, A. Kaczmarek, H. Forbert, M. Havenith and D. Marx, Science, 2009, 324, 1545–1548.
31. M. Letzner, S. Gruen, D. Habig, K. Hanke and T. Endres, et al. , J. Chem. Phys., 2013, 139, 154304.
32. H. Forbert, M. Masia, K. Kaczmarek-Kedziera, N. N. Nair and D. Marx, J. Am. Chem. Soc., 2011, 133, 4062–4072.
33. D. Mani, R. P. de Tudela, R. Schwan, N. Pal, S. Körning and H. Forbert, et al. , Sci. Adv., 2019, 5, 8179.
34. M. Calatayud, D. Courmier and C. Minot, Chem. Phys. Lett., 2003, 369, 287–292.
35. A. L. Sobolewski and W. Domcke, J. Phys. Chem., 2003, 107, 1557–1562.
36. H. Tachikawa, ACS Omega, 2023, 8, 10600–10606.
37. H. Tachikawa, Carbon. Res., 2022, 8, 36–57.
38. H. Tachikawa, J. Phys. Chem. A, 2022, 126, 119–124.
39. H. Tachikawa, Phys. Chem. Chem. Phys., 2022, 24, 3941–3950.
40. T. Yanai, T. D. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
41. M. Saitow, U. Becker, C. Riplinger, E. F. Valeev and F. Neese, J. Chem. Phys., 2017, 146, 164105.
42. 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, et al., Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2013.
43. H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. Daniels, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 2001, 114, 9758–9763.
44. S. S. Iyengar, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 2001, 115, 10291–10302.
45. H. B. Schlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 2002, 117, 8694–8704.
46. H. Tachikawa, J. Comput. Chem., 2017, 38, 1503–1508.
47. M. A. Kochman, A. Bilbc and C. A. Morrison, Phys. Chem. Chem. Phys., 2013, 15, 10803–10816.
48. Q. H. Meng, Y. C. Chi, L. D. Zhang, P. Zhang and L. S. Sheng, Phys. Chem. Chem. Phys., 2019, 21, 5232–5242.
49. R. Dushanan, S. Weerasinghe, D. P. Dissanayake and R. Senthilinithy, Mol. Sim., 2022, 48, 1464–1475.
50. H. Tachikawa, Surf. Sci., 2016, 647, 1–7.
51. H. Tachikawa, J. Phys. Chem. A, 2020, 124, 3048–3054.
52. H. Tachikawa, J. Phys. Chem. A, 2016, 120, 7301–7310.
53. D. M. Deaven and K. M. Ho, Phys. Rev. Lett., 1995, 75, 288–291.
54. S. K. Gregurick, M. H. Alexander and B. Hartke, J. Chem. Phys., 1996, 104, 2684–2691.
55. C. Roberts, R. L. Johnston and N. T. Wilson, Theor. Chem. Acc., 2000, 104, 123–130.
56. J. C. Marques and F. B. Pereira, Chem. Phys. Lett., 2010, 485, 211–216.
57. J. L. Llanio-Trujillo, J. M. C. Marques and F. B. Pereira, J. Phys. Chem. A, 2011, 115, 2130–2138.

---

Footnote

† Electronic supplementary information (ESI) available: ESI is available on the RCS website: the optimized geometrical parameters of HCl(H₂O)_n (n = 3–5) and IPs, sample trajectories for ssIP and cIP, effect of the functional on the reaction mechanism (wB97XD/6-311G(d,p) level), and the results of direct AIMD calculations using the wB97XD functional are provided. See DOI: https://doi.org/10.1039/d3cp05715a

---