Esterification of citric acid (CA) with the primary alcohols and the hydroxyl groups of cellulose chain (n = 1-2) in parched condition were investigated by using density functional theory (DFT) method and a two-layer ONIOM approach. Geometry and energy of reactants, products, and transition state (TS) structures were optimized at B3LYP/6-311g (d, p) level and ONIOM (B3LYP/6-311g (d, p):PM3MM) level. The computational results show that the esterification occurs in the two main steps: the first step is the dehydration reaction of CA to form anhydrides of 5-membered ring and 6-membered ring and the second step is the ring opening reaction with the hydroxyl –OH groups to form the ester products. The energy barrier of dehydration reaction step is much higher than that of ring opening reaction step. Effect of substituent R in primary alcohol R-CH2OH (R: CH=CH2, CH2NHCH3, CH2OCH3, CH2Cl) and cellulose chain (1G, 2G) on the reactivity, which has negative inductive effect –I, is significant. The combination of calculation data and experiment data were applied to make findings more rigorous. The activation energy of CA was determined by using differential scanning calorimetry (DSC) and thermal gravimetric (TG) analysis to be  = 47.8 kcal/mol; the experimental data favoured the dehydration reaction step of CA.

1. Introduction

Esterification is an important reaction in the chemical industries because the ester products are used in the many various applications such as solvents, flavors, plasticizers, pesticides, and emulsifiers [1]. Condensation reaction between a carboxylic acid and a primary alcohol in homogeneous acid catalysts such as sulfuric acid at 110°C–130°C or heterogeneous acid catalysts such as ion exchange resin at 50°C–80°C is often used to obtain organic ester because of simplicity, economy, and environmental benefit. For instance, the esterification reaction between ethyl alcohol and acetic acid in sulfuric acid catalyst has been investigated in a batch reactor; the order of forward reaction was found to be 1 [2], or a novel two-step process for the coproduction of ethyl and butyl acetates obtained higher purity using ion exchange resin (Dowex 50WX8) [3]. Furthermore, several innovative catalysts were also studied aiming to reach high selectivity, high efficiency, and mild temperature. For instance, the sulfated zirconia catalyst was used in the reaction of acetic acid and n-pentanol in batch reactor and the estimated enthalpy of the reaction was −2.03 kcal/mol and acetic acid adsorbed on the catalyst as the limiting step of the reaction [4], or a structured catalyst was developed from SiC foam materials with a ZSM-5 zeolite coating in the synthesis of n-amyl acetate [5].

On the other hand, some derivatives of carboxylic acid such as acyl halides and organic acid anhydride have higher reactivity than that of carboxylic acid; therefore, they can react with the primary alcohols in high efficiency and high selectivity without the presence of strong acid catalyst and solvent. The reactivity of derivatives with primary alcohols reduces in the sequence as follows: acyl halides > organic acid anhydride > carboxylic acid [6, 7]. CA, its salts, and esters are the important products in chemical industries such food, cosmetic, and pharmaceutical. CA is a special polycarboxylic acid; it can pyrolyze to form CA anhydrides; therefore, it can react with pentadecylic, dodecylic, and other alcohols at 120°C–180°C without the presence of strong acid catalysts and solvents [8].

Cellulose is a main component of the lignocellulosic materials. Its basic units are D-glucopyranose, which are linked together by β-(1 ⟶ 4)-glucosidic bonds. Structure of cellulose is a repeated two-sugar unit. Nucleophile sites are the hydroxyl –OH groups at positions of C6 and C2, especially the hydroxyl group of C6 position because of the twist structure of cellulose chain, hydrogen bond, and unfavourable rearrangement at positions of C2 and C3 when the reaction occurs. Surface degrees of substitution are in the range from 0.06 to 1.5 when using esterification reaction [911]. The use of CA as a cross-link agent was carried out to change properties and functions of the lignocellulosic materials in parched condition. For instance, Ramachandran et al. modified cotton fabric with CA to enhance the ability of crease resistant [12], or cotton textiles modified with CA were used as an efficient antibacterial agent for the prevention of nosocomial infections [13].

Furthermore, the modified lignocellulosic materials reported the great ability to exchange heavy metal cations in aqueous solution. The outlet water can meet the permissible standards of domestic water. Therefore, they can replace the commercial cationite resins in the treatment of groundwater by using ion exchange method. The modification work was appreciated because of environmentally friendly and saving water treatment cost [14, 15]. For instance, Marshall et al. connected CA to soybean hulls in parched condition to adsorb and exchange Cu (II) ions in aqueous solution [16], or orange peel was modified by CA to remove some heavy metal cations such as Cd (II), Zn (II), Co (II), and Ni (II) effectively [17]; Diaz-Munoz et al. modified avocado kernel seeds using CA to treat some heavy metal cations such as Cd (II), Cu (II), Ni (II), Pb (II), and Zn (II) in aqueous solution [18]. These heavy metal ions caused many adverse effects on living environment and ecosystem as well as neurotoxicity, genetic mutation, and cancer [19, 20].

To the best of our knowledge, the esterification mechanism of CA with the primary alcohols and the hydroxyl groups of cellulose in parched condition has not been reported before (seen in Figure 1). The esterification of CA with the hydroxyl –OH groups is predicted to occur in the two main steps: first, CA is dehydrated to form the anhydrides of 5-membered ring and 6-membered ring; second, the intermediate products join in the ring opening reaction with the hydroxyl groups to form the ester products [2123]. DFT method can reach fast converge with the acceptable accuracy level when comparing with other quantum mechanics methods; therefore, it helps to save much time and computer resource. However, DFT method is just applied effectively for small and medium chemical systems [24, 25]. Regarding large and complex chemical systems such as biomolecular system, scientists apply the computer power in which it is needed. At the heart of the system, calculations are based on quantum mechanics (QM). Beyond the scope of this area, they are based on molecular mechanics (MM) and, at the outermost layers, atoms and molecules are even lumped together to a homogeneous mass [2629]. Morokuma and coworkers developed a hybrid method, which can combine any computational method of QM and MM, called ONIOM method. ONIOM method can, in principle, be easily expanded to multiple layers [30, 31]. In this study, DFT method and the two-layer ONIOM approach were combined to probe the esterification mechanism.

2. Methodology

2.1. Computational Process

In the two-layer ONIOM method, the target total energy is the high level energy of the real system, which could be approximated by the following equation:

The chemical formula of reactants, products, and TS structures in the four reaction pathways (seen in Figure 2) were built with GaussView 03 software package and optimized with Gaussian 03 software package. In the scope of this study, B3LYP/6-311g (d, p) level is used for calculation in the case of primary alcohols and B3LYP/6-311g (d, p):PM3MM level is used for calculation in the case of cellulose chain (n =1-2) with the acceptable accuracy level in a reasonable computational effort.

The two keywords are mainly applied in the calculation process to be opt and freq. Vibration analysis was carried out to confirm TS structures and imaginary frequencies. Intrinsic reaction coordinate (IRC) method was also used to check again the energy pathways. Ea (kcal/mol) is the energy barrier of activation step.

2.2. Thermodynamic Functions, Rate Constant, and Solvent Model

The statistical thermodynamic method and Eyring transition state theory were used to calculate the thermodynamic functions and rate constants [32, 33].where k (T), T, , kb, h, and ѵ are reaction rate constant, Kelvin temperature, Wigner emendation factor, Boltzmann factor, Planck’s constant, and imaginary frequency of TS structure, respectively. ∆H, ∆S, and ∆G are a mole activation enthalpy, mole activation entropy, and mole activation Gibbs function of a reaction step.

Mole percent ratio of products A5 and A6 is determined based on the value of k (T) and the diagram of first-order parallel reaction [33].

where a is the mole amount of reactant CA at starting time; k1, k2 and X1, X2 are reaction rate constant and mole percent ratio to form products A5 and A6, respectively.

Onsager reaction field model is often applied to probe the effect of solvent because of simple model and fast converge; where solute occupies a fixed spherical cavity of radius a0 within the solvent field, solvent dipole and molecular dipole will interact with each other and leading to net stabilization (as seen in Figure 3) [34]. The two keywords are added into the calculation process to be volume and SCRF = dipole.

2.3. Experimental Process

The lignocellulosic materials such as cotton, Acacia auriculiformis sawdust, and coconut shell powder are modified as follows: a mixture of 0.2 mol/L NaOH/mixed solvent of 70% v/v ethanol and 30% v/v water was used to remove lignins in the initial raw materials. The delignified materials were immersed in 200 g/L CA solution for 24 h so that the solution could infiltrate deeply between the fibers and then were filtered out and dried at room conditions until their moisture content was lower than 15 wt (%); finally, the activated materials were rinsed thoroughly with distilled water, dried, and stored in room conditions after they were activated at 120°C for 3 h in parched condition (as seen in Figure 4).

The morphological characteristics of modified materials were observed by using HITACHI S4800 scanning electron microscope (SEM) machine; specific surface area of materials were measured by using Nova Station A instrument with Brunauer–Emmett–Teller (BET) method and the nitrogen gas. The functional groups were characterized by Fourier transform infrared spectroscopy (FT-IR) by using TENSOR 37 machine. Differential thermal analysis (DTA) and thermal gravimetric analysis (TG) of citric acid monohydrate were carried out in a temperature range from 50°C to 250°C with the speed of 2°C/min.

3. Results and Discussion

3.1. Esterification of Citric Acid with the Primary Alcohols in Parched Condition
3.1.1. Esterification of CA with Methanol

When the substituent R in R-CH2OH is H, the chemical formula of primary alcohol is methanol. The four possible parallel reaction pathways between CA and methanol in parched condition are presented in detail as follows:

Potential energy surface (PES) was established by setting the energy level of CA as the basepoint energy level; other energy levels of the reactant, products, and TS structures are determined comparatively towards the zero energy level. The computational results are summarized and shown in Figure 5. The energy level diagram has the two TS structures, its shape like a double-hill. The estimated enthalpy of the reaction is in a range from −2.1 kcal/mol to −2.5 kcal/mol; the esterification is a slightly exothermic and domino reaction.

Esterification mechanism of CA with methanol in parched condition could be described in short as follows. In the first reaction, CA is dehydrated by heat to form the very reactive anhydrides of 5-membered ring A5 and 6-membered ring A6 through TS structures TS_A5 and TS_A6; pathway (6) to form product C did not occur because the energy barrier of Ea3 is much higher than that of Ea1 and Ea2; note that the cyclic anhydrides are not intermediates of a stepwise reaction but the stable derivatives of CA. In the second reaction, anhydrides of A5 and A6 will take place at the ring opening reaction with methanol to form the ester products. The two ester products of A52 and A61 will compete each other in the esterification because the energy barrier of Ea4, Ea5, and Ea6 is not so much different.

3.1.2. Effect of Substituent R in the Primary Alcohols by Using DFT Method

Effect of substituent R on the reactivity of primary alcohols is very large, because they will modify the nucleophilic character of the corresponding alcohols with respect to that of methanol; therefore, it should be probed in depth. The change in the energy barrier of ring opening step was investigated when the substituent R was altered between electron-releasing (ER) groups such as H, CH3, and CH2CH3 based on positive inductive effect +I. The result reveals that the change in substituent R, which has +I effect, impacts negligibly the energy barriers of ring opening step (as seen in Table 1).

When the substituent R was altered between electron-withdrawing (EW) groups such as CH=CH2, CH2NHCH3, CH2OCH3, and CH2Cl based on negative inductive effect –I, the result showed that the energy barriers of ring opening step are affected significantly. It can be explained by the energy barrier of Ea which is related to negative charge density on O-atom in the hydroxyl group of alcohol; negative charge density on O-atom is less and alcohol molecule is more difficult to insert into the cyclic anhydride to form the ester product. Therefore, the presence of EW groups will make a rise in the energy barrier of Ea. The strength of the EW groups increases respectively with the strength of -I effect in the sequence as follows: CH=CH2 < CH2NHCH3 < CH2OCH3 < CH2Cl. However, the energy barrier of Ea still depends mainly on the geometrical structure of cyclic anhydrides in the sequence as follows: Ea6 < Ea4 < Ea5 (as seen in Figure 6).

3.2. Thermodynamic Functions, Rate Constant, and Onsager Model of the Dehydration Step

Repta and Higuchi were able to synthesize and isolate the monomolecular unsymmetrical anhydride of CA; the identification of this compound was investigated by elementary analysis, potentiometric titration, and cryoscopic and NMR measurements [37, 38]. They postulated that the unsymmetrical anhydride was mainly formed; it was based on the fact observation about reactivity and soluble level of CA anhydride [38]. According to our heretofore results, CA was dehydrated to form the cyclic anhydrides and the products were then dehydrated to form the aconitic acid anhydride. The energy barrier of C=C bond forming step is higher, 10 kcal/mol, than that of cyclic anhydrides forming step [39].

According to the prior literature, the DSC–TG spectra of CA monohydrate confirmed that CA released the two water molecules and did not decarboxylate in the decomposition process. However, the use of DSC-TG method obtained a mixture of main products such as aconitic acid, itaconic acid, and anhydrides. The different results depended on some main parameters such as temperature range, heating rate, and catalysts. Therefore, the authors, who use the DSC-TG spectra, often postulate the products of decomposition process in the following sequence: CA ⟶ aconitic acid ⟶ aconitic acid anhydride [40, 41].

The dehydration step of CA is the rate-determining reaction of this domino process. Rate constants k (T) and mole percent ratio X of the first reaction at 120°C are presented in Table 2. The main product of dehydration step is the anhydride of 6-membered ring with mole percent ratio of X2 = 99.5% and rate constant of k (393) = 2.0 x 10−11 s−1.

Energy barriers of Ea in gas-phase, activation Gibbs functions of at various temperatures, and energy barriers of in water solvent are presented in Table 3. The activation energy barrier from the calculated data is in good agreement with the activation energy from the experimental data [41]. The results reveal that the values of Ea, , and just have small energy gap when a range of reaction temperature is not so high. If we ignore the effect of hydrogen bond and the size of cyclic anhydride being formed, the effect of water solvent on the energy barrier of dehydration step using Onsager model is not significant in the case of self-catalysis esterification.

3.3. The Use of ONIOM Model for the Hydroxyl Groups of Cellulose Chain (n = 1-2)

The first highlighted line in the Table 4 belongs to the case of methanol above. When comparing with 1G model, the energy barrier of ring opening step has a small increase of 2-3 kcal/mol. The result could mainly be due to the -I effect of D-glucose. When changing integer n from 1 to 2, the energy barrier of ring opening step tends to increase. It could be explained by steric effect and hydrogen bond.

3.4. Evaluate and Refer to Experimental Data and Calculational Data

According to our experimental data, a chemical linkage was formed between the –COOH group of CA and the –CH2OH groups of cellulose chain, and the esterification occurs through CA anhydrides. DTA–TG spectra of citric acid monohydrate show two peaks at 146.6°C and 191.5°C, and the mass change in the temperature range is −17.89% wt (or 37.4 g/mol). The two water molecules were released from citric acid in this stage [39]. On the other hand, FT-IR spectra of the modified materials show the strong vibration peak of the C=O bond in the -COO ester group at wavenumber of 1741 cm−1 in the case of modified coconut shell powder [36] and at wavenumber of 1737 cm−1 in the case of modified cotton, but it does not show the speciality peak of the C=C double bond. In addition, specific surface area of initial cotton is measured by using BET method with nitrogen gas to be 724.955 m2/g, while the modified cotton has the value of 0.606 m2/g. SEM images of cotton before and after modification show that the diameter of fibers is unchanged and dispersed in the range from 10 to 15 μm. The surface of modified cotton is heterogeneous. Many white colour dots exposed the outside surface through the crevices. The negatively charged functional groups existed on the cellulose chain after modification process; these –COOH groups play the main role in the strong adsorbability of CA-modified lignocellulosic materials towards toxic metal cations in aqueous solution such as Fe3+, Pb2+, and Ni2+.

According to the prior literature, the use of CA to modify the lignocellulosic materials is more effective than that of other dicarboxylic acids in parched condition [42]; and ion exchange capacity of CA-modified lignocellulosic materials is reduced significantly when the reaction time prolongs [43]. The result could be explained by the two reasons: first, the –COOH groups of CA are distributed to comply with a rule in order to form the anhydride of 5-membered ring or the anhydride of 6-membered ring; second, the cross-link reaction of CA makes the amount of –COOH groups reduce when reaction time increases. Noordorver et al. assessed the reactions of CA with several primary alcohols and secondary diols by using NMR and FI-IR spectra with the aim to gain more insight into the mechanism of the chemical modification. The results show that CA derived anhydrides are already formed at temperatures between 140 and 160°C and the esterification reaction occurs through CA anhydrides. The temperature range is high enough to form reactive anhydrides but low enough to prevent the extensive formation of C=C unsaturated structures [44].

On the other hand, Yu et al. used DFT method to probe the direct esterification mechanism between succinic acid and ethylene glycol in the absence of foreign catalyst and solvent at 123°C. The reaction pathways were studied to be involve in the self-catalysis (acid, alcohol, and water) and non-self-catalysis. The result showed that the reaction pathways without self-catalysis had the activation energy barrier above 31.1 kcal/mol and rate constant closest to the experimental data [32]. Besides that, Lawal et al. investigated the esterification reaction of acetic acid with methanol in uncatalyzed condition using DFT method. The free energy of activation in uncatalyzed condition is 36 kcal/mol, and the TS structures of 4-membered ring and 6-membered ring formed during the esterification process [45]. The calculational data and experimental data also favoured the results of this study.

4. Conclusions

The esterification mechanism of CA with primary alcohols and cellulose chain (n = 1-2) in parched condition occurred in the two main steps: dehydration step and ring opening step. The effect of substituent R owning –I effect on the reactivity is significant. The CA anhydrides of 5-membered ring and 6-membered ring are formed at 120°C, which were supported by the experimental observation. These results contributed the initial basic knowledge to investigate the other effects such as the degree of polymerization (DP) and hydrogen bond.

Data Availability

All data used to support the findings of this study are available in the article and in the supplementary materials.

Conflicts of Interest

The authors declare no conflicts of interest.


This research was funded by Ho Chi Minh City University of Technology (HCMUT), VNU-HCMC, under Grant no. BK-SDH-2020-51205011. The authors acknowledge the support of time and facilities from Ho Chi Minh City University of Technology (HCMUT), VNU-HCM, for this study.

Supplementary Materials

Table S1: energy level of substances, products, and transition states of the dehydration reaction step. Table S2: energy level of substances, products, and transition states of the ring opening reaction step using 1G model. Figure S1: DTA-TG spectra of citric acid monohydrate. Figure S2: digital photographs of initial cotton and CA-modified cotton; SEM images of initial cotton and CA-modified cotton with different magnifications. Figure S3: FT-IR spectra of initial cotton and CA-modified cotton. Figure S4: FT-IR spectra of initial coconut shell powder and CA-modified coconut shell powder. Figure S5: specific surface area of initial cotton. Figure S6: specific surface area of CA-modified cotton. (Supplementary Materials)