Nuclear quantum effects on intramolecular hydrogen bonds and backbone structures in biuret analogues

Abstract

Biuret analogues, dithiobiuret and thiobiuret, form six-membered ring structures via intramolecular hydrogen bonds. The proton donor and acceptor atoms differ between these molecules, leading to varying energy barriers for proton transfer. We performed path integral molecular dynamics (PIMD) simulations for dithiobiuret and thiobiuret to investigate the correlation between proton transfer and changes in the bond alternation pattern in the backbone structure. In addition, we performed PIMD simulations for deuterated species. The results indicate that the fluctuations of the backbone originate not from secondary nuclear quantum effects (NQEs) of the transferring hydrogen nucleus, but from the NQEs of the heavy nuclei in the backbone structure.

---

1. Introduction

The hydrogen bond is one of the important interactions that exists in various systems ranging from organic and inorganic compounds to living things.^1–4 In particular, a low-barrier hydrogen bond (LBHB), which has a low energy barrier for proton transfer, is known to play important roles in expression of biomolecular functions and enzyme catalysis.^5–8 Hydrogen bonds such as X–H⋯Y (where X and Y are highly electronegative atoms, such as O or N) have been extensively studied, in the context of not only normal hydrogen bonds but also the LBHB. On the other hand, hydrogen bonds involving sulfur, such as S–H⋯Y and X–H⋯S, exhibit different bonding properties and chemical behavior, suggesting that they play a unique role in proton transfer reactions and in determining molecular structures.^9–13 For example, Zhang et al. demonstrated through theoretical calculations that replacing one of the oxygen atoms, which acts as either a proton donor or a proton acceptor, with a sulfur atom in 3-hydroxy-2-(naphthalen-2-yl)-4H-chromen-4-one alters the molecular absorption and fluorescence wavelengths. This is expected to be applicable to the modulation of luminescence properties through sulfur substitution, as the direction of wavelength shift depends on whether the donor or acceptor is replaced.¹¹ Additionally, in terms of molecular structure, hydrogen bonds involving sulfur atoms have been shown to play an important role in stabilizing the α-helix structure of proteins.¹²

Recent advances in computational methods and computer technology have expanded the scope of theoretical studies, particularly electronic structure analysis, toward understanding chemical phenomena and predicting the activity of biofunctional molecules.^14,15 It should be noted that nuclear motions are generally neglected in theoretical calculations; thus, nuclear quantum effects (NQEs) are not taken into account. The Born–Oppenheimer approximation is usually valid because even the lightest nucleus, the proton, is about 1800 times heavier than an electron. However, NQEs become important in certain systems,^16–20 with one example being systems involving the LBHB.

The path integral molecular dynamics (PIMD) method is one of the computational approaches that can account for NQEs. The PIMD method has been employed to analyze a wide range of NQEs, such as H/D isotope effects, in various systems and reactions.^21–29 In the PIMD method, nuclei are represented as a ring of beads connected by virtual springs, enabling the incorporation of NQEs. Moreover, PIMD simulations can capture not only NQEs but also thermal effects, as the molecular dynamics is used to obtain statistical averages.

Recently, we analyzed the NQEs in curcumin, a molecule featuring an intramolecular LBHB, using the PIMD method. Our analysis revealed that, due to NQEs, the hydrogen-bonded proton tends to reside near the midpoint between the two oxygen atoms, corresponding to the transition state (TS) on the potential energy surface, as a diffusive quantum particle.²⁹ Furthermore, although it is generally believed that intramolecular proton transfer is almost completely coupled with changes in the bond alternation pattern in the OCCCO backbone (Fig. 1), our PIMD analysis showed that the coupling between them is relatively weak.²⁹ More recently, we extended our PIMD analysis to biuret and biguanide, both of which have a backbone containing nitrogen atoms. Notably, in biguanide, in which the proton donor and acceptor are nitrogen atoms, we found that π-electrons in the backbone are delocalized, despite the high proton transfer barrier.³⁰

Fig. 1 Intramolecular proton transfer and π-electron reorganization in the backbone.

Since the ease of proton transfer varies depending on the atomic species of the proton donor and acceptor atoms (X and Y in X–H⋯Y), extending this analysis to systems containing heteroatoms other than nitrogen may provide further insight into the behavior of proton transfer and the associated changes in the backbone structure. Such studies would enable a more systematic investigation of the relationship between proton transfer and changes in the bond alternation pattern within the backbone in more detail.

In this study, we investigated the NQEs in dithiobiuret and thiobiuret (Fig. 2), where sulfur atoms are involved in the intramolecular hydrogen-bonded structures. Furthermore, by analyzing the H/D isotope effect, we aimed to clarify whether the observed NQEs originate from the quantum effects of hydrogen nuclei or from those of heavier atoms.

Fig. 2 Structure and definition of atomic labels and structural parameters for dithiobiuret (X₁ and X₂ = S) and thiobiuret (X₁ = O and X₂ = S).

2. Computational methods

Since PIMD simulations require a large number of atomic force calculations, it is crucial to select a method that can reproduce the results of high-accuracy calculations within a reasonable computational time. Therefore, in this study, we chose density functional theory (DFT) to compute atomic forces. Table 1 summarizes the activation energies for intramolecular proton transfers in dithiobiuret, biuret and biguanide, as obtained using each method. Note that intramolecular proton transfer does not occur in thiobiuret. The activation energy is defined as the electronic energy difference between the stable structure and its corresponding TS structure for proton transfer. In this study, the activation energy calculated using the MP2/cc-pVTZ method is regarded as the reference value. To evaluate the reliability of the MP2/cc-pVTZ level for activation energy estimation, CCSD(T)/aug-cc-pVTZ single-point calculations were carried out on the MP2/cc-pVTZ optimized geometries of dithiobiuret, which exhibits an extremely small activation barrier. The resulting barrier of 0.4 kcal mol⁻¹ supports the validity of the MP2/cc-pVTZ method for this purpose. We tested the performance of four functionals, B3LYP, CAM-B3LYP, M06-2X, and APFD, with the 6-31G** basis set for evaluating the activation barriers. Among them, the B3LYP/6-31G** method most accurately reproduced the MP2/cc-pVTZ results. Therefore, we adopted the B3LYP/6-31G** method to calculate atomic force in the PIMD simulations.

Table 1 Activation energies [kcal mol⁻¹] for intramolecular proton transfer in dithiobiuret, biuret and biguanide obtained by several methods

Method	Dithiobiuret	Biuret^a	Biguanide^a
a Ref. 30. b CCSD(T)/aug-cc-pVTZ//MP2/cc-pVTZ.
MP2/cc-pVTZ	0.1 (0.4)^b	1.0	5.6
B3LYP/6-31G**	0.2	0.9	5.1
CAM-B3LYP/6-31G**	0.6	0.8	5.0
M06-2X/6-31G**	0.2	0.7	4.5
APFD/6-31G**	0.0	0.6	4.3

---

In the present PIMD simulation, each nucleus was represented by 16 beads to account for the quantum nature of nuclei. A time step of 0.1 fs was used, and 50 000 steps were sampled after thermal equilibration of 10 000 steps. For the symmetric molecules, namely, dithiobiuret, biuret, and biguanide, the distributions obtained from the simulations were folded and symmetrized with respect to the molecular symmetry axis. Since thiobiuret is asymmetric, 100 000 steps were sampled after thermal equilibration of 10 000 steps. In the classical molecular dynamics (CLMD) simulation, each nucleus was represented using a single bead, and 250 000 steps were sampled after 50 000 equilibration steps for symmetric molecules, whereas 500 000 steps after 50 000 equilibration steps were sampled for thiobiuret. All steps, except for the preliminary calculations for thermal equilibration, were used in the analysis. Both PIMD and CLMD simulations were carried out in the canonical ensemble (NVT) using a Nosé–Hoover chain³¹ at 300 K. All electronic calculations were carried out using the GAUSSIAN16 program package, while an in-house code was used to carry out the PIMD simulations.^32,33

Here, we introduce the parameter δR_XH to describe the position of the intramolecular hydrogen-bonded proton, as defined below:

δR_XH = R_X1H8 − R_X2H8, (1)

In the definition, δR_XH = 0 corresponds to the TS structure, in which the hydrogen-bonded proton is located exactly at the midpoint between X1 and X2. This parameter has been used to analyze the position of the hydrogen-bonded proton in various previous studies.^{16,19,29,30,34}

To discuss the extent of π-electron delocalization in the X1C3N4C5X2 backbone, we also introduced a π-delocalization index, λ:³⁵

λ = 1/2[n_1,2 + n_3,4], (2)

n_1,2 = 1/2[(n₁ − 1) + (2 − n₂)], (3)

n_3,4 = 1/2[(n₃ − 1) + (2 − n₄)], (4)

where n₁, n₂, n₃, and n₄ are the bond orders of the X1–C3, C3–N4, N4–C5, and C5–X2 bonds, respectively. The bond orders were calculated based on Peter's definition (eqn (5)):³⁶

d_ABn = r_A1 + r_B1 − 10|Δχ_AB| − (C_A + C_B − 17|Δχ_AB|)log n, (5)

where d_ABn is the actual bond length, r_A1 and r_B1 are the single covalent bond radii, C_A and C_B are the multiple bond parameters, Δχ_AB is the difference in electronegativities between elements A and B, and n is the bond order. For the single covalent bond radii, multiple bond parameters, and electronegativity of the C, N, O, and S atoms, the values reported by Peter et al.³⁶ were used (Table 2). According to the definition, λ = 0.5 corresponds to the fully π-delocalized structure (TS), while λ = 0 or 1 corresponds to the fully π-localized structure.

Table 2 Bonding parameters for C, N, O and S atoms³⁶

Atoms	Single bond radius,^ar	Multiple bond parameter, C	Electronegativity,^bχ
a Picometers. b Allred–Rochow scale.
C	77	35	2.50
N	73	38	3.07
O	74	45	3.50
S	103	29	2.44

---

3. Results and discussion

3.1. Dithiobiuret

Table 3 lists the optimized interatomic distances (R^EQ/TS_S1H8, R^EQ/TS_S2H8, and R^EQ/TS_S1S2) and calculated δR^EQ/TS_SH values. The δR^EQ_SH value is less negative than the δR^EQ_OH values for biuret, indicating that the hydrogen-bonded proton is located closer to the midpoint between the two sulfur atoms even in the EQ structure. The TS exhibits C_2v symmetry, meaning that the proton is located exactly at the midpoint between S1 and S2.

Table 3 Optimized interatomic distances [Å] at EQ and TS in dithiobiuret and biuret obtained at the B3LYP/6-31G** level of calculations

		EQ	TS
a Ref. 30.
Dithiobiuret	R ^EQ/TSS1H8	1.460	1.621
	R ^EQ/TSS2H8	1.856	1.621
	R ^EQ/TSS1S2	3.241	3.183
	δR^EQ/TSSH	−0.397	0.000
Biuret^a	δR^EQ/TSOH	−0.456	0.000

---

Fig. 3 shows the one-dimensional distributions of (a) R_SH, (b) R_S1S2, and (c) δR_SH for dithiobiuret obtained from PIMD and CLMD simulations. Due to the C_2v symmetry, both R_S1H8 and R_S2H8 are plotted together in Fig. 3(a). The average values from the PIMD and CLMD results are also shown in Fig. 3(b). The statistical errors shown in parentheses were estimated using the blocking method.³⁷ To symmetrize the distribution in Fig. 3(c), the values with inverted signs were plotted together with the original values.

Fig. 3 One-dimensional distributions of (a) R_SH, (b) R_S1S2, and (c) δR_SH, (d) relative free energies, and (e) π-delocalization index (λ) for dithiobiuret obtained via PIMD and CLMD simulations.

The distribution of R_SH (Fig. 3(a)) obtained from the CLMD simulation exhibits two distinct peaks around the R^EQ_S1H8 and R^EQ_S2H8 values. Compared to the R_OH distribution for biuret obtained from the CLMD simulation,³⁰ the population near the TS region is clearly larger in dithiobiuret. This is because the energy barrier for proton transfer is significantly lower. In the CLMD simulation, the position of the proton fluctuates due to thermal effects, resulting in a slight shift of the second peak away from the EQ structure. In contrast, in the PIMD-H results, the population near the TS region is substantially increased. As a result, the second peak (R_SH ≅ 2 Å) disappears, yielding a single-peaked distribution, consistent with the behavior observed in the PIMD results for biuret.

In addition to the R_SH distributions, the one-dimensional distributions of R_S1S2 obtained from the CLMD and PIMD-H simulations are compared in Fig. 3(b). While the peak positions in both cases are located near the R^EQ_S1S2 value, the CLMD distribution is slightly more extended toward longer distances. The average R_S1S2 value decreases when NQEs are included. This trend is also observed in biuret,³⁰ and the underlying mechanism is the same: quantum fluctuations make the hydrogen-bonded proton more likely to interact with the sulfur atoms.

Next, we focused on the one-dimensional distribution of δR_SH (Fig. 3(c)), which represents the position of the hydrogen-bonded proton. The distribution obtained from the CLMD simulation shows clear peaks around δR_SH = ±0.6 Å, slightly shifted from the δR^EQ_SH values. This shift observed in the CLMD results can be attributed to vibrational excitation induced by thermal effects. The SS stretching mode was found at 206.57 cm⁻¹. The corresponding characteristic vibrational temperature can be calculated as follows:³⁸

where h, [small nu, Greek, tilde], c, and k_B are the Planck constant, the vibrational wavenumber, the speed of light, and the Boltzmann constant, respectively. The characteristic vibrational temperature for 206.57 cm⁻¹ is calculated to be 297 K. Therefore, the SS stretching mode is thermally activated at a simulation temperature of 300 K. This activation leads to an elongation of the R_S1S2 distance, which weakens the hydrogen-bond interaction and induces both a contraction of the covalent S–H bond and an elongation of the hydrogen-bonded S–H⋯S distances. As a result, the peaks in the δR_SH distribution shift toward larger |δR_SH| values. In contrast to the CLMD result, where a clear valley appears at δR_SH = 0 Å, the PIMD distribution exhibits a single peak centered at δR_SH = 0, indicating delocalization of the shared proton.

The relative free energy profile along the proton coordinate, δR_SH, provides valuable insights into the behavior of the proton. The relative free energy, ΔF, can be estimated using the following expression:

ΔF(δR_XH) = −k_BT{ln(P(δR_XH)) − ln(P(δR^max_XH))}, (7)

where k_B is the Boltzmann constant, T is the temperature, and P(δR_XH) is the probability distribution as a function of δR_XH. P(δR^max_XH) represents the maximum P(δR_XH) value of the probability distribution. The free energy at this point is set to 0.0 kcal mol⁻¹ as a reference. The relative free energy profiles of dithiobiuret obtained from the CLMD and PIMD simulations are shown in Fig. 3(d). The CLMD free energy profile of dithiobiuret exhibits a double-peaked shape with a shallow barrier at δR_SH = 0.0 Å. In contrast, the PIMD-H simulation yields a single-well free energy profile for dithiobiuret, with a minimum at δR_SH = 0.0 Å, suggesting that the hydrogen-bonded proton does not experience an energy barrier for the proton transfer and can freely delocalize between the two sulfur atoms due to NQEs.

Fig. 4 shows the two-dimensional distributions of δR_SH and R_S1S2 for dithiobiuret, overlaid with a contour map of the corresponding potential energy surface. In the CLMD results, the distributions are located in relatively low-energy regions and are skewed around the EQ structures. Notably, the CLMD distribution exhibits two distinct peaks, corresponding to the proton being localized near either sulfur atom. Furthermore, high-population regions extend toward longer R_S1S2 values, which can be attributed to the aforementioned thermal excitation of the SS stretching mode. In contrast, the PIMD-H results exhibit a slightly more contracted distribution along the R_S1S2 axis, with a single peak centered at δR_SH = 0 Å. This increase in population occurs particularly at R_S1S2 values larger than those in the TS structure, indicating that proton migration tends to occur in the elongated R_S1S2 region during the PIMD-H simulation.

Fig. 4 Two-dimensional distributions of the δR_SH and R_S1S2 for dithiobiuret obtained from (a) PIMD-H and (b) CLMD simulations. The black circles and triangles correspond to EQ and TS structures, respectively, and the dashed line corresponds to the IRC pathway.

Fig. 3(e) shows the one-dimensional distribution of the π-delocalization index, λ, for dithiobiuret. It should be noted that the λ values for the EQ structure are λ^EQ = 0.405 and 0.595, indicating that the S1C3N4C5S2 backbone structure in dithiobiuret is significantly delocalized even in the EQ structure. The CLMD distribution exhibits a double peak, with a modest valley remaining at λ = 0.5. The peaks observed in the CLMD results are located near the λ^EQ values. In contrast, the distribution becomes completely single-peaked with a peak at λ = 0.5, in the PIMD-H result, suggesting that NQEs enhance the formation of π-delocalized structures, as also observed in the PIMD simulations of biuret.³⁰

Fig. 5 shows the two-dimensional distribution of δR_SH and λ values obtained from the PIMD and CLMD simulations of dithiobiuret. In the CLMD results, clear high-density regions appear around the EQ structures (δR^EQ_SH = ±0.397 Å and λ^EQ = 0.405 and 0.595). In contrast, the distribution obtained from the PIMD-H simulation does not exhibit distinct high-density regions. Instead, the population around the TS structure (δR^TS_SH = 0.0 Å and λ^TS = 0.5) increases, indicating that the hydrogen-bonded proton tends to be located near the midpoint between two sulfur atoms and that the π-electrons in the S1C3N4C5S2 backbone tend to be delocalized due to the NQEs. The distribution is spread along the λ direction, implying that the correlation between the extent of proton migration and the bond alternation in the S1C3N4C5S2 backbone is not particularly strong, similar to the case observed in biuret.

Fig. 5 Two-dimensional distributions of δR_SH and λ for dithiobiuret by (a) PIMD-H, (b) PIMD-D, and (c) CLMD simulations.

3.2. Thiobiuret

Table 4 lists the optimized interatomic distances in the EQ structure of thiobiuret, obtained at the B3LYP/6-31G** level of theory. The proton-transferred structure does not correspond to an energy-minimum structure; therefore, a TS structure does not exist for thiobiuret.

Table 4 Optimized interatomic distances [Å] at EQ in thiobiuret obtained at the B3LYP/6-31G** level of calculations

	R ^EQO1H8	R ^EQS2H8	R ^EQO1S2	R ^EQO1H8–R^EQS2H8
EQ	1.023	1.944	2.888	−0.921

---

Fig. 6 shows the one-dimensional distributions of (a) R^EQ_O1H8, (b) R^EQ_S2H8, (c) R^EQ_O1S2, and (d) R^EQ_O1H8–R^EQ_S2H8 for thiobiuret obtained from PIMD and CLMD simulations. The average values are also depicted in Fig. 6(c). First, we focus on the R_O1H8 distributions (Fig. 6(a)). The peaks of both PIMD-H and CLMD distributions appear near R^EQ_O1H8, whereas the PIMD-H distribution is broader and slightly shifted toward longer R_O1H8 values, due to the anharmonicity of the potential energy curve along the covalent bond length direction. In contrast, for the R_S2H8 distribution (Fig. 6(b)), the PIMD-H distribution is shifted toward shorter R_S2H8 values, opposite to the trend observed in the R_O1H8 distribution. One possible explanation is that elongation of the covalent R_O1H8 bond length at the proton donor site leads to a shortening of the R_S2H8 distance. In addition to the elongation of the R_O1H8, changes in the bond angles within the O1C3N4C5S2 backbone caused by NQEs also contribute to the contraction of the R_S2H8 distance. The angles of H8O1C3, O1C3N4, C3N4C5, and N4C5S2 are found to be slightly smaller in the PIMD-H results than in the CLMD simulations (Table 5), indicating that NQEs induce a closing distortion of the H8O1C3N4C5S2 six-membered ring structure. This structural change further contributes to the shortening of the R_S2H8 distance.

Fig. 6 One-dimensional distributions of (a) R_O1H8, (b) R_S2H8, (c) R_O1S2, (d) R_O1H8–R_S2H8, (e) relative free energies, and (f) π-delocalization index (λ) for thiobiuret obtained from PIMD and CLMD simulations.

Table 5 Average values of angles (ϕ) [degrees] in the O1C3N4C5S2 backbone in thiobiuret obtained from PIMD-H and CLMD simulations

	PIMD-H	CLMD
ϕ(H8O1C3)	106.6	107.0
ϕ(O1C3N4)	127.8	128.3
ϕ(C3N4S5)	121.6	122.1
ϕ(N4C5S2)	128.2	128.4

---

The one-dimensional distributions of the R_O1S2 distance obtained from the PIMD-H and CLMD simulations are similar in both shape and peak height. The average value of the PIMD-H distribution is slightly smaller than that of the CLMD distribution, consistent with the trends observed in dithiobiuret, biuret, and biguanide.

Both the PIMD-H and CLMD distributions of R_O1H8–R_S2H8 exhibit a single peak around R^EQ_O1H8–R^EQ_S2H8 (−0.921 Å), indicating that the hydrogen-bonded proton is always localized on the oxygen side and never migrates to the sulfur side. However, the NQEs extend the distribution towards R_O1H8–R_S2H8 = 0 Å.

Fig. 6(e) shows the one-dimensional relative free energy curves for thiobiuret. Unlike those for dithiobiuret, both the PIMD-H and CLMD free energy curves exhibit a single minimum. The free energy increases as the value of R_O1H8–R_S2H8 deviates from R^EQ_O1H8–R^EQ_S2H8. Although the free energy curve obtained from the CLMD simulation is nearly harmonic, the curve from the PIMD-H exhibits anharmonic behavior, becoming less steep as R_O1H8–R_S2H8 approaches 0 Å from negative values. This behavior indicates that NQEs cause the hydrogen-bonded proton to localize near the midpoint between the oxygen and sulfur atoms.

Fig. 7 shows the two-dimensional distributions of R_O1H8–R_S2H8 and R_O1S2 obtained from the PIMD-H and CLMD simulations. In the CLMD results, the distribution is spread along the R_O1S2 axis in the low-energy region centered around the EQ structure. In contrast, the PIMD-H results show a broader distribution along the R_O1H8–R_S2H8 axis than the CLMD results because the primary NQE of the hydrogen nucleus, which directly affects the R_O1H8 covalent bond length, is more pronounced than its secondary NQE acting on the heavy atom (R_O1S2).

Fig. 7 Two-dimensional distributions of the R_O1H8–R_S2H8 and R_O1S2 for thiobiuret obtained by (a) PIMD-H and (b) CLMD simulations. The black circle corresponds to EQ.

One-dimensional distributions of π-delocalization index (λ) obtained from the PIMD and CLMD simulations for thiobiuret are shown in Fig. 6(f). Both PIMD-H and CLMD distributions exhibit a single peak around λ^EQ, indicating that the π-electrons in the O1C3N4C5S2 backbone of thiobiuret are significantly delocalized even in the CLMD simulation. The λ distribution obtained from the PIMD-H simulation is broader than that from the CLMD simulation, indicating that the O1C3N4C5S2 backbone undergoes larger fluctuations, although the peak positions of the two distributions are nearly identical.

To gain further insights into the fluctuation of the backbone structure, Fig. 8 shows the two-dimensional distributions of R_O1H8–R_S2H8 and λ. Both the PIMD-H and CLMD distributions (Fig. 8(a) and (c)) are centered around R_O1H8–R_S2H8 = −0.9 Å and λ = 0.5, while the PIMD-H distribution is slightly broader along the λ axis. These results indicate that the π-electrons in the O1C3N4C5S2 backbone remain significantly delocalized, even though the hydrogen-bonded proton is localized on the oxygen side.

Fig. 8 Two-dimensional distributions of R_O1H8–R_S2H8 and λ for thiobiuret by (a) PIMD-H, (b) PIMD-D, and (c) CLMD simulations.

3.3. Deuterium isotope effect

To investigate the NQEs on the backbone structure in detail, PIMD simulations were performed for the deuterium-substituted compounds, in which H8 was replaced with deuterium. The two-dimensional distribution of δR_OH and λ for deuterated dithiobiuret is shown in Fig. 5(b) as the PIMD-D distribution. The PIMD-D distribution is more delocalized than the CLMD distribution; however, it is less delocalized than the PIMD-H distribution, indicating that a modest energy barrier remains in the PIMD-D simulation.

In the case of thiobiuret (Fig. 8), where proton transfer does not occur, the distribution around the EQ region is slightly denser in the PIMD-D distribution compared to the PIMD-H result. However, the overall distributions are very similar.

To gain deeper insights into the deuterium isotope effects, PIMD simulations were also conducted for deuterated biuret and biguanide. Fig. 9 shows the two-dimensional distributions of λ and δR_OH or δR_NH for biuret and biguanide, respectively. For biuret, the distribution obtained from PIMD-D shows low density near the TS region, similar to the case of dithiobiuret. In contrast, biguanide, which does not undergo proton transfer due to its high barrier, shows almost no distribution around δR_NH = 0, with populations clearly localized in either the positive or negative regions of δR_NH.

Fig. 9 Two-dimensional distributions of δR_OH and λ for biuret from (a) PIMD-H and (b) PIMD-D simulations and those of δR_NH and λ for biguanide from (c) PIMD-H and (d) PIMD-D simulations.

Including the distributions of dithiobiuret and thiobiuret shown in Fig. 5 and 8, we found that a significant H/D isotope effect appears along the δR direction only when the proton transfer barrier is low, as in dithiobiuret and biuret. In contrast, for biguanide and thiobiuret, where proton transfer does not occur or occurs only rarely, the distributions along the δR coordinate are clearly separated, and no substantial H/D difference is observed. Importantly, in all systems, no appreciable H/D isotope effect is found along the λ coordinate. This indicates that while the H/D isotope effect has a substantial impact on proton transfer, it does not significantly influence the fluctuations of the backbone. Furthermore, since the PIMD distributions for dithiobiuret and thiobiuret are clearly broader along the λ coordinate than those obtained from the CLMD, this broadening is attributed not to secondary NQEs of hydrogen nuclei, but rather to the NQEs of the heavy nuclei themselves.

4. Conclusions

In this study, we analyzed the NQEs on the intramolecular hydrogen bond and the X1C3N4C5X2 backbone of dithiobiuret and thiobiuret using the PIMD method.

In dithiobiuret, due to the low activation barrier, proton transfer readily occurs, and the proton is predominantly distributed near the TS on the Born–Oppenheimer potential energy surface at room temperature. In contrast, although proton transfer does not occur in thiobiuret, the λ value was found to fluctuate significantly, indicating that π-electrons in the backbone structure are delocalized. These results reveal that fluctuations in the backbone structure can occur, to some extent, independent of proton transfer.

To further investigate the structural fluctuation of the backbone, we examined the deuterium isotope effects in dithiobiuret, thiobiuret, biuret, and biguanide. The results showed that although the distributions are significantly broader along the λ coordinate in the PIMD simulations than in the CLMD simulations, no notable difference was observed between the PIMD-H and PIMD-D results. This indicates that the structural fluctuations originate not from the secondary NQE of the transferring hydrogen nucleus, but from the NQEs of the heavy nuclei in the backbone structure.

Conflicts of interest

There are no conflicts to declare.

Data availability

The supplementary information includes the Cartesian coordinates of the DFT-optimized structures, their total electronic energies, and the vibrational frequencies of the optimized structures. See DOI: https://doi.org/10.1039/d5cp02587g

The data supporting the findings of this study are presented in the main text. Additional data supporting this article has been included in the Supplementary information.

Acknowledgements

This study was partly supported by JSPS KAKENHI, Grant No. 23K13827 (to K.K.), 25H00428 (to M.T. and T.U.), 21H00026 and 23K17905 (to M.T.), and 25K08559 (to T.U.). We acknowledge the computational time allotted by the Research Center for Computational Science, Okazaki, Japan (Project Nos: 24-IMS-C349 and 25-IMS-C299 to T.U.). This work was also financially supported by JST SPRING, Grant Number JPMJSP2125. Furthermore, H.T. acknowledges the Interdisciplinary Frontier Next-Generation Researcher Program of the Tokai Higher Education and Research System for the support.

References

1. An Introduction to Hydrogen Bonding, ed. G. A. Jeffery, Oxford University Press, New York, 1997.
2. A. M. Beatty, Coord. Chem. Rev., 2003, 246, 131–143.
3. P. A. Frey, S. A. Whitt and J. B. Tobin, Science, 1994, 264, 1927–1930.
4. R. P. Sijbesma, F. H. Beijer, L. Brunsveld, B. J. B. Folmer, J. H. K. Ky Hirschberg, R. F. M. Lange, J. K. L. Lowe and E. W. Meijer, Science, 1997, 278, 1601–1604.
5. W. W. Cleland and M. M. Kreevoy, Science, 1994, 264, 1887–1890.
6. S. Yamaguchi, H. Kamikubo, K. Kurihara, R. Kuroki, N. Niimura, N. Shimizu, Y. Yamazaki and M. Kataoka, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 440–444.
7. V. Borshchevskiy, K. Kovalev, E. Round, R. Efremov, R. Astashkin, G. Bourenkov, D. Bratanov, T. Balandin, I. Chizhov, C. Baeken, I. Gushchin, A. Kuzmin, A. Alekseev, A. Rogachev, D. Willbold, M. Engelhard, E. Bamberg, G. Büldt and V. Gordeliy, Nat. Struct. Mol. Biol., 2022, 29, 440–450.
8. C. S. Cassidy, J. Lin and P. A. Frey, Biochemistry, 1997, 36, 4576–4584.
9. J. A. Platts, S. T. Howard and B. R. F. Bracke, J. Am. Chem. Soc., 1996, 118, 2726–2733.
10. Monu, B. K. Oram and B. Bandyopadhyay, Comput. Theor. Chem., 2023, 1225, 114133.
11. Z. Zhang and H. Fang, J. Lumin., 2024, 265, 120209.
12. P. Zhou, F. Tian, F. Lv and Z. Shang, Proteins, 2008, 76, 151–163.
13. A. Amadeo, M. F. Torre, K. Mráziková, F. Saija, S. Trusso, J. Xie, M. Tommasini and G. Cassone, J. Phys. Chem. A, 2025, 129, 4077–4092.
14. P. E. M. Siegbahn and R. H. Crabtree, J. Am. Chem. Soc., 1997, 119, 3103–3113.
15. H. Chen, M. Ikeda-Saito and S. Shaik, J. Am. Chem. Soc., 2008, 130, 14778–14790.
16. H. Tanaka, K. Kuwahata, M. Tachikawa and T. Udagawa, ACS Omega, 2022, 7, 14244–14251.
17. T. Meier, S. Petitgirard, S. Khandarkhaeva and L. Dubrovinsky, Nat. Commun., 2018, 9, 2766.
18. Y. Kanematsu, H. Kamikubo, M. Kataoka and M. Tachikawa, Comput. Struct. Biotechnol. J., 2016, 14, 16–19.
19. G. Cassone, J. Phys. Chem. Lett., 2020, 11, 8983–8988.
20. S. Dasgupta, G. Cassone and F. Paesani, J. Phys. Chem. Lett., 2025, 16, 2996–3003.
21. D. Marx and M. Parrinrllo, Z. Phys. B: Condens. Matter, 1994, 95, 143–144.
22. D. Marx, M. E. Tuckerman, J. Hutter and M. Parrinello, Nature, 1999, 397, 601–604.
23. P. Durlak and Z. Latajka, Phys. Chem. Chem. Phys., 2014, 16, 23026–23037.
24. M. Dračínský, L. Čechová, P. Hodgkinson, E. Procházková and Z. Janeba, Chem. Commun., 2015, 51, 13986–13989.
25. R. Pohl, O. Socha, P. Slavíček, M. Šála, P. Hodgkinson and M. Dračínský, Faraday Discuss., 2018, 212, 331–344.
26. M. E. Tuckerman, D. Marx, M. L. Klein and M. Parrinello, Science, 1997, 275, 817–820.
27. W. Fang, J. Chen, M. Rossi, Y. Feng, X. Li and A. Michaelides, J. Phys. Chem. Lett., 2016, 7, 2125–2131.
28. M. Daido, Y. Kawashima and M. Tachikawa, J. Comput. Chem., 2013, 34, 2403–2411.
29. T. Udagawa, H. Yabushita, H. Tanaka, K. Kuwahata and M. Tachikawa, Phys. Chem. Chem. Phys., 2023, 25, 10917–10924.
30. K. Nishikawa, H. Tanaka, K. Kuwahata, M. Tachikawa and T. Udagawa, Phys. Chem. Chem. Phys., 2024, 26, 24364–24369.
31. G. J. Martyna and M. L. Klein, J. Chem. Phys., 1992, 97, 2635–2643.
32. M. Shiga, M. Tachikawa and S. Miura, Chem. Phys. Lett., 2000, 332, 396–402.
33. M. Shiga, M. Tachikawa and S. Miura, J. Chem. Phys., 2001, 115, 9149–9159.
34. S. D. Pino, E. D. Donkor, V. M. Sánchez, A. Rodriguez, G. Cassone, D. Scherlis and A. Hassanali, J. Phys. Chem. B, 2023, 127, 9822–9832.
35. P. Gilli, V. Bertolasi, L. Pretto, A. Lyčka and G. Gilli, J. Am. Chem. Soc., 2002, 124, 13554–13567.
36. L. Peter, J. Chem. Educ., 1986, 63, 123–124.
37. H. Flyvbjerg and H. G. Petersen, J. Chem. Phys., 1989, 91, 461–466.
38. K. Suzuki, M. Shiga and M. Tachikawa, J. Chem. Phys., 2008, 129, 144310.

---