Microsolvation of triazole in water: structural and energetic behaviors

Abstract

In this study, the hydration process of the triazole heterocycle is explored, considering different isomers of triazole–(H₂O)_n (n = 1–8) clusters. Among the used (28) density functional (DFT) exchange–correlation functionals, the M11 functional exhibits the lowest mean absolute deviation (MAD) with respect to the high-level DLPNO-CCSD(T)/CBS method used as reference. Subsequently, we analyzed the relative stabilities of these clusters over a temperature range of 20–400 K. Results indicate that triazole–water clusters are stabilized by strong hydrogen-bond interactions, and that the most stable configurations adopt compact, folded, cage-like structures. The hydration free energies per water molecule (ΔG) and hydration enthalpies per water molecule (ΔH) at room temperature and over the 20–400 K have been obtained in the framework of the cluster-continuum solvation approach (SMD). Results show that the hydration free energy follows a linear dependence on temperature, while the enthalpy remains nearly constant. For the TZ(H₂O)₈ cluster, the estimated hydration free energy is computed to be 0.56 kcal mol⁻¹, which is in good agreement with the available experimental counterpart (≈0.48 kcal mol⁻¹), supporting the reliability of the used computational approach.

---

1. Introduction

The five-ring aromatic compound triazole contains three nitrogen atoms that allow it to act as both a hydrogen bond donor and acceptor, forming directional interactions with polar solvents such as water. Among the various biologically active heterocycles, triazole has attracted increasing interest due to its multiple applications in chemical,^1,2 biological^3,4 and pharmaceutical sciences, including its antifungal,⁵ antiviral,^6,7 anticancer,^8,9 anti-inflammatory^9,10 and antidepressant^11,12 properties. For these reasons, understanding its solvation mechanism in aqueous solutions and, in particular, the knowledge of its interaction mechanisms with water molecules, is of considerable importance. For this reason, some theoretical^13–26 and experimental^17,27–32 investigations have appeared in the literature. These studies have focused on the possible hydrogen bonds, possible structures and solvation energies, in particular.

The evaluation of the solvation free energy (ΔG) has been previously reported both in the gas phase^{13,14,16–18,20,22,25,26,29} and in aqueous solutions.^15,19,21,30

Recent theoretical studies²⁵ on triazole–(H₂O)_n=1–5 clusters, in the gas phase and aqueous solutions (PCM, SMD, CPCM), at both MP2/aug-cc-pVDZ and DFT (B3LYP/6-311++G*) levels of theory, suggest that the B3LYP functional performs comparably well with the MP2 results, supporting the earlier DFT findings of Abboud et al.¹⁷ A comparative analysis of solvation free energies (ΔG), derived from gas-phase microsolvation with explicit water molecules, and those obtained via continuum solvation models (PCM, SMD, and CPCM) demonstrated that the latter approaches yielded results that are close to the experimental counterpart (4.1 kcal mol⁻¹ ref. 29 and 0.48 kcal mol⁻¹ in gas phase and aqueous solution,³⁰ respectively). In order to achieve good agreement with experimental data, these studies also underline the importance of including a sufficient number of water molecules around the triazole in the computation. This finding underscores the critical role played by the intermolecular interactions and, in particular, hydrogen bonding between the solvent and solute, which allow accurate solvation energy values to be obtained. Although there are studies in the literature in which few water molecules surround the triazole,^{13–17,21,22,25} detailed investigations on the role of the hydrogen-bond (HB) networks in the triazole–(H₂O)_n clusters are rare. Furthermore, in these small cluster studies, a comprehensive exploration of the global minimum on the potential energy surface (PES) has not been performed. Based on the previous studies concerning the reliability of density functional theory (DFT) in accurately predicting solvation free energies,^17,25 we have chosen to employ this method to study the hydration process of triazole, considering a higher number of water molecules. For this purpose, we were interested in verifying which exchange and correlation potential to use for this type of system and phenomenon. We therefore used 26 different exchange–correlation functionals available in the literature to predict the structural parameters and the binding energies of TZ(H₂O)_n (n = 1–8) clusters of hydrated triazole. In order to identify the best-performing functional, the results have been benchmarked against DLPNO-CCSD(T)/CBS calculations. The most suitable ones have been used to study the geometrical structures and energetic behaviors of the considered clusters. Additionally, a temperature dependence of the isomeric distribution of triazole–(H₂O)_n clusters has been studied theoretically in a temperature range of T = 20–400 K. Solvation enthalpies (ΔH) and solvation free energies (ΔG) for the most stable configurations of the triazole–water clusters in both gas and water solution have been computed and compared with the available experimental data.

2. Methodology

2.1 Geometry sampling

As a first step, we searched for the most stable conformations of the TZ(H₂O)_n, (n = 1–8) triazole–water clusters. For each cluster size, a systematic exploration of possible geometries was conducted using the ABCluster (artificial bee colony algorithm for cluster global optimization) algorithm developed by Zhang and Dolg.^33,34 This algorithm employs an energy model that combines the Lennard-Jones and electrostatic potentials, with parameters taken from the CHARMM force field,³⁵ and ranks the generated structures from the most to the least stable based on their classical interaction energies. More details on the used structural sampling procedure are in the original publications of Zhang and Dolg^33,34 and in some of our recent works.^36–43

2.2 Temperature dependence

To assess the relative stability of isomers corresponding to a given cluster size, at a given temperature, we computed their canonical probabilities, defined as follows:where β = 1/k_BT (k_B, is the Boltzmann constant, and G^(k)_n(T) is the free energy of the kth isomer of cluster size n at temperature T, as obtained from the Tempo code.⁴⁴ Only temperatures equal to or above 25 K have been considered, since quantum effects become non-negligible at lower temperatures, and the canonical ensemble remains valid above this threshold.

2.3 Solvation free energy and enthalpy

The solvation free energy (ΔG) and solvation enthalpy (ΔH) of triazole clusters were computed both in the gas phase and in aqueous solution using the SMD continuum solvation model that combines explicit solvation, where the triazole molecule is directly coordinated by a limited number of water molecules, with an implicit bulk solvent modeled as a polarizable dielectric continuum. This hybrid approach takes into account the hydrogen bonding interactions in the first hydration shell and the long-range electrostatic effects due to the surrounding solvent environment.

The solvation of triazole can be represented by the following equilibrium reaction:

Triazole (g) + (H₂O)_n (s) → triazole (H₂O)_n (s) (2)

The corresponding solvation free energy (ΔG) and enthalpy (ΔH) per water molecule are calculated as follows:

ΔG_s(triazole)_n = ΔG_s[triazole(H₂O)_n] − ΔG_s[(H₂O)_n] − ΔG_gas(triazole)/n (3)

ΔH_s(triazole)_n = ΔH_s[triazole(H₂O)_n] − ΔH_s[(H₂O)_n] − ΔH_gas(triazole)/n (4)

where the subscript s denotes the aqueous solvent phase and g the gas phase. The convergence of solvation energies is assumed once the calculated values of ΔG_s(triazole)_n and ΔH_s(triazole)_n stabilize with increasing cluster size n, and the saturation of the first solvation.

2.4 Computational details

As is well known, the performance of the DFT method is strictly related to the use of a specific exchange–correlation functional. Due to the lack of information in the literature on which of the proposed functionals is the most suitable to give reliable results for our studied properties, we made a benchmark using 28 different exchange–correlation functionals and compared the results obtained with those coming from the high-level DLPNO-CCSD(T)/CBS reference method. All the most stable isomers of the TZ(H₂O)_n (n = 1–8) clusters coming from the use of the ABCluster algorithm have been considered. In particular, the following 26 DFT functionals have been tested: APFD,⁴⁵ B3LYP-D3,⁴⁶ B3PW91-D3,⁴⁶ B97D3,⁴⁷ BLYP-D3,^48,49 BMK-D3,⁵⁰ BP86-D3,⁵¹ BPBE-D3,⁵² CAM-B3LYP-D3,⁵³ LC-ωPBE-D3,⁵⁴ M05-2X-D3,⁵⁵ M05-D3,⁵⁵ M06-2X-D3,⁵⁶ M06-HF-D3,⁵⁶ M06-2X-D3,⁵⁶ M06-D3,⁵⁶ M08-HX,⁵⁷ M11,⁵⁸ MN12-SX,⁵⁹ MN15,⁶⁰ N12-SX,⁵⁹ PBE0-D3,⁶¹ PBE-D3,⁶² PW6B95-D3,⁶³ TPSS-D3,⁶⁴ ωB97X-D.⁶⁵ The Def2-TZVPP basis set, originally described by Weigend and Ahlrichs,^66–68 was employed throughout this study due to its proven accuracy in describing non-covalent interactions and its suitability for benchmarking DFT functionals against high-level ab initio methods, as reported by Malloum et al.⁶⁹

Vibrational analysis was performed to ensure that the characterized geometries are real minima. These frequencies were also employed to derive the free energies and enthalpies of the micro-solvated triazole clusters. All geometry optimizations and frequency calculations were carried out using the Gaussian 16 suite of programs.⁷⁰ To enhance the accuracy of the results, the ‘tight’ optimization criterion and an ultrafine integration grid were consistently applied.

All optimizations were initially conducted in the gas phase. For each cluster size (n = 1–8), only the global minimum structure was subsequently re-optimized using the solvent model density (SMD) approach.⁷¹

As described in the previous section, the cluster–continuum solvation model is intrinsically temperature dependent. Among the various isomers located, only those with significant Boltzmann populations contribute meaningfully to the ensemble average and, hence, to the computed solvation free energies and enthalpies.

3. Results and discussion

We begin this section by evaluating the performance of 26 density functional theory (DFT) functionals, analyzing the mean absolute deviations (MADs) against the DLPNO-CCSD(T)/CBS reference data (see Section 3.1). The most accurate functional was subsequently employed to explore the potential energy surface of TZ(H₂O)_n for n = 1–8 clusters (Section 3.2). The effects of temperature on the relative stability of the various isomers were then examined for each cluster size over the temperature range of 20–400 K (Section 3.3). In the last part of the section, the solvation free energy and enthalpy in both gas phase and aqueous environment are reported at room temperature as well as across the considered temperature range (20–400 K).

3.1 DFT benchmarking

The computed binding energy data of solvated triazole for each cluster size (n = 1–8) are provided in Table S1. The mean absolute deviations (MADs) of the 26 functionals employed are shown in Fig. 1 and Table S2. The MAD values of the tested hybrid functionals range from 0.7 to 7.5 kcal mol⁻¹, and the three most suitable DFT functionals follow the order M11 > MN15 > B97D3. Therefore, the M11 functional was found to be the most appropriate for our study on the triazole–(H₂O)_n clusters.

Fig. 1 Mean absolute deviations (kcal mol⁻¹) of DFT binding energies with respect to CCSD(T)/CBS used as the reference.

Interestingly, the M11/def2-TZVPP functional (MAD = 0.7 kcal mol⁻¹) provides more accurate binding energies than B3LYP, which had previously been suggested as a good compromise between accuracy and computational cost for predicting molecular structures and solvation energies.²⁵ It is worth noting that the improved performance of the M11 functional may be mainly attributed to its range-separated hybrid meta-GGA design, which incorporates a higher percentage of exact exchange at long range, enhancing the description of non-covalent interactions, such as bonding networks and dispersion forces.

It is worth noting that the calculated binding energies (see Table S1 in the SI section) for all density functionals increase with the cluster size n, indicating that the stability of the solvated system improves as more water molecules are added. This trend reflects the cumulative stabilization arising from additional hydrogen bonding interactions and enhanced solute–solvent cooperativity within the cluster.

3.2 Structures of hydrated triazole

The initial configurations, generated using ABCluster, were subsequently optimized at the M11/def2-TZVPP level of theory. The optimized structures of the hydrated triazole monomer and dimer are shown in Fig. 2. For the solvated triazole monomer, TZ(H₂O)₁, the water molecule forms a strong hydrogen bond (HB) with the triazole, with the water acting as a proton acceptor. In the case of TZ(H₂O)₂, the most stable minimum, TZW2_1, showed four hydrogen bonds, forming a cyclic structure (see Fig. 2). Water molecules are directly connected with the nitrogen sites, acting as hydrogen bond donors and acceptors. This result is consistent with a previous study of neutral water clusters, which highlights the strong stability of cyclic configurations for small clusters.⁷² The TZW2_2 isomer is 6.9 kcal mol⁻¹ higher in energy than the most stable one. This structure is characterized by only two hydrogen bonds: one OH–N bond resulting from the interaction between the first water molecule and the nitrogen site of the triazole ring, and one OH–O bond due to the interaction between water molecules. The lower stability of TZW2_2 can be attributed to the reduced number of hydrogen bonds²⁴ within the network, previously identified as a key factor influencing the stability of solvated triazole molecules.

Fig. 2 M11/def2-TZVPP-calculated structures of the TZ–(H₂O)_n=1, TZ–(H₂O)_n=2 and TZ–(H₂O)_n=3 complexes. Zero-point-corrected relative energies (kcal mol⁻¹) with respect to the most stable isomer are indicated in blue.

Ten initial geometries of TZ(H₂O)₃ have been optimized, resulting in five distinct stable structures (Fig. 2). These configurations all lie in a range of 7.6 kcal mol⁻¹, and the most stable is the TZW3_1, characterized by a folded cyclic configuration and a network of six hydrogen bonds (HBs). The stability of this isomer can be attributed to both the cyclic structure and the higher number of hydrogen bonds. Notably, the presence of hydrogen bonds that encapsulate the triazole sensibly stabilizes the cluster topology. This finding is consistent with a recent study by Neha Chopra et al.,⁷³ which demonstrates that the stability of azole–water complexes increases as the azole molecules become more embedded within extended hydrogen-bonded networks. The second and third most stable isomers, TZW3_2 and TZW3_3, are located at 0.6 and 1.1 kcal mol⁻¹ above the global minimum, and both structures are characterized by cyclic configurations. In TZW3_2 and TZW3_3 structures, five and four hydrogen bonds (HBs) are present, respectively, and two water molecules directly interact with the triazole, establishing only three HBs. In the TZW3_4 and TZW3_5 isomers, located at 1.8 and 7.6 kcal mol⁻¹ above the global minimum, only one water molecule interacts with the nitrogen site of the triazole ring.

Concerning the solvated triazole tetramer, TZ(H₂O), eleven stable configurations have been identified on its potential energy surface (PES); the energies range from 0.1 to 3.6 kcal mol⁻¹ (see Fig. 2). As for TZW₂ and TZW₃, in this case, the most stable structure adopts a cyclic configuration, forming a chain of water that interacts with the triazole ring (see TZW4_1 to TZW4_6 in Fig. 2). The two lowest-energy isomers, TZW4_1 and TZW4_2, are nearly isoenergetic and are stabilized by a network of eight hydrogen bonds: four involve only the water molecules and the other act between water and the triazole ring. Moreover, four nearly degenerate isomers, TZW4_3, TZW4_4, TZW4_5, and TZW4_6, lie 0.2, 0.4, 0.5, and 0.6 kcal mol⁻¹ above the global minimum, respectively. Each of these structures exhibits seven hydrogen bonds. Two additional branched cyclic isomers, TZW4_7 and TZW4_8, were located at 0.6 and 1.1 kcal mol⁻¹ above the minimum, respectively. As shown in Fig. 2, higher-energy structures such as TZW4_9, TZW4_10, and TZW4_11 adopt more branched topologies and are stabilized by only six and five hydrogen bonds, respectively.

Starting from 20 geometries for the solvated triazole–water pentamer, the geometry optimization returns 17 stable ones (see Fig. 3) in an energy range of up to 4.6 kcal mol⁻¹. The most stable isomer, TZW5_1, adopts a double cyclic configuration with a robust hydrogen-bonding network involving water–water interactions, and interactions between water molecules and the triazole ring. A detailed analysis of the TZW5 series reveals that the most stable isomers (TZW5_1 to TZW5_15) are characterized by a compact, cage-like topology. In contrast, higher-energy structures (TZW5_16 and TZW5_17) adopt more open geometries.

Fig. 3 M11/def2-TZVPP structures of the TZ–(H₂O)_n=5 complexes. Zero-point-corrected relative energies (kcal mol⁻¹) with respect to the most stable isomer are indicated in blue.

For the TZ–(H₂O)_n=6, we retained 18 low-energy structures from the 25 initially optimized conformers (see Fig. 4). The most stable isomers (from TZW6_1 to TZW6_1) are stabilized by a three-dimensional hydrogen-bonded network, in which six water molecules form a distorted cubic cage-like geometry interconnected by hydrogen bonds with the triazole ring. Despite their structural similarities, some differences in their relative stability arise from the spatial arrangement and connectivity between water molecules, as well as from their interactions with the triazole ring. In fact, TZW6_1 and TZW6_2 exhibit a compact hydrogen-bonded cage in which all six water molecules are strongly interconnected and tightly attached to the triazole ring via directional hydrogen bonds. In contrast, isomers from TZW6_3 to TZW6_6 show increasing geometric distortion, with the triazole ring becoming partially excluded from the hydrogen-bonding network. The fewer direct interactions between the triazole and the surrounding water molecules can contribute to the slightly reduced stability. In the TZW6_6 configurations, which lie more than 1.0 kcal mol⁻¹ above the global minimum, we note that they exhibit more open and fewer or weaker hydrogen-bond networks. This behavior again underlines the critical role of both the hydrogen-bond network and the incorporation of the triazole within the water cage in determining the overall stability of the hydrated triazole clusters.

Fig. 4 M11/def2-TZVPP structures of the TZ–(H₂O)_n=6 complexes. Zero-point-corrected relative energies (kcal mol⁻¹) with respect to the most stable isomer are indicated in blue.

For the TZ–(H₂O)₇ cluster, 25 different isomers were located on the potential energy surface. Among them, 15 structures within a relative energy cutoff of 5.8 kcal mol⁻¹ are presented in Fig. 5. As shown in this figure, nearly all isomers adopt compact, folded cage structures. The structural differences mainly arise from the spatial orientation and arrangement of the surrounding water molecules. The less stable isomers feature fewer and weaker hydrogen-bond interactions between the water molecules and the triazole ring. This observation highlights a clear correlation: the greater the number and strength of hydrogen bonds involving the triazole ring, the higher the stability of the corresponding isomer. This behavior is consistent with previous works on microsolvated heterocycles, where the cooperative hydrogen bonding plays a crucial role in the geometry stabilization.^74–76

Fig. 5 M11/def2-TZVPP structures of the TZ–(H₂O)_n=7 complexes. Zero-point-corrected relative energies (kcal mol⁻¹) with respect to the most stable isomer are indicated in blue.

For the TZ–(H₂O)₈ cluster, only the five most stable isomers were reported. The complete list of located TZ–(H₂O)₈ isomers is provided in the SI (see Fig. S1). Also, in these systems, the overall stability of the hydrated triazole clusters is primarily governed by the number and strength of hydrogen-bond interactions between the water molecules and the triazole ring. Furthermore, the spatial arrangement of the triazole within the three-dimensional water framework plays a crucial role in stabilizing folded cage structures (Fig. 6).

Fig. 6 M11/def2-TZVPP structures of the TZ–(H₂O)_n=8 complexes. Zero-point-corrected relative energies (kcal mol⁻¹) with respect to the most stable isomer are indicated in blue.

3.3 Temperature effects on the relative stability

The effects of temperature on the stability of the treated systems have been evaluated by using the Boltzmann distribution, as expressed in eqn (1). The probabilities of the different isomers of TZ–(H₂O)_n have been computed over a temperature range from 20 to 400 K. Results for the cluster dimensions ranging from n = 2 to n = 6 are reported in Fig. S2, while those for n = 7 and n = 8 are shown in Fig. 7. For the di-hydrated complex, the relative population at all considered temperatures is exclusively due to the most stable TZW2_1 isomer. The computations for the tri-hydrated complex (n = 3) (Fig. S2) indicate that the probability of the most stable triazole–water trimer, TZW3_1, decreases slightly with increasing temperature, while it remains the most favored isomer across all temperature range investigated. Isomers TZW3_2 and TZW3_5 contribute only marginally to the population at elevated temperatures (≥150 K). Overall, the population analysis shows that isomers TZW3_2 and TZW3_5 make a negligible contribution (less than 15% to the total population of the trimer). For the n = 4 water molecules, competition was observed between TZW4_1 and TZW4_2. TZW4_1 dominates at T < 50 K, while TZW4_2 becomes the most populated isomer from 50 K to 340 K, peaking at approximately 50% around 70 K before gradually decreasing (Fig. S2). Above 340 K, TZW4_4 becomes the most probable one, with a population of around 25%. At higher temperatures, the other isomers (TZW4_3, TZW4_5, and TZW4_6) contribute negligibly to the overall population of the cluster. We note that all the isomers contributing to the relative population of solvated triazole possess cyclic configurations, suggesting that the compact structures are energetically favored at lower temperatures. This agrees well with previous studies.⁷⁷

Fig. 7 Calculated relative populations for the solvated triazole clusters, TZ(H₂O)₇ and TZ(H₂O)₈ in water at the M11/def2-TZVPP level of theory. Populations are determined from Boltzmann distributions according to eqn (1).

For the penta-hydrated system TZW5, the most stable folded cage-like isomers, TZW5_1 and TZW5_2, dominate the cluster population. In particular, TZW5_1 prevails at temperatures below 115 K, while TZW5_2 becomes the most populated above this threshold. TZW5_3 and TZW5_4 clusters, with relative energies within approximately 1.0 kcal mol⁻¹, make a smaller contribution to the population of the triazole–water pentamer. As illustrated in Fig. S2(c), TZW5_3 appears at temperatures below 60 K with a population not exceeding 15%, whereas TZW5_4 starts contributing above 120 K and maintains a probability of up to 20% at high temperatures.

As for the population of penta-hydrated clusters, the temperature dependence of the TZW6, TZW7, and TZW8 clusters follows similar trends (Fig. S2 and S7). The most stable folded cage-like structures, with relative energy differences lower than 1.0 kcal mol⁻¹, compete with each other, with some being favored at low temperatures while others are favored at high temperatures. In this case, we can also underline that closed cage-like hydrogen-bonded interactions are a key factor influencing the stability of solvated triazole systems.

In conclusion, the stability of triazole–water clusters depends on several factors: firstly, increasing the number of water molecules enhances cluster stability due to the formation of additional hydrogen bonds; secondly, compact structures, which maximize interactions between water molecules and the triazole ring, tend to be more stable, especially at lower temperatures. Furthermore, the orientation of water molecules relative to the positions of hydrogen bond donor and acceptor sites influences the nature and strength of these interactions, thereby affecting the overall stability of the system.

3.4 Hydration free energy and enthalpy

The knowledge of the free energy and enthalpy provides precious insight into the solvation mechanisms, intermolecular interactions, and the overall stability of solute–solvent systems. Previous studies^25,78,79 have shown that the micro-solvation phenomena, such as the solvation free energies, can be reproduced theoretically, in good agreement with the experiment, if a sufficient number of water molecules are explicitly considered in the computations. For this reason, we have computed the hydration free energies including up to eight water molecules around the triazole molecule (see eqn (2)). The absolute values obtained for cluster sizes ranging from n = 1 to 8 at 298.15 K are summarized in Table 1.

Table 1 Solvation free energy of hydrated triazole in the gas phase as a function of the water molecules in the clusters. ΔG_T=298.15K

Solvation free energy ΔGT=298.15K
Cluster formula	Gas phase (kcal mol^−1)	Experimental data^29
TZ(H2O)1	−1.26631
TZ(H2O)2	3.50212
TZ(H2O)3	4.04365
TZ(H2O)4	1.7156	4.1
TZ(H2O)5	1.36483
TZ(H2O)6	3.67343
TZ(H2O)7	4.00098
TZ(H2O)8	2.46798

---

As shown in this table, the ΔG values fluctuate, indicating the limited accuracy of the model in which only solute–water interactions are considered. On adding more water molecules,²⁵ no significant improvement in the accuracy of hydration free energy predictions was observed. To overcome this problem, and following previous suggestions,²⁵ we also considered the bulk contributions in our computation, employing the solvation density model (SMD). Following this hybrid approach (more details are given in the Methods section), the hydration free energies and enthalpies were computed using eqn (3) and (4). For each cluster size n, only the global minimum of the TZ(H₂O)_n complex has been considered. Geometry optimizations have been redone at the M11/def2-TZVPP level of theory, and thermodynamic quantities have been evaluated at room temperature and in the range of 20–400 K.

The variations in hydration free energy (ΔG) and enthalpy (ΔH) at 298.15 K per water molecule for the TZ(H₂O)₇ and TZ(H₂O)₈ clusters are reported in Fig. 8, while the behaviors for TZ(H₂O)_n (n = 2–6) clusters are given in Fig. S3. Both properties exhibit oscillations around their respective mean values. The estimated average ΔG and ΔH for triazole are approximately 0.6 kcal mol⁻¹ and 1.0 kcal mol⁻¹, respectively. Notably, at n = 8, ΔG converges toward the experimental value³⁰ (0.48 kcal mol⁻¹) versus (0.56 kcal mol⁻¹). This suggests that a significant number of explicit water molecules are needed to closely approximate experimental data. These findings underscore the importance of accounting for both local (explicit) and long-range (continuum) electrostatic interactions when modeling solvation.

Fig. 8 Variations in the hydration free energy and enthalpy as a function of the cluster size at room temperature (298.15 K).

The impact of temperature on hydration thermodynamics can be obtained by studying the temperature dependence of hydration free energy and enthalpy over the 20–400 K range. The results for ΔG and ΔH across cluster sizes from n = 1 to 8 are shown in Fig. 9. Our analysis reveals that the hydration enthalpy remains nearly constant across the temperature range for a given cluster size, suggesting a weak temperature dependence. In contrast, the hydration free energy exhibits an almost linear increase with temperature (see Fig. 9), indicating that this behavior is mainly driven by entropy, as reflected in the slope of the free energy curves. This trend is consistent with our previous findings on the microsolvation of phenol in water,⁸⁰ where similar temperature behavior was observed: a linear increase in hydration free energy, and a nearly constant enthalpy. These consistent observations reinforce the conclusion that the entropy contribution plays a dominant role in the temperature dependence of hydration thermodynamics, while the enthalpic contribution remains largely invariant.

Fig. 9 Variations in the hydration enthalpy and solvation free energies as a function of the number of water molecules for the triazole micro-solvation.

4. Conclusion

In this work, we have explored the potential energy surfaces of triazole–(H₂O)_n clusters (n = 1–8) to investigate the solvation of triazole in water. To assess the performance of various density functional theory (DFT) methods in predicting the structural and energetic parameters of hydrated triazole clusters, we have tested 26 DFT exchange–correlation functionals. We have found that the M11 functional yields the smallest mean absolute deviation (MAD) compared to the other functionals tested, when benchmarked against DLPNO-CCSD(T)/CBS results. Thus, the M11 method was adopted for the detailed study of the triazole–water clusters.

The analysis of the relative stability of the triazole–(H₂O)_n isomers shows that a more extensive hydrogen-bonded network around the triazole, together with a more compact structure and an optimal orientation of hydrogen bonds, enhanced the overall stability of the system. The isomeric distribution as a function of temperature (20–400 K) shows that folded cage-like conformers are favored at lower temperatures.

An accurate reproduction of the experimental hydration free energy requires a cluster with a sufficiently large number of explicit water molecules and the consideration of the continuum solvation contribution.

The hydration free energy (ΔG) and the hydration enthalpy (ΔH) per water molecule in the range of 20–400 K indicate that ΔG decreases linearly with increasing temperature, while ΔH remains nearly constant.

For the TZ(HO)₈ cluster, the estimated hydration free energy is approximately 0.56 kcal mol⁻¹, which agrees well with the reported experimental value (≈0.48 kcal mol⁻¹).

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp03255e.

The authors declare that all the data are available.

Acknowledgements

This work was supported by the Tunisian Ministry of Higher Education and Scientific Research. The authors gratefully acknowledge the Center for High Performance Computing (CHPC, Grant Number. CHEM0947), South Africa, for providing the computational resources. We also thank the Department of Chemistry and Chemical Technologies at the University of Calabria, Italy, for their valuable collaboration and support.

References

1. A. R. Katritzky, C. A. Ramsden, J. A. Joule and V. V. Zhdankin, Handbook of Heterocyclic Chemistry, Elsevier, 2010.
2. S. Rachwal and A. I. Katritzky, Compr. Heterocycl. Chem. III, 2008, 1.
3. M. J. Vaishnani, S. Bijani, M. Rahamathulla, L. Baldaniya, V. Jain, K. Y. Thajudeen, M. M. Ahmed, S. A. Farhana and I. Pasha, Green Chem. Lett. Rev., 2024, 17, 2307989.
4. A. Lauria, R. Delisi, F. Mingoia, A. Terenzi, A. Martorana, G. Barone and A. M. Almerico, Eur. J. Org. Chem., 2014, 3289–3306.
5. N. Abdou, S. Soliman and S. Abou, Bull. Fac. Pharm., 1990, 28, 29–31.
6. J. H. Cho, D. L. Bernard, R. W. Sidwell, E. R. Kern and C. K. Chu, J. Med. Chem., 2006, 49, 1140–1148.
7. A. J. Srivastava, S. Swarup, V. Saxena and B. Chowdhury, Chem. Inf., 1992, 23, 1217–1227.
8. M. M. Alam, Arch. Pharm., 2022, 355, 2100158.
9. Z. Xu, S.-J. Zhao and Y. Liu, Eur. J. Med. Chem., 2019, 183, 111700.
10. P. S. Rao, C. Kurumurthy, B. Veeraswamy, G. S. Kumar, Y. Poornachandra, C. G. Kumar, S. B. Vasamsetti, S. Kotamraju and B. Narsaiah, Eur. J. Med. Chem., 2014, 80, 184–191.
11. F. Donato, M. G. de Gomes, A. T. R. Goes, N. Seus, D. Alves, C. R. Jesse and L. Savegnago, Life Sci., 2013, 93, 393–400.
12. M.-X. Song, B.-Q. Rao, B.-B. Cheng, H. Zeng, Y.-G. Luo and X.-Q. Deng, CNS Neurol. Disord.: Drug Targets, 2017, 16, 187–198.
13. J. Catalán, M. Sánchez-Cabezudo, J. L. G. De Paz, J. Elguero, R. W. Taft and F. Anvia, J. Comput. Chem., 1989, 10, 426–433.
14. F. Tomas, J. L. M. Abboud, J. Laynez, R. Notario, L. Santos, S. O. Nilsson, J. Catalan, R. M. Claramunt and J. Elguero, J. Am. Chem. Soc., 1989, 111, 7348–7353.
15. J. R. Cox, S. Woodcock, I. H. Hillier and M. A. Vincent, J. Phys. Chem., 1990, 94, 5499–5501.
16. C. Törnkvist, J. Bergman and B. Liedberg, J. Phys. Chem., 1991, 95, 3123–3128.
17. J. L. M. Abboud, C. Foces-Foces, R. Notario, R. E. Trifonov, A. P. Volovodenko, V. A. Ostrovskii, I. Alkorta and J. Elguero, Eur. J. Org. Chem., 2001, 3013–3024.
18. S. E. Murdock, R. Lynden-Bell, J. Kohanoff, C. Margulis and G. J. Sexton, Phys. Chem. Chem. Phys., 2002, 4, 5281–5288.
19. P. I. Nagy, F. R. Tejada and W. S. Messer, J. Phys. Chem. B, 2005, 109, 22588–22602.
20. W. P. Ozimiński, J. C. Dobrowolski and A. P. Mazurek, J. Mol. Struct., 2003, 651, 697–704.
21. H. A. Dabbagh, E. Rasti and A. N. Chermahini, J. Mol. Struct.: THEOCHEM, 2010, 947, 92–100.
22. R. M. Balabin, J. Chem. Phys., 2009, 131, 154105.
23. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618.
24. M. Khodiev, U. Holikulov, A. Jumabaev, N. Issaoui, L. N. Lvovich, O. M. Al-Dossary and L. G. Bousiakoug, J. Mol. Liq., 2023, 382, 121960.
25. W. P. Oziminski, Tetrahedron, 2013, 69, 3197–3205.
26. V. Jiménez and J. B. Alderete, J. Mol. Struct.: THEOCHEM, 2006, 775, 1–7.
27. J. Elguero, C. Marzin, A. Katritzky and P. Linda, Suppl. J. Heterocycl. Chem., 1976, 1, 124.
28. E. D. Raczyńska, W. Kosińska, B. Ośmiałowski and R. Gawinecki, Chem. Rev., 2005, 105, 3561–3612.
29. M. Bertrup, C. Nielsen, L. Nygaard, S. Samdal, C. Sjøgren and G. Sørensen, Acta Chem. Scand., 1988, 42, 500–514.
30. A. Albert and P. J. Taylor, J. Chem. Soc., Perkin Trans. 2, 1989, 1903–1905.
31. L. Xue, F. Zhao, X. Xing, Z. Zhou, K. Wang, S. Xu, J. Yi and R. Hu, J. Solution Chem., 2012, 41, 17–24.
32. M. Bellagamba, L. Bencivenni, L. Gontrani, L. Guidoni and C. Sadun, Struct. Chem., 2013, 24, 933–943.
33. J. Zhang and M. Dolg, Phys. Chem. Chem. Phys., 2015, 17, 24173–24181.
34. J. Zhang and M. Dolg, Phys. Chem. Chem. Phys., 2016, 18, 3003–3010.
35. K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes and I. Vorobyov, J. Comput. Chem., 2010, 31, 671–690.
36. A. Hattab, Z. Dhaouadi, S. Lahmar, A. Malloum, J. J. Fifen, N. Russo and M. Toscano, J. Mol. Liq., 2022, 360, 119414.
37. A. Hattab, Z. Dhaouadi, A. Malloum, J. J. Fifen, S. Lahmar, N. Russo and E. Sicilia, J. Comput. Chem., 2019, 40, 1707–1717.
38. A. Hattab, A. Malloum, Z. Dhaouadi and N. Russo, J. Comput. Chem., 2025, 46, e70066.
39. A. Hattab, Z. Dhaouadi, A. Malloum, J. J. Fifen, S. Lahmar, N. Russo and E. Sicilia, Theor. Chem. Acc., 2019, 138, 1–7.
40. A. Malloum and J. Conradie, New J. Chem., 2020, 44, 17558–17569.
41. A. Malloum, J. J. Fifen and J. Conradie, J. Chem. Phys., 2019, 150, 124306.
42. A. Malloum, J. J. Fifen and J. Conradie, J. Comput. Chem., 2020, 41, 21–30.
43. A. Malloum, J. J. Fifen, Z. Dhaouadi, S. G. N. Engo and J. Conradie, New J. Chem., 2019, 43, 13020–13037.
44. J. J. Fifen, M. Nsangou, Z. Dhaouadi, O. Motapon and N.-E. Jaidane, J. Chem. Phys., 2013, 138, 184304.
45. A. Austin, G. A. Petersson, M. J. Frisch, F. J. Dobek, G. Scalmani and K. Throssell, J. Chem. Theory Comput., 2012, 8, 4989–5007.
46. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
47. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
48. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098.
49. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
50. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184–32215.
51. T. Weymuth, E. P. Couzijn, P. Chen and M. Reiher, J. Chem. Theory Comput., 2014, 10, 3092–3103.
52. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
53. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
54. O. A. Vydrov, J. Heyd, A. V. Krukau and G. E. Scuseria, J. Chem. Phys., 2006, 125, 074106.
55. Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Phys., 2005, 123, 161103.
56. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
57. R. Li, J. Zheng and D. G. Truhlar, Phys. Chem. Chem. Phys., 2010, 12, 12697–12701.
58. R. Peverati and D. G. Truhlar, J. Phys. Chem. Lett., 2011, 2, 2810–2817.
59. R. Peverati and D. G. Truhlar, Phys. Chem. Chem. Phys., 2012, 14, 16187–16191.
60. H. S. Yu, X. He, S. L. Li and D. G. Truhlar, Chem. Sci., 2016, 7, 5032–5051.
61. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
62. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
63. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2005, 109, 5656–5667.
64. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
65. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
66. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
67. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
68. A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577.
69. A. Malloum and J. Conradie, J. Phys. Chem. A, 2021, 125, 2345–2356.
70. M. Frisch, G. Trucks, H. B. Schlegel, G. Scuseria, M. Robb, J. Cheeseman, G. Scalmani, V. Barone, G. Petersson and H. Nakatsuji, Gaussian, Gaussian Inc., Wallingford CT, 2016.
71. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
72. A. Malloum and J. Conradie, Polyhedron, 2021, 193, 114856.
73. N. Chopra, D. Kaur and G. Chopra, ACS Omega, 2018, 3, 12688–12702.
74. G. Qiu, D. Diao, L. Chaussy, S. Michaud-Chevallier, A. J. Simaan, P. Nava, A. Martinez and C. Colomban, Dalton Trans., 2022, 51, 10702–10706.
75. F. N. Keutsch and R. J. Saykally, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 10533–10540.
76. D. J. Wales and M. P. Hodges, Chem. Phys. Lett., 1998, 286, 65–72.
77. K. H. Lemke, Phys. Chem. Chem. Phys., 2021, 23, 18990–18998.
78. D. Riccardi, H.-B. Guo, J. M. Parks, B. Gu, L. Liang and J. C. Smith, J. Chem. Theory Comput., 2013, 9, 555–569.
79. L. Tomaník, E. Muchová and P. Slavíček, Phys. Chem. Chem. Phys., 2020, 22, 22357–22368.
80. A. Malloum and J. Conradie, Mol. Simul., 2023, 49, 403–414.

---