A DFT study on lignin dissolution in imidazolium-based ionic liquids

Abstract

Density functional theory (DFT), atoms in molecules (AIM) theory, natural bond orbital (NBO) analysis, and reduced density gradient (RDG) analysis were employed to investigate the mechanism of lignin dissolution in imidazolium-based ionic liquids (ILs). Lignin was modeled with guaiacyl glycerol-β-guaiacyl ether (GG), which is one type of β-O-4 linked dimers. Hydrogen bonds (H-bonds) are studied specifically and characterized by different methods to evaluate the strength of interaction between ILs and the lignin model compound. From the theoretical results, it is observed that H-bonds between anions and the GG model are stronger than those between cations and the GG model. Also, anions have the strongest interaction at the α-OH position of GG, while cations have the strongest interaction at the γ-OH position of GG. In addition, anions Cl, OAc and MeSO₄ have much stronger H-bonding ability than PF₆, and the length of the alkyl chain does not have a significant influence on the cation–GG interaction. This work also simulates the interaction between the GG model and ion pairs, with the results suggesting that anions in ion pairs play a key role in forming H-bonds, and cations have a π-stacking interaction with GG. The calculation data provide the interaction mechanism of lignin dissolution in ILs to some extent.

---

1. Introduction

Second to cellulose, lignin is the most abundant bio-renewable and bio-degradable polymer in nature.^1,2 Its complicated structure is derived from coniferyl (G unit), sinapyl (S unit), and p-coumaryl (H unit) alcohols, which are linked by various types of C–O and C–C bonds, including β-O-4, α-O-4, 4-O-5, β-5 and β–β linkages. β-O-4 is the dominant linkage among these, constituting more than half of the linkage units of lignin in both softwoods and hardwoods.³ G units constitute approximately 90% of softwood, whereas roughly equal proportions of G unit and S unit appear in hardwood.⁴ It has been proved that H-bonds between neighboring carboxylic groups, hydroxyl groups and ether groups, and π–π interactions between aromatic moieties are responsible for the associated and complicated structures of lignin.^5–8 Based on the above features, lignin may be used as a potential source of phenolic compounds in order to substitute for petroleum-derived chemicals, and its prospects for creating high value-added chemicals are obvious. However, in order to separate lignin from biomass, traditional processes using acid, alkali, and organic solvents have the disadvantages of pollution, corrosion, odor, high temperature, high pressures and high water use.^9,10 In general, such reagents are incompatible with green chemistry.

Recently, room temperature ionic liquids (RTILs)^11–14 composed entirely of cations and anions have drawn the attention of chemists as promising green media for synthesis, catalysis and separation.^15–17 Rogers et al. have shown that some of the RTILs can dissolve cellulose,¹⁸ and Gomes et al.^19–21 extended many works both experiment and simulation to understand the solvation and dissolution of cellulose in BmimOAc in the presence of cosolvent. It is worth mentioning that some gases, such as CO₂, SO₂ and CH₄, show great solubility in functional ILs.^22–27 Their use may be extended to lignocellulose composites.^28,29 These experiments proved that RTILs have effectively promoted the dissolution and catalytic conversion of biomass materials. Since then, many experimental efforts have been devoted to developing novel ILs in order to dissolve lignin. Fort²⁸ found BmimCl was capable of dissolving both cellulose and lignin simultaneously, and cellulose was easily regenerated by adding co-solvent. Lee³⁰ reported that EmimOAc was able to selectively extract lignin from wood flour, with a highly concentrated solution of chemically unmodified lignin obtained in the recovery of the solvent. Pu et al.⁷³ reported that the solubility of residual softwood kraft pulp lignin in BmimCl reached 13.9 g L⁻¹, while reaching 344 g L⁻¹ and 312 g L⁻¹ in MmimMeSO₄ and BmimMeSO₄ respectively; it was also reported in this work that lignin was insoluble in BmimPF₆. It seems that ILs with Cl, OAc, MeSO₄ anions are better solvents for lignin dissolution, while PF₆-based ILs are unfavorable. However, it is unclear why different anions result in the dramatically different solubility.

Besides reliable results from experiments, molecular simulations provide specific interaction information that reveals deeper mechanism while reducing manual work and experimental costs.³¹ Using the COSMO-RS model, Balaji³² computed Hildebrand solubility parameters for 34 cations and 34 anions, and found that high Hildebrand parameters of imidazolium, pyridinium, pyrrolidinium and ammonium-based ILs suggested higher potential for dissolving lignin than other ILs. The energetic contribution of entropic and enthalpy were used to estimate the dissolution mechanism of gases in ILs,²⁷ and free energy was closely related to solubility of polymers¹⁹ and crystals,³³ thus the thermodynamics of lignin dissolution are speculated to be correctly estimated with these terms. In addition, excess enthalpy of IL/lignin mixtures affects the affinity of lignin with different ILs. Casas³⁴ proved that strong exothermic excess enthalpy was beneficial for lignin dissolution, and that FTIR spectroscopy of regenerated lignin solids further demonstrated good performance of the selected ILs. To study the mechanism of lignin dissolution in ILs, a new prospective different from macroscopic thermodynamic properties is used to reveal the interactions between ILs and polymers at the molecular level, typically Density Functional Theory (DFT) calculations.^35,36 Ji et al.³⁷ concluded that stronger H-bonding interaction in AmimCl–LigOH than in LigOH–LigOH was responsible for lignin dissolution, with the regenerability of lignin in water explained by the reduced H-bonds in AmimCl–LigOH when adding water to this system. Recently, Janesko⁵ performed comprehensive work to elucidate the noncovalent interactions between various ILs and lignocellulose. They highlighted the importance of IL cations in tuning the relative solubilities of lignin and cellulose. Furtherly, they found that ILs with π-conjugated cations are favorable for enhancing lignin solubility. The above studies have shown that the conformers' binding energy is effective to reveal the dissolution mechanism, and demonstrated that the components of the cation and anion in the ILs can be fine-tuned. However, they were mainly focused on the multiple structures and energies of ILs-lignin conformers. To the best of our knowledge, there are few works concerning NBO, AIM and RDG analysis to study the mechanism of lignin dissolution, and very few computational studies shed light on the question why different anion-typed ILs give lignin dramatically different solubility,³⁸ as well as the influences of alkyl substituents on the performance of ILs.

Thus, our work aims at exploring the detailed interactions, especially H-bonding interactions between ILs and the lignin model compound. Due to limited computational capacity, lignin is modeled by the dimer guaiacyl glycerol-β-guaiacyl ether (GG, Fig. S1†), which consists of G units connected with β-O-4 linkages, and whose structure is thought to represent the building blocks of natural lignin. Furthermore, efficiently cleaving the β-O-4 linkage could be of critical importance for the degradation and decomposition of lignin.³ In the present work, the interaction between the GG model and a series of imidazolium-based cations (Mmim, Emim, Bmim and Hmim), anions (Cl, OAc, MeSO₄, PF₆) and ion pairs are studied to reveal the dissolution mechanism by employing theoretical calculations. Cations and anions are used to study the effects of alkyl chain length³⁹ and coordination ability on the strength of the interaction. Interaction sites are studied to figure out which position is helpful for cleavage of the β-O-4 linkage. Interaction energy is used to estimate approximately the most efficient solvent among the studied ILs at the molecular level. Natural bond orbital (NBO) analysis, atoms in molecules (AIM) theory, and reduced density gradient (RDG) are used in combination to specifically investigate the hydrogen bonds (H-bonds) and π-stacking interactions.

2. Computational details

All calculations are performed with the Gaussian 09 package.⁴⁰ The DFT method with the Becke's three-parameter functional and the nonlocal correlation of Lee, Yang, and Parr (B3LYP)⁴¹ together with the 6-31++g** is used for all calculations. Here, we also highlight many studies got reasonable results in correspondence with experiments by using the same level of theory.^14,42–44 All the obtained structures are confirmed by frequency analysis to ensure no imaginary frequencies. It is known that the B3LYP/6-31++g** level is an excellent compromise between computational cost and accuracy of the computational results; further, the 6-31++g** basis set was thought to be effective for considering the double-ξ quality for valence electrons with the consideration of certain diffuse functions.⁴⁵ Therefore, the basis set is suitable for the ILs and lignin system.^46–48 To get comparable results with Janesko,⁵ dispersion-corrected^49–52 DFT is applied to determine the significant π-stacking interactions at the B3LYP-D3/6-311+g**,⁴¹ and several reasonable geometries are also performed using M06-2X/6-311+g** level of theory.^27,53 In these calculations, the electrostatic potential (ESP) method is used to predict the most plausible sites for electrophilic and nucleophilic attack;⁵⁴ the initial configurations of the cation–GG, anion–GG, and ion pairs with GG are structured based mainly on the ESP (Fig. S2†) on the isosurface of cations, anions and GG. The most stable conformer in this work is obtained by comparing the electron energies of different initial structures designed and optimized in the calculations. The interaction energy (ΔE) is defined as the energy difference between the conformer and the corresponding isolated ions.

ΔE = 2625.5 × (E_ab − E_a − E_b), kJ mol⁻¹

Natural bond orbital (NBO) analysis, atoms in molecules (AIM) theory and reduced density gradient (RDG) analysis are performed with the software Multiwfn,⁵⁵ VMD⁵⁶ and CYL view⁵⁷ at the same level to confirm the existence of H-bonds. In the NBO analysis, E(2) is used to describe the intensity of orbital interaction between electron donors and electron acceptors well within the equation.

The higher the value of E(2), the more electrons tend to transfer from the donor.⁵⁸ Where ε_i and ε_j are diagonal elements, q_i is the donor orbital occupancy and F(i,j) is the off-diagonal NBO Fock matrix element.

According to the AIM theory, the electron density (ρ_BCP) and its Laplacian value (∇²ρ_BCP)^59,60 are used to indicate the strength and nature, respectively, of H-bonds at the bond critical point (BCP). The H-bonds and ionic bonds are characterized by ∇²ρ_BCP > 0 in the closed shell interactions, while the covalent bonds are characterized by ∇²ρ_BCP < 0; the extremely strong H-bonds can also have negative Laplacian value. RDG scatter plots and isosurfaces visualize^55,61,62 all non-covalent interaction in the system, as indicated in Fig. S3.†

3. Results and discussion

3.1 Interaction between cations/anions and GG

3.1.1 Geometrics and interaction energy of cations/anions–GG. The structures of isolated cations, anions and the GG model compound are optimized at the B3LYP/6-31++g** level from different initial guesses (Fig. S4–S8†). The lowest energy conformers are used to construct the initial structures for the interaction between cations and GG, and between anions and GG. All of the anion–GG and cation–GG conformers can be found in the ESI;† the strongest interactions of conformers Cl–GG, OAc–GG, MeSO₄–GG, PF₆–GG, Mmim–GG, Emim–GG, Bmim–GG and Hmim–GG are depicted in Fig. 1. The strongest interaction between anions and GG appears when H-bonds are formed between the electronegative atoms (Cl, O, F) of anions and the α-OH of GG. Alternately, the strongest interaction between cations and GG appears when H bonds are formed between C2–H of the imidazolium ring and the oxygen atom of γ-OH of GG. It is noteworthy that, due to a greater number of electronegative oxygen atoms, anions such as MeSO₄ tend to form multiple H-bonds compared to Cl and OAc, and that cations have two H-bonds with the exposed oxygen atoms of the hydroxyl group. The formation of C–H⋯Cl H-bonds is determined if the distance between H and Cl is <2.95 Å,⁶³ which is the sum of the van der Waals radii of H and Cl atoms. The corresponding distances for C–H⋯O, C–H⋯F H-bonds should be shorter than 2.72 Å and 2.67 Å,⁶³ respectively. The corresponding H-bond distances are labeled in Fig. 1 and the interaction energies are summarized in Table 1.

Fig. 1 Optimized anion–GG conformers (a–d), and cation–GG conformers (e–h) at the level of 6-31++g**. H-Bonds are indicated by dashed lines.

Table 1 H bonds and interaction energy of anions–GG, cations–GG conformer

Structure		H-Bond	Length (Å)	Angle (°)	ΔE (kJ mol^−1)
Fig. 1(a)	Cl–GG	Cl⋯H43–O42	2.154	166.92	−95.21
		Cl⋯H9–C6	2.604	139.29
Fig. 1(b)	OAc–GG	O50⋯H43–O42	1.625	175.97	−121.46
		O49⋯H9–C6	2.169	153.97
		O49⋯H15–C13	2.439	141.06
Fig. 1(c)	MeSO4–GG	O48⋯H43–O42	1.873	172.46	−81.98
		O45⋯H15–C13	2.341	139.02
		O46⋯H9–C6	2.305	149.38
Fig. 1(d)	PF6–GG	F50⋯H43–O42	1.913	172.13	−61.69
		F49⋯H15–C13	2.319	169.75
		F48⋯H9–C6	2.366	144.49
		F48⋯H18–C17	2.613	135.13
Fig. 1(e)	Mmim–GG	O35⋯H50–C46	2.001	155.75	−72.93
		O35⋯H53–C52	2.593	144.70
Fig. 1(f)	Emim–GG	O54⋯H7–C3	2.007	155.72	−70.76
		O54⋯H18–C16	2.560	144.86
Fig. 1(g)	Bmim–GG	O60⋯H7–C3	2.027	152.66	−68.46
		O60⋯H17–C15	2.509	146.41
Fig. 1(h)	Hmim–GG	O66⋯H7–C3	2.013	154.36	−67.94
		O66⋯H17–C15	2.559	145.43

---

Fig. 1 shows two important features of these geometries: (i) α-OH-involved H-bonds (e.g. 2.154 Å in structure (a)) have shorter bond lengths than other H-bonds (e.g. 2.604 Å in structure (a)) formed between GG model compound and anions; and (ii) the C2–H of cations gets closer to γ-OH than the methyl group of cations (e.g. 2.001 Å and 2.593 Å of the two H-bonds in the structure (e)). All of these optimized conformers show that the optimal interaction site for anion–GG interaction is the α-OH group of GG, and the γ-OH of GG for cation–GG interaction. Comparing the interaction energies of these structures, the order of anions to form H-bonds with GG follows: OAc > Cl > MeSO₄ > PF₆, with the sequence of cations being: Mmim > Emim > Bmim > Hmim. It can be found that increasing alkyl chain length leads to a decrease of interaction energy between cations and GG, but the change is not obvious. Hart⁶⁴ demonstrated that the proportion of lignin in saturated solutions of lignin in EmimMeSO₄ was nine times that of EmimPF₆, where the low H-bonding basicity of PF₆ means that it cannot dissolve lignin (the theoretical results are consistent with those obtained by the experiments). Therefore, we speculate that these intermolecular H-bonds are expected to be a key factor resulting in the solubility of lignin in ILs. Furthermore, anions may play more important roles than cations in the dissolution of lignin because of stronger interaction with the lignin model compound. These results also may be used to show that lignin cannot be dissolved in PF₆-based ILs because of weak interaction in the PF₆–GG conformers.

3.1.2 NBO analysis of cations/anions and GG. Natural bond orbital (NBO) method⁶⁵ characterizes H-bonds in terms of hyper-conjugative donor–acceptor interactions and has been used to study the bonding properties in theses conformers. Table 2 shows the main donor–acceptor interactions between anions/cations with GG, as well as their second-order perturbation stabilization energies, E(2). The values of E(2) denote the strength of donor–acceptor interaction, and the larger the value is, the stronger the interaction will be. From the table, the obvious and efficient overlaps are found between lone-pair orbitals of anions and anti-bonding orbitals of GG, and between anti-bonding orbitals of cations and lone-pair orbitals of GG. The E(2) values show that some of the H-bonds are strong and some are weak. Among these, the strongest H-bonds in anion–GG conformers are formed between the electronegative atoms of anions and α-OH of GG, while the strongest H-bonds in cation–GG conformers are formed between the C2–H of cations and γ-OH of GG. It seems that OAc and Mmim are the most favorable anion and cation for the interaction, since the E(2) values of Lp O49 → σ* O42–H43 and Lp O35 → σ* C46–H50 are 244.92 and 52.91 kJ mol⁻¹, respectively. It can be concluded that Cl, OAc, and MeSO₄ have stronger interaction with GG than PF₆, which may be attributed to the weak acceptors of hydrogen atoms in PF₆. The length of the alkyl chain has little effect on cation–GG interaction, as can be seen from the small differences in E(2) energies. It has been confirmed that anions rather than cations play a key role when lignin is dissolved in ILs.

Table 2 The main electron donor–acceptor interactions in the anion–GG, cation–GG conformers and their second-order perturbation stabilization energies (E(2)) at B3LYP/6-31++g** level

Structure	Donor (i)	Acceptor (j)	E(2), kJ mol^−1
Cl–GG	Lp Cl44	σ* C6–H9	20.43
	Lp Cl44	σ* O42–H43	91.67
OAc–GG	Lp O49	σ* C6–H9	29.97
	Lp O49	σ* C13–H15	8.62
	Lp O50	σ* O42–H43	244.92
MeSO4–GG	Lp O45	σ* C13–H15	6.82
	Lp O46	σ* C6–H9	16.62
	Lp O46	σ* O42–H43	72.88
PF6–GG	Lp F48	σ* C6–H9	8.92
	Lp F48	σ* C17–H18	3.14
	Lp F49	σ* C13–H15	11.64
	Lp F50	σ* O42–H43	41.27
Mmim–GG	Lp O35	σ* C46–H50	52.91
	Lp O35	σ* C52–H53	3.85
Emim–GG	Lp O54	σ* C3–H7	51.61
	Lp O54	σ* C16–H18	3.77
Bmim–GG	Lp O60	σ* C3–H7	47.05
	Lp O60	σ* C15–H17	5.65
Hmim–GG	Lp O66	σ* C3–H7	50.11
	Lp O66	σ* C15–H17	4.60

---

3.1.3 AIM analysis of cations/anions and GG. To obtain more information about the intermolecular interaction, the AIM theory is used to analyze bonding characteristics. The values of ρ_BCP and ∇²ρ_BCP for the intermolecular H-bonds in anion–GG conformers and cation–GG conformers are summarized in Table 3. There are a set of criteria for ρ_BCP and ∇²ρ_BCP proposed at bond critical points (BCPs) for the conventional H-bonds. Both parameters for the closed-shell interactions as H-bonds are positive within the ranges of 0.002–0.035 a.u. for the electron density and 0.024–0.139 a.u. for its Laplacian value.⁶⁶ For most H-bonds considered here, the values of ρ_BCP and ∇²ρ_BCP lie in the relative proposed ranges, and the extreme strong H-bond-like O50⋯H43–O42 in the OAc–GG conformer has very large values of ρ_BCP (0.0528 a.u.) and ∇²ρ_BCP (0.1439 a.u.). Therefore, for the observed conformer, the ρ_BCP and ∇²ρ_BCP of H-bonds fall within 0.0060–0.0528 a.u. and 0.0261–0.1439 a.u., respectively. The Laplacian of electron density at BCPs have positive values that show why the nature of these H-bonding interactions is noncovalent. Moreover, the negative values of energy density at BCPs indicate these H-bonds have the characteristic of covalent bonds, and the positive values of H_BCP indicate these H-bonds have electrostatic properties. It is observed that the strength of H-bonds follows the order of OAc–GG > Cl–GG > MeSO₄–GG > PF₆–GG in the anion–GG interactions, and Mmim–GG > Emim–GG > Bmim–GG > Hmim–GG in the cation–GG interactions. This result again confirms that OAc is the most effective anion and Mmim is the most effective cation in the interaction with lignin model compound, which is consistent with the previous results. Comparing the electron properties of anion–GG and cation–GG conformers, H-bonding interactions between anions and GG are stronger than those between cations and GG. It is reasonable that α-OH (O42–H43) of GG plays a key role in anion–GG interactions and O35/O54/O60/O66 in γ-OH of GG are the main sites when interacting with the C2–H group of cations, which lead to these stronger H-bonds in the same conformer.

Table 3 Properties of electron density at bond critical points (BCPs) for the interaction of anion–GG and cation–GG (a.u.)

Structure	H-Bond	ρBCP	∇^2ρBCP	HBCP (10^−3 a.u.)
Cl–GG	Cl⋯H43–O42	0.0285	0.0587	−1.6143
	Cl⋯H9–C6	0.0128	0.0392	1.3348
OAc–GG	O50⋯H43–O42	0.0528	0.1439	−2.7072
	O49⋯H9–C6	0.0183	0.0504	−0.3137
	O49⋯H15–C13	0.0101	0.0334	0.7030
MeSO4–GG	O46⋯H43–O42	0.0315	0.1010	−0.9943
	O45⋯H9–C6	0.0132	0.0412	0.3625
	O45⋯H15–C13	0.0092	0.0310	0.7458
PF6–GG	F50⋯H43–O42	0.0219	0.0711	−0.8172
	F48⋯H9–C6	0.0098	0.0387	0.9949
	F48⋯H18–C17	0.0060	0.0261	1.2580
Mmim–GG	O35⋯H50–C46	0.0252	0.0671	−0.9019
	O35⋯H53–C52	0.0074	0.0266	0.9087
Emim–GG	O54⋯H7–C3	0.0248	0.0663	−0.8775
	O54⋯H18–C16	0.0074	0.0263	0.9106
Bmim–GG	O60⋯H7–C3	0.0238	0.0641	−0.8031
	O60⋯H17–C15	0.0089	0.0303	0.7937
Hmim–GG	O66⋯H7–C3	0.0235	0.0658	−0.8407
	O66⋯H17–C15	0.0080	0.0280	0.8634

---

3.1.4 RDG analysis of cations/anions and GG. RDG analysis is employed as another useful method to further study noncovalent interaction in this work. This method can supply more reliable data compared to the AIM method. Reduced density gradient (RDG) versus the electron density multiplied by the sign of second Hessian eigenvalue (sign(λ₂)ρ)⁶² is plotted by scatter graph; the visualization of these inter- and intra-molecular weak interactions is shown in Fig. 2. It can be seen that many RDG spikes in scatter graph, as well as the spikes corresponding to H-bonding interactions, van der Waals interactions and steric effects, are shown as a function of sign(λ₂)ρ from negative to positive values. The spikes and disc-shaped blocks that represent the strongest H-bonds in OAc–GG and Mmim–GG conformers are marked by the blue and red circles, respectively. As can be seen, the spike corresponding to O50⋯H43–O42 is located at 0.0528 a.u. and the spike corresponding to O35⋯H50–C46 is located at 0.0252 a.u. Distinctly, the value of sign(λ₂)ρ in Fig. 2(A) is larger than that in Fig. 2(C), and the corresponding region color of the disc-shaped block in structure B is dark blue in contrast to the light blue in structure D between H and O atoms. It can be interpreted that the interaction between GG and OAc is stronger than Mmim. The results are consistent with geometrical, NBO and AIM analysis. The same analysis is suitable for other conformers formed between GG and anions/cations in this study, with the results shown in Fig. S10–S11.† The much stronger attraction force between anions and GG than between cations and GG suggests that anions of ILs more favorably interact with the model compound when dissolving lignin. Reasoning that anions with strong coordinating ability have strong H-bonding ability which is associated with the H-bond basicity at the same time,⁶⁷ and it is previously reported that the H-bond basicity and the size of the anion play an important role in disrupting the structure of the lignin.⁶⁸ Meanwhile, the length of alkyl chain do not affect the interaction significantly because it has little influence on charge distribution and the positive charges are always distributed around imidazolium ring.

Fig. 2 RDG scatter plots (isovalue = 0.5 a.u.) and surface plots (s = 0.7 a.u.) of (A and B) OAc–GG interaction and (C and D) Mmim–GG interaction. The isosurfaces are colored on a blue-green-red scale according to values of sign(λ₂)ρ, ranging from −0.03 to 0.02 a.u. Blue indicates strong attractive interactions and green indicates the π-stacking interaction.

3.2 Interaction between ion pairs and GG

Since the interaction between GG and single anion/cation has been calculated, and because anions and cations appear in the form of ion pairs in the real structure of ILs, it is therefore of great importance to study the structural behavior between model compound and ion pairs which is essential to reveal the dissolution mechanism. The conformers formed between ion pairs and GG model compound were optimized at B3LYP/6-31++g** level. In this part, the interaction energies of 48 possible conformers are compared in Fig. 3(a), wherein we see that increasing the length of the alkyl chain results in a decrease of interaction energies between ion pairs and GG; these changes are not obvious, but are more apparent for anions in ion pairs. Therefore, four Bmim-based ion pairs with different anions are selected to analyze the interaction between GG and ion pairs at different sites in Fig. 3(b). These have much stronger interaction at α-OH of GG, which suggests that anions in ion pairs play a key role in lignin dissolution. The following analysis is performed for the conformers BmimCl–GG, BmimOAc–GG, BmimMeSO₄–GG and BmimPF₆–GG at the α-OH position.

Fig. 3 Interaction energy comparison of all possible structures for ion pairs–GG interaction.

3.2.1 Geometrics and interaction energy. The interaction between ion pairs and GG is concentrated in the four ion pairs. In Fig. 4, structures A, B and C are the possible conformers of BmimCl–GG at α-OH, p-OH and γ-OH positions, respectively; structures D, E and F are the conformers of BmimOAc–GG, BmimMeSO₄–GG and BmimPF₆–GG at the α-OH position of GG. As can be seen, both anions and cations can form H-bonds with GG simultaneously; these are denoted by dashed lines. The optimized bond length, bond angles and interaction energies for the six conformers are summarized in Table 4. Usually, the stronger the interaction between hydrogen donor and acceptor, the shorter the bond length is. Angles of H-bonds are within the range of 90° to 180°. The interaction energies of A, B and C are −437.35 kJ mol⁻¹, −420.47 kJ mol⁻¹ and −417.19 kJ mol⁻¹; thus, the favorable position of GG follows α-OH > p-OH> γ-OH position, a sequence that also applies for other ion pairs–GG interactions. The interaction energies of structures BmimCl–GG, BmimOAc–GG, BmimMeSO₄–GG and BmimPF₆–GG at α-OH position of GG are −437.35 kJ mol⁻¹, −474.53 kJ mol⁻¹, −417.87 kJ mol⁻¹, and −365.09 kJ mol⁻¹, respectively; thus, the order for the H-bonding interactions is BmimOAc > BmimCl > BmimMeSO₄> BmimPF₆. All structures show that the α-OH group is the key site when ion pairs interact with GG model compound, and the p-OH is the secondary site favorable for H-bonds. The interaction between ion pairs and GG at γ-OH is not as strong as the first two sites because of the intramolecular H-bond formed between γ-OH and the adjacent methoxy group, making it difficult for ILs to approach γ-OH of GG. In addition, strong interaction between BmimOAc/BmimCl/BmimMeSO₄ and GG favors lignin dissolution, while weak interaction between PF₆-based ILs and GG is a key factor that explains the insolubility of lignin in these kinds of ILs.

Fig. 4 Optimal structures of interaction between GG and Bmim based ILs at B3LYP/6-31++g** level and H bonds are indicated by dashed line. Structures (A), (B), and (C) are BmimCl–GG interaction at α, p, and γ position respectively. Structures (D), (E) and (F) are BmimOAc–GG, BmimMeSO₄–GG, and BmimPF6–GG interaction at α position.

Table 4 H bonds and interaction energy of ion pairs with GG at level of 6-31++g**

Structure		H-Bond	Length (Å)	Angle (°)	ΔE (kJ mol^−1)
Fig. 4(A)	BmimCl–GG	Cl69⋯O42–H43	2.197	161.59	−437.35
		O21⋯C58–H59	2.431	138.58
Fig. 4(B)	BmimCl–GG	Cl69⋯O10–H11	2.188	152.34	−420.47
		O10⋯C58–H59	2.240	171.65
Fig. 4(C)	BmimCl–GG	Cl69⋯O35–H36	2.107	167.44	−417.19
		O35⋯C52–H54	2.224	164.17
Fig. 4(D)	BmimOAc–GG	O74⋯O42–H43	1.719	168.21	−474.53
		O42⋯C52–H54	2.362	167.41
		O75⋯C13–H15	2.412	166.41
Fig. 4(E)	BmimMeSO4–GG	O71⋯O42–H43	1.890	165.18	−417.87
		O70⋯C13–H15	2.538	144.92
		O21⋯C52–H54	2.425	149.71
Fig. 4(F)	BmimPF6–GG	F73⋯O42–H43	1.906	159.03	−365.09
		O42⋯C52–H54	2.427	170.77
		F71⋯C13–H15	2.690	136.34

---

As described before, the interaction of ion pairs–GG is more intensive at α-OH than at γ-OH and p-OH. Also, the binding energies for cation–GG and anion–GG at the B3LYP/6-31++g** level can well explain the H-bonding interactions in this system. The above results reveal that the electrostatic interactions play a key role in the system when there exists only one ion (Section 3.1). The electronegative position of anion tends to attack the –OH groups of GG, and C2–H of imidazolium cations prefer forming H-bonds with O atoms of GG. All of these results are a bit different from Janesko's work.⁵ Nevertheless, we also confirm there exist π-stacking effects in the conformers of ion pairs–GG as they reported previously. As tested by Grimme's,⁴⁹ B3LYP-D3 is considered to be an ideal choice to investigate large molecular system, taking the dispersion interaction into account to enhance its capability for noncovalent interaction. The improvement of the B3LYP-D3/6-311+g** level is indicated by the changes in H-bond distances and angles (Table S12 in ESI†). In Table 5, 16 kinds of conformers of ion pair–GG are studied at B3LYP/6-31++g** and B3LYP-D3/6-311+g** levels of theory. As can be observed, the binding energies for ion pairs are −382.12 kJ mol⁻¹, −378.22 kJ mol⁻¹, −374.99 kJ mol⁻¹ and −373.85 kJ mol⁻¹ at the B3LYP/6-31++g** level for MmimCl, EmimCl, BmimCl, and HmimCl, respectively. That increase the alkyl chain of substituent results in a decrease of the interaction energies between cations and anions. The rule is applied for ΔE_1,CA of other anion-typed ILs, and the similar results are obtained (Table 5). The elongation of alkyl chain leads to the increase of the ion volume, promoting charge distribution to a greater extent, which leads to a decline of binding energies interaction between cations and anions, i.e. ΔE_1,CA and ΔE_2,CA. Refer to the data of ΔE_1,CA–GG, BmimMeSO₄–GG system has relatively higher binding energy, −64.15 kJ mol⁻¹, while the other three's binding energies are −62.36 kJ mol⁻¹ (BmimCl–GG), −61.72 kJ mol⁻¹ (BmimOAc–GG), and −44.01 kJ mol⁻¹ (BmimPF₆–GG), respectively. This result can explain the solubilities of lignin in different ILs in some degree. MeSO₄ anion based ILs have good H-bond basicity and multiple electron donors, possessing superior ability to dissolve lignin, while PF₆ anion based ILs is too sluggish to dissolve lignin.

Table 5 Comparison of interaction energies at the B3LYP/6-31++g** and B3LYP-D3/6-311+g** levels of theory (C – cation, A – anion, CA – ion pairs)

Structures	B3LYP/6-31++g**, kJ mol^−1			B3LYP-D3/6-311+g**, kJ mol^−1
	ΔE1,CA	ΔE1,CA–GG	ΔE1,C–A–GG	ΔE2,CA	ΔE2,CA–GG	ΔE2,C–A–GG
MmimCl–GG	−382.12	−62.37	−444.49	−391.21	−116.72	−507.92
EmimCl–GG	−378.22	−61.98	−440.20	−389.77	−117.87	−507.64
BmimCl–GG	−374.99	−62.36	−437.35	−386.98	−121.53	−508.51
HmimCl–GG	−373.85	−62.56	−436.41	−385.79	−124.13	−509.92
MmimOAc–GG	−421.44	−60.55	−481.99	−434.10	−121.75	−555.85
EmimOAc–GG	−416.53	−61.77	−478.30	−429.77	−127.03	−556.79
BmimOAc–GG	−412.81	−61.72	−474.53	−428.94	−129.28	−558.22
HmimOAc–GG	−411.52	−61.61	−473.13	−427.60	−132.34	−559.93
MmimMeSO4–GG	−363.00	−65.06	−428.06	−392.37	−123.91	−516.28
EmimMeSO4–GG	−356.81	−66.40	−423.21	−386.54	−131.81	−518.35
BmimMeSO4–GG	−353.72	−64.15	−417.87	−384.77	−135.21	−519.98
HmimMeSO4–GG	−352.53	−62.94	−415.47	−383.58	−138.85	−522.43
MmimPF6–GG	−324.77	−45.05	−369.82	−354.59	−100.18	−454.77
EmimPF6–GG	−323.32	−44.92	−368.24	−356.29	−97.85	−454.14
BmimPF6–GG	−321.08	−44.01	−365.09	−357.15	−99.44	−456.59
HmimPF6–GG	−319.98	−44.24	−364.22	−356.61	−102.33	−458.94

---

The interaction energy sequence of ΔE_1,CA–GG follows MmimMeSO₄–GG > MmimCl–GG MmimOAc–GG > MmimPF₆, which is inconsistent with the H-bond ability.⁶⁷ According to the previous work,^69,70 we further consider the substituent and π–π stacking effects on ILs–GG system by reoptimizing those obtained geometries at B3LYP-D3/6-311+g** level, which are expected to have great contribution to the total interaction. It is observed that interaction energies of ion pairs (ΔE_2,CA) are nearly 10–30 kJ mol⁻¹ higher than ΔE_1,CA, and ΔE_2,CA–GG and ΔE_2,C–A–GG are increased by 54–76 kJ mol⁻¹ and 63–107 kJ mol⁻¹, respectively, showing D3 term can fairly well consider the dispersion effects in ion pairs–GG system. Owing to the dominant electrostatic interactions between cations and anions, the values of ΔE_1,CA and ΔE_2,CA do not show such great differences, while the large increased values of ΔE_1,C–A–GG and ΔE_2,C–A–GG show a big proportion of electrostatic interaction between them. Refer to ΔE_2,CA and ΔE_2,CA–GG of the conformers (MmimMeSO₄–GG, EmimMeSO₄–GG, BmimMeSO₄–GG and HmimMeSO₄–GG), there is a weakening of cation–anion interaction with the elongation of alkyl substituent (ΔE_2,CA), accompanied by increased binding energies (ΔE_2,CA–GG) because ILs can provide larger spaces to contact with lignin GG model. However, it is clear from the data that increasing alkyl chain is not significant enough to influence the interaction energies of the four conformers. As reported by Janesko,⁵ cations are more important for lignin dissolution than cellulose, but anions still play key roles in governing the dissolution process, which also has been proved by experiment.⁷¹ Moreover, other factors, such as the steric hindrance and viscosity⁷² of ILs would affect the dissolution of lignocellulose. In addition, the previous basis B3LYP/6-31++g** may not be as accurate as B3LYP-D3/6-311+g** for thoughtless of dispersion interaction, which still gives a reasonable tendency in consistent with experimental results. The interaction energy is a metric that should be as an approximation to solubility in the lack of a better metric such as the free energy. We believe that the descriptor of interaction energy can provide an interesting prospective at the molecular level to understand the dissolution mechanism.

3.2.2 NBO analysis of ILs–GG. In addition to interaction energy, NBO analysis also gives the electron-related properties. The most important acceptor–donor interactions and their second-order perturbation stabilization energies are listed in Table 6 and their atomic labels can be found in Fig. S12 (more data can be found in Table S7†). Cl, O and F can act as electron donors and GG can be considered as electron acceptor in the ion pairs–GG conformers. As studied by Janesko,⁵ the extended π-systems of imidazolium ILs yield stronger interactions with lignin, the π-stacking interactions in our work are also observed. The dominant effect is not only the strong H-bond between the anions and GG, but also comparable π–π and CH–π interactions between the imidazolium cations and benzene ring of GG model. The maximum values of E(2) for H-bonds in each ion pair–GG conformer are Lp O74 → σ* O42–H43 (94.08 kJ mol⁻¹), Lp Cl69 → σ* O42–H43 (71.57 kJ mol⁻¹), Lp O62 → σ* O42–H43 (69.97 kJ mol⁻¹) and Lp F69 → σ* O42–H43 (26.25 kJ mol⁻¹). Fig. 5 shows that the strongest H-bonds in the conformers are formed due to the orbital overlapping between lone pairs of anions and the anti-bonding orbital of α-OH (O42–H43) in GG. The E(2) of BmimOAc–GG conformer with the donor–acceptor orbitals Lp O74 → σ* O42–H43 (94.08 kJ mol⁻¹) is larger than those of configurations BmimCl–GG, BmimMeSO₄–GG and BmimPF₆–GG. In BmimOAc–GG conformer, the much stronger interaction of Lp O74 → σ* O42–H43 compared to Lp O42 → σ* C52–H54 suggests that H-bonds formed between GG and anions of ion pairs are much stronger than those between GG and cations, indicating that anions contribute more to the total interaction energy. What's more, the sequence of E(2) is found to correlate well with the calculated interaction energy, which follows the order: BmimOAc–GG > BmimCl–GG > BmimMeSO₄–GG > BmimPF₆–GG. In order to understand the CH–π and π–π interaction effect, the E(2) energies indicating π-stacking are analyzed, for the BmimOAc–GG conformer, there are three pairs of E(2) related to the π-stacking interaction, namely, π C1–C2 → π* C46–N47, π C44–C45 → π* C3–C4 and π C1–C2 → σ* C52–H54. The stabilization energies are 2.77 kJ mol⁻¹, 0.63 kJ mol⁻¹ and 0.84 kJ mol⁻¹, respectively. Though π–π interaction is more significant than CH–π interaction, the total 4.24 kJ mol⁻¹ of π-stacking does not play a major part in this system as compared to E(2) of H-bonds. The same phenomenon occurs in the BmimCl–GG, BmimMeSO₄–GG and BmimPF₆–GG conformers.

Table 6 Electron donor orbitals, electron acceptor orbitals, and the corresponding second-order interaction energy E(2) of ion pair–GG conformers

Structure	π-Stacking interaction, kJ mol^−1			H-Bond interaction, kJ mol^−1
	Donor	Acceptor	E(2)	Donor	Acceptor	E(2)
Fig. 5(a)	π C1–C2	π* C46–N47	2.77	Lp O74	σ* O42–H43	94.08
	π C1–C2	σ* C52–H54	0.84	Lp O75	σ* C13–H15	7.27
	π C44–C45	π* C3–C4	0.63	Lp O42	σ* C52–H54	6.59
Fig. 5(b)	π C1–C2	π* C46–N47	3.77	Lp Cl69	σ* O42–H43	71.57
	π C44–C45	π* C3–C4	0.81	Lp Cl69	σ* C13–H15	5.08
	σ C62–H64	π* C22–C24	0.60	Lp O42	σ* C52–H54	2.86
Fig. 5(c)	π C1–C2	π* C46–N47	1.47	Lp O62	σ* O42–H43	69.97
	π C44–C45	π* C3–C4	0.92	Lp O61	σ* C13–H15	9.70
	π C1–C2	σ* C52–H54	0.67	Lp O60	σ* C6–H9	7.52
Fig. 5(d)	σ C1–C2	π* C46–N47	3.36	Lp F69	σ* O42–H43	26.25
	π C1–C2	σ* C52–H54	0.67	Lp O42	σ* C52–H54	14.66
	π C44–C45	π* C3–C4	0.67	Lp F73	σ* C13–H15	4.28

---

Fig. 5 Natural bond orbital interaction in (a) BmimOAc–GG, (b) BmimCl–GG, (c) BmimMeSO₄–GG, and (d) BmimPF₆–GG conformer calculated at B3LYP-D3/6-311+g** level (pink: isovalue = 0.05, blue: isovalue = −0.05).

3.2.3 AIM analysis of ILs–GG. In order to reveal the chemical origin of the interaction between ion pairs and GG model compound, AIM analysis is employed to validate the intermolecular interactions. The values of ρ_BCP, ∇²ρ_BCP and H_BCP are obtained from AIM calculations, as shown in Table 7. The positive electron density values confirmed that there are good electrostatic interactions between GG and ion pairs in each conformer. For most H-bonds considered here, the ρ_BCP and ∇²ρ_BCP values lie in the relative proposed ranges. For the BmimOAc conformer, ρ_BCP and ∇²ρ_BCP at BCPs confirm that H-bond O42–H43⋯O74 is stronger than H-bonds C52–H54⋯O42 and C13–H15⋯O75. The negative value of H_BCP suggests that O42–H43⋯O74 has the covalent characteristic and positive values of the other two H-bonds, and that C52–H54⋯O42 and C13–H15⋯O75 have the electrostatic characteristic. It is also found that electron densities in the BmimOAc–GG conformer are larger than those of three H-bonds in the BmimPF₆–GG conformer. The results show that H-bond interactions in BmimOAc–GG are stronger than those in the other three conformers, which agree with the previous interaction energy, and that anions of ion pairs are found to play an important role in ion pairs–GG interactions. Experimentally, ionic liquid anions are found to play an important role in dissolving lignin and lignocellulose,^30,73 which can be explained by the strong H-bond interactions, whereas non-coordinating anions such as PF₆ are unfavorable for H-bonds, as they cannot dissolve lignin.

Table 7 Properties of electron density at bond critical point (BCP) for the interaction of ILs–GG (a.u.)

Structure	H-Bond	ρBCP	∇^2ρBCP	HBCP (10^−3 a.u.)
BmimOAc–GG	O42–H43⋯O74	0.045	0.141	−3.531
	C13–H15⋯O75	0.013	0.044	1.506
	C52–H54⋯O42	0.012	0.037	1.002
BmimCl–GG	O42–H43⋯Cl69	0.029	0.068	−1.025
	C13–H15⋯Cl69	0.006	0.016	0.684
	C52–H54⋯O42	0.011	0.033	0.962
BmimMeSO4–GG	O42–H43⋯O62	0.036	0.127	−0.311
	C13–H15⋯O61	0.012	0.037	1.146
	C6–H9⋯O60	0.012	0.038	1.118
BmimPF6–GG	O42–H43⋯F69	0.022	0.091	1.596
	C52–H54⋯O42	0.016	0.052	1.479
	C13–H15⋯F73	0.010	0.037	1.129

---

3.2.4 RDG analysis of ILs–GG. Studied with the reduced density gradient along with the sign of the second Hessian eigenvalue, the noncovalent interactions in this part reveal the strong attraction, van der Waals effect and steric repulsion in the system. The weak interactions are shown by the spikes in RDG scatter plots and RDG surface plots (Fig. 6). Because of the similarity among the four conformers (Fig. S13†), we only take BmimMeSO₄–GG as an example. The sign(λ₂)ρ ranges from −0.05 to 0.05, and the RDG ranges from 0 to 2.0. Spikes at the position −0.036 a.u. marked by blue circles indicate the strongest O42–H43⋯O62 H-bond interaction in the BmimMeSO₄–GG conformer, other spikes near zero are attributed to π-stacking interactions, and these spikes are more condensed. Surface plots in Fig. 6 show the noncovalent interactions, wherein it may be observed that strong H-bonds exist between anions and C2–H of cations. Meanwhile, anions can form strong H-bonds with α-OH of GG, as highlighted by the red circles. Also, steric repulsions appear in the benzene ring and imidazolium ring indicated by red ellipsoid. Other large green surfaces can be classified as π–π stacking effects between the two rings. The RDG isosurfaces with yellow or green colors indicate that such interactions are not so strong.⁷⁰ The face-to-face orientation is the prototypical system for studying π-stacking effect, and it is generally agreed that dispersion forces play the major role.^69,74,75 Hence, we can safely draw the conclusion that there exists both π-stacking and H-bonds interactions in the BmimMeSO₄–GG conformers, and this type of interactions are believed to exist in other conformers (Fig. S13†). Therefore, these insights illustrate that dissolution of lignin in ILs may be synergistically affected by H-bonding interactions of lignin with the ILs anions and π–π stacking interactions of lignin with the ILs cations.

Fig. 6 (A) RDG scatter plots (isovalue = 0.5 a.u.) and (B) isosurfaces (s = 0.7 a.u.) of BmimMeSO₄–GG conformer. The isosurfaces are colored on a blue-green-red scale according to values of sign(λ₂)ρ, ranging from −0.03 to 0.02 a.u. Blue indicates strong attractive interactions and green indicates the π-stacking interaction.

4. Conclusion

Structures and electron properties in the interactions between ILs and lignin model compound are studied using a density functional theory at B3LYP/6-31++g** level and dispersion corrected DFT at B3LYP/6-311+g** level. Lignin is modeled with guaiacyl glycerol-β-guaiacyl ether (GG), and the interactions between single anion/cation and GG, and ion pairs and GG, are taken into consideration. B3LYP/6-31++g** method generally provides reasonable results for the interaction between ILs and GG model, and B3LYP-D3/6-31++g** method gives more desirable calculations in consideration of the dispersion effects and electrostatic interactions. H-Bond interactions are particularly investigated and characterized by interaction energy (ΔE), NBO analysis, AIM analysis and RDG method to judge the strength of interactions. On the basis of geometries, NBO, AIM and RDG analyses, we find that both cations (Mmim, Emim, Bmim, Hmim) and anions (Cl, OAc, MeSO₄, PF₆) can form H-bonds with GG model. The analyses also show that anions have the strongest interactions at the α-OH position of GG and cations have the strongest interactions at γ-OH of GG. Furthermore, anions–GG interactions are stronger than cations–GG interactions, and the changing lengths of alkyl chains of cations do not significantly influence H-bond interactions in all of the cations–GG interactions. Anions with electronegative atoms or conjugated electrons are favorable for H-bonds, and non-coordinating anions do not have strong interactions with lignin. In the ion pairs–GG interactions, anions still play a key role in providing H-bond interactions, and cations can afford a π-stacking interaction with model compound GG and from H-bonds with anions simultaneously. NBO and RDG analysis show the importance of π-stacking interactions in lignin dissolution, but its effect is not as significant as anions do. In all, the calculation data explain the dissolution mechanism to some extent, and H-bond is of great importance to facilitate the dissolution of lignin in ILs, which is in consistent with the experimental evidence. Thus, to design an efficient IL solvent for lignin, anions with good H-bond basicity, smaller volume and forming more H-bonds are favorable for lignin dissolution. The long alkyl substituent may result in higher viscosity and steric effect, thus causing reduction of dissolution ability. Cations with too short alkyl chain also have strong association with anions, thus we need come to a balance between chain length and solubility in practices.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 21406230), the Key Program of National Natural Science Foundation of China (Grant No. 21336002), the Open Research Fund of State Key Laboratory of Multiphase Complex Systems, IPE, CAS (Grant No. MPCS-2014-D-04), and by the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (second phase).

References

1. F. H. Isikgor and C. R. Becer, Polym. Chem., 2015, 6, 4497–4559.
2. V. K. Thakur and M. K. Thakur, Int. J. Biol. Macromol., 2015, 72, 834–847.
3. F. S. Chakar and A. J. Ragauskas, Ind. Crops Prod., 2004, 20, 131–141.
4. B. Kamm, P. R. Gruber and M. Kamm, Biorefineries–Industrial Processes and Products, Wiley Online Library, 2007.
5. B. G. Janesko, Phys. Chem. Chem. Phys., 2011, 13, 11393–11401.
6. G. Cheng, M. S. Kent and L. He, et al., Langmuir, 2012, 28, 11850–11857.
7. X. Qiu, Q. Kong and M. Zhou, et al., J. Phys. Chem. B, 2010, 114, 15857–15861.
8. Y. Deng, X. Feng and M. Zhou, et al., Biomacromolecules, 2011, 12, 1116–1125.
9. S. S. Y. Tan, D. R. MacFarlane and J. Upfal, et al., Green Chem., 2009, 11, 339.
10. N. Sun, H. Rodriguez and M. Rahman, et al., Chem. Commun., 2011, 47, 1405–1421.
11. J. S. Wilkes and M. J. Zaworotko, J. Chem. Soc., Chem. Commun., 1992, 965–967.
12. P. Walden, Bull. Russ. Acad. Sci., 1914, 8, 405–422.
13. A. A. Fannin Jr, D. A. Floreani and L. A. King, et al., J. Phys. Chem., 1984, 88, 2614–2621.
14. H. He, H. Chen and Y. Zheng, et al., Chem. Eng. Sci., 2015, 121, 169–179.
15. L. J. Conceiçao, E. Bogel-Łukasik and R. Bogel-Łukasik, RSC Adv., 2012, 2, 1846–1855.
16. J. G. Huddleston, H. D. Willauer and R. P. Swatloski, et al., Chem. Commun., 1998, 1765–1766.
17. T. Welton, Chem. Rev., 1999, 99, 2071–2084.
18. R. P. Swatloski, S. K. Spear and J. D. Holbrey, et al., J. Am. Chem. Soc., 2002, 124, 4974–4975.
19. J. M. Andanson, A. A. H. Padua and M. F. Costa Gomes, Chem. Commun., 2015, 51, 4485–4487.
20. J.-M. Andanson, E. Bordes and J. Devemy, et al., Green Chem., 2014, 16, 2528–2538.
21. S. Velioglu, X. Yao and J. Devémy, et al., J. Phys. Chem. B, 2014, 118, 14860–14869.
22. H. Wu and E. J. Maginn, Fluid Phase Equilib., 2014, 368, 72–79.
23. S. Budhathoki, J. K. Shah and E. J. Maginn, Ind. Eng. Chem. Res., 2015, 54, 8821–8828.
24. J. L. Anthony, J. L. Anderson and E. J. Maginn, et al., J. Phys. Chem. B, 2005, 109, 6366–6374.
25. J. L. Anthony, E. J. Maginn and J. F. Brennecke, J. Phys. Chem. B, 2002, 106, 7315–7320.
26. C. Cadena, J. L. Anthony and J. K. Shah, et al., J. Am. Chem. Soc., 2004, 126, 5300–5308.
27. G. B. Damas, A. B. A. Dias and L. T. Costa, J. Phys. Chem. B, 2014, 118, 9046–9064.
28. D. A. Fort, R. C. Remsing and R. P. Swatloski, et al., Green Chem., 2007, 9, 63–69.
29. N. Sun, M. Rahman and Y. Qin, et al., Green Chem., 2009, 11, 646–655.
30. S. H. Lee, T. V. Doherty and R. J. Linhardt, et al., Biotechnol. Bioeng., 2009, 102, 1368–1376.
31. A. K. Sangha, L. Petridis and J. C. Smith, et al., Environ. Prog. Sustainable Energy, 2012, 31, 47–54.
32. C. Balaji, T. Banerjee and V. V. Goud, J. Solution Chem., 2012, 41, 1610–1630.
33. A. Lüttge, J. Electron Spectrosc. Relat. Phenom., 2006, 150, 248–259.
34. A. Casas, J. Palomar and M. V. Alonso, et al., Ind. Crops Prod., 2012, 37, 155–163.
35. R. S. Payal, R. Bharath and G. Periyasamy, et al., J. Phys. Chem. B, 2012, 116, 833–840.
36. B. B. Cao, J. Y. Du and D. M. Du, et al., Carbohydr. Polym., 2016, 149, 348–356.
37. W. Ji, Z. Ding and J. Liu, et al., Energy Fuels, 2012, 26, 6393–6403.
38. J. Zakzeski, P. C. Bruijnincx and A. L. Jongerius, et al., Chem. Rev., 2010, 110, 3552–3599.
39. P. Sista, B. Xue and M. Wilson, et al., Macromolecules, 2012, 45, 772–780.
40. M. J. Frisch, G. W. Trucks and H. B. Schlegel, et al., Gaussian 09, Revision D. 01, Gaussian, Inc., Wallingford, CT, 2013.
41. R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728.
42. K. Dong, Y. Cao and Q. Yang, et al., Ind. Eng. Chem. Res., 2012, 51, 5299–5308.
43. K. Dong, Y. Song and X. Liu, et al., J. Phys. Chem. B, 2012, 116, 1007–1017.
44. K. Dong and S. Zhang, Chem.–Eur. J., 2012, 18, 2748–2761.
45. W. G. Cheng, F. Xu and J. Sun, et al., Synth. Commun., 2016, 46, 497–508.
46. Z. D. Ding, Z. Chi and W. X. Gu, et al., Carbohydr. Polym., 2012, 89, 7–16.
47. J. Guo, D. Zhang and C. Liu, J. Theor. Comput. Chem., 2010, 9, 611–624.
48. X. Zhang, W. Yang and W. Blasiak, Energy Fuels, 2011, 25, 4786–4795.
49. S. Grimme, J. Antony and S. Ehrlich, et al., J. Chem. Phys., 2010, 132, 1–19.
50. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
51. T. Lu and F. Chen, J. Mol. Model., 2013, 19, 5387–5395.
52. L. Goerigk and S. Grimme, Phys. Chem. Chem. Phys., 2011, 13, 6670–6688.
53. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
54. R. Fu, T. Lu and F. Chen, Acta Phys.–Chim. Sin., 2014, 30, 628–639.
55. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
56. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
57. H. Guernon and C. Y. Legault, Organometallics, 2013, 32, 1988–1994.
58. B. Miehlich, A. Savin and H. Stoll, et al., Chem. Phys. Lett., 1989, 157, 200–206.
59. R. F. W. Bader, in Encyclopedia of Computational Chemistry, John Wiley & Sons, Ltd, 2002.
60. J. Cioslowski, Science, 1991, 252, 1566–1568.
61. J. Contreras-Garcia, E. R. Johnson and S. Keinan, et al., J. Chem. Theory Comput., 2011, 7, 625–632.
62. E. R. Johnson, S. Keinan and P. Mori-Sanchez, et al., J. Am. Chem. Soc., 2010, 132, 6498–6506.
63. A. Bondi, J. Phys. Chem., 1964, 68, 441–451.
64. W. E. S. Hart, J. B. Harper and L. Aldous, Green Chem., 2015, 17, 214–218.
65. E. D. Glendening, C. R. Landis and F. Weinhold, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 1–42.
66. R. G. Parr, in Horizons of Quantum Chemistry, Springer, 1980, pp. 5–15.
67. P. A. Hunt, C. R. Ashworth and R. P. Matthews, Chem. Soc. Rev., 2015, 44, 1257–1288.
68. B. J. Cox, S. Jia and Z. C. Zhang, et al., Polym. Degrad. Stab., 2011, 96, 426–431.
69. A. L. Ringer and C. D. Sherrill, J. Am. Chem. Soc., 2009, 131, 4574–4575.
70. W. Gao, Y. Tian and X. Xuan, J. Mol. Graphics Modell., 2015, 60, 118–123.
71. J. Zakzeski, P. C. Bruijnincx and A. L. Jongerius, et al., Chem. Rev., 2010, 110, 3552–3599.
72. K. C. Badgujar and B. M. Bhanage, Bioresour. Technol., 2015, 178, 2–18.
73. Y. Pu, N. Jiang and A. J. Ragauskas, J. Wood Chem. Technol., 2007, 27, 23–33.
74. S. E. Wheeler, Acc. Chem. Res., 2013, 46, 1029–1038.
75. C. R. Martinez and B. L. Iverson, Chem. Sci., 2012, 3, 2191–2201.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra27059j

---