Theoretical analysis of H/D isotope effect in K₃H(SO₄)₂ and its influence on phase transition temperature

Abstract

K₃H(SO₄)₂ (KHS) is a typical zero-dimensional hydrogen-bonded dielectric material. KHS remains in a paraelectric state as the temperature decreases from room temperature. In contrast, its deuterium-substituted counterpart, K₃D(SO₄)₂ (DKHS), undergoes a phase transition to an antiferroelectric state at a critical temperature of 84 K. In this study, to clarify the mechanism underlying this isotope effect, path integral molecular dynamics simulations are performed on both KHS and DKHS systems, taking into account both thermal and nuclear quantum effects. The results show that quantum effects lower the free energy barrier, facilitating the distribution of hydrogen atoms between the two stable structures. In other words, quantum effects reduce the degree of order in the arrangement of hydrogen atoms and favor the paraelectric state, consistent with experimental results. Furthermore, quantum effects shift the hydrogen atom distribution toward the central position between two oxygen atoms, thereby drawing the oxygen atoms closer and shortening the oxygen–oxygen distance. This trend is consistent with experimental results, which indicate that the oxygen–oxygen distance in KHS is shorter than that in DKHS.

---

1. Introduction

Ferroelectric materials are widely used in various domains, such as ceramic capacitors, ignition devices, and non-volatile memory applications such as ferroelectric random-access memory.¹ Similarly, antiferroelectric materials have a broad range of applications, including capacitors, actuators, and energy storage devices.² In recent years, ferroelectricity has been discovered in organic crystals, leading to increased interest in hydrogen-bonded dielectrics owing to their potential for use in low-cost and environmentally friendly organic electronics. Notably, hydrogen-bonded dielectrics exhibit significant changes in phase transition temperature when hydrogen atoms are substituted with deuterium atoms. This isotope effect has been attributed to tunneling^3–5 and geometrical isotope effects.^6,7 However, a consensus regarding the mechanism underlying this effect remains to be reached, underscoring the need for further experimental and theoretical investigations. K₃H(SO₄)₂ (KHS) is a typical zero-dimensional hydrogen-bonded dielectric material (Fig. 1(a)). KHS remains in a paraelectric state even as the temperature decreases from room temperature. In contrast, its deuterium-substituted counterpart, K₃D(SO₄)₂ (DKHS), undergoes a phase transition to an antiferroelectric state at a critical temperature T_C of 84 K.⁸ In KHS, hydrogen atoms tend to be positioned closer to one of the two oxygen atoms in each hydrogen bond. This offset leads to unique electric properties owing to the distortion of surrounding atoms. A paraelectric state emerges when hydrogen atoms are either randomly distributed between the two oxygen atoms throughout the crystal or symmetrically located at the bond center. In contrast, an antiferroelectric state emerges when the displacement of hydrogen atoms in adjacent hydrogen bonds is anti-parallel. Therefore, the distribution of hydrogen atoms plays a crucial role in the H/D isotope effect in hydrogen-bonded dielectrics and its influence on the phase transition temperature.

Fig. 1 (a) Structure of K₃H(SO₄)₂ and (b) definition of δ_OH.

Experiment-based structural analysis^9–11 and theoretical studies^12,13 have extensively investigated this hydrogen distribution. X-ray structural analysis has indicated that at room temperature, the electron density distribution of hydrogen exhibits a double-peak profile, which changes to a single peak below 100 K. However, direct experimental observation of the behavior of hydrogen and deuterium atoms remains challenging. Theoretical studies have explored the influence of this isotope effect on hydrogen distribution by calculating the zero-point energy via harmonic approximation. These studies have shown that compared to deuterium atoms, hydrogen atoms are more likely to occupy central positions. However, the relevant models did not account for temperature effects, which are essential for investigating the phase transition temperatures.

Path integral molecular dynamics (PIMD) rigorously incorporates both temperature effects and quantum effects, making it a powerful tool for theoretical calculations of hydrogen-bonded dielectrics.^14,15 Therefore, this study aims to apply PIMD simulations to K₃H(SO₄)₂ at various temperatures to investigate influence of temperature on the behavior of hydrogen (or deuterium) in K₃H(SO₄)₂.

2. Computational details

2.1. PIMD simulations

To investigate the influence of nuclear quantum effects on the phase transition, we performed PIMD simulations for KHS and DKHS at 300 K, 200 K, 100 K, and 50 K. For simplicity, the PIMD simulations for KHS and DKHS are referred to as QM(H) and QM(D), respectively. These simulations incorporated the quantum effects of the atomic nucleus through bead expansion, using 16, 24, 48, and 96 beads for the simulations at 300 K, 200 K, 100 K, and 50 K, respectively. The simulation temperatures were controlled using the massive Nosé–Hoover chain method¹⁶ and all simulations were conducted using a program developed in our laboratory.^17,18 The mass of each beads ware determined as λ^(2s−1) = λ^(2s−2) = 2P[1 − cos{2π(s − 1)/P}]. A thermostat was applied to all degrees of freedom for all beads.

Each simulation was run for up to 53 000 steps with a time step of 0.2 fs. The first 3000 steps were used for equilibration, and the subsequent 50 000 steps were used for analysis. For comparison, we also performed PIMD simulations without quantum effects by setting the number of beads to 1 (denoted as CL) under the same conditions. These CL simulations ran for up to 210 000 steps with a time step of 0.5 fs. The first 10 000 steps were used for equilibration, and the subsequent 200 000 steps were used for analysis.

2.2. Electronic structure calculation

The forces and total energies at each step of the PIMD simulations were obtained using density functional theory (DFT), as implemented in the Vienna Ab initio Simulation Package (VASP 6.3.1).^19,20 The exchange–correlation functional was the revised Perdew–Burke–Ernzerhof (RPBE) functional,²¹ with a plane-wave cutoff energy of 500 eV and a k-point mesh of 1 × 2 × 2. Preliminary calculations indicated that the activation energy for hydrogen atom transfer obtained using the RPBE functional showed excellent agreement with that obtained using the SCAN+rVV10 functional. The SCAN functional combined with a dispersion correction has been reported to reproduce experimental structural parameters of hydrogen-bonded systems with high accuracy.²² Before performing the PIMD simulations, the lattice constants were determined from conventional molecular dynamics (CMD) simulations implemented in VASP. According to CMD simulations in the NPT ensemble at 1000 hPa, the average cell volumes were 492.2 ų, 495.9 ų, 500.1 ų, and 511.5 ų at 50 K, 100 K, 200 K, and 300 K, respectively. Detailed lattice constants are presented in Table S4 in the SI.

3. Results and discussion

The phase transition of KHS involves changes in the ordering of hydrogen atoms. To analyze the relative positions of hydrogen atoms in hydrogen bonds, we introduce the parameter δ_OH, representing the difference in the interatomic distances between the hydrogen and two oxygen atoms (Fig. 1(b)). Additionally, the parameter R_OO represents the distance between the two oxygen atoms in a hydrogen bond. Section 3.1, describes the distribution of hydrogen atoms. Section 3.2 discusses the results from Section 3.1 and phase transitions from the perspective of free energy surfaces. Sections 3.3 and 3.4 outline the relationship between the oxygen–oxygen distance R_OO and behavior of the hydrogen atom in terms of δ_OH.

3.1. One-dimensional distribution of δ_OH

Fig. 2 shows the one-dimensional distributions of δ_OH obtained from the PIMD simulations. The histograms were obtained from the structural data of all beads at all simulation steps (for example, 50 000 steps × 16 beads × 2 O–H bonds for the QM(H) simulation at 300 K). At 300 K (Fig. 2(a)), the CL simulation (black line) exhibits a double peak, indicating that in the absence of the nuclear quantum effects, the hydrogen atom is localized near one of the two oxygen atoms. In contrast, the distributions in the QM(H) (blue line) and QM(D) (red line) simulations extend toward δ_OH = 0, suggesting that nuclear quantum effects allow the hydrogen atom to overcome the activation barrier at the center of the hydrogen bond.

Fig. 2 One-dimensional distributions of δ_OH from QM(H) (blue), QM(D) (red), and CL (black) simulations at (a) 300 K, (b) 200 K, (c), 100 K, and (d) 50 K.

Next, we focus on the effects of temperature. In the CL simulation, the distribution height at δ_OH = 0 is 0.2 Å⁻¹ at 300 K. As the temperature decreases, the distribution at δ_OH = 0 diminishes, becoming nearly zero at 100 K. Simultaneously, the peak values around δ_OH = ±0.5 Å increase with decreasing temperature in the CL simulations. This indicates that in the absence of nuclear quantum effects, the hydrogen atoms are unable to overcome the activation barrier at low temperatures and remain at stable configurations. A similar trend is observed in the QM(D) simulations, although the change in the distribution is less pronounced. Furthermore, in the QM(H) simulations, the distribution remains nearly flat around δ_OH = 0, with minimal temperature dependence. This suggests that under the influence of nuclear quantum effects, hydrogen atoms can cross the activation barrier even at 50 K. At low temperatures, nuclear quantum effects become more prominent as kinetic energy decreases. Structural analysis experiments have suggested that the distribution of hydrogen atoms changes from a double-peak profile to a single peak at temperatures below 100 K.¹⁰ However, such a transition is not observed in our simulations. As discussed in Sections 3.3 and 3.4, the oxygen–oxygen distances in our simulation are slightly larger than those observed in experiments. This elongation likely contributes to the separation of the hydrogen distribution into two peaks instead of the emergence of a single peak.

3.2. Free energy surface

To determine whether the atomic motion can overcome the energy barrier, the barrier height must be compared with the kinetic energy. Therefore, we compared the free energy barrier with the kinetic energy per degree of freedom of a single particle . The free energy of hydrogen atoms is computed from the one-dimensional histogram of δ_OH shown in Fig. 2 using the following equation:²³

ΔF = −k_BT{ln(P(δ_OH)) − ln(P(δ^max_OH))},

where k_B, T, P, and δ^max_OH denote the Boltzmann constant, temperature, probability density, and δ_OH at maximum probability density, respectively. Fig. 3 shows the free energy surfaces as a function of δ_OH. At 300 K, all simulations exhibit a double-well shape with a barrier at δ_OH = 0. The free energy barrier in the CL simulation is approximately 1.17 kcal mol⁻¹, considerably larger than the thermal kinetic energy at 300 K (approximately 0.3 kcal mol⁻¹). In contrast, the free energy barriers in the QM(D) and QM(H) simulations are approximately 0.204 kcal mol⁻¹ and 0.104 kcal mol⁻¹, respectively, lower than the kinetic energy (approximately 0.3 kcal mol⁻¹). This indicates that quantum effects allow hydrogen atoms to overcome the barrier, contributing to their deviation from a single stable point. In other words, at 300 K, hydrogen and deuterium atoms can move freely between the two stable points and are randomly arranged throughout the crystal, leading to the paraelectric state. The results of the simulation at 200 K are similar to those at 300 K. At 100 K, the free energy barrier in the QM(H) simulation (0.041 kcal mol⁻¹) remains smaller than the kinetic energy (0.1 kcal mol⁻¹), whereas the free energy barrier in the QM(D) simulation (0.114 kcal mol⁻¹) exceeds the kinetic energy. This suggests that hydrogen atoms can move freely between the two stable points, forming a paraelectric state, whereas deuterium atoms remain localized at one of the stable points, enabling the emergence of the antiferroelectric state. The results at 50 K show the same trend as those at 100 K. These findings confirm that a significant isotope effect appears at temperatures below 100 K, which is qualitatively consistent with the experimentally determined phase transition temperature of 85 K.

Fig. 3 Free energy surfaces as a function of δ_OH at (a) 300 K, (b) 200 K, (c), 100 K, and (d) 50 K. The inset shows an enlarged view of the region where the free energy is lower than 0.2 kcal mol⁻¹. Green crosses represent the activation energy.

3.3. One-dimensional distribution of R_OO

Fig. 4 shows the one-dimensional distributions of the oxygen–oxygen distance (R_OO) and their expectation values, along with experimental values. Across all simulations, the expectation value of R_OO increases with increasing temperature. This trend is consistent with the results obtained from experiments,^9–11 as shown in (Fig. 4(e)). The CL simulations yield the largest R_OO values, followed by the QM(D) and QM(H) simulations, except at 300 K. At 300 K, the R_OO value from QM(D) is smaller than that from the QM(H). However, the difference between the two values falls within the range of the error bars, suggesting that it lies within the range of numerical error. These results indicate that nuclear quantum effects shorten the oxygen–oxygen distance. Notably, the experimentally observed value at room temperature for KHS (2.493 Å) is shorter than that for DKHS (2.519 Å),^9–11 consistent with the simulation results that indicate that quantum effects reduce the R_OO. However, the theoretical R_OO values from the QM(H) and QM(D) simulations at 300 K are 2.540 Å and 2.534 Å, respectively, slightly larger than the experimental values. Notably, in this work, the lattice constants are determined by CMD simulations in the NPT ensemble, which are slightly larger than those observed in experiments (Table S4 in SI), leading to an overestimation of R_OO in the simulations. Additionally, the same lattice constants are applied to both the QM(H) and QM(D) simulations. Consequently, the difference in oxygen–oxygen distances between QM(H) and QM(D) (0.006 Å) is smaller than the experimentally observed difference between KHS and DKHS (0.026 Å).

Fig. 4 One-dimensional distributions of R_OO and expectation values at (a) 300 K, (b) 200 K, (c) 100 K, and (d) 50 K. (e) Expectation values obtained from the simulations and experimental values of R_OO^9–11 as a function of temperature. Error bars obtained using the blocking method.²⁴

Overall, this analysis demonstrates that quantum effects shorten the oxygen–oxygen distance. However, as the quantum effects on the oxygen atoms are minimal, the decrease in R_OO in the QM simulations is primarily attributable to the influence of hydrogen atoms. Therefore, as discussed in the next section, we explore the relationship between δ_OH and R_OO.

3.4. Two-dimensional distributions of δ_OH and R_OO

Fig. 5 shows the two-dimensional distributions of δ_OH and R_OO obtained from the PIMD simulations. Across all simulations, a positive correlation is observed between |δ_OH| and R_OO. This correlation weakens as nuclear quantum effects increase from the CL to QM(H) simulations or as the temperature increases. This positive correlation is attributable to the hydrogen atom drawing both oxygen atoms inward when located centrally between them.

Fig. 5 Two-dimensional distributions of δ_OH and R_OO obtained from PIMD simulations.

Focusing on δ_OH = 0, we observe that the QM(H) and QM(D) simulations exhibit larger distributions at this position than the CL simulation. This is because quantum effects allow the hydrogen atoms to overcome the energy barrier. This larger distribution in δ_OH = 0 and the positive correlation between |δ_OH| and R_OO explain the shorter bond length of R_OO seen in Fig. 4. Specifically, nuclear quantum effects shift the hydrogen atom toward the center, and the hydrogen atom located at the center pulls the oxygen atoms at both ends, resulting in a shorter R_OO in the QM simulations.

4. Conclusion

We performed PIMD simulations on K₃H(SO₄)₂ to theoretically clarify the influence of the H/D isotope effect on the phase transition temperature. One-dimensional histograms of δ_OH show a prominent double peak in the CL simulation and broader distribution in the QM simulation, indicating that quantum effects promote hydrogen atoms to cross the activation barrier.

To further investigate the motion of hydrogen atoms, the free energy surfaces of hydrogen atoms were derived from the distribution of δ_OH.

The comparison of the barrier height and kinetic energy in the QM simulations at 300 K and 200 K shows that both QM(H) and QM(D) exhibit thermal kinetic energy higher than the free energy barrier. This means that hydrogen and deuterium atoms can move freely between the two stable points, leading to a disordered arrangement throughout the crystal, which corresponds to a paraelectric state. At temperatures below 100 K, differences between QM(H) and QM(D) become apparent. In the QM(H) simulations at 100 K and 50 K, the kinetic energy remains higher than the free energy barrier, resulting in the emergence of a paraelectric state. In contrast, in the QM(D) simulations, the free energy barrier exceeds the kinetic energy, suggesting that deuterium atoms remain localized at one of the stable points. This localization leads to an ordered arrangement throughout the crystal, giving rise to the antiferroelectric state. These findings confirm that a significant isotope effect appears below 100 K, consistent with the experimentally observed phase transition temperature of 84 K. This analysis highlights that comparing the free energy barrier and kinetic energy allows for a quantitative discussion of the phase transition temperature.

Furthermore, we investigated the influence of quantum effects on R_OO and its correlation with the distribution of hydrogen atoms. The results show that quantum effects decrease R_OO and a positive correlation exists between |δ_OH| and R_OO. This behavior is attributable to the quantum effects shifting the hydrogen atom toward a position closer to the transition state. The shifted hydrogen atom draws the oxygen atoms at both ends inward, reducing R_OO. This behavior is qualitatively consistent with the experimental results, which show that R_OO in KHS is shorter than that in DKHS.

Author contributions

Yuu Ishii: formal analysis; investigation; visualization; writing – original draft. Kazuaki Kuwahata: conceptualization; funding acquisition; investigation; methodology; software; validation; visualization; writing – review & editing. Tomomi Shimazaki: conceptualization; funding acquisition; resources; supervision. Masanori Tachikawa: conceptualization; funding acquisition; methodology; project administration; resources; software; supervision; writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available upon request.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5cp02885j.

Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers 22KJ2564 and 23K13827 (to K. K.), and 25H00428, 25H01271, and 25K22246 (to M. T.), and by JST, CREST Grant Number JPMJCR2234, Japan. Generous allocation of computational time from the Research Center for Computational Science, Okazaki, Japan (Project Number 23-IMS-C256 to M. T.) is also gratefully acknowledged.

Notes and references

1. R. Blinc, Ferroelectrics, 2002, 267, 3–22.
2. X. Hao, J. Zhai, L. B. Kong and Z. Xu, Prog. Mater. Sci., 2014, 63, 1–57.
3. R. Blinc, J. Phys. Chem. Solid, 1960, 13, 204–211.
4. M. Tokunaga and T. Matsubara, Prog. Theor. Phys., 1966, 35, 581–599.
5. M. Tokunaga and T. Matsubara, Ferroelectrics, 1987, 72, 175–191.
6. M. Ichikawa and K. Motida, J. Phys. Soc. Jpn., 1988, 57, 2217–2218.
7. E. Matsushita and T. Matsubara, Prog. Theor. Phys., 1982, 67, 1–19.
8. K. Gesi, J. Phys. Soc. Jpn., 1980, 48, 886–889.
9. Y. Noda, S. Uchiyama, K. Kafuku, H. Kasatani and H. Terauchi, J. Phys. Soc. Jpn., 1990, 59, 2804–2810.
10. Y. Noda, K. Kafuku, Y. Watanabe, H. Terauchi and K. Gesi, J. Phys. Soc. Jpn., 1990, 59, 3249–3253.
11. Y. Noda, H. Kasatani, Y. Watanabe and H. Terauchi, J. Phys. Soc. Jpn., 1992, 61, 905–915.
12. Y. Suwa, J. Yamauchi, H. Kageshima and S. Tsuneyuki, Mater. Sci. Eng. B, 2001, 79, 31–44.
13. Y. Suwa, J. Yamauchi, H. Kageshima and S. Tsuneyuki, Mater. Sci. Eng. B, 2001, B79, 98–112.
14. K. T. Wikfeldt and A. Michaelides, J. Chem. Phys., 2014, 140, 041103.
15. K. T. Wikfeldt, J. Phys.: Conf. Ser., 2014, 571, 012012.
16. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
17. M. Shiga, M. Tachikawa and S. Miura, Chem. Phys. Lett., 2000, 332, 396–402.
18. M. Shiga, M. Tachikawa and S. Miura, J. Chem. Phys., 2001, 115, 9149–9159.
19. G. Kresse and J. Furthmuller, Comput. Mater. Sci., 1996, 6, 15–50.
20. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
21. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413–7421.
22. J. G. Brandenburg, J. E. Bates, J. Sun and J. P. Perdew, Phys. Rev. B, 2016, 94, 115144.
23. H. Tanaka, K. Kuwahata, M. Tachikawa and T. Udagawa, ACS Omega, 2022, 7, 14244–14251.
24. M. Jonsson, Phys. Rev. E, 2018, 98, 043304.

---