Microscopic origin of quantum plasticity in small H₃⁺(H₂)_n (n = 1–3) clusters revealed by path integral molecular dynamics simulations

Abstract

We investigate the microscopic origin of “quantum plasticity,” i.e., flexible molecular geometries induced by nuclear quantum effects, in small H₃⁺(H₂)_n (n = 1–3) clusters using first-principles path integral molecular dynamics simulations. Nuclear quantum effects soften the rotational barriers of the H₂ ligands, stabilizing structurally flexible configurations absent in classical descriptions.

---

Neutral clusters such as (H₂)_n, (D₂)_n, and ⁴He_n have long served as prototypical systems for investigating quantum effects at low temperatures and have therefore been extensively studied both theoretically and experimentally.^1–5 In particular, small (p-H₂)_n clusters (n ≤ 20) exhibit pronounced quantum behavior and are known to retain liquid-like properties even near the ground state. In contrast, the corresponding (o-D₂)_n clusters tend to be more rigid and freeze into solid-like configurations.^2,3 For p-H₂ clusters, superfluid behavior analogous to that observed in ⁴He_n clusters has been demonstrated, with nuclear quantum effects (NQEs) playing an essential role in the emergence of this behavior.^6–12 In addition, it has been reported that H₂ molecules and helium atoms confined within carbon nanotubes exhibit striking quantum delocalization.^13,14 Furthermore, in metal-doped quantum clusters, the position of the metal atom—whether it is trapped within the cluster or resides on its surface—strongly depends on the nature of the metal.^15,16

Beyond neutral systems, cationic clusters have also attracted considerable attention over decades of experimental investigations.^17–22 Among them, in H₃⁺(H₂)_n clusters, the central H₃⁺ core is sequentially solvated by H₂ molecules, leading to solvation shells that close at the magic numbers, such as n = 3, 6, and 8. Moreover, infrared predissociation spectroscopy has provided compelling evidence that H₃⁺(H₂)₃ adopts a structure in which three H₂ ligands bind to the vertices of the triangular H₃⁺ core.^21,22 Numerous experimental and theoretical studies on these H_n⁺ clusters have been reported to date.^23–38 Several studies have pointed out that, since the proton is the lightest nucleus, NQEs become significant for these H_n⁺ clusters, affecting their structures and binding properties.^{31,32,35–37} Therefore, incorporating NQEs is indispensable for conducting precise theoretical studies on the H_n⁺ clusters.

There have been several theoretical studies on H_n⁺ clusters. Among them, Štich et al. reported path integral molecular dynamics (PIMD) studies,^31,32 in which NQEs were represented by expanding each nucleus into a necklace-like series of beads, for H_n⁺ clusters at the generalized gradient approximation (GGA) level of density functional theory (DFT). They introduced the interesting concept of “quantum plasticity”, referring to flexible changes in molecular geometries arising from the NQEs on the ligand H₂ molecules, which drastically alter the dynamics of the H_n⁺ clusters by enabling facile rotations of the ligand H₂ molecules. This phenomenon is not limited to the H_n⁺ clusters but can be a key factor in various hydrogen-bonded or proton-transfer systems where NQEs play a dominant role, as it represents a fundamental aspect of how nuclear quantum fluctuations can reshape potential energy surfaces and molecular dynamics.

The spatial distributions and orientational configurations of the surrounding H₂ ligands are the main focus of this study. Although Štich et al. introduced the concept of quantum plasticity, referring to the NQEs that facilitate ligand rotation, its origin has not yet been analyzed in detail. In this study, we aim to reveal the mechanism behind quantum plasticity through highly accurate PIMD simulations.

However, PIMD simulations require a large number of force evaluations, which makes highly accurate ab initio treatments impractical for achieving sufficient statistical sampling. Therefore, it is important to employ an electronic-structure method that can estimate atomic forces with low computational cost while retaining high accuracy. In this study, the B3LYP/6-31G method with Grimme's D3 dispersion correction (GD3) was used for PIMD simulations because this method can reproduce not only the dissociation energies for H₃⁺(H₂)_n → H₃⁺(H₂)_n−1 + H₂ (n = 1–3) calculated using the highly accurate CCSD(T)/CBS method but also the optimized geometrical parameters obtained by CCSD(T)/aug-cc-pVQZ calculations. The CCSD(T)/CBS energies for reference were obtained by extrapolating CCSD(T) energies calculated with the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets according to the following expression:³⁴

E(x) = E^CBS + Be^(−(x−1)) + Ce(^−(x−1)2), (1)

where x is the cardinal number of the basis set (x = 2, 3, and 4 for D, T, and Q, respectively). All PIMD simulations were performed using a Nosé–Hoover chain at T = 50 K, and each nucleus was represented by P = 96 beads. In addition to the PIMD simulations, classical molecular dynamics (CLMD) simulations with P = 1 (i.e., without NQEs) were also carried out for comparison. The time step was set to Δt = 0.1 fs. For each cluster, PIMD (CLMD) simulations were performed for 40 000 (450 000) steps after an equilibration period of 10 000 (50 000) steps. For the H₅⁺ cluster, an additional 50 000 steps of PIMD were carried out, resulting in a total of 90 000 steps. The PIMD simulations were conducted using an in-house code.

The optimized structures of the H₃⁺(H₂)_n clusters obtained by B3LYP+GD3/6-31G calculations are displayed in Fig. 1. All clusters have a structure in which the H₃⁺ core is surrounded by ligand H₂ molecule(s). The ligand H₂ molecule(s) are oriented perpendicular to the plane defined by the H₃⁺ core. The core–ligand distance in the H₅⁺ cluster (1.360 Å) is clearly shorter than those in the larger clusters (1.586 Å for H₇⁺, and 1.682 Å for H₉⁺). This is due to the low activation barrier for the migration of the inner hydrogen atom in the H₅⁺ cluster from one H₂ side to the other. The transition state (TS) structure is presented in Fig. 1(b), which is only 0.14 kcal mol⁻¹ higher in energy than the equilibrium structure (Fig. 1(a)). For comparison, the corresponding relative energy at the CCSD(T)/CBS level is 0.20 kcal mol⁻¹. The static electronic barrier for proton migration is very small and comparable to typical quantum contributions such as zero-point vibrational energy. When ZPE corrections are included, the energy of the TS is found to lie below that of the equilibrium structure, indicating that the apparent barrier becomes submerged. Specifically, the ZPE-corrected relative TS energies are −0.87 kcal mol⁻¹ at the B3LYP+GD3/6-31G level and −0.81 kcal mol⁻¹ at the CCSD(T)/CBS level. This clearly demonstrates that NQEs can qualitatively alter the relative energetics and the associated dynamics. Therefore, a treatment that explicitly accounts for NQEs, such as PIMD, is essential for a reliable description of this system.

Fig. 1 Optimized structures and representative interatomic distances [Å] of (a) H₅⁺ (C_2v), (b) H₅⁺ (D_2d, transition state), (c) H₇⁺ (C_2v), and (d) H₉⁺ (D_3h), obtained by static B3LYP+GD3/6-31G calculations.

Fig. 2 illustrates the definition of the angle ϕ between the plane of the H₃⁺ core and the axis of the ligand H₂ molecule. Before analyzing the PIMD and CLMD trajectories, we first estimated the rotational barriers for the ligand H₂ rotation in the H₅⁺, H₇⁺, and H₉⁺ clusters. The equilibrium and TS structures correspond to ϕ = 90° and 0°, respectively. The energy barriers are notably small: 0.21, 0.08, and 0.05 kcal mol⁻¹ for the H₅⁺, H₇⁺, and H₉⁺ clusters, respectively. The aforementioned shorter core–ligand distance in the H₅⁺ cluster may be responsible for its slightly higher rotational barrier.

Fig. 2 Schematic illustration of the definition of angle ϕ between the plane of the H₃⁺ core and the axis of the ligand H₂ molecule.

Since the rotational freedom around the H₃⁺ plane decreases as the angle ϕ approaches 90°, the observed angular distribution was normalized by a geometrical factor 2πr cos ϕ, which corresponds to the arc length available to the ligand at a given angle. The normalized distributions are shown in Fig. 3.

Fig. 3 Normalized distributions of the ligand H₂ rotation in (a) H₅⁺, (b) H₇⁺, and (c) H₉⁺ clusters obtained from PIMD and CLMD simulations.

The CLMD distributions exhibit a peak around ϕ = 90°, which is consistent with the expectation based on the aforementioned rotational potential energy curves (PECs). In contrast, the peak heights decrease in the PIMD distributions, resulting in broader profiles. These changes correspond precisely to the quantum plasticity reported by Štich and co-workers.^31,32

To reveal the mechanism underlying quantum plasticity, we conducted additional analyses. First, we estimated the zero-point vibrational energies of the normal modes for the ligand H₂ rotation. Fig. 4 displays the normal modes related to the ligand H₂ rotation of the H₇⁺ cluster. Two of them consist purely of ligand H₂ rotational modes, whereas the remaining two involve both ligand H₂ rotation and distortion of the H₃⁺ core structure. The sum of the zero-point vibrational energies of pure rotational modes is 0.84 kcal mol⁻¹, and the sum of those of these four modes is 1.72 kcal mol⁻¹, which are clearly larger than the rotational barriers for the ligand H₂ rotation (0.21, 0.08, and 0.05 kcal mol⁻¹ for the H₅⁺, H₇⁺, and H₉⁺ clusters, respectively). Therefore, the energy barrier for the H₂ rotation could be overcome by the NQEs (zero-point vibrational energies) associated with the rotation.

Fig. 4 Selected normal modes related to the ligand H₂ rotational motion. Normal mode vectors for these modes are listed in Table S2. Modes 1 and 5 consist purely of ligand H₂ rotation, whereas modes 2 and 6 involve both ligand H₂ rotation and distortion of the H₃⁺ core.

To further substantiate our interpretation of quantum plasticity as arising from the NQEs, we next analyzed the rotational PECs in more detail. First, the dependencies of the rotational barrier on the distance between the H₃⁺ core and the midpoint of the ligand H₂ (R_c-lMid), and on the covalent H₂ bond length of the ligand (R_ligand) were investigated, since interatomic distances fluctuate during the simulation and these fluctuations are enhanced by the NQEs. The rotational PECs were calculated while fixing either (a) R_c-lMid or (b) R_ligand at specific values (Fig. 5). These PECs demonstrated that elongation of R_c-lMid and contraction of R_ligand slightly reduce the rotational barrier. This result aligns with chemical intuition, as such structural changes should reduce repulsive interactions. From the above analyses, it can be deduced that the distributions of R_c-lMid and R_ligand may unveil the reason behind quantum plasticity. The one-dimensional distributions of R_c-lMid and R_ligand (Fig. 6) are highly informative in this regard. The one-dimensional distributions (Fig. 6(a) and (b)) clearly exhibit the NQEs on the R_c-lMid and R_ligand distances. The distributions are broadened by the NQEs, and the average values for R_c-lMid and R_ligand distances obtained from the PIMD simulation are larger than those obtained from the CLMD simulation. According to the perspective from Fig. 5, it can be speculated that the NQEs on the R_c-lMid distance facilitate the rotation of the ligand H₂, while the NQEs on the R_ligand distance inhibit the rotation.

Fig. 5 Potential energy curves for the rotational barrier as a function of (a) the distance between the H₃⁺ core and the midpoint of the ligand H₂, and (b) R_ligand calculated for the H₇⁺ cluster.

Fig. 6 One-dimensional distributions of (a) R_c-lMid and (b) R_ligand obtained from PIMD and CLMD simulations.

To visualize the net effect of the NQEs, Fig. 7 presents the two-dimensional distribution of R_c-lMid and R_ligand distances obtained from (a) CLMD and (b) PIMD simulations. At each R_c-lMid and R_ligand geometry, the corresponding rotational barrier was estimated as the energy difference between the structures at ϕ = 90° and 0°. The bold contour line corresponds to an energy difference of 0.0 kcal mol⁻¹ relative to the true rotational barrier, while the solid and dashed lines represent positive (unstable) and negative (stable) deviations, respectively. The contour map indicates that the barrier-lowering effect of the NQEs is mainly driven by the R_c-lMid distance.

Fig. 7 Two-dimensional distributions of R_c-lMid and R_ligand distances obtained from (a) CLMD and (b) PIMD simulations. The bold contour line corresponds to a zero energy difference relative to the true rotational barrier, while the solid and dashed lines represent positive (unstable) and negative (stable) deviations, respectively. The contour interval is 0.02 kcal mol⁻¹.

The PIMD distribution is significantly broadened by the NQEs, and the center of the distribution shifts toward the longer R_c-lMid region. In this region, the facilitation from R_c-lMid is greater than the inhibition from R_ligand. Our analysis clearly reveals that the rotation of the ligand H₂ is facilitated by the elongation of the R_c-lMid distance, thereby demonstrating that we successfully unveiled the origin of quantum plasticity for the first time.

To summarize, we analyzed the NQEs in H₃⁺(H₂)_n (n = 1–3) using the PIMD method. The NQEs induced an intriguing behavior, referred to as “quantum plasticity”, in which they facilitate the rotation of the H₂ ligands. In this study, we successfully unveiled the underlying mechanism of this phenomenon in the H_n⁺ clusters for the first time. The elongation of the distance between the H₃⁺ core and the midpoint of the ligand H₂ (R_c-lMid) caused by the NQEs is primarily responsible for quantum plasticity, because the decrease in the activation barrier for ligand H₂ rotation due to the elongation of R_c-lMid is greater than the inhibition arising from the elongation of R_ligand.

Thus, this study demonstrates that PIMD calculations employing DFT methods, which can reproduce high-level ab initio potential energy surfaces very well through appropriate choices of functionals and basis sets, are extremely useful and powerful tools for deepening our understanding of the NQEs in various systems, including H_n⁺ clusters.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article are included as part of the supplementary information (SI). Supplementary information: normal mode vectors for selected vibrational modes. See DOI: https://doi.org/10.1039/d6cp00022c.

Acknowledgements

This study was partly supported by JSPS KAKENHI, Grant No. 25K08559 (to T. U.), 23K13827 (to K. K.), 23K17905 (to M. T.), and 25H00428 (to T. U. and M. T.). We acknowledge the computational time allotted by the Research Center for Computational Science, Okazaki, Japan (Project No: 24-IMS-C349 and 25-IMS-C299 to T. U.).

Notes and references

1. R. N. Barnett and K. B. Whaley, Z. Phys. D, 1994, 31, 75–84.
2. M. V. Rama Krishna and K. B. Whaley, Z. Phys. D, 1991, 20, 223–226.
3. D. Scharf, G. J. b Martyna and M. L. Klein, Chem. Phys. Lett., 1992, 197, 231–235.
4. C. Chakravarty, Phys. Rev. Lett., 1995, 75, 1727–1730.
5. A. Castillo-García, A. W. Hauser, M. P. de Lasa-Castells and P. Villarreal, Molecules, 2021, 26, 5783.
6. P. Sindzingre, D. M. Ceperley and M. L. Klein, Phys. Rev. Lett., 1991, 67, 1871–1874.
7. P. Sindzingre, M. L. Klein and D. M. Ceperley, Phys. Rev. Lett., 1989, 63, 1601–1604.
8. Y. Kwon and K. B. Whaley, Phys. Rev. Lett., 2002, 89, 273401.
9. M. P. de Lasa-Castells and A. O. Mitrushchenkov, J. Phys. Chem. Lett., 2011, 2, 2145–2151.
10. H. Li, R. J. Le Roy, P. N. Roy and A. R. W. McKellar, Phys. Rev. Lett., 2010, 105, 133401.
11. S. Grebenev, B. G. Sartakov, J. P. Toennies and A. F. Vilesov, J. Chem. Phys., 2010, 132, 064501.
12. S. Grebenev, B. Sartakov, J. P. Toennies and A. F. Vilesov, Science, 2000, 289, 1532–1535.
13. M. P. de Lasa-Castells and A. O. Mitrushchenkov, Phys. Chem. Chem. Phys., 2021, 23, 7908–7918.
14. M. P. de Lasa-Castells and A. O. Mitrushchenkov, J. Phys. Chem. Lett., 2020, 11, 5081–5086.
15. D. Scharf, G. J. Martyna and M. L. Klein, J. Chem. Phys., 1993, 99, 8997–9012.
16. S. Broude and R. B. Gerber, Chem. Phys. Lett., 1996, 258, 416–420.
17. R. Clampitt and L. Gowland, Nature, 1969, 223, 815–816.
18. K. Hiraoka and P. Kebarle, J. Chem. Phys., 1975, 62, 2267–2270.
19. K. Hiraoka, J. Chem. Phys., 1987, 87, 4048–4055.
20. D. Mathur and J. B. Hasted, Nature, 1979, 280, 573–574.
21. M. Okumura, L. I. Yeh and Y. T. Lee, J. Chem. Phys., 1988, 88, 79–91.
22. M. Okumura, L. I. Yeh and Y. T. Lee, J. Chem. Phys., 1988, 83, 3705–3706.
23. S. L. Bennett and F. H. Field, J. Am. Chem. Soc., 1972, 94, 8669–8672.
24. R. Johnsen, C. Huang and M. A. Biondi, J. Chem. Phys., 1976, 65, 1539–1541.
25. M. T. Elford, J. Chem. Phys., 1983, 79, 5951–5959.
26. R. J. Beuhler, S. Ehrenson and L. Friedman, J. Chem. Phys., 1983, 79, 5982–5990.
27. T. Oka, Chem. Rev., 2013, 113, 8738–8761.
28. Y. Yamaguchi, J. F. Gaw and H. F. Schaefer, J. Chem. Phys., 1983, 78, 4074–4085.
29. U. Nagashima, K. Morokuma and H. Tanaka, J. Phys. Chem., 1992, 96, 4294–4300.
30. D. Marx and M. Parrinello, Z. Phys. B, 1994, 95, 143–144.
31. I. Štich, D. Marx, M. Parrinello and K. Terakura, Phys. Rev. Lett., 1997, 78, 3669–3672.
32. I. Štich, D. Marx, M. Parrinello and K. Terakura, J. Chem. Phys., 1997, 107, 9482–9492.
33. M. Barbatti and M. A. C. Nascimento, J. Chem. Phys., 2003, 119, 5444–5448.
34. R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Phys. Chem. A, 2003, 107, 4768–4772.
35. Y. Ohta, K. Ohta and K. Kinugawa, J. Chem. Phys., 2004, 121, 10991–10999.
36. R. P. de Tudela, P. Barragán, R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Phys. Chem. A, 2011, 115, 2483–2488.
37. M. Sugimoto, M. Shiga and M. Tachikawa, Comput. Theor. Chem., 2011, 975, 31–37.
38. P. Barragán, R. P. de Tudela, C. Qu, R. Prosmiti and J. M. Bowman, J. Chem. Phys., 2013, 139, 024308.

---