Insights on the cellulose pretreatment at room temperature by choline-chloride-based deep eutectic solvents: an atomistic study

The pretreatment or disruption of a cellulose I β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} crystallite by four deep eutectic solvents (DES): choline-chloride ethylene glycol, choline-chloride oxalic acid, choline-chloride urea, and choline-chloride levulinic acid, was described from the atomistic interactions observed in molecular dynamics and ab-initio (NCI-AIM) studies. The cellulose I β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} disruption was studied considering plausible correlations between the Kamlet-Taft (α and β) solvent parameters, and a series of thermodynamic, structural, and energetic properties. It was found that the Kamlet-Taft parameters correlated with the thermodynamic properties of the DES, as well as their variations after the addition of the cellulose crystallite. Structural analysis revealed that the weaker the interactions within the molecules of the solvent, the stronger the interactions between the hydroxyl group from cellulose with the chloride anion and with the hydrogen bond donor. Further analysis indicated that the R-CO-R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R-CO-R$$\end{document} moieties in the hydrogen bond donor within the solvent, displayed best interplaying with the cellulose. The hydrogen bond occupancies within the cellulose crystallite, evidenced that the main O6-H6⋯O2/O3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O6-H6\cdots O2/O3$$\end{document} and O2-H2⋯O6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O2-H2\cdots O6$$\end{document} interchain hydrogen bonds in the glucan located at the edge of the solute, were replaced by weak O6-H6⋯O4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O6-H6\cdots O4$$\end{document} hydrogen bonds in all solvents. This effect was related to the O/ClDES⋯H-donor\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${O/Cl}_{DES}\cdots H-donor$$\end{document} and ODES-HDES⋯acceptor\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${O}_{DES}-{H}_{DES}\cdots acceptor$$\end{document} HBs between cellulose and DES molecules, and it was confirmed by the non-covalent interactions obtained through DFT calculations. The energetic interactions and the atomistic degree of disruption of the cellulose crystallite, were not completely described by the Kamlet-Taft β or α parameters when considered separately. Surprisingly, by using the net basicity (β-α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta -\alpha $$\end{document}) definition such correlations were improved, suggesting that both parameters must be considered together to develop new, green, and sustainable solvents for cellulose pretreatment process.

variations after the addition of the cellulose crystallite. Structural analysis revealed that the weaker the interactions within the molecules of the solvent, the stronger the interactions between the hydroxyl group from cellulose with the chloride anion and with the hydrogen bond donor. Further analysis indicated that the R − CO − R moieties in the hydrogen bond donor within the solvent, displayed best interplaying with the cellulose. The hydrogen bond occupancies within the cellulose crystallite, evidenced that the main O6 − H6 ⋯ O2∕O3 and O2 − H2 ⋯ O6 interchain hydrogen bonds in the glucan located at the edge of the solute, were replaced by weak O6 − H6 ⋯ O4 hydrogen bonds in all solvents. This effect was related to the O∕Cl DES ⋯ H − donor and O DES − H DES ⋯ acceptor HBs between cellulose and DES molecules, and it was confirmed by the noncovalent interactions obtained through DFT calculations. The energetic interactions and the atomistic degree of disruption of the cellulose crystallite, were not completely described by the Kamlet-Taft β or α parameters when considered separately. Surprisingly, by using the net basicity ( − ) definition such correlations were improved, suggesting that both parameters must be considered together to develop new, green, and sustainable solvents for cellulose pretreatment process.

Introduction
The global demand of energy is constantly increasing (Zhang et al. 2017), in consequence new efficient, ecological, and economical energy resources and processes are urgently needed (van Putten et al. 2013;Clarke et al. 2018;Körner et al. 2019). Lignocellulosic biomass (LCB) is the largest and greenest source of carbon on Earth (Tang et al. 2017;Körner et al. 2019). The use of the LCB as renewable source has been proposed to produce chemicals, materials, fuels, and energy, currently based on fossil oils (Gunny et al. 2015;Clarke et al. 2018;Tiong et al. 2018;Chen and Mu 2019;Zhou et al. 2020). The efficient use of LCB requires solvents capable to pretreating or dissolving its components (van Putten et al. 2013;Rocha et al. 2017;Tiong et al. 2018), namely, cellulose (40-80%), hemicellulose (15-30%), and lignin (10-25%) (van Osch et al. 2017;Kumar et al. 2020). These processes facilitate accessibility of either incoming third particles (such as reactants, catalysts, etc.) or enzymes (Procentese et al. 2015), to produce valuable products. Unfortunately, the common solvents used for these purposes are organic solvents, as well as strong mineral solvents H 2 SO 4 , H 3 PO 4 , and NaOH, which are widely recognized by their potential environment concerns (Ren et al. 2016a;Tang et al. 2017;Sert et al. 2018;Kumar et al. 2020).
The experimental ability of common DESs, including those with choline-chloride (ChCl) as HBA, to dissolve LCB is low compared to that achieved by ILs (Francisco et al. 2012;Zhang et al. 2012;Sharma et al. 2013;van Osch et al. 2017;Morais et al. 2020;Sosa et al. 2020). However, it has been reported that DESs are able to pretreat the LCB (Xia et al. 2014;Procentese et al. 2015; Kalhor and Ghandi 2019; Chen and Mu 2019; Morais et al. 2020), providing a green and novel way to synthesize valuable derivatives (Sirviö et al. 2015(Sirviö et al. , 2016. On this regard, the LCB pretreament with DESs can enhance subsequent steps on the biomass valorization (Xia et al. 2014;Procentese et al. 2015;Lu et al. 2016), such as the enzymatic hydrolysis (Xia et al. 2014;Gunny et al. 2015) or acting as an effective media for the cellulose modification (Sirviö et al. 2015(Sirviö et al. , 2016Tang et al. 2017;Smirnov et al. 2020).
Recently, the solvent behavior of DESs to dissolve or pretreat the LCB has been related to the ability of the DESs to form hydrogen bonds (HBs) with LCB components (Sharma et al. 2013;Ren et al. 2016a;Lu et al. 2016;Häkkinen and Abbott 2019). Experimentally, Ren et al. (2016aRen et al. ( , 2016b, Xia et al. (2018), and Liang et al. (2021) have found that the solubility, pretreatment, or product yield of the LCB can be correlated to the Kamlet-Taft α and β parameters of the DESs. Furthermore, experimental results on LCB solvation, reported in the literature (Procentese et al. 2015;Sirviö et al. 2016;Liu et al. 2019;Tong et al. 2021;Sun et al. 2022), have the possibility to be also correlated with the available Kamlet-Taft parameters (Kundu et al. 2020).
The Kamlet-Taft α and β parameters describe the solvent polarity and the hydrogen bonding ability of the solvent (Abe et al. 2010(Abe et al. , 2015Cláudio et al. 2014;Ren et al. 2016b). Particularly, the ability of the solvent to donate protons in the hydrogen bonds formation (also known as acidity) towards the solute, can be measured by the Kamlet-Taft α parameter. On the other hand, the ability of the solvent to accept protons to form hydrogen bonds (also known as basicity) is measured by the Kamlet-Taft β parameter (van Osch et al. 2017). The Kamlet-Taft α and β parameters can also be used as pre-screening tool to propose, select, and develop effective task-specific solvents (Islam et al. 2020), including DESs for LCB pretreatment or dissolution (Ren et al. 2016b;Tong et al. 2021). Unfortunately, the experimental values for Kamlet-Taft parameters are scarce (Kundu et al. 2020), not to mention the difficulties related to the experimental procedures to obtain them (Cláudio et al. 2014;Kundu et al. 2020). On this regard, several authors have made significant efforts in the search for new factors that may influence the dissolution or pretreatment of the LCB and provide a reliable pre-screening of the effectiveness of DESs pretreatment. Ren et al. (2016a) evaluated the ability of a series of DESs to dissolve cellulose. The authors reported that dissolution of cellulose was correlated with the Kamlet-Taft β parameter values, but also with thermodynamic properties such as viscosity and conductivity of the DESs. Recently, Liang et al. (2021) correlated the experimental degree of LCB pretreatment in series of lactic acid-based DESs with their Kamlet-Taft α parameter values. Surprisingly, the authors also found that a notable correlation was obtained by using the DESs density. Interestingly, the use of the density and viscosity has been exploited by Abe et al. (2010Abe et al. ( , 2015 to design and synthetize efficient ionic liquids solvents to dissolve cellulose, including the Kamlet-Taft parameter values. Thermodynamic properties of the DES are related to the intrinsic interactions between HBD and HBA in the bulk of the DES (Tolmachev et al. 2022). These solvent-solvent interactions govern the solute-solvent interactions that may emerge in the system (Sosa et al. 2020;Tolmachev et al. 2022), as it was suggested by Xia et al. (2018) through a series of experiments, and further investigated by means of ab-initio calculations carried out by Fu et al. (2020). Unfortunatelly, the interaction mechanism between DESs and solutes such as cellulose and the relation with its dissolution or pretreatment process is still unclear (Chen and Mu 2019;Smirnov et al. 2020;Alizadeh and Kirchner 2021). Theoretically, molecular dynamics (MD) simulations have provided an efficient route to reveal, understand, and predict the microscopic solvent-solvent and solute-solvent interactions that drive the LCB dissolution (Wang et al. 2012;Mohan et al. 2017;van Osch et al. 2017;Xia et al. 2018;Fu et al. 2020;Smirnov et al. 2020;Zhou et al. 2020). Since the cellulose is the main component within the LCB (40-80%) (van Osch et al. 2017;Kumar et al. 2020), most of the LCB-DES theoretical works reported in the literature are focused on both cellulose-DES (Smirnov et al. 2020;Alizadeh and Kirchner 2021) and glucose-DES (Mohan et al. 2017;Fu et al. 2020) systems. Mohan et al. (2017) and Smirnov et al. (2020) have used classical MD simulations to reveal the glucose-DES and the cellulose surface-DES interactions, respectively. The authors found that the anion from the HBA and the hydroxyl groups from the HBD are the main interactions towards the hydroxyl groups from the carbohydrates, and therefore responsible of the glucose dissolution and the surface deconstruction process. Alizadeh and Kirchner (2021) carried out MD simulations in order to understand the interaction of cellulose in two DESs: choline-chloride urea (ChCl/ urea) and choline-acetate urea. Their results suggested that choride anion in ChCl/urea play a sustancial role in the initial stage of the cellulose dissolution, as a consequence of its proximity to glucose units.
Recently, Fu et al. (2020) performed quantum chemical calculations in order to evaluate the strength of the interactions between cellobiose and DES molecules. The authors found that HBs between the HBA and the HBD are weakened when they interact with the cellobiose. This interaction reduces the cellobiose hydroxyl bond strength; consequently, the cellobiose HB network becomes weaker, favoring its dissolution. Accordingly, HBs between solute and solvent and between DES molecules might be the driving forces of such processes.
As indicated, the fundamental interactions within the DESs as well as between solute and DESs that conduct the LCB pretreatment, or its dissolution could be revealed through MD simulations (Mohan et al. 2017;Fu et al. 2020;Smirnov et al. 2020;Alizadeh and Kirchner 2021). Additionally, such processes can be further described in terms of the Kamlet-Taft α and β parameters (Ren et al. 2016a(Ren et al. , 2016bXia et al. 2018;Liang et al. 2021). Furthermore, as reported Ren et al. (2016a) and Liang et al. (2021), such processes might be related with the thermodynamic properties of the DESs. On this regard, there could be a bridge between solubility or pretreatment degree of the LCB (described through the Kamlet-Taft parameters) and the pure DES solvent properties, explained from the fundamental solvent-solvent and solute-solvent interactions. The aforementioned is an invaluable opportunity and our motivation to provide novel knowledge in order to pre-screen, design and develop efficient DES solvents for LCB dissolution or pre-treatment.
The objective of the present work was to describe the pretreatment or disruption of a cellulose I crystallite by four deep eutectic solvents (DES): cholinechloride ethylene glycol, choline-chloride oxalic acid, choline-chloride urea, and choline-chloride levulinic acid, from the atomistic interactions observed in molecular dynamics and ab-initio studies. The cellulose I disruption was studied considering plausible correlations between the Kamlet-Taft (α and β) solvent parameters, and a series of thermodynamic, structural, and energetic properties. Naturally, the HB formation is the key concept in this work. Besides, the conjunction of Kamlet-Taft solvent parameters was included in the correlation studies.
The present study comprised two stages: (i) the DESs thermodynamic properties calculations, their changes after the addition of a cellulose I crystallite (cellulose-DES mixture), and their further relation with the Kamlet-Taft α and β parameters values of the DESs, i.e. the DES solvent behavior. (ii) A description on the ability of DESs to disrupt the cellulose I crystallite and its relationship with the Kamlet-Taft α and β parameters of the DESs, i.e. the cellulose crystallite disruption from an atomistic perspective. For such purposes, density ( ), molar volume ( V m ), molar volume of HBD ( V m.HBD ), enthalpy of vaporization ( ΔH vap ), cohesive energy density (c), Hildebrand parameter ( ), radial distribution functions (RDFs), cellulose-DES HBs, HB occupancies, non-bonded potential energy of cellulose crystallite ( E non−bonded,cellulose ), and short-range non-bonded pair energies ( E SR ) have been calculated for ChCl/ethylene glycol, ChCl/oxalic acid, ChCl/urea, and ChCl/levulinic acid DESs, as well as for their corresponding cellulose I crystallite-DES mixture. Lv et al. (2012) found that the addition of cellulose in ILs changes the pure thermodynamic properties of the solvent. Accordingly, thermodynamic properties of the cellulose-DES mixtures were also evaluated.
All the MD simulations were conducted at room temperature since desired green and sustainable applications require low energy consumptions (Gunny et al. 2015;Welton 2015;Lu et al. 2016;Clarke et al. 2018;Körner et al. 2019).
To the best of our knowledge, there is scarce literature information concerning the interactions between a cellulose crystallite I model and ChClbased DESs, as well as thermodynamic, atomistic, and structural properties of DESs and cellulose-DES mixtures regarding the cellulose pretreatment.

Force fields
All the ChCl-based DESs were modeled with the force fields (FFs) developed by Doherty and Acevedo (2018), also-called OPLS-DES FFs. These FFs use a ± 0.8e scaled charge for the cations and anions present in the DES, providing a remarkable agreement between calculated and experimental properties including densities, viscosities, heat capacities, surface tensions, and diffusion coefficients (Doherty and Acevedo 2018). The corresponding HBA:HBD molar ratios for the evaluated ChCl/ethylene glycol, ChCl/oxalic acid, ChCl/urea, and ChCl/levulinic acid DESs in this work, were 1:2, 1:1, 1:2, and 1:2, respectively, as is also indicated in Table S1 in the Supplementary Information section and shown in Fig. 1a-d. On the other hand, the cellulose crystallite ( Fig. 1e-f) was modeled with the OPLS 2005 FF (Banks et al. 2005), which has demonstrated to be a reliable FF for the modelling of carbohydrates, (Stortz et al. 2009;Rabideau et al. 2013Rabideau et al. , 2014Rabideau and Ismail 2015;Velioglu et al. 2014;Kubicki et al. 2018;Stalker et al. 2019) and recently validated in our previous work (Sánchez-Badillo et al. 2021). The geometric combining rule was applied for all mixed ij and ij interactions along with a scaling factor of 0.5 for the 1-4 Lennard-Jones and the 1-4 electrostatic interactions.

Molecular simulations details
The MD simulations were carried out using the GROMACS (2020.4) package (Abraham et al. 2015).
The v-rescale (Bussi et al. 2007) thermostat and Berendsen (Berendsen et al. 1984) barostat were used to keep constant the temperature and the pressure at 300 K and 1 bar, respectively. The v-rescale and Berendsen algorithms have been successfully applied in the literature for the calculation of thermodynamic and transport properties of either DES or ILs through MD simulations (Doherty and Acevedo 2018;Doherty et al. 2018). A coupling time of 1 ps was used for both temperature and pressure. Newton's equations of motion were integrated with a time step ( Δt ) value of 1 fs by using the leap-frog algorithm. Periodic boundary conditions were used along with a cut-off value of 16 Ȧ and 30 Ȧ for the calculation of the non-bonded interactions for condensed and gas phase systems, respectively. For the condensedphase systems, the long-range non-bonded interactions beyond the cut-off were calculated by using the smooth Particle-Mesh-Ewald (PME) (Essmann et al. 1995). On the contrary, for the gas phase systems, the cut-off method was used instead. The validation of the cut-off method for gas phase was evaluated against the PME method, these results are indicated in Table S2 in the Supplementary Information section. The LINCS algorithm (Hess et al. 1997) was employed to constraint all the bonds with hydrogen atoms. Equilibration and production runs were carried out for all the systems as indicated below.
Quantum chemical calculations were also carried out in order to reveal the electron density magnitude related to the HBs between cellulose unit cell and DES molecules by using the graphics processing units for atoms and molecules (GPUAM) package  Atomic labels and molecular structures for the deep eutectic solvents and cellulose-crystallite models. a Choline-chloride ethylene glycol (ChCl/ethylene glycol) with 1:2 HBA:HBD molar ratio. b Choline-chloride oxalic acid (ChCl/oxalic acid, 1:1). c Choline-chloride urea (ChCl/urea, 1:2). d Choline-chloride levulinic acid (ChCl/levulinic acid, 1:2). e Glucan atom labeling, including the commonly used atom numbering (in parenthesis); the OD atom label refers to the oxygen atom in the hydroxyl group, the OC atom label includes the oxygen atom within the six-member ring, the OE atom label refers to the oxygen atom in glycosidic bond, and the HS atom label denotes the hydrogen atom attached to OD atom. f Cellulose crystallite composed by 39 glucan of 12 glucose units (GLU residues) each as indicated, arranged in seven layers; left: projection onto the a-b base plane, right: isometric projection. Glucan in these pictures have been green-scale colored for clarity. Color/atom code: blue/nitrogen, green/chlorine, grey/carbon, red/oxygen, and white/hydrogen. All images have been created with the visual molecular dynamics (VMD) software (Humphrey et al. 1996) electron density, Laplacian of electron density, electrostatic potential, non-covalent interactions (NCI), among others. The versatility of the GPUAM package relies on its ability to associate one thread of the GPU unit with one point of the selected grid to evaluate the requested field . This feature allows the evaluation of several electronic properties of large systems in reasonable time (Cruz et al. 2019).
The NCI from quantum chemical calculations were evaluated for the cellulose unit cell and DES molecules systems. For such purpose, the cellulose unit cell in conjunction with the nearest DES molecules were extracted from the cellulose-DES systems. Extracted molecules were optimized in the TeraChem software (Ufimtsev and Martínez 2008). Then, the NCI were obtained in the GPUAM package. Finally, the (r) magnitude and geometric conformation of the NCI (including HBs) were extracted from the GPUAM log output file.
Thermodynamic and structural properties for pure DES and cellulose-DES systems Cubic boxes containing the corresponding ChCl and HBD pairs as indicated in Table S1, were constructed using the Packmol package (Martínez et al. 2009). A minimization stage was performed by using the steepest-descent method for 10,000 steps in order to remove any possible steric clash. Then, the DESs systems were equilibrated for 10 ns in the NPT ensemble at 300 K and 1 bar, followed by a production run of 40 ns at same conditions. The cellulose I crystallite-DESs systems were constructed as follows: a cellulose I crystallite comprising 39 glucan with a degree of polymerization (DP) of 12 each, see Fig. 1f, was built using the cellulose-builder toolkit (Gomes and Skaf 2012). The constructed model implies the current crystallographic standard for cellulose from Nishiyama et al. (2002). The cellulose crystallite was placed in the center of a cubic box and then surrounded by the same number of DES molecules as in their corresponding pure systems, see Table S1. A minimization stage of 10,000 steps followed by 100 ns of equilibration, and by a production run of 50 ns were carried out for all cellulose-DES systems in the NPT ensemble at the same conditions as the pure state of DESs. Structural properties for both DES and cellulose-DES systems were calculated from the production stages by using the TRAVIS software (Brehm et al. 2020).

Enthalpy of vaporization (ΔH vap )
The enthalpy of vaporization for both DES and cellulose-DES systems was calculated by using the proposed methodology of Salehi et al. (2019) as follows: With N as the total number of molecules in the condensed phase. In Eq. 1, the E liq (N) and E liq (N − 1) correspond to the total energy of the condensed phase before and after removing one molecule of the entity that evaporates, ΔE vap (T, P) , while R and T are the universal gas constant and temperature, respectively. Despite there is a lack of evidence regarding the composition of the gas phase for the studied DESs in this work (Ferreira et al. 2016;Salehi et al. 2019), results from MD calculations (Ferreira et al. 2016;Salehi et al. 2019Salehi et al. , 2021 suggest that DES gas phase is largely composed by HBDs. In order to evaluate the ability of the selected OPLS-DES FF to reproduce such behavior, i.e. the vapor-liquid equilibrium, molecular simulations at 300, 400, and 500 K for the ChCl/ethylene glycol system were carried out in the NVT ensemble. Detailed methodology and results can be found in the Supplementary Information section. From Fig S2, it can be observed that gas phase is dominated by the ethylene glycol molecules without presence of ChCl molecules at 400 and 500 K, in agreement with the recent finding of Salehi et al. (2021). Accordingly, the E gas (1) required in Eq. 1, can be computed as the total energy of one HBD molecule in the ideal gas phase at the same temperature than corresponding condensed phase. The E gas (1) was evaluated as follows: the HBD was placed in the center of a cubic box of length 100 Ȧ to avoid interactions with its periodic images. Then, a run of 1000 steps for a minimization stage was developed, followed by an NVT equilibration run of 10 ns at 300 K, and by a NVT production run of 15 ns at 300 K. The E liq (N) values for both DES and cellulose-DES systems, were obtained from their corresponding production runs. The output coordinates from the E liq (N) simulation, were modified in order to eliminate one HBD molecule, and then used as input coordinates for the E liq (N − 1) calculation. The N − 1 DESs were equilibrated for 25 ns followed by 25 ns of a production period in the NPT ensemble at 300 K. In the case of the N − 1 cellulose-DES systems, they were equilibrated for 10 ns followed by a production of 10 ns length in the NPT ensemble to avoid further deconstruction of the cellulose crystallite.
Hildebrand solubility parameter (δ) The solubility of one material in another is facilitated when their non-bond attractive forces are similar. Therefore, the dissolution capacity of a solvent regarding to a specific solute can be predicted if solvent ∼ solute (Rinaldi and Reece 2017). To explore this behavior in the pretreatment of cellulose by DESs, the parameters for the DESs were calculated through of the following equations (Salehi et al. 2019): Here c is the cohesive energy density and V m,HDB is the molar volume of vaporizing entity in the condensed phase (HBD). The V m,HDB was obtained by using the volume and the number of HBD molecules in the simulation box before one molecule of HBD vaporizes.

NCI-AIM analysis for HBs between cellulose and DES molecules through DFT calculations
The cellulose unit cell, composed by five cellobiose (two glucose molecules linked by the β(1→4') glycosidic bond) molecules was extracted from the center top and from the edge of the cellulose crystallite along with the DES molecules within 8 Ȧ from the unit cell. The missing hydroxyl groups in the C − OH bonds in the glucose units were included by using GaussView6 package (Dennington et al. 2016). The molecular configuration of the system (containing around 350 atoms) was minimized in the Tera-Chem software (Ufimtsev and Martínez 2008) at the PBE0 (Perdew et al. 1996(Perdew et al. , 1997; Adamo and Barone 1999)/6-311 g(d,p) level of theory with empirical D3 dispersion corrections (Grimme et al. 2010), maintaining the heavy atoms fixed, in order to directionally optimize all the HBs. After minimization, a single point energy calculation was evaluated at the same level of theory in TeraChem. The generated molden output file was then used to generate the electronic wavefunction through MOLDEN2AIM software (Zou 2021). The NCI were obtained through the analysis of the bond critical points of the electronic density (ρ(r)) from wavefunction in the GPUAM package (

DESs solvent behavior
In the search for new factors that may influence the pretreatment of the cellulose and provide a reliable pre-screening of the effectiveness of DESs pretreatment, a series of correlations between thermodynamic properties and the Kamlet-Taft parameters were investigated. Experimentally, the Kamlet-Taft α and β parameters have been correlated with the solubility degree (Ren et al. 2016a(Ren et al. , 2016bTong et al. 2021), extraction rate (Xia et al. 2018) and the pretreatment (Liang et al. 2021  Kamlet-Taft parameters with a deviation lower than 3.14% compared to the available experimental values. For this reason, the set of Kamlet-Taft parameters values from ab-initio calculations by Kundu et al. (2020) were included in this work (Table 1).

Density and molar volume
The calculated density ( sim ) and molar volume ( V m,sim ) for pure DES and cellulose-DES systems are shown in Table 1. It can be noticed that calculated density of the DESs ( sim,DES ) agree with the available experimental density values, presenting an overall average deviation around 2.0%. The addition of the cellulose crystallite ( The effect of DES density in the cellulose pretreatment was explored through the relationship between calculated pure DES density and the Kamlet-Taft β parameter, as is shown in Fig. 2a. From Fig. 2a, a notable correlation emerged for the relationship between DESs density values and β values, presenting a squared correlation coefficient ( R 2 ) of 0.91. Similarly, from Fig. S4a the relation between DESs density values and α values also showed an important correlation ( R 2 = 0.84 ). Experimental evidence indicates that larger values of α and β parameters are related to larger LCB dissolution (Ren et al. 2016a;Liang et al. 2021). On this regard, the correlations observed in Fig. 2a and S4a suggest that low-density DESs such as ChCl/ethylene glycol and ChCl/levulinic acid could provide favorable environment for cellulose dissolution or pretreatment than high density DESs such as ChCl/oxalic acid. In the literature, the relation between DES density and LCB dissolution is divided. Jablonsky et al. (2019) have reported that DESs with low densities favor the cellulose dissolution. On the other hand, Liang et al. (2021) have indicated that DESs with high densities enhanced the LCB treatment. Results from the present work agrees with the conclusions reported by Jablonsky et al. (2019).
Häkkinen and Abbott (2019) indicated that the dissolution of cellulose is favored in some ILs due to their crystalline structure, since the addition of cellulose produces an increase in the entropy of the system. Therefore, larger entropy changes could conduct to large changes in densities from solvent upon the addition of solute. On this basis, the calculated density variations between the pure DES and the cellulose-DES mixture ( Δ = sim,cellulose−DES − sim,DES ) are indicated in Table 1, and their relationship with the Kamlet-Taft β parameter is shown in Fig. 2b. It was observed that the larger Δ values (~ 0.053 g/ cm 3 ) correspond to the low-density DESs, including Similarly, to the previous relations, the molar volume of the system ( V m,sim,DES ) and its change were also evaluated. As observed from Fig. 2cd, the larger V m values are associated to the larger ΔV m ( V m,sim,cellulose−DES − V m,sim,DES ) values, as indicated in Table 1. However, both V m and ΔV m displayed slight correlations with the Kamlet-Taft β values ( R 2 = 0.66 , R 2 = 0.62 , respectively). Interestingly, the correlations between the V m and ΔV m and the Kamlet-Taft α values were slightly better correlated ( R 2 = 0.75 , R 2 = 0.71 , respectively), as shown in Fig. S4c-d. The smallest ΔV m (and V m ) values corresponded to the ChCl/oxalic acid DES for both cases. This behavior could suggest that the small variations in the V m of the system could produce repulsive interactions between solute and solvent that hinder the cellulose pretreatment in the DES bulk.
The molar volume of the HBD ( V m,sim,HBD ) was also calculated for the DES and cellulose-DES systems and reported in Table 1 (indicted in parenthesis). Unfortunately, there was no correlation between such values and the Kamlet-Taft parameters for the evaluated solvents ( Fig. 2e and Fig. S4e) contrary to the V m,sim,.DES behavior ( Fig. 2c and Fig. S4c). This could suggest that both HBD and HBA play a relevant role in the pretreatment or dissolution of cellulose.
Results from Fig. 2a-d and Fig. S4a-d, indicate that bulk properties such as density and molar volume of the system can be used as effective parameters to identify or develop new and sustainable solvents for cellulose pretreatment. Particularly, the DES density has been reported as a thermodynamic property that can reflect the potential applicability of new DESs for LCB pretreatment or dissolution (Jablonsky et al. 2019).

Enthalpy of vaporization
Results for the energetic thermodynamic properties of the bulk, including ΔH vap , c , and δ obtained for both DES and cellulose-DES systems through Eqs. 1-3, are displayed in Table 2. The calculated values of the enthalpy of vaporization ( ΔH vap,sim ), through the methodology proposed by Salehi et al. (2019), agree with both available data sets: MD estimated ( ΔH vap,sim,lit ) (Salehi et al. 2019 and experimental ( ΔH vap,exp ) (Shahbaz et al. 2016;Ravula et al. 2019).
The relationship between the ΔH vap,sim and the β and α parameters revealed low correlation coefficient values, R 2 = 0.50 , and R 2 = 0.59 respectively, as displayed in Fig. 2f and Fig. S4f. Interestingly, the addition of the cellulose crystallite into the bulk of the DESs produced positive enthalpy of vaporization changes ( ΔΔH vap~4 .85 kcal/mol), as consequence of the E pot,liq reduction (repulsion interactions) of the mixture compared to the pure system, as it can be noticed in Table 2. The increase of the ΔH vap in the cellulose-solvent mixture was also observed by Liu et al. (2010) for a series of cellulose and 1-ethyl-3-methylimidazolium acetate ((C 2 mim][OAc]) IL mixtures through MD simulations. Maximo et al. (2010) reported the elevation of the normal boiling point as consequence of the increase in glucose concentration in water, which leads to an increment in the necessary heat for the vaporization process, in agreement with the obtained results in this work.
Surprisingly, after the cellulose crystallite addition, the correlation between the ΔΔH vap and the Kamlet-Taft β values raised up to R 2 = 0.74 (Fig. 2g), and R 2 = 0.83 for the Kamlet-Taft α parameter (Fig.  S4g). This is probably caused by the attractive interactions between cellulose crystallite and DES molecules that contribute to the total energy required to reach the vaporization of the system. In other words, strong solute-solvent interactions (described by larger Kamlet-Taft parameters values) result in larger values of ΔΔH vap , whereas weak or repulsive solute-solvent interactions (related to small Kamlet-Taft parameter values) favors the gas phase formation due to a small ΔΔH vap change. Accordingly, the larger the Kamlet-Taft α and β parameters the larger ΔΔH vap for the evaluated DESs.
Cohesive energy density Tang et al. (2017) indicated that the reduced activity of the Cl − ion within the ChCl-based DESs was probably related to the restricted HB network compared to 1-butyl-3-methylimidazolium chloride ([C 4 mim][Cl]) IL. The aforementioned despite that the Kamlet-Taft β parameter (proton accepting) of the [C 4 mim][Cl] IL is larger (β = 0.94 (Lungwitz and Spange 2008)) than the ChCl/levulinic acid DES (β = 0.61 (Kundu et al. 2020)). Thus, it seems reasonable that low c values might provide a favorable media to dissolve or pretreat the cellulose. The calculated values for the c of both DES and cellulose-DES systems are indicated in Table 2, and the corresponding relationships between calculated c of the pure DES ( c sim,DES ) and the Kamlet-Taft β and α values are shown in Fig. 2h and Fig. S4h, respectively. It was not found a correlation between cohesiveness and the Kamlet-Taft parameters ( R 2 = 0.00 and 0.02 for β and α, respectively). On the contrary, the relationship between Δc values and Kamlet-Taft (β or α) parameters, produce a feasible correlation ( R 2 = 0.84 for β in Fig. 2i and R 2 = 0.90 for α in Fig. S4i). This was to be expected since Δc involved the reduction of the E pot,liq (shown in Table 2) in the cellulose-DES mixture. From Fig. 2i and Fig. S4i, it can be noted that large Δc values corresponded to the DESs with larger both Kamlet-Taft β and α values. This suggests the formation of attractive interactions between solute and solvent emerge at same time the solvent maintains its own interactions. This is in accordance with the findings of Häkkinen and Abbott (2019) who evaluated the interactions between glucose and ChCl/ethylene glycol DES. They found that such DES was able to maintain its interactions with itself at the same time the ChCl and ethylene glycol HBD form new HBs with glucose, increasing the cohesiveness of the mixture, and thus the augmented energy to vaporize it ( Fig. 2g and Fig. S4g).

Hildebrand parameter
The Hildebrand parameter ( ) is a key property to select solvents for separation process (Salehi et al. 2019). As it can be noted from Table 2 Rinaldi and Reece (2017) have summarized the available Hildebrand parameter data for cellulose ( cellulose ). Unfortunately, the cellulose experimental values presented large variation between them, since they range from 25.7 up to 39.9 MPa 0.5 , therefore a direct comparison with sim,DES could provide unclear conclusions. Despite the evaluated DESs are not able to dissolve the cellulose crystallite under the evaluated temperature (Morais et al. 2020;Sosa et al. 2020;Alizadeh and Kirchner 2021), it could suggest a possible relationship between sim and Kamlet-Taft parameters to describe the pretreatment or the solvation power of DESs, as it was pointed out by Mäki-Arvela et al. (2010) to correlate cellulose dissolution in ILs through Hildebrand parameters and recently by Lang et al. (2021) to correlate cellulose dissolution in DESs through Hansen parameters. The relationship between the sim,DES and Kamlet-Taft parameters is shown in Fig. 2j for β parameter and in Fig. S4j for α parameter. Unfortunately, the sim,DES values did not correlate with the Kamlet-Taft parameters. It seems that the cellulose pretreatment ability of the DESs, described by the Kamlet-Taft parameters, can be correlated directly with pure both density and molar volume of the DESs (Fig. 2 and Fig. S4), providing a valuable solvent pre-screen for LCB pretreament. However, there are no direct correlation with the remaining properties. In order to investigate the origin of such correlations, atomistic structure and energetics involving in the disruption of the cellulose crystallite were also investigated, as well as their correlations with the Kamlet-Taft parameter values.

Cellulose crystallite disruption
The final configurations of the cellulose crystallite in the evaluated DESs in this work are shown in Fig. 3. From the a-b base plane projection, a slightly swelling of the cellulose crystallite in the ChCl/ethylene glycol DES can be observed (Fig. 3a). Additionally, a reduced swelling effect is also noticed in the ChCl/ urea and ChCl/levulinic acid DESs in Fig. 3c and Fig. 3d, respectively. On the other hand, an apparent null swelling was observed in the ChCl/oxalic acid DES (Fig. 3b). This swelling behavior of the cellulose crystallite has been recently reported by Alizadeh and Kirchner (2021) in urea-based DESs through MD simulations. Accordingly, DESs could perturb the crystallinity of the cellulose even at low temperature. In order to identify the relevant interactions that cause this behavior as well as that observed in the analysis of thermodynamic properties in the previous section, structural analyzes through radial distribution functions (RDFs) were carried out.

Cellulose-DES radial distribution functions
Briefly, the RDF indicates the probability g(r) that one particle finds another particle as a function of distance, r . The calculated cellulose-DES RDFs involved the interactions between OD, OC, and HS atoms from cellulose crystallite and the atoms from the DES: H CHOL (aliphatic hydrogen atoms in the choline cation), HY (hydrogen atom in the hydroxyl group in choline cation), HO (hydrogen atom in the hydroxyl group of HBD), H HBD (aliphatic hydrogen atoms hydrogen atoms in HBD as well as hydrogen atoms from urea), Cl (chloride anion), OY (oxygen atom in the hydroxyl group in choline cation), OH (oxygen atom in the hydroxyl group in HBD), and O HBD (oxygen atoms from HBD, except those in the hydroxyl group), are shown in Fig. 4, Fig. S5, Fig. S6,   Fig. 3 The a-b base plane and the isometric projection of the cellulose crystallite at the end of the NPT simulation in: a ChCl/ethylene glycol. b ChCl/oxalic acid. c ChCl/urea. d ChCl/levulinic acid (DES molecules have been omitted). Glucan in these pictures have been green-scale colored for clarity. All images have been created with the VMD software (Humphrey et al. 1996) Table 3 First peak distance position (A) and corresponding coordination number (Ncoord, in parenthesis) for pair interactions between cellulose and DESs, and between HBA and HBD in pure DES and cellulose-DES systemsa a The "-" indicates no defined peaks were found for such interaction.

Interaction
Cellulose-ChCl/ethylene glycol and Fig. S7 for the cellulose-ChCl/ethylene glycol, cellulose-ChCl/oxalic acid, cellulose-ChCl/urea, and cellulose-ChCl/levulinic acid systems, respectively. For all the interactions, the first peak position and its corresponding total coordination (coordination number, N coord ), calculated by integrating the g(r) until the first minimum, have been summarized in Table 3. From Table 3 and Fig. 4a, S5a, S6a, and S7a the OD − HY first peak is larger or better positioned than OD − H CHOL interaction for all the systems. From the RDF analysis, a strong HB is generally related to a short interaction distance with a large and welldefined g(r) peak (Doherty et al. 2018). Accordingly, the OD − HY interactions are stronger and more relevant than the OD − H CHOL interactions in all the systems.
The shortest OD − H CHOL and OD − HY interactions were found in the cellulose-ChCl/oxalic system (Table 3). An important feature from the RDF analysis is the number of molecules or atoms surrounding the selected molecules or atoms in the first solvation shell, given by the coordination number. Table 3 reveals that the OD atom from cellulose in ChCl/ oxalic system is surrounded by the larger number of hydrogen atoms from choline cation than other systems, indicating that such interactions play a relevant role in that solvent. This behavior could hinder other relevant and stronger cellulose-DES interactions such as those involving Cl − and the HBD, in agreement with the smallest Kamlet-Taft α and β values (1.1921 and 0.2707, respectively) (Kundu et al. 2020). Interestingly, the remaining cellulose-DES systems (which also have larger Kamlet-Taft β values) showed weaker OD − H CHOL and OD − HY interactions. These results may explain the ΔΔH vap as well as the change in the cohesiveness of the cellulose-DES mixture reported in Table 2 and Fig. 2i.
The RDFs involving the OC oxygen atom and aliphatic hydrogen atoms from choline ( OC − H CHOL ), continuous black lines in Fig. 4a, Fig. S5a, S6a, S7a, shows no-defined peaks for most of the systems. The OC − HY interactions (dotted black lines in same figures) show similar behavior than OD − HY with smaller g(r) probability values, as it was reported by Alizadeh and Kirchner (2021) for the cellulose in ChCl/urea DES. This is probably caused by a steric effect surrounded the OC oxygen atom in the glucose ring (see Fig. 1e), thus playing a minor role in the LCB pretreatment.
Interactions involving the oxygen atoms from cellulose and hydrogen atoms from HBD are shown in Fig. 4b, S5b, S6b, and S7b. As it can be noted, the OD − HO interactions (continuous blue lines) show the shortest interaction distances (Table 3) and higher g(r) peaks than those related to the H CHOL ( OD − H CHOL ) and HY ( OD − HY ) hydrogen atoms. The interactions involving the OD − H HBD ( H HBD involves aliphatic hydrogen atoms from HBD as well as hydrogen atoms from urea HBD) were also evaluated, however they presented larger interaction distances than OD − HO interplay. This ability to form stronger hydrogen bonds from the acidic OH and HY hydrogen atoms of the HBD towards the OD atoms from cellulose could be explained from the intrinsic pure DES interactions. On this regard, small Cl − HO and Cl − HY interaction distances (strong HBs) within the pure DES might hinder the formation of HBs with solute Zhang et al. 2020). This effect can be found in the ChCl/oxalic acid DES, presenting the shortest both Cl − HO (2.017 Ȧ ) and Cl − HY (2.15 Ȧ ) distances among the evaluated DESs as well as larger coordination numbers for such interactions (0.734 and 0.831, respectively). Accordingly, large interaction distances between cellulose and HO and HY hydrogen atoms (weak HBs) can be found for this system. In other words, small Cl − HY and Cl − HO distance in conjunction with the larger N coord in pure DESs, could be an indicative of the HB network strength within the DES as well as their limiting role to deconstruct the cellulose crystallite.
The OC − HO interactions presented larger interaction distances than the OD − HO with small g(r) values (continuous lines in Fig. 4b, S5b, S6b, and S7b). Considering that the OC oxygen atom from cellulose is covered by steric repulsions, the acidic (HO) hydrogen atoms from solvent need to promote the deconstruction towards the OC atom in cellulose. On this regard, the OC − HO interaction becomes stronger in the following order: ChCl/oxalic acid < ChCl/ethylene glycol < ChCl/levulinic acid, because its corresponding distance decreases as 2.01 > 1.95 > 1.81 Ȧ . Interestingly, this behavior agrees with the corresponding Kamlet-Taft α values (1.19, 1.47, and 1.55) (Kundu et al. 2020).
From Table 3 and Fig. 4c, S5c, S6c, and S7c, the HS − Cl interaction was located about 2.02 Ȧ for all cellulose-DES systems. This interaction was the most prominent and thus frequent (compared to HS − OY ) due to its large g(r) value, in agreement with the findings of Häkkinen and Abbott (2019), Smirnov et al. (2020), and Alizadeh and Kirchner (2021). However, other interactions involving the HS hydrogen atom are positioned at shorter distances. For example, the HS − OH in the cellulose-ChCl/ethylene glycol system (Fig. 4d), is located at 1.83 Ȧ , and the HS − O HBD in the cellulose-ChCl/urea system at 1.97 Ȧ (Fig.  S6d), as it was reported also in the literature (Smirnov et al. 2020). Analizing the interactions involving the chloride atom, it can be noted that the first peak for the Cl − HY interaction decreases its length after the cellulose-crystallite addition, for the most of the evaluated DESs excepting the ChCl/oxalic acid DES, which increased its length (Table 3). A similar behavior can be found for the Cl − HO distances. The reduction of the Cl − HY and Cl − HO distances might explain the arise of the cohesiveness in mixture, providing the large calculated Δc values for the ChCl/ethylene glycol, ChCl/urea, and ChCl/levulinic acid DESs shown in Fig. 2i. Meanwhile, the smallest Δc value for the ChCl/oxalic acid DES can be associated to the increment (weakening) of the Cl − HY and the Cl − HO distances in mixture.
Finally, the interactions between HS and oxygen atoms from HBD are displayed in Fig. 4d, S5d, S6d, and S7d. As it can be observed, the stronger interactions are found in the cellulose-ChCl/ethylene glycol and cellulose-ChCl/urea systems. Interestingly, the hydroxyl oxygen atom (OH) in the -COOH moiety present in ChCl/oxalic acid (Fig. S5d) and ChCl/ levulinic acid (Fig. S7d) DESs, formed the weakest interactions with the HS atom, compared to the HS − O HBD HB found in the cellulose-ChCl/oxalic acid system (Fig. S5d), cellulose-ChCl/urea (Fig.  S6d), and in the cellulose-ChCl/levulinic acid (Fig.  S7d) systems. This suggests that stronger HS-oxygen interactions between cellulose and HBD, might be enhanced for R-CO-R structures than -COOR moieties within the HBD. In order to validate such an assumption, a larger number of HBD must be evaluated.
The largest g(r) values for HS − Cl (hydrogen bond basicity) and the relevant OD − HO (hydrogen bond acidity) interactions, indicated that such HBs might play a significant role in the decrystallization of the cellulose crystallite. This process is related to the weakening or breaking of the inter-and intrachain HBs within the solute (Procentese et al. 2015;van Osch et al. 2017;. In order to explore this effect, the total (inter-and intrachain) HBs have been analyzed for all cellulose-DES systems.

HBs within the cellulose crystallite
The disruption of the inter-and intrachain HBs must be a key step towards the understanding of the cellulose pretreatment with DESs (Procentese et al. 2015). On this regard, the total interchain and intrachain number of HBs within the cellulose crystallite has been plotted as function of the last 50 ns of simulation (production time) and displayed in Fig. S8. The number of HBs between the cellulose crystallite and choline cation ( HBs cellulose−choline ), chloride anion ( HBs cellulose−chloride ), and HBD ( HBs cellulose−HBD ), was also calculated and displayed in Table 4. All the HBs in this analysis have been calculated through the Hydrogen Bonds plugin available in the VMD program (Humphrey et al. 1996), which involves a donor-acceptor (D-A) distance less than 3.5 Ȧ and an acceptor-donor-hydrogen (A-D-H) angle less than 30° (Li et al. 2015). From Fig. 5a-b and Fig. S9a-b, there were not found correlations for the intrachain or interchain HBs against Kamlet-Taft α and β parameters. However, from Fig. S8a and Fig. S8d, a reduction in the number of intrachain HBs can be observed for the ChCl/ethylene glycol and ChCl/levulinic acid DESs, contrary to the behavior shown in ChCl/oxalic acid and even ChCl/urea DESs (Fig. S8b and Fig. S8c). On the other hand, the number of the interchain HBs remained almost constant for all cellulose-DES systems along the production time. Furthermore, the interchain HBs number was almost similar to each other as indicated in the Table 4, suggesting that the cellulose crystallite deconstruction (pretreatment) process starts by reducing the intrachain HBs by action of the DES. The slight decrease in the number of intrachain HBs in ChCl/ethylene glycol and ChCl/levulinic acid DESs along time, suggested that the rate of change (the slope in Fig. S8) could provide useful information. We calculated this slope ( r intrachainHBs ), obtaining -0.403, -0.216, -0.277, and -0.514 intrachain HBs/ns for cellulose-ChCl/ethylene glycol, cellulose-ChCl/oxalic acid, cellulose-ChCl/urea, and cellulose-ChCl/levulinic acid systems, respectively. The calculated negative slope indicates the decrystallization of the cellulose crystallite in real time, which is larger for ChCl/levulinic acid DES. Surprisingly, a meaningful correlation between these values and the Kamlet-Taft β and α values is displayed in Fig. 5c and Fig.  S9c. The HBs stability was also elucidated through their occupancies, to further understand such decrystallization process.

HB occupancies within the cellulose crystallite
The HB occupancy indicates the percentage of occurrence (or fraction, and thus the stability) of a particular HB along the simulation. Intrachain and interchain HB occupancies were calculated in two different sections of the cellulose crystallite: chain 0 (located at the edge of the cellulose crystallite) and chain 17 (located at the middle of the cellulose crystallite), for all systems (Fig. S10 to S13).
As it can be noted from Fig. S10a-b, S11a-b, S12ab and S13a-b, the intrachain HB occupancies present in both the chain 17 and chain 0 were related to the typical OD − HS ⋯ OD and OD − HS ⋯ OC HBs, specifically the O2 − H2 ⋯ O6 , O3 − H3 ⋯ O5 , O6 − H6 ⋯ O3 , and O6 − H6 ⋯ O2 HBs, as reported for cellulose I (Shen and Gnanakaran 2009). Despite the chain 0 was surrounded by only two glucan chains, larger intrachain HB occupancies were absent or less frequent than in the chain 17. This suggests also that the cellulose disruption process started by diminishing the HB occupancies, and weakening the chains located at the edge of the cellulose crystallite by action of the DES molecules.
From the obtained interchain HB occupancies in chain 17 and chain 0, it can be noted that the strong interchain HBs for the glucan at chain 0 in contact with the solvent, were replaced by the O6 − H6 ⋯ O4 HBs, i.e., HBs involving the oxygen atom in the glycosidic bond (OE). The emerging occurrence of the O6 − H6 ⋯ O4 HBs by action of the DESs, promotes the releasing of the O2 and O3 oxygen atoms. The aforementioned favors the functionalization of the hydroxyl group in the synthesis of cellulose derivatives, for example, as required in the acetylation of wood (Sun et al. 2019). This effect suggest the use of the DES as potential sustainable media to synthetize cellulose derivatives due to promotion of the interactions between cellulose and third incorming particles. In order to explore the interplay between cellulose and DES molecules, HBs between them were also evaluated.

HBs between cellulose crystallite and DES molecules
From the RDF analysis, it was found that the HBs between the cellulose crystallite and the choline cation may hinder the role of the HBD and chloride anion towards the cellulose disruption. This effect can be observed in Fig. 5d and Fig S9d, since the lower the Kamlet-Taft β (or α) parameter the larger the number of HBs cellulose−choline . The large number of HBs cellulose−choline was found in the ChCl/oxalic acid DES ( β = 0.27, = 1.19 ), as displayed in Table 4, and could be related to the strong OD − H CHOL and OD − HY correlations as well as the large coordination number for the OC − HY interaction, indicated in Table 2. The HBs cellulose−choline and the RDF analysis suggest that interaction between cellulose and the oxalic acid HBD, is limited, however, it is able to replace strong interchain HBs to the O6 − H6 ⋯ O4 HBs, as shown in Fig. S11d. On the other hand, the small number of HBs cellulose−choline in the remaining solvents (Table 3) might be a consequence of the largest Cl − HY , Cl − HO , Cl − H HBD distances (Table 2). Accordingly, choline cation and HBD will compete to form interactions with cellulose crystallite. Nevertheless, as RDF indicates, the interaction between cellulose crystallite and HBD are stronger than those with choline cation, and thus larger number of HBs cellulose−HBD could be found, as shown in Table 4. However, there was no correlation between HBs cellulose−HBD and the Kamlet-Taft parameters of the evaluated ChCl-based DESs ( R 2 = 0.13 in Fig. 5f for β parameter and R 2 = 0.08 in Fig. S9f for α parameter).
The HBs cellulose−chloride shown in Fig. 5e and Fig.  S9e, indicate that the largest number of HBs was found for the cellulose-ChCl/ethylene glycol system (161 ± 8), in concordance with the largest distance found for the Cl − HY HB in pure DES (i.e., small interactions between anion and cation). The smallest number of HBs cellulose−chloride (and small number of HBs cellulose−HBD ) was found in the cellulose-ChCl/ levulinic acid system. This agrees with the stronger Cl − HO HBs before and after addition of the cellulose crystallite, as indicated in Table 3, which might hinder better interactions from the levulinic acid HBD towards the cellulose crystallite, and probably caused the low presence of the O6 − H6 ⋯ O4 HBs (Fig. S13d). It is interesting to note that ChCl/levulinic acid DES seems to be limited by the interactions present between its HBA and HBD.

NCI-AIM analysis for HBs between cellulose and DES molecules
Calculated HBs from AIM analysis are shown in Table 5 for the cellulose unit cell in ChCl/ethylene glycol DES, and in Table S3 and Table S4 for the remaining systems. After hydrogen atoms optimization, it can be noted that the number of interchain HBs for a single cellobiose chain does no decrease, except for ChCl/ethylene glycol DES, regardless the cellulose unit cell location (top or edge). In contrast, the number of intrachain HBs for the same cellobiose chain decreased for the unit cell at the edge of the cellulose crystallite in all systems. The aforementioned is in agreement with the observations from classical MD results in the present work.
The solute-solvent interplays are indicated in Table 5, Table S3, and Table S4 for two cellulose unit cell location. From these tables, it can be noticed that the HBs between unit cell at the edge of the cellulose crystallite and DES molecules (Table 5 and  Table S4), shows the largest (r) values compared to those related to the unit cell at the top of the cellulose crystallite (Table 5 and Table S3). Consequently, large deconstruction of the cellulose crystallite could be obtained from the edge. The total number of the HBs (based on (r) criterion) between cellulose crystallite at top (Table 5 and Table S3) and edge (Table 5  and Table S4) positions and DES molecules, displayed in Table 4, did not correlated with the Kamlet-Taft α and β parameters, as shown in Fig. S9g and Fig. 5g, respectively.
Interestingly, results from Table 5, Table S3, and Table S4, suggest a synergetic role of the DES molecules towards the cellulose disruption. This is because the D DES − H DES ⋯ A (D: donor, A: acceptor) HBs (related to the hydrogen bond acidity) play a larger or similar role than the Cl∕O DES ⋯ H − D at the top of the cellulose crystallite; on the contrary, the Cl∕O DES ⋯ H − D HBs (related to the hydrogen bond basicity) are largely present at the edge of the cellulose crystallite. This synergetic effect from the top to the edge on the cellulose crystallite, probably promotes the breaking of the both intra-and interchain HBs in the solute, suggesting that both (α and β) Kamlet-Taft parameters must be considered simultaneously (together). Hauru et al. (2012) suggested an effective dissolution window defined as "net basicity" ( − ) , which compares Kamlet-Taft α and β values, to correlate the dissolution/regeneration of cellulose in water-IL mixtures. Liu et al. (2019) and Hong et al. (2020) have also proposed similar correlations for the extraction of LCB components by using DES. Surprisingly, by calculating the net basicity of the evaluated DESs, a notable correlation arises with previously computed properties (Fig. 6) including: cohesive energy density (Fig. 6b, R 2 = 0.65 ), Hildebrand parameter (Fig. 6c, R 2 = 0.65 ), intrachain HBs (Fig. 6d, R 2 = 0.71 ), HBs cellulose−chloride (Fig. 6e, R 2 = 0.95 ), and the number of HBs between cellulose unit cell and DES molecules from DFT calculations (Fig. 6f, R 2 = 0.88 ), compared from those obtained with the Kamlet-Taft β or α parameters when considered separately. Interestingly, the enthalpy of vaporization showed a modest correlation with the α parameter ( R 2 = 0.59 in Fig. S4f), β parameter ( R 2 = 0.51 in Fig. 2f), and net basicity ( R 2 = 0.51 in Fig. 6a). Regarding the Hildebrand solubility parameter, Fig. 6c indicated that cellulose should be larger than 25 MPa 0.5 since ChCl/ urea and ChCl/ethylene glycol DESs have shown the larger deconstruction of the cellulose crystallite. Table 5 Obtained interactions between the unit cell located at the top and the edge of the cellulose crystallite and DES molecules in ChCl/ethylene glycol, calculated at PBE0/6-311 g(d,p) level of theory a,b,c a All the interactions are displayed as atoms from DES ⋯ atoms from glucose units. Atom labels are indicating in Fig. 1a-e A: acceptor, D: donor b In bold, identified HBs based on the electronic density values from 0.005 to 0.05 a.u. (Johnson et al. 2010) c In bold italics, identified HBs based on the geometric definition (Hassan et al. 2021) Interchain HBs for single chain: 9 6 Intrachain HBs for single chain: 2 1 Interaction Interestingly, this value is in close agreement with the lower limit of the Hildebrand solubility values for cellulose reported by Rinaldi and Reece (2017), which varies from 25.7 up to 39.9 MPa 0.5 . Furthermore, by evaluating the net basicity ( − ) notable correlations can be obtained towards the experimental LCB dissolution or pretreatment for the evaluated DESs (Procentese et al. 2015;Ren et al. 2016a;Sun et al. 2022). On this regard, larger dissolution, pretreatment, or product yield for LCB can be obtained with those DESs with large net basicity. This highlight suggests that the conjunction of Kamlet-Taft parameters is important to understand the role of the DESs in the cellulose pretreatment. Results and observations of the net basicity ( − ) on the deconstruction of the cellulose crystallite, can be elucidated in terms of the cellulose crystallite non-bonded potential energy and short range (SR) van der Waals and coulombic non-bonded pair energies ( E SR ) for all cellulose-DES systems.

Cellulose crystallite potential energy and non-bonded pair energies
The overall effect on the cellulose crystallite decrystallization caused by DESs can be described through the potential energy of the cellulose crystallite ( E non−bonded,cellulose ). The E non−bonded,cellulose was evaluated by computing the single point (molecular mechanics) energy of the final configuration of the cellulose crystallite (shown in Fig. 3). The E non−bonded,cellulose values are displayed in Table 4. As it observed in Fig. 5h and Fig. S9h, there was no correlation between E non−bonded,cellulose and Kamlet-Taft parameters separately. Interestingly, the E non−bonded,cellulose correlates with the net basicity of the DESs ( R 2 = 0.88 ), as shown in Fig. 6g. In order to explore the effect of the individual DES molecules (cation, anion, and HBD) towards the cellulose deconstruction, the E SR was calculated.
The E SR was evaluated among the molecules belonging to the groups CEL (i.e. the cellulose crystallite), CHOL (choline cations), Cl (chloride anions), and HBD in an additional NPT MD simulation of 25 ns for all cellulose-DES systems. Remaining MD simulation details corresponded to those used in the NPT condensed phase production stage, as indicated in the methodology section. It is important to indicate that the E SR corresponded only to the LJ(SR) (van der Waals) and Coul(SR) (electrostatic) terms from the GROMACS energy output file. The pair energy interactions between cellulose crystallite and DES molecules (cellulose crystallite ( E SR,CEL−CEL ), cellulose-choline ( E SR,CEL−CHOL ), cellulose-chloride ( E SR,CEL−Cl ), and cellulose-HBD ( E SR,CEL−HBD )) are displayed in Fig. S14. As observed, intermolecular energies within the cellulose crystallite are bigger than those between cellulose and solvent, thus no dissolution was found for this simulation. The largest E SR,CEL−CEL value (more stable cellulose crystallite) was found within the ChCl/levulinic acid DES (-47,035 kcal/mol in Fig.  S14a) in agreement with the large intra-and inter HBs, presented in Table 4, as well as the smallest net basicity value for the solvent. This behavior is probably caused by the small energy contributions from E SR,CEL−CHOL , E SR,CEL−Cl , and E SR,CEL−HBD in the ChCl/levulinic acid DES, and associated to the small number of HBs from unit cell reported in Table S3,  and Table S4. On the other hand, larger interactions towards cellulose crystallite were found in the ChCl/ ethylene glycol and ChCl/urea DESs (Fig. S14 b-d). This behavior was probably caused by the smallest HS − Cl distance (Table 3) as is also suggested by the stronger Cl DES ⋯ H − D HBs from Table 5 in the case of ChCl/ethylene glycol (Fig. S14c); and with the small HS − O HBD (Table 3, Fig. S14d, and Table S4) and the larger OD − H CHOL , OD − HY and Cl − H HBD distances (Table 3, Fig. S6) in the case of ChCl/urea DES; both in agreement with their larger net basicity values.
Finally, from Fig. S14b, it can be observed that E SR,CEL−CHOL was more stable in ChCl/oxalic acid DES, in agreement with the short OD − H CHOL and OD − HY interaction distances and the large N coord shown in Table 3, as well as the stronger O DES − H DES ⋯ A HBs reported in Table S3 and  Table S4. These large interactions do not seem to reduce the interactions with Cl − anion (Fig.  S14c) when compared to the remaining cellulose-DES systems; however, a largest reduction in the E SR,CEL−HBD can be noticed (Fig. S14d). Interestingly, by combining these interaction energy contributions regarding the cellulose crystallite (i.e. E SR,CEL−CHOL + E SR,CEL−Cl + E SR,CEL−HBD ) from the DESs with net basicity ( − ) definition, an important correlation is obtained, as shown in Fig. 6h ( R 2 = 0.73 ), unlike the null correlation observed with the Kamlet-Taft parameters β and α separately.
Despite the ChCl/levulinic acid DES had the largest Kamlet-Taft β parameter it was not able to disrupt further the cellulose crystallite as indicated by the E SR,CEL−CEL in Fig. S14a and E non−bonded,cellulose . Additionally, the cellulose pretreatment with this solvent produced small number of O6 − H6 ⋯ O4 HB occupancies compared to the remaining solvents. On this regard, it seems that there is a limited role of this DES towards the cellulose crystallite, which it cannot be explained directly by its Kamlet-Taft β parameter. In order to explore this point, the Cl − CHOL , Cl − HBD , CHOL − HBD , CHOL − CHOL and HBD − HBD E SR were also calculated and displayed in Fig. 7.
As it can be noticed from Fig. 7d-e, the interactions between the choline cation and HBD and between HBDs in ChCl/levulinic acid DES showed the largest negative values (-16.15 and -63.67 kcal/ mol, respectively) among all the cellulose-DES systems. These interactions within the DES might be the origin of the low cellulose disruption performance as is reported in Table 4 and Fig. S14a. This finding suggested that both choline cation and the levulinic acid HBD are both interacting strongly, limiting the role of the DES molecules in the cellulose pretreatment. This effect is also correlated with its smallest net basicity value ( − = −0.9427 ) among the evaluated DESs, as shown in Fig. 6i. This competing effect has been reported by Zgrzeba et al. (2015) and by Crowhurst et al. (2003) in ILs, indicating the importance to take into account both Kamlet-Taft parameters.
From Fig. 7a, it can be noted that the largest E SR,Cl−CHOL was found for the ChCl/oxalic acid DES, in agreement with the observations reported in Table 3 for the small Cl − HY interaction distance. Regarding to the E SR,Cl−HBD in Fig. 7b, the ChCl/oxalic acid DES shows the most positive value (-14.95 kcal/mol); this potential advantage over the Cl − anion to interact better with cellulose crystallite, turned into a large E SR,Cl−CHOL interaction. On , and e Hydrogen bond donors ( E SR,HBD−HBD ) in cellulose-ChCl/ethylene glycol system (black bar), cellulose-ChCl/oxalic acid (blue bar), cellulose-ChCl/urea (yellow bar), and cellulose-ChCl/ levulinic acid (green bar) systems. Since the number of clusters were different for all cellulose-DES systems, the interactions were scaled by the number of the corresponding DES clusters present in each system (as indicated in Table S1). Mean values are indicated for each interaction, error bars indicate one standard deviation the contrary, the large E SR,Cl−HBD value was obtained for the ChCl/ethylene glycol DES (-25.58 kcal/mol in Fig. 7b), however, it seems that such interaction did not limit the performance of the chloride anion towards the cellulose disruption as reported in Fig.  S14a and Fig. S14c. This is probably related to the small energy value of the E SR,Cl−CHOL interaction, which releases the chloride anion, in agreement with the large Cl − HY distance shown in Table 3. Interestingly, the ChCl/urea and ChCl/ethylene glycol DESs showed similar E SR,Cl−HBD values (probably related to their similar net basicity values). The ChCl/ urea DES also showed the most positive E SR,HBD−HBD value (52.82 kcal/mol in Fig. 7e). Accordingly, there was not competing between urea HBDs to interact with the cellulose crystallite, as is indicated by the correlation between net basicity values and E SR,HBD−HBD energies, as shown in Fig. 6i. Regarding the role of the interactions between choline cations ( E SR,CHOL−CHOL , Fig. 7c), it can be noticed that all of them were similar and showed repulsive interactions between them.
Finally, from Fig. 7e, the largest repulsion energies between HBDs molecules within the DES could promotes a better performance of HBD towards the cellulose crystallite. This is the case of ChCl/ethylene glycol and ChCl/urea DESs, which also presented the largest E SR,CEL−HBD as indicated in Fig. S14d, and similar net basicity. Interestingly, such behavior can be screening from both α and β Kamlet-Taft parameters. Thus, special attention must be paid to the net basicity values and not only to the Kamlet-Taft α and β parameters separately, in order to pre-screen potential new, green, and sustainable DESs for an effective cellulose pretreatment and lignocellulosic biomass valorization applications.

Conclusions
A series of correlations between Kamlet-Taft α and β parameters and calculated thermodynamic properties of choline-chloride deep eutectic solvents have been investigated in order to provide novel indicators, factors, and insights to pre-screen, design, and develop efficient DESs solvents for dissolution or pretreatment of the LCB, particularly cellulose. Thermodynamic properties such as density and molar volume displayed meaningful correlations towards the Kamlet-Taft α and β parameters. Enthalpy of vaporization showed a modest correlation towards these parameters.
Further insights regarding the interplay between a cellulose crystallite and solvent molecules as well as the decrystallization of cellulose were found through atomistic analysis including radial distribution functions, intra and interchain hydrogen bonds, hydrogen bond occupancies, solute-solvent hydrogen bonds, and interaction energies. These results indicate the slow breaking of the cellulose crystallite as effect of the DES molecules. Radial distribution analysis indicates that Cl − HY interactions in conjunction with its corresponding N coord within the solvent could be an indicative of the hydrogen bond network strength, thus restraining the capability of the solvent to pretreat the cellulose. Additionally, it was found that the HS − Cl and the OD − HO interactions between cellulose and solvent molecules play a relevant role in the deformation or decrystallization of the cellulose crystallite. Interestingly, the OD − HO interplay between glucose units and hydrogen bond donor presented shorter distances than the HS − Cl interaction. Furthermore, it was found that strong HS − oxygen interactions between cellulose and hydrogen bond donor could be enhanced for R − CO − R structures within the hydrogen bond donor.
Analysis of the hydrogen bonds within the cellulose crystallite indicated that the intrachain hydrogen bonds were reduced first compared to the interchain hydrogen bonds. These values, particularly the change of the intrachain hydrogen bonds along time indicate that the deep eutectic solvents can slowly break up the cellulose crystallite. Additionally, hydrogen bond occupancies indicated that the stronger O6 − H6 ⋯ O2∕O3 and O2 − H2 ⋯ O6 interchain hydrogen bonds present in inner chains, are slightly reduced and replaced by weak O6 − H6 ⋯ O4 hydrogen bonds in the glucan located at the edge, i.e. glucan in direct contact with the solvent.
The analysis of the atomistic properties of the cellulose crystallite including hydrogen bonds, did not show a complete correlation with the Kamlet-Taft β or α parameters separately. Accordingly, both parameters need to be considered to further correlate with such properties. The aforementioned was found through the analysis of the hydrogen bonds at quantum chemistry level through density functional theory calculations. Accordingly, the net basicity ( − ) definition of the solvent is required to understand the interactions between cellulose crystallite and eutectic solvents. Interestingly, thermodynamic properties such as cohesive energy and Hildebrand parameter correlated with the net basicity of the evaluated solvents. Furthermore, atomistic properties correlated also with such definition. Energetic analysis such as the pair energies contributions suggested that even though the choline-chloride levulinic acid has the largest Kamlet-Taft β parameter (0.61) among the evaluated solvents, it also has a large hydrogen bond donor-hydrogen bond donor pair energy interaction, which strongly limits the role of the levulinic acid towards cellulose crystallite. Revealing that both parameters are important in order to developed new deep eutectic solvents.
Summarizing, for the evaluated choline-chloridebased deep eutectic solvents, an effective solvent for the cellulose pretreatment could be related (but not limited) to: low density, large molar volume, cohesive energy density of at least 140 cal/cm 3 , Hildebrand parameter of at least 25 MPa 0.5 , R − CO − R moieties within the HBD, large net basicity value, and small HBD-HBD interaction. Fortunately, all these features can be obtained through in-silico evaluations. Therefore, solvent engineering can be used before solvent synthesis. We hope the aforementioned insights provides relevant information for the development of new, sustainable, and more efficient DESs solvents for the cellulose pretreatment.