In silico investigation of the microscopic structural features of a glucose-based deep eutectic solvent

Abstract

Atomistically detailed force field is used to investigate the microscopic structure of a naturally abundant deep eutectic solvent (NADES) composed of glucose, urea, and water in a 6 : 4 : 1 mass ratio at 328 K. The spatial distribution of the constituent molecules is analysed using pair correlation functions and by classifying the types of triangular spaces formed by three randomly chosen molecules. Key structural features such as water molecule clustering, relative orientation of urea molecular planes, and statistics of multi-hydrogen bonded electronegative atoms are also reported. The analysis reveals predominance of asymmetrically spaced molecules, with very few instances of symmetrical arrangements, such as molecular positioning at the vertices of equilateral or isosceles triangles. We also report significant deviations in the water structure from tetrahedrality, with water mostly existing in clusters of 4 to 5 molecules. Urea molecules predominantly adopt orthogonal relative orientations, a behaviour primarily driven by steric hindrance resulting from molecular crowding. Finally, strong hydrogen bond interactions between water and glucose are observed, with water oxygen atoms typically forming 2 to 3 hydrogen bonds, while glucose oxygen atoms form 1 to 2. In contrast, urea–urea interactions are favoured over interactions with glucose or water, with urea–water hydrogen bonding being the least preferred.

---

1. Introduction

The widespread use of conventional organic solvents in industrial processes presents serious environmental and safety concerns due to high volatility, toxicity, and persistence in the environment.^1,2 In response to these challenges, the search for greener alternatives has prompted the development of novel solvent classes, including gas-expanded liquids, room-temperature supercritical fluids, deep eutectic solvents (DESs), and ionic liquids.^3–7 Among these, DESs have emerged as particularly promising due to their low toxicity, cost-effective preparation, biocompatibility, and favourable physicochemical properties.^8,9 Deep eutectic solvents are formed from eutectic mixtures of Lewis or Brønsted acids and bases, which can involve a wide range of anionic and/or cationic species. A key feature of DESs is their extensive hydrogen-bonding network, which induces a significant depression of the melting point compared to the individual components.⁶ This phenomenon contributes to their ability to remain in the liquid phase under physiological conditions, thereby facilitating practical applications. In addition to this, DESs allow precise control over the viscosity and polarity of the reaction medium, a trait especially advantageous for maintaining the necessary conformational stability of flexible macromolecules during chemical reactions.^10,11 Hence, an in-depth understanding of the microscopic interactions and dynamic behaviour of DESs is critical for optimizing their performance in various applications.

A particularly eco-friendly subclass of DESs, known as naturally abundant deep eutectic solvents (NADESs), are derived from primary metabolites such as sugars, organic acids, amino acids, and urea.^8,12,13 NADESs cause significant interest owing to their inherent sustainability and green credentials.^14,15 For example, NADESs have been successfully employed as media for therapeutic systems and bio-transformation hosts.^16,17 Among the components used to synthesize NADESs, carbohydrates have been identified as particularly effective in creating biodegradable and non-ionic mixtures.^18–20 The physicochemical characteristics of carbohydrate-based DESs have been investigated and utilized in diverse sectors, notably as extraction media for important natural products like saponins, alkaloids, anthraquinones, flavonoids, etc.²¹ Computational studies on glucose-based DESs have revealed that a 1 : 1 molar ratio of glucose and choline chloride results in the most effective melting point depression. This can be attributed to the formation of a strong and extensive hydrogen-bonded network between the solutes.²² The strong interactions between choline chloride and glucose primarily arise from the chloride ion's simultaneous interactions with both choline and glucose molecules. This behaviour is reported by Sailau et al. in a combined density functional theory (DFT) and classical molecular dynamics study.²³ Based on this concept, Biswas and co-workers recently report a novel, fully biodegradable, non-ionic DES comprising glucose [C₆H₁₂O₆], urea [NH₂CONH₂], and water in a 6 : 4 : 1 mass ratio.²⁴ This mixture remains in the liquid state at 300 K and exhibits a glass transition temperature at around 236 K. A temperature-dependent dynamic fluorescence anisotropy study demonstrates that the relaxation time (τ) related to the rotational dynamics of hydrophobic and hydrophilic solutes in this NADES obeys a fractional dependence on viscosity (η), characterized by the relation . This behaviour indicates pronounced temporal heterogeneity and a deviation from classical hydrodynamic predictions. Fluorescence emission spectroscopy further reveals mild spatial heterogeneity in the system. Further insights into the microscopic dynamics of this NADES are provided by Baksi et al., who perform computational analysis to investigate how the structural and dynamic behaviour of water molecules is altered in this crowded environment.²⁵ Their results show that water molecules experience significantly restricted jumping–caging motions due to spatial confinement. Both centre-of-mass translational dynamics and hydrogen bond relaxation processes are observed to be manifold slower than that in bulk water. Recently, Müller-Plathe and co-workers elucidate the dependence of hydrogen bond relaxation times on the number of hydrogen bonds formed. The authors report that the relaxation dynamics are strongly connected to the hydrogen bond count for water and urea, while the slower translational and rotational motions of glucose molecules primarily govern the hydrogen-bonding behaviour.²⁶

In this simulation study, we investigate the glucose-based natural deep eutectic solvent system reported by Tarif et al. As discussed, previous experimental work on this system has primarily focused on the temperature-dependent behaviour of hydrophobic and hydrophilic solutes to probe microscopic-level heterogeneity.²⁴ In contrast, prior in silico investigations have mainly examined the dynamical properties of the constituent components.^25,26 Given that urea and glucose together account for over 90% of the system's mass, obtaining detailed insights into these two components is essential for a thorough understanding of the NADESs. To this end, we set up a simulation and computation protocol designed to characterize the system's microscopic structural features. Specifically, we investigate the extended hydrogen-bonded network, cluster populations, and the relative orientations of the constituents. This manuscript is structured as follows: Section 2 details the technical methodology of the simulations employed. Section 3 presents and discusses the primary findings of this computational investigation. Finally, Section 4 offers concluding insights and outlines directions for future research.

2. Model and simulation details

We investigate a NADES system composed of glucose (C₆H₁₂O₆), urea (NH₂CONH₂), and water in a 6 : 4 : 1 mass ratio, as reported in the experimental study conducted by Tarif et al.²⁴ Following the in silico study conducted by Baksi et al.,²⁵ we construct a system containing 66 glucose, 133 urea, and 111 water molecules, thereby maintaining the experimentally reported mass ratio. The cubic simulation box is constructed using PACKMOL²⁷ and populated with the specified number of molecules. A schematic representation of the molecular structures of the components, including relevant atom names, is provided in Fig. S1 of the SI.

We carry out all-atom molecular dynamics simulations using the GROMACS package (version 2021.1).^28,29 The force-field parameters for glucose and urea are obtained from the OPLS-AA³⁰ library, while water molecules are modeled using the SPC/E³¹ potential. The atomic partial charges for glucose and urea are adopted from Baksi et al.,²⁵ and the LigParGen³² web server is used to generate the required topology files. The system configuration created in PACKMOL is initially energy-minimized using the steepest-descent³³ algorithm to eliminate steric clashes. For the sake of better sampling, the minimized structure is randomized at 500 K under the canonical ensemble (NVT), followed by ten successive 1 ns NVT steps that gradually reduce the temperature to 328 K. These simulations employ the V-rescale³⁴ thermostat with a coupling time constant of 0.5 ps. Consequently, a 20 ns equilibration run is performed under the isothermal–isobaric (NpT) ensemble at 328 K, using the Berendsen³⁵ barostat with a 2 ps time constant to maintain an average pressure of 1 atm. The simulation temperature is selected based on the availability of both experimental²⁴ and computational²⁵ data. The simulated density at 328 K is found to be 1.38 g cm⁻³. Densities calculated across a range of temperatures are presented in Table S1 of the SI. At all temperatures, the calculated densities deviate by less than 2% from the experimental²⁴ value. A spherical cutoff of 1.4 nm is employed for real-space electrostatic and van der Waals interactions. Periodic boundary conditions are applied in all three spatial directions to eliminate edge effects. Short-range Lennard-Jones interactions are treated using the minimum image convention,³⁶ and long-range electrostatics are computed via the particle mesh Ewald (PME)³⁷ method. All bond lengths are constrained using the LINCS algorithm.³⁸ The production run spans 200 ns, utilizing a 1 fs integration time step with the leap-frog³⁶ integration algorithm. The trajectory data are saved at every 2 ps. System properties are analysed using built-in GROMACS modules²⁸ along with custom analysis scripts. Molecular visualizations are performed using VMD version 1.9.3.³⁹

3. Results and discussion

3.1. Distribution of the constituents

3.1.1. Calculation of the radial distribution function. We present the radial distribution function, g(r), calculated between different atomic pairs of the constituents in Fig. 1. The coordination number obtained by integrating the first peak of g(r) is provided in Table S2 of the SI. In Fig. 1a, the distribution of water oxygen atoms (OW) around a central OW atom is shown. A prominent peak is observed at approximately 0.28 nm, corresponding to the first solvation shell of a central water oxygen atom. The structured nature of this peak indicates a closely packed arrangement of water molecules, suggesting potential clustering behaviour. Both the position and intensity of this peak are consistent with previous in silico findings.²⁵ The coordination number analysis indicates that a central water oxygen atom (OW), within the distance corresponding to the first minimum of the radial distribution function, is surrounded by approximately 1–2 neighbouring OW atoms, as reported in Table S2. In Fig. 1b, the intermolecular interactions within urea are analysed. The nitrogen–oxygen pair (NU–OU, black curve) exhibits a peak at approximately 0.31 nm, indicating a relatively strong interaction. In contrast, the nitrogen–nitrogen interaction (NU–NU, red curve) yields a broader and weaker peak centred around 0.35 nm, suggesting a lower probability of close approach between NU atoms. The comparatively higher peak for NU–OU implies a preferential spatial association of NU atoms with OU atoms. This can be attributed to complementary hydrogen bonding characteristics, where NU functions as a hydrogen bond donor and OU as an acceptor. Integration of the first peak yields coordination numbers of approximately 3 for the NU–OU pair and approximately 2 for the NU–NU pair. The oxygen–oxygen distribution (OU–OU, blue curve) displays a minor peak at ∼0.39 nm and a more prominent peak at ∼0.5 nm, indicating limited and less favourable self-association among OU atoms. The first small peak corresponds to a coordination number of ∼1.38. The spatial distribution of the urea molecules' centre-of-mass is represented by the carbon atoms (CU–CU, green curve), which shows a distinct peak around 0.44 nm. This peak exhibits the highest intensity among all atom–atom correlations within urea, indicating a significant spatial ordering of CU atoms with the respective coordination number of ∼6. Fig. 1c presents the radial distribution of the same oxygen atoms of different glucose molecules. The exocyclic oxygen atoms (OG1, OG2, OG3, OG4, and OG5) exhibit a pronounced peak near 0.3 nm, reflecting favourable proximity and possible hydrogen bonding due to their exposure to a less sterically hindered environment. Conversely, the ring oxygen atom (OG6) shows no discernible structured peak, likely due to its limited accessibility and lack of hydrogen bonding capability. The diffuse nature of its g(r) curve indicates a low degree of positional correlation. Integration of the first peak of the radial distribution function yields coordination numbers in the range of approximately 0.1–0.4 for these glucose atom pairs, with the lowest value corresponding to the cyclic oxygen atom. Additional radial distribution functions for other atomic pairings within glucose are provided in Fig. S2(a–c) of the SI.

Fig. 1 Radial distribution function calculated between different atoms of (a) water, (b) urea, (c) glucose, (d) urea–water, (e) glucose–urea, and (f) glucose–water molecules.

The spatial distribution of electronegative atoms from water and urea molecules is presented in Fig. 1d. As shown, the water oxygen atoms (OW) are similarly distributed around both the nitrogen (NU) and oxygen (OU) atoms of urea, with g(r) peaks centred at approximately 0.3 nm. Depending on the choice of the central atom, the calculated coordination numbers range from approximately 1.0 to 1.4, as shown in Table S2. This observation suggests that water, possessing both hydrogen bond donor and acceptor capabilities, interacts comparably with NU and CU. The position and intensity of the g(r) peak are consistent with the findings of Baksi et al.²⁵ In Fig. 1e, the radial distribution function between glucose oxygen atoms (OG) and urea nitrogen atoms (NU) is presented. A distinct peak at 0.3 nm is observed across all possible OG–NU combinations. Although the ring oxygen (OG6) can function as a hydrogen bond acceptor, its location in a sterically crowded environment results in the lowest NU population around OG6 compared to other exocyclic OG atoms. The coordination number analysis reveals that each glucose oxygen atom (OG) is surrounded by approximately 0.7–1 nitrogen atoms of urea (NU), whereas a central NU atom is, on average, coordinated by exactly half as many OG atoms. The g(r) profiles for OG–OU interactions are included in Fig. S2d of the SI. Lastly, the g(r) between OG and OW atoms is computed and shown in Fig. 1f. This distribution exhibits a prominent peak at approximately 0.28 nm, closely aligning with the OW–OW g(r) peak, indicating a similar spatial proximity. Among all g(r) peaks discussed so far, the g(r) peak of OG–OW is the second most intensified, following the OW–OW interaction with the coordination number approximately equal to or below 1. A close inspection of all the g(r) peaks reveals the presence of well-defined short-range ordering and a lack of long-range structural correlations. This pattern is characteristic of a liquid-like structural organization within the system. To examine the size-independence of the results presented in this study, an additional 10 ns simulation is performed on a system twice the original size, comprising 132 glucose molecules, 266 urea molecules, and 222 water molecules. The comparative results, provided in Fig. S3 of the SI, demonstrate a strong consistency with those obtained from the smaller system.

3.1.2. Occupancy of triangular space. To investigate local structural organization, we randomly select three constituent molecules and classify the geometry of the triangular space they span based on their mutual distances. This triplet-based analysis provides insight into the relative spatial arrangement of three molecules, which is not captured by conventional pair-correlation functions. The analysis of such triangular configurations is particularly significant as it yields information about locally dense packing structures, which are closely linked to the thermodynamic behaviour of complex fluids.^40,41 For instance, Dhabal et al. employed triplet-correlation functions in systems comprising room-temperature ionic liquids and high-temperature ionic melts to characterize angle-resolved atomic or ionic organization.⁴² Their findings reveal distance-dependent correlations among specific ionic triplets, with a focus on those forming equilateral triangles. Following this rationale, we categorize the various types of triangular configurations formed by molecular triplets and quantify their occurrence probabilities. In our system, ten distinct molecular triplets can be formed from water (W), urea (U), and glucose (G): W–W–W, W–W–U, W–W–G, U–U–U, U–U–W, U–U–G, G–G–G, G–G–W, G–G–U, and W–U–G. For each of these, we assess the frequency of four possible structural types: equilateral, isosceles, scalene, and collinear. No instances of collinear configurations are observed, and thus, we report the probabilities of occurrence for the remaining three categories in Fig. 2a. To determine side length equality in equilateral and isosceles triangles, we apply a tolerance threshold of 0.01 nm.

Fig. 2 Population distribution of (a) occurrences of different triangular spaces formed by three randomly chosen molecules and (b) the perimeter length of equilateral triangles formed by (I) minimum two water molecules, (II) minimum two urea molecules, (III) minimum two glucose molecules, and (IV) simultaneously by water–urea–glucose molecules.

As shown in Fig. 2a, all molecular triplets exhibit similar distributions among the triangle types. Scalene triangles dominate with an occurrence probability slightly exceeding 97.5%, while isosceles triangles occur in approximately 2.47% of cases. Equilateral triangles are the least common, with probabilities on the order of 0.01%. Despite their rarity, equilateral configurations are of particular interest, as they might suggest a high degree of symmetry and structural order. In Fig. 2b(I), we present the population distribution of equilateral triangles composed of at least two water molecules, categorized by the perimeter length. The distribution follows a Gaussian profile, spanning from 1 nm to 10 nm, with a peak near 6 nm and an associated probability of occurrence between 10% and 15%. This Gaussian behaviour implies that perimeter variations arise from random, independent fluctuations. Notably, molecules separated by ∼2 nm most frequently form equilateral triangles. When compared to the peak positions in the radial distribution functions (see Fig. 1), it becomes evident that closely spaced molecules rarely form equilateral triangles, likely due to spatial constraints resulting from molecular crowding. Fig. 2b(II) shows the distribution for equilateral triangles formed by at least two urea molecules. The overall trend mirrors that of the water-based configurations, except for the U–U–U triplets. Here, a notable population (∼8%) appears at perimeter lengths below 2 nm, comparable to the peak near 6 nm. A perimeter less than 2 nm corresponds to a molecular separation of ∼0.6 nm, which aligns with the first minimum in the urea–urea radial distribution function (Fig. 1b), suggesting a higher tendency for urea–urea association. In Fig. 2b(III), we show the population distribution of equilateral triangles containing at least two glucose molecules. The distribution again follows a Gaussian pattern, but the peak (population: ∼10%) shifts to a lower range (4–5 nm) compared to water or urea. This shift can be attributed to the larger molecular size of glucose, which results in greater spatial occupancy; consequently, the probability of encountering a glucose molecule at a relatively shorter distance increases significantly. Finally, Fig. 2b(IV) illustrates the distribution for equilateral triangles formed by W–U–G triplets. The resulting distribution is Gaussian having a peak with a population of ∼10% near 6 nm.

In contrast to previous studies on ionic liquid systems, which reported equilateral triangle formation at molecular separations below 1 nm, marking highly packed symmetric arrangements,⁴² our findings indicate that, in densely crowded systems, a minimum separation of ∼2 nm is required to form equilateral triangles. Therefore, in our system, equilateral configurations at these larger separations do not signify closely packed structures. Rather, configurations at shorter distances exhibit greater asymmetry, reflecting enhanced spatial disorder.

3.2. Structural characteristics of water molecules

3.2.1. Quantification of water clusters. A set of n water molecules is defined as a cluster if the oxygen atoms (OW) of the constituent molecules are within a cut-off distance of 0.35 nm, corresponding to the first minimum in the OW–OW radial distribution function g(r), with n ≥ 3. In the present calculation, it should be noted that the results are strongly dependent on the choice of the distance cut-off. As evident from Fig. 1a, the peak intensity at 0.28 nm is nearly five times greater than that of the second peak observed at ∼0.45 nm. This indicates that any water clustering, if present, can only be captured within such a short-distance cut-off. In contrast, selecting a distance cut-off in the range of 0.4–0.6 nm, corresponding to the second peak, yields either a negligible number of water clusters or none, which is consistent with the very low intensity of the second peak. Fig. 3 presents the analysis of water clusters, confirming their presence in the system.

Fig. 3 Population distribution of (a) the number of water clusters present and (b) the number of water molecules in each cluster.

Fig. 3a shows the population distribution of the number of simultaneous clusters, which follows a Gaussian profile ranging from 3 to 17 clusters, peaking at approximately 10 clusters with a ∼17% probability. This indicates that the most thermodynamically favourable configuration under the given system conditions (e.g., temperature, system size, and density) consists of 9–10 water clusters, although fluctuations around this value are common. The Gaussian nature of the distribution suggests that cluster formation and dissociation occur randomly but within a bounded range, reflecting liquid-like behaviour where clustering dynamics are spatially and temporally uncorrelated. The breadth of the distribution reflects the degree of fluctuation in the clustering, which might be related to density fluctuations. Notably, for larger systems, on the absolute scale, we might expect a distribution with a different mean and data spread without altering the qualitative trend. The absence of collapse into one or two large clusters indicates that the water molecules do not form a percolating or system-spanning hydrogen bond network, but instead a dynamic, fragmented network.

Fig. 3b illustrates the distribution of the number of molecules per cluster, which follows a decay from approximately 5 to 30 molecules per cluster. The peak at 4–5 molecules per cluster implies that the system is dominated by small, short-lived clusters, likely formed via transient local hydrogen bonding. The lack of a long-tailed distribution further supports the absence of large, system-spanning clusters, highlighting the localized and ephemeral nature of water molecule associations in the system.

3.2.2. Calculation of the tetrahedral order parameter (TOP). The tetrahedral order parameter (TOP), also referred to as the orientational order parameter, quantitatively characterizes the local structural arrangement of water molecules by evaluating the extent to which a central water molecule and its four nearest neighbours adopt an ideal tetrahedral geometry. Numerous previous studies have employed the distribution of the TOP for binary and ternary mixtures to gain insights into the structural organization of water.^43,44 The tetrahedral order parameter, q_i, is defined as,^45,46where Ψ_jk represents the angle between the vectors r⃑_ij and r⃑_jk, which connect the central oxygen atom of the ith water molecule to its jth and kth nearest neighbours (NNs), respectively. This angle Ψ_jk effectively captures the angular deviation from the ideal tetrahedral angle. The value of q_i ranges from 0 to 1. A value of q_i = 1 corresponds to a perfect tetrahedral configuration, while q_i = 0 indicates a completely disordered, non-tetrahedral environment. It is important to note that in a crowded medium or near macromolecular surfaces, water molecules may complete their local tetrahedral geometry by forming interactions with non-hydrogen atoms from the surrounding molecules.^47–49 Accordingly, the calculation of the tetrahedral order parameter for a given water molecule can be performed using two distinct approaches: (i) considering only other water molecules as the four nearest neighbours or (ii) including non-hydrogen atoms (N, O, etc.) from other molecules in addition to water. While following the second approach, we choose all electro-negative atoms, namely, OW, NU, OU, OG1, OG2, OG3, OG4, OG5, and OG6, as the partner of a central water oxygen. Calculations using two different approaches are necessary, as a true “tetrahedral order” is only observed when all possible hydrogen-bonded neighbours of the central water oxygen are considered. Restricting the analysis to only water oxygen atoms (OW) introduce deviations from the ideal tetrahedral arrangement.

We provide the population distribution of q_i in Fig. 4a. The black and red curves correspond to calculations of q_i using the first and second approaches, respectively. The red curve exhibits a peak centred around 0.5, whereas the black curve shows a broad and flat distribution, indicative of the absence of a preferred local structure and significant deviation from tetrahedrality. The distribution obtained by considering OW as the neighbouring partners aligns well with the findings reported by Baksi et al.²⁵ A close inspection of the distribution reveals that the probability of occurrence of q_i values greater than 0.5 is merely 7% when only OW atoms are considered as partners for a central OW. However, this probability increases markedly to approximately 47% when all electronegative atoms are included as potential partners, suggesting a substantial recovery of tetrahedrality. At this point, this observation may appear inconsistent with the prominent peak observed in the OW–OW radial distribution function (Fig. 1a). This apparent contradiction arises because the presence of cosolvent molecules induces a more compact packing of water molecules, although restricting their ability to adopt the optimal orientations required for forming a perfect tetrahedral arrangement. A similar effect is reported by Elola et al. in aqueous formamide solutions at various concentrations.⁵⁰ We also compute the distance distribution of the four nearest neighbours of a central OW based on the chosen partnering scheme. The corresponding results are provided in Fig. S4 of the SI.

Fig. 4 Population distribution of (a) the tetrahedral order parameter (q_i) of a central water oxygen atom, and (b) contribution of different electronegative atoms from water, urea, and glucose molecules as the nearest neighbours of a central water oxygen atom.

Fig. 4b illustrates the relative contributions of electronegative atoms from water (represented in black bars), urea (red bars), and glucose (blue bars) toward the calculation of q_i at different neighbouring positions. Among the first nearest neighbours (1st NN), oxygen atoms from water molecules (OW) contribute the most significantly, followed by glucose and urea. However, by the fourth nearest neighbours (4th NN), glucose atoms become the dominant contributors, while water atoms contribute the least. A steady decreasing trend in the contribution of OW atoms is observed as the neighbour shell increases from the 1st to the 4th NN. Specifically, the contribution of water decreases from approximately 50% at the 1st NN to around 20% at the 4th NN. This 30% reduction in water's contribution is largely offset by a corresponding increase in the contribution from electronegative atoms of urea, which rises from about 10% at the 1st NN to nearly 40% at the 4th NN. This trend suggests a progressive replacement of water molecules by urea molecules with increasing distance from a central water oxygen atom (OW). In contrast, the contribution from glucose atoms remains relatively constant at approximately 40% across all neighbour levels. This consistent contribution is likely due to the presence of multiple electronegative sites within a glucose molecule, allowing it to simultaneously influence multiple neighbouring positions. Overall, the qualitative trends in atomic contributions to q_i observed in Fig. 4b are consistent with the variations in peak positions and intensities of the radial distribution functions shown in Fig. 1. It should be noted that urea is typically classified as a water structure–destabilizing osmolyte, whereas glucose, like other polyols, is generally considered a structure–stabilizing osmolyte.⁵¹ Contrary to this conventional view, our observations indicate that both urea and glucose actively participate in determining the local water structure. They can displace other water molecules from the first solvation shell of a central water molecule, thereby functioning as agents that disrupt the water network.

3.2.3. Hydrogen bond analysis. We present a detailed analysis of hydrogen bonding centred on a water molecule, considering all possible donor–acceptor combinations. A hydrogen bond is deemed successful based on established geometric criteria for hydrogen bonding,⁵² which are used to identify valid interactions involving the central water molecule. This methodological approach has been previously used with success in elucidating the complex microscopic structures of binary and ternary mixtures, as reported by Chakrabarti and co-workers.^47,48 The total number of possible donor–acceptor combinations for various central atoms is listed in Table S3 of the SI.

As shown in Fig. 5a, we quantify the number of hydrogen bonds simultaneously formed by a central water molecule acting as a donor, acceptor, or both. This calculation is performed under two conditions: (i) considering only water molecules as hydrogen-bonding partners (black curve) and (ii) including contributions from urea and glucose as well (red curve). This framework is analogous to that used in the computation of the tetrahedral order parameter. Our findings indicate that when only water is considered, the mono-hydrogen-bonded state is the most prevalent. This can be attributed to the inability of multiple surrounding water molecules to simultaneously fulfil the geometric criteria required for hydrogen bond formation with the central water molecule. In contrast, when all potential donor–acceptor combinations (including urea and glucose) are considered, the central water molecule commonly forms two to three hydrogen bonds concurrently. In this scenario, the bi- and tri-hydrogen-bonded states emerge as the most populated, each with an occurrence probability of approximately 35%. These observations are consistent with prior experimental⁵³ and simulation⁴⁸ studies. Formation of an ideal tetrahedral water structure requires four simultaneous hydrogen bonds. However, analysis of the distribution of hydrogen bonding states reveals that mono- and tetra-hydrogen-bonded states are relatively infrequent, each occurring with a probability between 10% and 15%. The low population of the tetra-hydrogen-bonded state contributes to the overall structural distortion of water molecules.

Fig. 5 Population distribution of (a) different multi-hydrogen-bonded water oxygen atoms, (b) contribution of different electronegative atoms from water, urea, and glucose molecules for (I) bi-hydrogen bonded, (II) tri-hydrogen bonded water oxygen atoms.

Furthermore, we have characterized the constituent composition of water molecules in bi- and tri-hydrogen-bonded states. The hydrogen bonding modes are represented as (aW–bU–cG), where a, b, and c denote the number of atoms from water (W), urea (U), and glucose (G), respectively, involved in hydrogen bonding. Fig. 5b(I) presents all six possible modes for bi-hydrogen-bonded water molecules in bar-plot form. The dominant modes involve glucose either as the sole hydrogen-bonding partner or in combination with water. Specifically, the (1W–0U–1G) mode contributes ∼35% and the (0W–0U–2G) mode contributes ∼25% to the bi-hydrogen-bonded population. In contrast, hydrogen bonding modes involving urea are less frequent, each contributing approximately 10% or less. These trends are consistent with the hydrogen bond dynamics data reported by Sarkar et al.²⁶ and the qualitative interaction patterns observed in the radial distribution functions and tetrahedral order parameter calculations (Fig. 1 and 4, respectively). In Fig. 5b(II) we present all ten possible tri-hydrogen-bonded configurations. Among these, the dominant hydrogen bond modes involve atoms from glucose molecules. The most frequent mode is (1W–0U–2G), accounting for 25% of the tri-hydrogen-bonded population, followed by two equally populated modes: (0W–0U–3G) and (1W–1U–1G).

3.3. Structural characteristics of urea molecules

3.3.1. Relative orientation of the molecular plane. To gain deeper insight into the spatial organization of urea molecules, we examine the relative orientation of their molecular planes as a function of intermolecular distance. This analysis is particularly relevant given that urea molecules are known to exhibit stacking tendencies, especially with aromatic side chains in protein structures.⁵⁴ Furthermore, guanidinium ions, which are structurally related to urea, are reported to display parallel stacking in aqueous environments.⁵⁵ In this study, we construct a two-dimensional colour map based on two parameters: the distance ‘r’ between the carbon atoms (CU–CU) of two randomly selected urea molecules, and the angle ‘θ’ between their respective molecular planes. A schematic depiction of ‘r’ and ‘θ’ is provided in Fig. 6a.

Fig. 6 (a) Schematic representation of two urea molecules inclined at an angle ‘θ’ and separated by a distance ‘r’. The unit vector perpendicular to the urea plane is denoted by û. (b) Colour map constructed using (r, θ). Yellow and violet colours correspond to probabilities of 0 and 1, respectively. For convenience, the g_CU–CU(r) curve is also presented in green.

As Fig. 1b is evident, the radial distribution function for the CU–CU pair exhibits a pronounced peak at approximately 0.44 nm, with the first coordination shell extending up to about 0.6 nm. This peak indicates a significant accumulation of urea molecules within this spatial range. The corresponding colour map, presented in Fig. 6b, reveals that the most probable configuration within 0.6 nm involves an interplanar angle of approximately 90° between the molecular planes of interacting urea pairs. This pronounced orientational preference, centred around 0.44 nm, suggests that urea molecules favour a perpendicular arrangement when in closest possible proximity. This orthogonal orientation can be attributed to steric constraints arising from molecular crowding, which prevent the molecules from adopting a parallel configuration. Instead, the ∼90° arrangement likely contributes to the overall structural stability of the system. Moreover, this perpendicular alignment facilitates dense and symmetrical short-range packing of urea molecules, as illustrated in Fig. 2b(II). Notably, the absence of significant populations at 0° or 180° indicates a lack of parallel stacking among urea molecular planes. In systems characterized by complex hydrogen bonding networks, like the one we are investigating here, such relative orientations can substantially influence urea–urea hydrogen bonding interactions. A detailed analysis of these interactions is presented in the subsequent section.

3.3.2. Hydrogen bond analysis. We present a comprehensive analysis of hydrogen bonding centred on the NU and OU atoms of urea molecules, incorporating all possible donor–acceptor pairings. As presented in Fig. 7a, we quantify the number of hydrogen bonds simultaneously formed by a urea atom functioning as a donor, an acceptor, or both. The results reveal that NU (black curve) and OU (red curve) exhibit similar distributions across different hydrogen-bonding states, despite that OU can only act as an acceptor. This similarity likely arises from a preferential urea–urea interaction, wherein NU and OU complement each other in terms of hydrogen bond donation and acceptance. We find that the mono- and bi-hydrogen-bonded configurations are dominating, with approximate population probabilities of 40% and 35%, respectively.

Fig. 7 Population distribution of (a) different multi-hydrogen bonded nitrogen and oxygen atoms of urea, (b) contribution of different electronegative atoms from water, urea, and glucose molecules for (I) mono-hydrogen bonded, (II) bi-hydrogen bonded nitrogen and oxygen atoms of urea.

To further elucidate these findings, we have characterized the compositional environment surrounding a central urea atom in both mono- and bi-hydrogen-bonded states. Fig. 7b(I) presents the three possible mono-hydrogen-bonded configurations in bar plot format. For both NU and OU atoms, urea is the dominant interacting species (40–50%), followed by glucose, with water contributing the least. A closer examination of Fig. 7b(I) suggests that OU atoms serve as more favourable interaction sites than NU atoms in urea–urea interactions, a conclusion also supported by the peak intensities shown in Fig. 1b. In Fig. 7b(II), we depict all six bi-hydrogen-bonded configurations. Among these, urea-involved hydrogen bonding remains dominant, followed by interactions involving glucose molecules. Specifically, for NU atoms, the (0W–1U–1G) configuration is most prevalent, accounting for approximately 30% of the bi-hydrogen-bonded population. In contrast, for OU atoms, the (0W–2U–0G) mode dominates with an occurrence of ∼40%. The second most populated modes are exactly inversely distributed: (0W–2U–0G) for NU (∼20%) and (0W–1U–1G) for OU (∼25%). A plausible explanation for the enhanced urea–urea interaction is their relatively perpendicular molecular orientation, which may facilitate optimal hydrogen-bonding configurations between the hydrogen donor groups (–NH₂) and the carbonyl acceptor sites (C=O) of neighbouring molecules.

3.4. Structural characteristics of glucose molecules

Fig. 8 presents the distribution of the number of hydrogen bonds formed simultaneously by a central glucose molecule. While previous analyses focused on the atomic-level characteristics of hydrogen bonding, this molecular-level perspective is crucial due to the relatively large size and complex structure of glucose. Specifically, glucose possesses multiple hydrogen bond donor and acceptor sites dispersed over a broad molecular region. As a result, the hydrogen bonding network around a central glucose molecule is more spatially extended compared to smaller molecules such as water or urea.

Fig. 8 Population distribution of (a) different multi-hydrogen bonded glucose molecules and (b) contribution to form hydrogen bonds by water, urea, and glucose molecules to a central glucose molecule.

As shown in Fig. 8a, the distribution of the number of hydrogen bonds exhibits a Gaussian-like profile. This indicates that the hydrogen bonding behaviour of glucose arises from numerous small, random interactions fluctuating around a most probable value, with no inherent bias toward forming either very few or an excessive number of hydrogen bonds. These fluctuations likely reflect dynamic changes in the local molecular environment and configurations over time. At any given moment, a central glucose molecule can form as few as one hydrogen bond and up to approximately thirteen simultaneous hydrogen bonds with neighbouring molecules. This broad range highlights the conformational flexibility of glucose. The most frequently observed state corresponds to the formation of 7 to 8 hydrogen bonds, each occurring with an approximate probability of 20%. This is consistent with the molecular structure of glucose, which contains five hydroxyl groups and a ring oxygen atom, each capable of participating in hydrogen bonding.

Fig. 8b displays, in bar chart format, the probability of hydrogen bonding between a central glucose molecule and its surrounding molecules. Each bar represents a specific combination of hydrogen bonding partners, denoted by a three-letter code where each letter (W: water, U: urea, G: glucose) is prefixed with either ‘y’ (Yes, indicating bonding) or ‘n’ (No, indicating no bonding). The most probable scenario, with a 72% occurrence, involves simultaneous hydrogen bonding of the central glucose molecule with all three types of neighbours (i.e., water, urea, and glucose). This is followed by cases where the central glucose forms hydrogen bonds with other glucose molecules and either water or urea, each contributing approximately 10% to the overall probability distribution.

We next present a detailed analysis of hydrogen bonding centred on the oxygen atoms of glucose molecules, accounting for all possible donor–acceptor pairings. As illustrated in Fig. 9a, we quantify the number of hydrogen bonds simultaneously formed by individual oxygen atoms of glucose acting as donors, acceptors, or both. The analysis reveals that the exocyclic oxygen atoms OG1, OG3, OG4, and OG5 display comparable distributions across different hydrogen-bonding states. In particular, mono- and bi-hydrogen-bonded configurations dominate, each exhibiting population probabilities of approximately 40%. This similarity likely arises from the shared exocyclic nature of these atoms and their positioning within similarly congested molecular environments. Also, these four oxygen atoms exhibit an ∼10% probability of existing in a non-hydrogen-bonded state. In contrast, although OG2 is also an exocyclic oxygen atom, its proximity to the ring oxygen OG6 restricts its accessibility, thereby hindering the formation of two simultaneous hydrogen bonds. Consequently, OG2 exhibits the highest probability (∼50%) of being mono-hydrogen bonded, with both bi-hydrogen bonded and non-hydrogen bonded configurations occurring with approximately 20% probability each. The ring oxygen OG6, lacking hydrogen bond donor capability and being embedded within the cyclic framework, is subject to pronounced steric hindrance. As a result, OG6 is the least favourable site for hydrogen bonding, with the non-hydrogen-bonded state being the most prevalent, followed by the mono-hydrogen-bonded configuration, which occurs with a probability of ∼40%.

Fig. 9 Population distribution of (a) different multi-hydrogen-bonded glucose oxygen atoms and (b) contribution of different electronegative atoms from water, urea, and glucose molecules for (I) mono-hydrogen bonded and (II) bi-hydrogen bonded glucose oxygen atoms.

To further contextualize these observations, we characterize the local chemical environment surrounding each glucose oxygen atom in both mono- and bi-hydrogen-bonded states. Fig. 9b(I) presents the distribution of species involved in the three possible mono-hydrogen-bonded configurations using bar plots. For all exocyclic oxygen atoms, water and glucose molecules dominate the interactions, each contributing approximately 30–40%, while urea appears to be the least preferred hydrogen bond partner. In contrast, the ring oxygen OG6 exhibits a different pattern, with urea accounting for the majority of interactions (∼40%), followed by water and glucose. Fig. 9b(II) depicts all six possible bi-hydrogen-bonded configurations. The qualitative trends observed in the mono-hydrogen-bonded cases persist here as well. For the exocyclic oxygen atoms, simultaneous hydrogen bonding by water and glucose molecules is the most probable configuration (∼25%), followed by the (0W–1U–1G) configuration, which has an ∼20% probability. However, OG6 demonstrates a distinct bonding pattern, with the (0W–2U–0G) configuration dominating (∼40%) and the (1W–1U–0G) configuration following at ∼20%. This deviation is attributed to the preferential interaction of exocyclic oxygen atoms with glucose and water molecules, whereas the ring oxygen OG6 is left behind to predominantly interact with urea.

These trends in both mono- and bi-hydrogen-bonded states are further corroborated by radial distribution function (RDF) analysis, as shown in Fig. 1. The oxygen atom of water (OW) exhibits over threefold higher spatial density around OG1–OG5 compared to OG6. In contrast, the spatial distribution of urea nitrogen (NU) atoms around all glucose oxygen sites is relatively uniform, supporting the conclusion that OG6 is largely excluded from preferred hydrogen bonding interactions with water and glucose, instead predominantly interacting with urea.

4. Conclusions

We utilize classical all-atom molecular dynamics simulations to investigate, at a microscopic level, the structural organization of all constituent molecules in a naturally abundant deep eutectic solvent (NADES) comprising water, urea, and glucose. Detailed structural characterization that includes pair-correlation functions of atomic species, clustering behaviour of water molecules, relative orientations of urea molecular planes, and statistical occurrences of multi-hydrogen-bonded electronegative atoms is reported.

Radial distribution function analysis reveals a pronounced population of water oxygen atoms surrounding glucose oxygen atoms, indicative of favourable water–glucose interactions. In contrast, urea–urea interactions dominate over urea–water and urea–glucose pairwise combinations, with urea oxygen atoms (OU) serving as preferential hydrogen bond sites. These interaction trends are corroborated by hydrogen bonding analysis, which shows a higher frequency of hydrogen bonds between various atom pairs originating from water and glucose, while urea molecules tend to form hydrogen bonds predominantly with other urea molecules. This observation is likely a consequence of the close packing behaviour of urea, as supported by triangular spatial occupancy analysis, where urea–urea–urea (U–U–U) triads exhibit equilateral configurations at short intermolecular distances, signifying symmetric and tightly packed molecular arrangements. A direct consequence of this preferential hydrogen-bonded structure is expected to manifest in the dynamical behaviour of the system's constituents, consistent with observations in other mixtures.⁵⁶ A comparison between the hydrogen-bond population reported here and the relaxation timescales reported by Sarkar et al.²⁶ reveals that urea–water hydrogen-bond relaxation occurs significantly faster than urea–urea or water–water relaxation. This observation further substantiates that the number of hydrogen bonds formed exerts a pronounced influence on the overall hydrogen-bond relaxation dynamics. Further analysis of the relative orientations of urea molecular planes demonstrates a preferred orthogonal arrangement, which may facilitate hydrogen bonding between urea molecules. Water clustering behaviour reveals a preferential formation of 9–10 water clusters in the system, each consisting typically of 4–5 water molecules. Examination of the water's local tetrahedral structure indicates a substantial deviation from ideal tetrahedrality. This deviation is attributed to the low prevalence (less than 20%) of tetra-hydrogen bonded water configurations and the dominance of bi- and tri-hydrogen bonded species.

Over recent decades, both experimental and computational investigations have aimed to elucidate the microscopic structure and dynamics of complex fluids such as ionic liquids and deep eutectic solvents.^{24,26,48,57–59} One of the notable industrial applications of deep eutectic solvents is in selective gas absorption and gas separation technologies.^60,61 While extensive experimental data exist on this subject, establishing direct correlations between gas absorption selectivity and molecular-level interactions remains challenging. Atomistic molecular dynamics simulations offer a powerful approach to bridge this gap. Identifying correlations between descriptors such as the free energy of solvation and enthalpy of absorption with gas uptake selectivity can yield mechanistic insights and inform the rational design of DESs for specific gas separation tasks. The structural features reported here, particularly hydrogen bonding patterns and water clustering phenomena, may contribute to a molecular-level understanding of gas solubilization in NADESs. Moreover, the simulation methodology and protocols employed in this study are readily extendable to other solvent systems and temperature regimes relevant to industrial applications.

Author contributions

Soham Sarkar: conceptualization, data curation, formal analysis, investigation, methodology, validation, visualization, writing – original draft, writing – review and editing.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

All the data and the codes required to reproduce the results are freely available in GitHub (https://github.com/soham9038/NADES_structure).

Supplementary information: Schematic representation of the structure of water, urea, and glucose; comparison of the simulated densities with experimental densities at different temperatures; radial distribution function calculated between different atoms of glucose and glucose–urea molecules; population distribution of the distance between the central water oxygen atom and the nearest neighbour (NN) atoms; the number of possible donor–acceptor combinations for different central electronegative atoms. See DOI: https://doi.org/10.1039/d5nj02678d.

Acknowledgements

S. S. acknowledges Dr Atanu Maity for technical help.

References

1. D. Majumdar (neé Som), C. Dutta, A. K. Mukherjee and S. Sen, Transp. Res. D: Transp. Environ., 2008, 13, 524–530.
2. B. W. Lee, K. T. Kelsey, D. Hashimoto, B. Yakes, T. Seitz and D. C. Christiani, J. Occup. Environ. Med., 1997, 39, 960–969.
3. J. E. Lewis, R. Biswas, A. G. Robinson and M. Maroncelli, J. Phys. Chem. B, 2001, 105, 3306–3318.
4. A. M. Scurto, K. W. Hutchenson and B. Subramaniam, Gas-Expanded Liquids and Near-Critical Media, American Chemical Society, 2009, vol. 1006, pp. 3–37.
5. M. J. Earle and K. R. Seddon, Pure Appl. Chem., 2000, 72, 1391–1398.
6. E. L. Smith, A. P. Abbott and K. S. Ryder, Chem. Rev., 2014, 114, 11060–11082.
7. Q. Zhang, K. D. O. Vigier, S. Royer and F. Jérôme, Chem. Soc. Rev., 2012, 41, 7108–7146.
8. M. Espino, M. de los Ángeles Fernández, F. J. V. Gomez and M. F. Silva, Trends Anal. Chem., 2016, 76, 126–136.
9. D. V. Wagle, H. Zhao and G. A. Baker, Acc. Chem. Res., 2014, 47, 2299–2308.
10. M. Kondo, X. Li and M. Maroncelli, J. Phys. Chem. B, 2013, 117, 12224–12233.
11. A. A. Neufeld, A. I. Burshtein, G. Angulo and G. Grampp, J. Chem. Phys., 2002, 116, 2472–2479.
12. W.-C. Huang, D. Zhao, N. Guo, C. Xue and X. Mao, J. Agric. Food Chem., 2018, 66, 11897–11901.
13. Y. H. Choi, J. van Spronsen, Y. Dai, M. Verberne, F. Hollmann, I. W. C. E. Arends, G.-J. Witkamp and R. Verpoorte, Plant Physiol., 2011, 156, 1701–1705.
14. A. P. Abbott, D. Boothby, G. Capper, D. L. Davies and R. K. Rasheed, J. Am. Chem. Soc., 2004, 126, 9142–9147.
15. A. Paiva, R. Craveiro, I. Aroso, M. Martins, R. L. Reis and A. R. C. Duarte, ACS Sustainable Chem. Eng., 2014, 2, 1063–1071.
16. I. Juneidi, M. Hayyan and M. A. Hashim, Process Biochem., 2018, 66, 33–60.
17. I. M. Aroso, R. Craveiro, Â. Rocha, M. Dionísio, S. Barreiros, R. L. Reis, A. Paiva and A. R. C. Duarte, Int. J. Pharm., 2015, 492, 73–79.
18. A. Hayyan, F. S. Mjalli, I. M. AlNashef, Y. M. Al-Wahaibi, T. Al-Wahaibi and M. A. Hashim, J. Mol. Liq., 2013, 178, 137–141.
19. Y. Dai, J. van Spronsen, G.-J. Witkamp, R. Verpoorte and Y. H. Choi, Anal. Chim. Acta, 2013, 766, 61–68.
20. M. Faggian, S. Sut, B. Perissutti, V. Baldan, I. Grabnar and S. Dall’Acqua, Molecules, 2016, 21, 1531.
21. L. Duan, L.-L. Dou, L. Guo, P. Li and E.-H. Liu, ACS Sustainable Chem. Eng., 2016, 4, 2405–2411.
22. S. Barani Pour, J. Jahanbin Sardroodi and A. Rastkar Ebrahimzadeh, J. Mol. Liq., 2021, 334, 115956.
23. Z. Sailau, N. Almas, A. Aldongarov and K. Toshtay, J. Mol. Model., 2022, 28, 235.
24. E. Tarif, J. Mondal and R. Biswas, J. Phys. Chem. B, 2019, 123, 9378–9387.
25. A. Baksi, J. Rajbangshi and R. Biswas, Phys. Chem. Chem. Phys., 2021, 23, 12191–12203.
26. S. Sarkar and F. Müller-Plathe, J. Chem. Phys., 2024, 161, 224505.
27. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
28. A. Lindahl, B. Hess and D. Van Der Spoel, GROMACS 2021. 1 Manual 2021.
29. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
30. W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 1988, 110, 1657–1666.
31. P. Mark and L. Nilsson, J. Phys. Chem. A, 2001, 105, 9954–9960.
32. L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, Nucleic Acids Res., 2017, 45, W331–W336.
33. S. S. Petrova and A. D. Solov’ev, Hist. Math., 1997, 24, 361–375.
34. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
35. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
36. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 1989.
37. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
38. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
39. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
40. D. Dhabal, A. H. Nguyen, M. Singh, P. Khatua, V. Molinero, S. Bandyopadhyay and C. Chakravarty, J. Chem. Phys., 2015, 143, 164512.
41. M. Singh, D. Dhabal, A. H. Nguyen, V. Molinero and C. Chakravarty, Phys. Rev. Lett., 2014, 112, 147801.
42. D. Dhabal, A. Gupta and H. K. Kashyap, J. Chem. Phys., 2017, 146, 094503.
43. S. Daschakraborty, J. Chem. Phys., 2018, 148, 134501.
44. D. Bandyopadhyay, S. Mohan, S. K. Ghosh and N. Choudhury, J. Phys. Chem. B, 2014, 118, 11757–11768.
45. J. R. Errington and P. G. Debenedetti, Nature, 2001, 409, 318–321.
46. P.-L. Chau and A. J. Hardwick, Mol. Phys., 1998, 93, 511–518.
47. S. Sarkar, A. Maity and R. Chakrabarti, J. Phys. Chem. B, 2021, 125, 13212–13228.
48. S. Sarkar, A. Maity and R. Chakrabarti, Phys. Chem. Chem. Phys., 2021, 23, 3779–3793.
49. L. Aggarwal and P. Biswas, J. Phys. Chem. B, 2018, 122, 4206–4218.
50. M. D. Elola and B. M. Ladanyi, J. Chem. Phys., 2006, 125, 184506.
51. A. Panuszko, P. Bruździak, E. Kaczkowska and J. Stangret, J. Phys. Chem. B, 2016, 120, 11159–11169.
52. A. Luzar and D. Chandler, Nature, 1996, 379, 55–57.
53. Z. Lu, M. Lanagan, E. Manias and D. D. Macdonald, J. Phys. Chem. B, 2009, 113, 13551–13559.
54. S. Goyal, A. Chattopadhyay, K. Kasavajhala and U. D. Priyakumar, J. Am. Chem. Soc., 2017, 139, 14931–14946.
55. D. Bandyopadhyay, K. Bhanja, S. Mohan, S. K. Ghosh and N. Choudhury, J. Phys. Chem. B, 2015, 119, 11262–11274.
56. S. Kumar, S. Sarkar and B. Bagchi, J. Chem. Phys., 2019, 151, 194505.
57. D. Maji and R. Biswas, J. Mol. Liq., 2025, 423, 127042.
58. S. Chatterjee, T. Chowdhury and S. Bagchi, J. Phys. Chem. B, 2024, 128, 12669–12684.
59. S. M. Jha, K. K. Athira, R. L. Gardas, N. Malek and S. P. Ijardar, J. Chem. Eng. Data, 2025, 70, 2631–2644.
60. S. Sarmad, Y. Xie, J.-P. Mikkola and X. Ji, New J. Chem., 2016, 41, 290–301.
61. R. B. Leron and M.-H. Li, Thermochim. Acta, 2013, 551, 14–19.

---