Structural insights into carboxylic-acid based DES across H-bond donor ratios: impact of CL&Pol refinement

Abstract

Deep eutectic solvents (DES) are of significant interest due to their eco-friendly nature and vast applications. Carboxylic-acid-based choline chloride (ChCl) DES are notable for their roles in electrochemical, drug delivery, and biomass processing applications, with efficiency influenced by the ChCl : carboxylic acid ratio. Understanding these mechanisms requires detailed knowledge of their structure. This study investigates the choline chloride–lactic acid (ChCl:LA) DES structure using ab initio molecular dynamics simulations to assess the accuracy of the transferable and polarizable CL&Pol force field. We observe that the CL&Pol force field qualitatively captures primary interactions within the system, despite numerical discrepancies due to its transferable nature. To refine the original force field, we incorporate two improvements: tuning the σ parameter of the strongest hydrogen-bond interactions and incorporating the Tang–Toennies damping function to correct chloride ion overpolarization. The first adjustment enhances the targeted interactions and significantly improves the short-range structure of the entire hydrogen-bond network. The second refinement, although minimally impacting the structure at low LA ratios, proves critical at higher ratios by correcting the oversegregation of ionic molecules in the original force field. Consequently, it becomes essential for reliably depicting the medium and long-range structure of the system, highlighting that the specific parameter of the force field to be refined depends on the structural scale under investigation. Notably, the long-range structure results from the competition between choline and carboxylic acid for chloride, rebalanced by the suggested modifications, especially the overpolarization correction.

---

1 Introduction

Deep eutectic solvents (DES) are intricate mixtures typically composed of a hydrogen bond donor (HBD) and a hydrogen bond acceptor (HBA), exhibiting a significant reduction in melting temperature at specific compositions compared to their individual constituents. DES are recognized for their ease of synthesis and environmentally friendly origins.¹

Among these solvents, choline chloride-based DES incorporating carboxylic acid-based HBDs have gained significant attention for several reasons: they are naturally abundant, highly cost-effective due to the plentiful availability of raw materials, and possess environmentally friendly properties. These DES are currently used in advanced applications such as therapeutic applications,^2,3 electrodeposition processes,⁴ spent battery recycling,⁵ and are considered promising alternatives for electrolytes in batteries and fuel cells.⁶

The literature has revealed that ChCl–carboxylic acid DES, display intricate structural microheterogeneities influenced by composition,⁷ which are closely tied to their potential applications.^8–13 Particularly, the microheterogeneous structure of ChCl-based DES materials has been shown to vary with factors such as the cation's alkyl chain length, the HBD ratio, and water content. Realizing the full potential of these materials hinges on elucidating the intricate relationship between their properties and atomic-level structure, as emphasized by Stephens et al.¹⁴

In this context, computational simulations have shown to be an ideal tool providing a cost-effective alternative to experimental methods but also offering detailed, atomic-level insights into the structure of the system. However, DES, along with ILs, present distinctive challenges for theoretical chemistry, stemming from their complex nature. On one hand, the dynamic study of large model systems is imperative due to the inherent fluctuations in DES structures and the presence of strong, long-range Coulomb interactions. These features often require the use of force-field-based calculations, which allow large size of model systems to be studied. On the other hand, the role of dispersion and induction interactions in these solvents is paramount, making non-polarizable force fields less reliable for addressing structural properties. Finally, the large variety of possible combinations hinders the availability of ad hoc parameterized force-fields.

Recent advancements, such as the work by Goloviznina et al.,¹⁵ have introduced a polarizable and transferable force field (CL&Pol) that seeks to improve the treatment of induction interactions by explicitly incorporating polarization effects. While this innovative approach has significantly advanced our ability to study these systems theoretically, it has also revealed certain limitations, particularly when applied to chlorine-containing DES solvents.¹⁶

One notable limitation pertains to what is commonly referred to as “naked hydrogens” within the force field. Specifically, this limitation arises from the absence of Lennard-Jones parameters for the hydrogen atoms of hydroxyl groups, which leads to the so-called polarization catastrophe. In the case of a chloride containing solvent, excessively short distances and strong interactions between the naked hydrogens and chlorides are obtained, and consequently, an overestimation of the global density of the liquid. This issue has been solved by adding a Tang–Toennies dumping function into all the induced dipole–charge interactions, and by increasing the σ parameter of the Lennard-Jones pair potential between the oxygen atom of the hydroxyl group and the chloride anion to 370 pm. Nevertheless, recent findings in the literature¹⁶ have raised concerns about the adequacy of this fixed σ value. It has been reported that this value may not be suitable, as it tends to underestimate choline-chloride and chloride–neutral molecule interactions in ChCl:ethylene glycol (EG) systems.

Furthermore, another significant challenge arises from the overpolarizability of chlorine atoms in this force field, as highlighted in studies by Szabadi et al.¹⁷ and Maglia de Souza et al.¹⁶ This overpolarizability results in excessively large nanoheterogeneities when applied to ChCl:EG systems, impacting the accuracy of simulations. Different solutions have been suggested in the literature to solve this issue, as reducing chlorine polarizability¹⁷ or applying the Tang–Toennies damping factor to chlorine,¹⁶ the latter of which has provided satisfactory results.

These issues have primarily been analyzed in DES with alcohol group-containing HBDs; however, their effects have not been explored in DES with carboxylic acid-based HBDs, despite their significant interest. Addressing the influence of these limitations will pave the road towards a more reliable representation of carboxylic group-based DES. In order to do so, in this work, we have chosen to study choline chloride with lactic acid (LA) DES, the latter being a HBD that displays an alcohol group and a carboxylic group. This will permit us to study the extent of the limitations of the force field previously described for alcoholic HBDs (1) when combined with other functional groups, and (2), in carboxylic HBDs. We will also address the effects of the mentioned issues in the structure of this DES as the presence of the HBD varies.

Beyond the interest in ChCl:LA DES due to the distinctive chemical groups, this DES is highly relevant due its wide range of applications during recent years: its notable conductivity properties have made it useful as an electrolyte,^18–20 while it has also been explored for drug delivery,²¹ the extraction of organic compounds,²² and non-linear optical phenomena.²³ Furthermore, it has consistently shown exceptional performance in lignin extraction from biomass.^24–26 Additionally, in certain applications such as non-linear optical properties, lignin extraction efficiency, and conductivity,²⁷ adjusting the ChCl : LA ratio has been found to influence the overall performance of the solvent. Despite the broad applicability of the system, no experimental studies of its structure have been published in the literature for the best of our knowledge, and most studied parameters as conductivity, viscosity, diffusion, and density, have been measured only for binary mixtures of water and ChCl:LA DES.^20,27

In summary, in this study, our primary objective is to elucidate the structure of ChCl:LA DES and investigate the effects of increasing the HBD ratio on the system. For that, we will initially address and improve the limitations found in the CL&Pol force field for ChCl-based materials, and assess its performance at the short-range as compared to DFT-based molecular dynamics simulations (AIMD). While AIMD calculations provide accurate short-range structural information, they are constrained to small systems and short simulation times. Therefore, we will initially focus on the 1 : 1 ratio. Subsequently, our investigation will be extended to a range of ChCl : LA ratios by means of classical molecular dynamics (CMD) simulations. Our overarching goal is to provide insights into the significant variations in experimental performance with increasing LA excess.

2 Computational details

CMD simulations of periodic cubic systems are carried out using the LAMMPS package.²⁸ Interparticle interactions are modeled using the extended version of the CL&Pol polarizable force field by Goloviznina et al.^15,29 based on the OPLS-AA framework. Within this force field the drude oscillator model, implemented in the USER-DRUDE package,³⁰ is employed to describe all atoms of the system excluding hydrogens. Force field parameters for choline and chloride molecules are taken from the fftool package,³¹ while parameters for LA molecules are obtained from the OPLS-AA based LigParGen web server.³² Subsequently, they are polarized analogously to the work by Goloviznina et al.,²⁹ treating the LA molecule as a single fragment. The predicted dispersion scaling factors between different type of fragments are shown in Table 1. The induced short-range electrostatic charge–dipole and dipole–dipole interactions resulting from the inclusion of polarization are dumped by the Tang–Toennies and Thole functions, respectively, as detailed in ref. 15 and 29. Electrostatic interactions are computed by the particle–particle particle-mesh scheme with an accuracy of 1 × 10⁻⁵. A cutoff of 1200 pm was considered for the Lennard-Jones potential including tail corrections.³³

Table 1 Dispersion scaling factors obtained for each pair of types of fragment considered in this work

Fragment i	Fragment j	k ij
Ch^+	Cl^−	0.48
Ch^+	LA	0.47
Cl^−	LA	0.51
LA	LA	0.66

---

The CMD simulations are performed within the NPT ensemble using the Nose–Hoover thermostat, as detailed by Dequidt et al.³⁰ The temperature is fixed at 333 K, pressure at 1 bar, and timestep to 1 fs. Initial random molecular configurations are generated using the Packmol package.³⁴ Input files are constructed with the fftool package, following instructions detailed by Goloviznina et al.³⁵ Unless otherwise stated, the equilibration period for the simulations are 112 ns, followed by a 30 ns production run. The amount of molecules included in the simulation box varies depending on the ChCl : LA DES ratio: 300 ChCl and 300 LA molecules for the 1 : 1 ratio, 200 ChCl and 400 LA molecules for the 1 : 2 ratio, 100 ChCl and 500 LA molecules for the 1 : 5 ratio, and 55 ChCl and 550 LA molecules for the 1 : 10 ratio. This choice of number of molecules results in similar averaged equilibrated cubic box dimensions for the various employed ratios and force-fields, ranging from 4482 pm to 4665 pm.

The structural analysis (radial distribution functions, combined distribution functions and Voronoi analysis) shown in this work is carried out on the production run with the TRAVIS software^36,37 while direct images of the simulated molecules are captured with the VMD³⁸ software.

AIMD trajectories are produced at 333 K by means of the CP2K software,^39–43 starting from a cubic box with dimensions of 1752 pm × 1752 pm × 1752 pm that contains 16 ChCl ion pairs and 16 LA molecules produced by CMD simulations. NVT ensemble calculations are performed using the Nose–Hoover thermostat.^44,45 A production of 109 ps is obtained, with a time step of 0.5 fs, at the BLYP-D3/DZVP-MOLOPT-SR-GTH^46–53 level. Due to the higher computational cost associated with AIMD calculations, several CMD simulations are carried out containing 16 molecules of each molecular compound. This allows for a direct comparison of results between the two levels of theory. In this case, the CMD simulations involve 2 ns of equilibration time followed by 10 ns of production time. The size and time convergence aspects of these calculations are further discussed in the ESI† (Section S1).

3 Results

3.1 Force field refinement: the short-range structure

Our first task is to address the limitations found in the force field, based on the corrections previously suggested in the literature and taking AIMD calculations as reference. The strongest interactions among those which form the hydrogen bond network established in our system are the charge–dipole interactions, thus, the hydrogen bonds between the chlorine anion and the hydroxyl groups. In particular, in our system there are three different hydroxyl groups (see Fig. 1): the hydroxil in the carboxylic group of the LA, the alcohol group in LA, and the hydroxil group in choline, referred to as LACOOH, LAOH and ChOH respectively hereafter.

Fig. 1 Chemical structure of the choline chloride–lactic acid DES.

Fig. 2 displays the RDFs obtained for Cl⁻–ChOH, Cl⁻–LACOOH, and Cl⁻–LAOH interactions, calculated using various levels of theory in this study (details are provided later in the text). These simulations were conducted with an equal number of molecules (16 of each molecular compound) in the simulation box to prevent any size-related artifacts that could affect the fitting process. The first-neighbour maxima of the RDFs obtained with AIMD calculations are located at 212 pm, 198 pm and 215 pm, respectively. The first two interactions exhibit significantly larger values for the first peak maxima: 10.7 and 11.5, respectively, whereas for the Cl⁻–LAOH RDF is 5.1.

Fig. 2 Radial distribution functions (RDF) between (a) Cl⁻–ChOH (b) Cl⁻–LACOOH and (c) Cl⁻–LAOH atoms for AIMD (black), FF1 (blue), FF2 (red) and FF3 (green) in a box of 16 molecules of each molecular compound.

CMD RDFs obtained with the extended version of the CL&Pol polarizable force field by Goloviznina et al.^15,29 (named FF1 hereafter) agree in the hydrogen bond distances (blue in Fig. 2), but with discrepancies in the g(r) value, as previously observed by Maglia de Souza et al. for the hydrogen bonds formed in ChCl:EG DES.¹⁶ In their case, the RDFs of the hydrogen bonds formed by chlorides employing the CL&Pol force field show lower first-neighbor peaks than the ones obtained with AIMD simulations. In our case the FF1 shows distinct deviations in the RDFs with respect to AIMD calculations depending on the molecule type with which chloride is forming a hydrogen bond: an overestimation with choline (Fig. 2(a)) and an underestimation with LA (Fig. 2(b and c)). This is specially striking for Cl⁻–LAOH RDF, where at the first peak, the g(r) value reaches 0.58, considerably lower than what a non-correlated pair of particles would exhibit. While keeping the rest of the interparticle σ parameters from the FF1 force field unchanged, we adjust the specifical σ parameters related to the Cl⁻–ChOH, Cl⁻–LACOOH, and Cl⁻–LAOH interactions, using the RDFs obtained from AIMD as a reference. The modified force field that includes these tuned parameters is named FF2 hereafter. The original and modified values of the σ parameters are shown in Table 2.¹⁵ The agreement between AIMD and FF2 (see RDFs in red in Fig. 2) is very satisfactory.

Table 2 Original and modified values of interparticle Lennard-Jones σ parameter for each interparticle interaction studied in this work

	σ FF1 (pm)	σ FF2 (pm)	σ FF3 (pm)
Cl^−–ChOH	370	380	376
Cl^−–LACOOH	370	366	364
Cl^−–LAOH	370	347	338

---

Regarding the overpolarization of the chlorine, its influence was monitorized in the literature by analyzing the Cl⁻–Cl⁻ RDF,^16,17 reporting the underestimation of Cl⁻–Cl⁻ distances as compared to AIMD values. Our AIMD simulations show quite a wide peak at around 600 pm for the Cl⁻–Cl⁻ RDF (see Fig. S4, ESI†), closer than related systems, as the ChCl:EG system by Alizadeh et al.,⁵⁴ which predict the maximum at around 770 pm. FF1 and FF2, on the contrary, predict this maximum at 500 pm approximately, similar to ChCl:EG CMD simulations.¹⁶ Nevertheless, FF-based RDFs in our system predict a higher and narrower peak, suggesting that the H-bond network with Cl⁻ is stronger in this case. The overpolarization of the force field was corrected in the literature, by introducing the Tang–Toennies damping function to Cl⁻,¹⁷ and in fact was seen to shift the maximum of the Cl⁻–Cl⁻ RDF to around 700 pm for ChCl:EG,¹⁶ reaching a better agreement with the corresponding AIMD data. This was done by applying the Tang–Toennies function to the previously σ-corrected force field. In our case, following an approach similar to that for FF2, we obtain a new set of σ values (see Table 2) by fitting the first peaks of the RDFs from AIMD calculations (shown in green in Fig. 2) including the Tang–Toennies dumping function between the induced dipoles that contain Cl⁻ ions. This procedure is analogous to that of Maglia de Souza et al., where we also set the same value for the parameter b = 4.5 in the Tang–Toennies potential.^15,16 This modified version of the force field is referred to as FF3 hereafter. Similar to the FF2 version, in FF3 only the σ parameters for the interactions listed in Table 2 are adjusted, while the σ parameters for all other interparticle interactions from FF1 remain unchanged.

In our case, although the value of the first peak in the Cl⁻–Cl⁻ RDF with the Tang–Toennies damping function approaches that obtained by AIMD, it holds its position at around 500 pm (see S2, ESI†). Consequently, the FF3 force field shows no improvement in the positioning of the first neighbor peak, contrarily to what observed by Maglia de Souza et al. for ChCl:EG.¹⁶ This suggest that our structure relies mostly on Cl⁻–ChOH and Cl⁻–LACOOH interactions, which makes the overall structure less dependant on the polarization effects than other, non carboxylic acid-based systems.

Although the RDFs obtained with the modified force fields FF2 and FF3 show satisfactory agreement with those from AIMD, some structural differences remain. Specifically, the second RDF peak in Fig. 2(a) present in all three force fields, disappears in the AIMD results. This difference arises from the tendency of chloride anions to occupy a secondary position closer to the ammonium group of choline in all CMD simulations (see Section S2 of the ESI†).

After establishing the three distinct force fields, it would be appropriate to evaluate their performance against available experimental data. Unfortunately, experimental data for ChCl:LA DES under comparable temperature and pressure conditions are scarce in the literature. To the best of our knowledge, the only reported data so far is the density for a 1 : 1 ChCl : LA ratio.⁵⁵Table 3 presents the experimentally obtained density value alongside those obtained from our CMD simulations employing the FF1, FF2, and FF3 force fields. The density predicted by the FF1 force field largely agrees with the experimental value, albeit slightly underestimating it. The modified FF2 and FF3 force fields rectify this deviation, yielding values within the range of experimental estimations. Overall, there is significant alignment between these experimental and computational density values, supporting the transferability of the CL&Pol force field. This conclusion aligns with the findings of Goloviznina et al.,^15,29 who conducted a comparative analysis of experimental and computational results across diverse ILs and DES types. The calculated diffusion coefficients and both partial and total X-ray structure factors are included in the ESI† (Sections S3 and S4) to facilitate direct comparisons with additional experimental data in the future.

Table 3 Density of ChCl : LA DES at 1 : 1 ratio, temperature T = 333 K and pressure P = 1 bar obtained through experiments⁵⁵ and by CMD simulations employing the force fields considered in this work

	ρ (g cm^−3)
Exp.^55	1.13–1.14
FF1	1.121
FF2	1.130
FF3	1.132

---

We now analyze the impact of the modified force fields on the overall hydrogen bond network of the liquid. We do it through the combined distribution functions (CDF) where the RDF between atoms forming the H-bond (x-axis) and the angular distribution function that contains these two atoms (y-axis) are shown in combination. Fig. 3 shows the CDFs of the H-bonds for which the σ parameter has been optimized. In terms of occurrence values, these three type of H-bonds emerge as the strongest ones throughout the simulations (see Fig. S5 (ESI†) for comparison with the remaining type of H-bonds). Notably, FF1 underestimates the interaction with LACOOH and LAOH and slightly overestimates the interaction with choline (see Fig. 3). The agreement between AIMD with the FF2 and FF3 force fields regarding the H-bonds shown in Fig. 3 is highly satisfactory, despite the large differences on the size of the system and simulation timescale employed at each level of theory. Importantly, the observed similarities extend beyond RDFs, which were the ones employed during the fitting process. The different ADFs also exhibit large agreement. This supports the use of small size systems containing 16 molecules of each chemical compound for the fitting process.

Fig. 3 CDFs of the Cl⁻–OH interactions taking place in the system: (a) Cl⁻–ChOH (b) Cl⁻–LACOOH and (c) Cl⁻–LAOH, at the different theory levels considered.

Concerning the remaining hydrogen bonds in the system, the corresponding CDFs are provided in Section S5 in the ESI.† Overall, FF1 tends to underestimate the strength of the hydrogen bond network in comparison to DFT, with a unique exception (Fig. 3(a)).

Regarding the influence of the modified force fields in the H-bond network, two distinct scenarios emerge: those in which LA acts as a hydrogen bond donor, interacting with either choline or another LA. For these specific H-bonds, the three distinct force fields give similar deviations with respect to AIMD (see ESI†). Conversely, when choline serves as a donor interacting with either LAs or other cholines, the underestimation by the FF1 force field is significantly rectified in FF2 and FF3, resulting in an improved agreement with AIMD results. This difference can be attributed to the σ correction, which weakens the Ch–Cl interaction and consequently, more choline molecules are available to form H-bonds with other cholines or LAs. This analysis suggests that the tuning of the selected σ parameters enhances not only the description of their respective H-bond interactions but also indirectly improves, on a global scale, the representation of the entire H-bond network.

3.2 ChCl:LA structure as the LA excess increases

In this section, we describe the structural changes of the system as the ChCl : LA ratio increases, with a particular focus on understanding how different versions of the force field influence the predicted structure of the material.

Fig. 4 shows the RDFs and coordination number between the Cl⁻–ChOH (a), Cl⁻–LACOOH (b) and Cl⁻–LAOH (c) pairs of atoms for each ChCl : LA ratio: 1 : 1, 1 : 2, 1 : 5 and 1 : 10. It is worth noting that the Cl⁻–LAOH population is substantially lower than that of the other two. This can be deduced by comparing the maximum values of the RDFs or the coordination numbers in Fig. 4. Even if g(r) increases upon the σ correction in FF2 and FF3, the Cl⁻–LAOH bond is still significantly less populated than the other two, more evidently as the ChCl : LA ratio increases. Two distinct trends emerge in these interactions as the ChCl : LA ratio varies, depending on whether the employed force field version exhibits overpolarization or not.

Fig. 4 (a–c) Radial distribution functions and coordination numbers of (a) Cl⁻–ChOH, (b) Cl⁻–LACOOH and (c) Cl⁻–LAOH interactions obtained with the FF1 (blue), FF2 (red) and FF3 (green) force fields for the different ChCl : LA ratios considered: the darker the color, the larger the LA molecule ratio. (d) Values of all these RDFs at r₀ = 800 pm for each type of interaction: from left to right, Cl⁻–ChOH, Cl⁻–LACOOH and Cl⁻–LAOH, respectively.

In the FF1 force field, the population of Cl⁻–ChOH remains consistently close to 1, despite changes in the LA ratio (see panel (a) in Fig. 4). This shows that although the relative LA molecule number increases with the ratio, they are not affecting the Cl⁻–ChOH population. In the FF2 force field, a similar behavior is observed, with the coordination number experiencing a minor decrease as the LA ratio increases. Conversely, the RDFs obtained with the inclusion of overpolarization correction in FF3 show a gradual decrease in the coordination number as the ratio shifts from 1 : 1 to 1 : 2, followed by a significant drop at a ratio of 1 : 5 and 1 : 10. This clearly indicates that, with this force field, the increase in LA molecules affects the Cl⁻–ChOH population by replacing some of the chloride ions with LA molecules.

In contrast, the Cl⁻–LACOOH coordination number (panel (b) in Fig. 4) experiences a consistent reduction across all considered force fields as the LA excess increases. Regardless of the force field used, this behavior is anticipated, as the relative number of Cl⁻ ions in the mixture diminishes with an increase in the LA ratio. FF1 and FF2 predict a coordination number lower than 0.5, which further decreases with increasing excess. FF3, however, predicts a slightly higher coordination number, and the reduction is not as pronounced. The Cl⁻–LAOH pair interaction (see Fig. 4(c)), exhibits significant differences between the original and modified force fields. These differences primarily stem from the notable changes in the σ parameter associated with this interaction. As a result, the modified force fields display larger coordination numbers.

Longer-range correlations can be inferred by analyzing the RDFs at greater distances. Fig. 4(d) presents the values of each RDF from Fig. 4(a–c) at a distance of r₀ = 800 pm, approximated to contain multiple solvation shells. For ease of comparison, a horizontal line representing the g = 1 value, indicative of an ideally uncorrelated pair of atoms, is included. Notably, in the Cl⁻–ChOH interaction, the RDFs from FF1 and FF2 show significant deviations from unity with increasing LA ratios, hinting at the formation of ChCl clusters in the sample. Conversely, the FF3 force field exhibits a qualitatively different behavior, yielding RDF values closer to an ideally uncorrelated system, especially at higher LA ratios (1 : 5 and 1 : 10). The RDF values for Cl⁻–LACOOH and Cl⁻–LAOH interactions corroborate this idea: FF3 consistently yields values closest to g = 1, regardless of the LA ratio. This trend is particularly noticeable, again, in the 1 : 5 and 1 : 10 LA ratios. In summary, these observation further suggests that, among the three force fields, FF3 tends to significantly achieve a more homogeneous mixing of the different molecular compounds in the system as the ratio of LA molecules increases.

This difference is expected to exert a substantial influence on the overall structure of the DES. For this purpose, we classify the molecules using Voronoi tessellation, distinguishing between ChCl molecules on one side and LA molecules on the other. The time evolution of ChCl cluster counts during the MD simulation is shown in the upper panel of Fig. 5. The behavior of this parameter exhibits noteworthy distinctions between different FF versions considered. For FF1, convergence is rapidly achieved, showing a consistent formation of approximately four clusters across various ChCl : LA ratios. With FF2, a similar amount of ChCl cluster formation is observed, albeit with a slight increase as the LA molecule ratio rises. Notably, FF3 presents the most distinct scenario: the number of clusters varies significantly among different ChCl : LA ratios, reaching its largest value at approximately 36 ChCl clusters for the 1 : 10 DES ratio.

Fig. 5 (a) Number of ChCl clusters and (b) ChCl clustering degree evolution through simulation time for the different force field considered in this work: FF1 (blue), FF2 (red) and FF3 (green). The ChCl : LA ratios examined (1 : 1, 1 : 2, 1 : 5, and 1 : 10) are indicated by the color intensity: darker color represent higher LA ratios. The values depicted are the averages of each 1000 snapshots. (c) Neighbour count at different ChCl : LA ratios and force fields (same color pattern) between pairs of molecular species: from left to right, Cl⁻–Ch⁺, Cl⁻–LA, and Ch⁺–LA.

The clustering degree⁵⁶ (central panel in Fig. 5) offers further insights into the characteristics of the formed clusters. Defined as the ratio between the number of ChCl molecules in the largest cluster and the total number of ChCl molecules in the system, it provides the proportion of ChCl molecules that form part of the largest cluster. In the context of FF1, this ratio is remarkably close to unity across different ChCl : LA ratios, indicating that nearly all the ChCl molecules belong to a single, large cluster. A similar pattern emerges with FF2, where all considered systems display clustering degrees approaching one, albeit with minor distinction for 1 : 10 ChCl : LA ratio, where a clustering degree of around 0.9 is obtained. Finally, for FF3, the clustering degree is also close to one for the 1 : 1 and 1 : 2 ratios. However, a significant difference from this trend is observed for the 1 : 5 (0.87) and 1 : 10 (0.30) ratios, indicating a lower level of clustering among ChCl molecules. It can be concluded that, in these instances, ChCl molecules are distributed more evenly across multiple clusters.

Based on the Voronoi tessellation scheme, we calculate the neighbour count which gives the number of times that a given molecular species shares a Voronoi contact area with another molecular species. Fig. 5(c) shows the neighbour count between Cl⁻–Ch⁺, Cl⁻–LA, and Ch⁺–LA molecules at different ratios with the three distinct force fields. The obtained results align with what is obtained in the number of clusters and clustering degree, FF3 tends to evenly distribute the molecular compounds conforming the DES, while FF1 tends to separate ChCl and LA molecules. The neighbor count also indicates that the structural differences given by each force field are enhanced as the LA ratio increases, showing that the overpolarization correction included in the FF3 force field does not show large structural differences with respect to FF2 at low LA ratios, but it significantly affects the structure of the liquid at large LA ratios.

These structural differences are also discernible when observing the last snapshots of the simulation box (see Fig. 6). At a DES ratio of 1 : 1, FF1 exhibits a clear segregation between the ions and neutral molecules, a phenomenon that becomes increasingly pronounced with rising LA ratios. FF2 begins to display segregation at a ratio of 1 : 2, a trend that persists at higher ratios. FF3, on the other hand, presents a markedly more homogeneous mixture across all examined ratios, indicating a qualitative improvement with respect to FF1 and FF2.

Fig. 6 Final snapshot of the simulation box captured for each ChCl : LA ratio (1 : 1, 1 : 2, 1 : 5, and 1 : 10) and force field (FF1, FF2, and FF3). ChCl molecules are represented in blue, while LA molecules in red.

Regarding the dynamical properties of the liquid, the calculated diffusivities (see ESI†) indicate that FF3 produces the highest ionic diffusivities and the lowest ones for lactic acid. This behavior is consistent with the fact that FF3 yields the most homogeneous liquid structure among the three force field versions, leading to fewer direct ionic interactions and thus increasing the mobility of the ions. At the same time, lactic acid forms stronger interactions with the ions, which reduces its diffusivity. Furthermore, we observe a general trend, independent of the force field employed, where the diffusivities of all molecular species increase as the lactic acid ratio increases.

Finally, note that the time required to obtain a converged value for these large-scale structural parameters is often around 100 ns. For instance, the clustering degree of the 1 : 10 ChCl : LA ratio with the FF2 force field takes approximately 80 ns to converge, while the collapse degree exhibited by all force fields with a 1 : 1 ratio is achieved in about 100–120 ns. This highlights that the dynamics of these clusters demand substantial simulation times to achieve an equilibrated structure.

4 Discussion

In this work, two different aspects of the CL&Pol force field are refined for the ChCl:LA system. The first of them, the σ parameter of the three strongest H-bond interactions are tuned based on AIMD 1 : 1 simulations. This parameterization is highly system-dependant, and affects significantly the accuracy of the force field in capturing the short-range structure of the DES. In particular, focusing on the main interactions in the system, the original force field overestimates the Cl⁻–ChOH interaction and underestimates Cl⁻–LA interactions, particularly Cl⁻–LAOH, as compared to AIMD. The agreement is optimized by tuning the σ value to fit the AIMD data, obtaining large agreement between both levels of theory. Furthermore, it is observed that as a consequence, the improvement is extended to the whole H-bond network, optimizing the overall agreement with AIMD without falling in the overparameterization.

The second aspect tackled is the overpolarization of Cl⁻, which is corrected by applying the Tang Toennies damping function between the induced dipoles interactions. Only minor modifications are observed in the 1 : 1 ratio, mainly in the Cl⁻–LAOH interaction, pointing out to the fact that this is the interaction that depends most on the charge–induced dipole interactions of the system. The description of the Cl⁻–Cl⁻ RDF, is also improved by the overpolarization correction (see Fig. S2, ESI†) as compared to the AIMD calculations, but only in terms of the g(r) value. The position of the first peak does not shift across the three employed force fields, which stays with a deviation of 100 pm from the AIMD spectra. The Tang Toennies correction did induce a shift of 200 pm the Cl⁻–Cl⁻ peak of the ChCl–EG system, leading to a disagreement of 70 pm with AIMD, not far from the difference in our system for all the force field versions. However, the Cl⁻ anions were predicted to be further away from each other by AIMD calculation in that system, pointing out to a less structured environment with significant influence of the polarization effects. We can conclude that the Tang Toennies correction does not crucially alter the structure of the system at this point. This already suggests the lack of relevance of the Cl⁻–LAOH interaction in the system, the most susceptible interaction among the three charge–dipole interactions. This is probably due to the fact that the main interaction between ChCl and LA is the one taking place with the carboxylic moiety, on the contrary to other systems in the literature having solely an alcohol group in the neutral molecule.

Interestingly, the relevance of both corrections is reversed as the ChCl–LA ratio increases. We observe that both FF1 and FF2 reported similar results in terms of coordination numbers in different ratios (see Fig. 4); increasing (decreasing) the σ value for the Cl⁻–ChOH (Cl⁻–LACOOH) interaction, weakening (strengthening) it, alters its coordination numbers to Cl⁻ only weakly. The influence of the overpolarization, nevertheless, is crucial as the ratio increases. The coordination number in the Cl⁻–ChOH interaction decreases upon the increasing the ratio in FF3, stating, therefore, that the overpolarization remarkably affects the population of the Cl⁻–ChOH interaction. While in ratio 1 : 1 this has no remarkable influence, as the ratio is increased, the higher availability of the Cl⁻ in FF3 induces more probable interactions with the LA, in line with several works that show ChCl-based DES hydrogen bond network to rely on the competition between cation and HBD for the anion.^8,57,58

It is the balance of these interactions that determines the structure of the material. It is worth noting that correcting the underestimation of the Cl⁻–LAOH interaction by FF1 does not imply a substantial increase of its population, and remains a second level interaction that does not govern the overall structure. Given that Cl⁻–ChOH is the strongest interaction taking place between the ions, it tunes the internal cohesion of the HBA, the ChCl. This internal cohesion is overestimated by FF1 and FF2 due to the overpolarization; when this is corrected, the interaction between the chlorine and LA is improved, therefore ameliorating the mixing between ChCl and LA. Evidently, this has important implications for the overall structure of the material that are reflected in distinct behaviours of the clustering and collapse degrees in each system.

The collapse degree and neighbour count analysis complement each other in providing a general picture of the structure at the different ratios and the different force field versions considered. Both FF1 and FF2 show similar structural results: the clustering degree is near to one, meaning that almost all the ChCl molecules form part of the same cluster, and therefore, a very strong segregation takes place in the system. The neighbour count slightly differs between both force fields, with FF2 tending towards a more homogeneous liquid.

In an homologous situation, Maglia de Souza et al.¹⁶ detected an artificial phase-separation problem in the 1 : 2 ChCl : EG system that was corrected by decreasing the σ value of Cl⁻–OH interactions and therefore strengthening them. In our case, σ fitting does not affect significantly the long range structural behaviour of our system, which is still oversegregated, even though the σ-corrected FF2 force field exhibits a modest decrease of the segregation between ions and neutral molecules compared to the original force field.

On the contrary, once the overpolarization correction is included in FF3, a significant structural change is predicted for the system as the LA excess increases. FF3 markedly enhances the structure of ionic and neutral clusters, achieving a more homogeneous mixture, especially at the highest DES ratios examined in this study: 1 : 5 and 1 : 10. This improvement is evidenced by the low clustering degree and large differences in neighbour count given by FF3 with respect to FF1 and FF2.

5 Conclusions

The artifacts resulting from the limitations detected in the CL&Pol force field for ChCl-based DES have been seen to affect also the carboxylic acid-based ChCl:LA both in the Cl⁻–LAOH and Cl⁻–LACOOH interactions. Particularly, fixing the σ value for all the Cl⁻–OH interactions affects mostly the short-range structure of the system, while the overpolarization results in profound differences when it comes to the long-range structure.

Indeed, the former leads to a better agreement with AIMD, but does not significantly alter its overall structure in different ratios, on the contrary to what was observed in the literature for other systems.¹⁶ It is interesting to observe that despite the Cl⁻–LAOH interaction being the most affected interaction in the system, correcting its underestimation by means of the σ parameterization does not result in a significant role on the overall structure of the liquid, which is mostly tuned by Cl⁻–ChOH and Cl⁻–LACOOH. On the other hand, it is worth noting that the optimal σ values obtained for the Cl⁻–LAOH interaction is very similar as the one obtained in the literature for the Cl⁻–EG interaction, but not that of Cl⁻–ChOH, probably due to the different competing interactions governing the overall structure.

Regarding the overpolarization correction, it does not significantly alter the short-range structure of the 1 : 1 system. The Cl⁻–LAOH interactions seem to be more affected by the contribution to the interaction than the Cl⁻–LACOOH interactions; nevertheless, we observe that the overpolarization correction alters very significantly the medium-long range structure by tuning the Cl⁻–ChOH and Cl⁻–LACOOH competition. This results in a better mixing of the HBD and HBA molecules, resulting in structural changes that become more relevant as the ChCl : LA ratio increases. The maleability of the Cl⁻–ChOH interaction is evidenced and its coordination number becomes lower as the LA content increases. This is due to the screening of the charge–induced dipole interactions, that have as a consequence the better availability of Cl⁻ and therefore, the competition and interplay between the intra HBA and inter HBA–HBD interactions.

Consequently, an evolution of the structure of the liquid as the ChCl:LA increases is observed, achieving, a better mixed, more homogeneous structure in the highest ratios (1 : 5, 1 : 10). This results in a higher contact area between HBD and HBA and a better availability of the main interacting moieties (namely Cl⁻ and OH) that might alter the performance of the DES in practical applications.

In addition to these observations, the evaluation of the three force fields studied in this work highlights several key aspects of their performance in predicting the structural properties of ChCl:LA. While FF1 produces a density estimation close to the available experimental data, the modified force field FF3 offers an improved estimate, falling within the experimental range. Additionally, the issue of chloride overpolarization, previously observed in similar systems such as ChCl:EG,¹⁶ has been addressed through modifications in FF3. This improvement is consistent with the modifications that corrected unphysical phase separation for ChCl:EG.

Given the limited experimental data currently available for ChCl:LA, particularly beyond density measurements, DFT calculations serve as the most reliable reference for validating force fields in the short-range structure. Both FF2 and FF3 demonstrate a significant improvement over FF1 in capturing the hydrogen-bond network, which is critical for accurately modeling the system. Based on these factors, FF3 emerges as the most suitable force field for classical molecular dynamics simulations of ChCl:LA. However, additional experimental data, such as diffusivity measurements or the characterization of ionic domain sizes, would be crucial for further validation and refinement of the force fields. The availability of such data would allow for a more definitive selection of the most appropriate force field for this system.

Author contributions

E. R. performed the AIMD calculations and supervised the project. J. Z. conducted the CMD calculations and data analysis. All authors contributed to the conceptualization of the work, and were involved in writing, reviewing, and editing the manuscript.

Data availability

Data for this article, including LAMMPS input files are available at github at https://github.com/jony-zz/ChCl-LA.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge the Basque Government - Eusko Jaurlaritza (IT1254-19, IT1584-22), for financial support and the SGi/IZO-SGIker UPV/EHU for generous allocation of computational resources. We thank Dr Óscar Pozo and Dr David de Sancho for useful discussion.

Notes and references

1. B. B. Hansen, S. Spittle, B. Chen, D. Poe, Y. Zhang, J. M. Klein, A. Horton, L. Adhikari, T. Zelovich, B. W. Doherty, B. Gurkan, E. J. Maginn, A. Ragauskas, M. Dadmun, T. A. Zawodzinski, G. A. Baker, M. E. Tuckerman, R. F. Savinell and J. R. Sangoro, Chem. Rev., 2021, 121(3), 1232–1285.
2. H. Shekaari, M. T. Zafarani-Moattar, M. Mokhtarpour and S. Faraji, Sci. Rep., 2023, 13, 1–11.
3. M. B. Bianchi, C. Zhang, E. Catlin, G. Sandri, M. Calderón, E. Larrañeta, R. F. Donnelly, M. L. Picchio and A. J. Paredes, Mater. Today Bio, 2022, 17.
4. D. Uschapovskiy, V. Vorobyova, G. Vasyliev and O. Linucheva, J. Electrochem. Sci. Eng., 2022, 12, 1025–1039.
5. A. B. Botelho Junior, G. Pavoski, M. D. C. R. da Silva, W. L. da Silva, D. A. Bertuol and D. C. R. Espinosa, Nano Technology for Battery Recycling, Remanufacturing, and Reusing, 2022, pp. 79–103.
6. Y. Hariyanto, Y. K. Ng, Z. Z. Siew, C. Y. Soon, A. C. Fisher, L. Kloyer, C. W. Wong and E. W. C. Chan, Energy Fuels, 2023, 37, 18395–18407.
7. K. Chai, Y. Zhou, X. Lu, T. Yamaguchi, K. Ohara, H. Liu and F. Zhu, Phys. Chem. Chem. Phys., 2023, 25, 10481–10494.
8. S. Spittle, D. Poe, B. Doherty, C. Kolodziej, L. Heroux, M. A. Haque, H. Squire, T. Cosby, Y. Zhang, C. Fraenza, S. Bhattacharyya, M. Tyagi, J. Peng, R. A. Elgammal, T. Zawodzinski, M. Tuckerman, S. Greenbaum, B. Gurkan, M. Dadmun, E. J. Maginn and J. Sangoro, Nat. Commun., 2022, 13(1), 1–14.
9. R. Elfgen, O. Hollóczki and B. Kirchner, Acc. Chem. Res., 2017, 50, 2949–2957.
10. S. Sahu, S. Banu, A. K. Sahu, B. V. Phani Kumar and A. K. Mishra, J. Mol. Liq., 2022, 350, 118478.
11. A. Malik and H. K. Kashyap, Phys. Chem. Chem. Phys., 2021, 23, 3915–3924.
12. S. Chatterjee, T. Haldar, D. Ghosh and S. Bagchi, J. Phys. Chem. B, 2020, 124, 3709–3715.
13. S. Kaur, A. Gupta, H. K. Kashyap and H. K. Kashyap, J. Phys. Chem. B, 2020, 124, 2230–2237.
14. N. M. Stephens and E. A. Smith, Langmuir, 2022, 38, 14017–14024.
15. K. Goloviznina, Z. Gong, M. F. Costa Gomes and A. A. H. Pádua, J. Chem. Theory Comput., 2021, 17(3), 1606–1617.
16. R. Maglia De Souza, M. Karttunen and M. C. C. Ribeiro, J. Chem. Inf. Model., 2021, 61, 5938–5947.
17. A. Szabadi, R. Elfgen, R. Macchieraldo, F. L. Kearns, H. Lee Woodcock, B. Kirchner and C. Schröder, J. Mol. Liq., 2021, 337 year.
18. V. K. Vorobiov, M. A. Smirnov, N. V. Bobrova and M. P. Sokolova, Mater. Lett., 2021, 283, 128889.
19. M. A. Smirnov, A. L. Nikolaeva, V. K. Vorobiov, N. V. Bobrova, I. V. Abalov, A. V. Smirnov and M. P. Sokolova, Polymers, 2020, 12, 350.
20. R. Ninayan, A. S. Levshakova, E. M. Khairullina, O. S. Vezo, I. I. Tumkin, A. Ostendorf, L. S. Logunov, A. A. Manshina and A. Y. Shishov, Colloids Surf., A, 2023, 679, 132543.
21. H. Shekaari, M. T. Zafarani-Moattar, M. Mokhtarpour and S. Faraji, Sci. Rep., 2023, 13, 11276.
22. P. Strižincová, I. Šurina, M. Jablonsky, V. Majová, A. Ház, K. Hroboňová and A. Špačková, Processes, 2024, 12, 653.
23. S. Karimarji, A. Khorsandi, G. Azimi and Z. Mardani, Opt. Mater., 2024, 148, 114912.
24. A. K. Kumar, S. Sharma, E. Shah and A. Patel, J. Mol. Liq., 2018, 260, 313–322.
25. Y. T. Tan, G. C. Ngoh and A. S. M. Chua, Ind. Crops Prod., 2018, 123, 271–277.
26. A. M. Da Costa Lopes, Acta Innovations, 2021, 40, 64–78.
27. R. Alcalde, A. Gutiérrez, M. Atilhan and S. Aparicio, J. Mol. Liq., 2019, 290, 110916.
28. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Compos. Phys. Commun., 2022, 271, 108171.
29. K. Goloviznina, J. N. Canongia Lopes, M. Costa Gomes and A. A. H. Pádua, J. Chem. Theory Comput., 2019, 15, 5858–5871.
30. A. Dequidt, J. Devemy and A. A. Padua, J. Chem. Inf. Model., 2016, 56, 260–268.
31. A. A. Pádua, fftool, 2013, https://github.com/paduagroup/fftool.
32. L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, Nucleic Acids Res., 2017, 45, W331–W336.
33. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, 2017.
34. L. Martnez, R. Andrade, E. G. Birgin and J. M. Martnez, J. Comput. Chem., 2009, 30, 2157–2164.
35. K. Goloviznina and A. A. Pádua, CL Pol, 2019, https://github.com/paduagroup/clandpol.
36. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51(8), 2007–2023.
37. M. Brehm, M. Thomas, S. Gehrke and B. Kirchner, J. Chem. Phys., 2020, 152(16), 164105.
38. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
39. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152 year.
40. J. Vandevondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
41. M. Frigo and S. G. Johnson, Proc. IEEE, 2005, 93, 216–231.
42. J. Kolafa, J. Comput. Chem., 2004, 25, 335–342.
43. J. VandeVondele and J. Hutter, J. Chem. Phys., 2003, 118, 4365–4369.
44. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
45. S. Nose, Mol. Phys., 1984, 52, 255–268.
46. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
47. C. Lee, W. Yang and R. G. Parr, Development, 1988, 37 year.
48. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
49. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132.
50. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
51. M. Krack, Theor. Chem. Acc., 2005, 114, 145–152.
52. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
53. S. Goedecker and M. Teter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
54. V. Alizadeh, F. Malberg, A. A. Pádua and B. Kirchner, J. Phys. Chem. B, 2020, 124, 7433–7443.
55. D. I. Al-Risheq, M. Nasser, H. Qiblawey, I. A. Hussein and A. Benamor, Sep. Purif. Technol., 2021, 255, 117737.
56. Y. Tang, S. Bera, Y. Yao, J. Zeng, Z. Lao, X. Dong, E. Gazit and G. Wei, Cell Rep. Phys. Sci., 2021, 2, 100579.
57. V. C. Ferreira, G. Marin, J. Dupont and R. R. Correia, J. Phys.: Condens. Matter, 2020, 33, 1–7.
58. Y. Zhang, D. Poe, L. Heroux, H. Squire, B. W. Doherty, Z. Long, M. Dadmun, B. Gurkan, M. E. Tuckerman and E. J. Maginn, J. Phys. Chem. B, 2020, 124, 5251–5264.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp03233k

---