Volume 9 - Year 2021 - Pages 14-19
DOI: 10.11159/ijtan.2021.003

Hydrolysis of Cellulose in Supercritical Water: Quantum Simulation

Taketo Oku1, Mitsuhiro Matsumoto1

1Graduate School of Engineering, Kyoto University
Kyoto-Daigaku-Katsura, Kyoto 615-8540, Japan
oku.taketo.57n@st.kyoto-u.ac.jp; matsumoto@kues.kyoto-u.ac.jp

Abstract - Cellulose nanofiber (CNF) is a high-strength nanomaterial made from cellulose fibers. Among several fabrication processes of CNF, we focus on the hydrolysis of cellulose in supercritical water and analyze the reaction mechanism by numerical simulation. In order to deal with the detailed chemical reaction, a series of quantum molecular dynamics simulations were performed based on the density functional theory coupled with the tight binding model. After locating the vapor-liquid critical point with a 100 water molecule system, we explored the hydrolysis reaction of cellulose using a simplified system consisting of a single cellobiose surrounded by 100 water molecules. We observed a cleavage of the 1,4-𝛽-glycosidic bond in some cases. Electric charge analysis suggests that the carbon atom at the cleavage site gives the electron to a water molecule approaching to the bond with sufficiently large velocity.

Keywords: Cellulose nanofiber, hydrolysis, supercritical water, MD simulation, DFTB.

© Copyright 2021 Authors - This is an Open Access article published under the Creative Commons Attribution License terms. Unrestricted use, distribution, and reproduction in any medium are permitted, provided the original work is properly cited.

Date Received: 2020-12-20
Date Accepted: 2020-12-22
Date Published: 2021-03-24

1. Introduction

Cellulose nanofiber (CNF) is a fairly new nanomaterial made from various resources of cellulose. It has recently attracted much attention in many fields [1], such as electronic devices, medical material, and food additives. One of the significant features is its high aspect ratio; a typical width is 4-100 nm while the length is often more than 5 μm. It is so light and tough that many applications as structural material are proposed, such as automobile bodies. Cellulose is a renewable resource since it is mainly made from wood pulp.

To fabricate CNF from conventional cellulose resources (e.g., wood pulp and grass), some treatment is required to smash macromolecules of cellulose. Currently exist two major fabrication methods: (i) Chemical treatment [2][3], which uses hydrolysis reactions with strong acid or catalysts. This is not environment-friendly nor cost-effective. (ii) Mechanical approach [4], in which the pulp is mechanically smashed or ground under high shear force. It generally gives fiber size larger than 15 nm. Less popular treatment is hydrolysis using supercritical water [5], which is environment-friendly and has a possibility of precise size control, but its detailed mechanism is not fully understood.

In this study, we focus on the cellulose hydrolysis in supercritical water. To analyze the chemical reaction in atomic scale, we executed a series of molecular dynamics (MD) simulations with quantum effects taken into account as the density functional tight binding (DFTB) model.

2. Methods

In the density functional (DF) theory, the wave function of single electron  and the eigen value (orbital energy) is obtained self-consistently from the Kohn-Sham equation,


where V(r) is the potential defining the interaction between an electron and the collection of atomic nuclei, VH(r) the potential for the Coulomb repulsion between the electron and the total electron density, and VXC(r) the potential for everything else.

In addition, we utilize the tight binding (TB) model in which every single-electron wave function is expressed as a linear combination of atomic orbitals,


where ϕμ(r) is the atomic orbital.

We have utilized the DFTB+ software [6] to perform the ab initio calculation. A self-consistent charge (SCC) DFTB calculation is executed with the third order correction [7]. The 3ob-3-1 parameter set [8] is adopted, which is known to be well applicable to a wide range of organic and bio molecules. Use of appropriate hydrogen bond correction is relevant to our investigation, and we adopted a damping method for modification of short-range potential energy U [7], as



where rab is the distance between atom a and b, Δqi the electric charge on atom i, Ui the Hubbard parameter of atom i, and S a short-ranged damping function [7].

Assuming the Born-Oppenheimer approximation, the equation of motion for each atom is integrated using the velocity Verlet dynamics with the Nosé-Hoover chain thermostat [9] algorithm. We used the 1 × 1 × 1 Monkhorst-Pack grid [10] for k-point sampling (the Γ point calculation).

3. Evaluation of critical point

Before starting the main simulations of hydrolysis, the critical point of the model water should be determined. Using 100 water molecules confined in a cubic simulation cell of various sizes with periodic boundary conditions, we performed the DFTB-MD simulations under constant temperature conditions. Three temperature conditions were investigated, 400, 650 and 800 K. The time step was chosen to be 0.5 fs. Typical snapshots are shown in Figure. 1.

(a) T= 650 K, ρ = 0.886 g/cm3 (cell size: 15 ×15 × 15 A3)

(b) T= 650 K, ρ= 0.191 g/cm3 (cell size: 25 ×25 × 25 A3)

Figure 1. Example of the pure water system with 100 molecules to determine the critical point.

After equilibrating the system, the mean pressure was evaluated. The obtained P- V curve is shown in Figure. 2, which suggests that the liquid-vapor critical point of this model water exists around (Tc,Vc ,Pc ) = (650 K , 80 A3, 21.9 MPa); 80 A3 per molecule corresponds to the density ρ ∼ 0.374 g/cm3. A reported experimental value is (647 K, 0.322 g/cm3, 22.1 MPa) [11], and the agreement is reasonable.

Figure 2. p-V relation for the pure water system.

To look at the fluid structure near the critical point, the oxygen-oxygen radial distribution function  is calculated at V= 80 A3. Based on the results shown in Figure. 3, coexisting of liquid and vapor states is suggested at 400 K, while uniform gas phase is realized at T = 800 K. A snapshot at T = 400 K (Figure. 4a) indicates the formation of a huge cluster. On the other hand, water molecules spread homogeneously at T = 800 K (Figure. 4b). Thus we suppose that the state in between (650 K) is close to the critical point. Based on these results, we performed the main simulation of hydrolysis at T≥ 650 K.

Figure 3. O-O radial distribution function of fluid water with  = 80 A3 per molecule.

(a) T = 400 K.

(b) T = 800 K.

Figure 4. Snapshots of the pure water system; V = 80 A3 per molecule, corresponding to the density ρ = 0.374 g/cm3

4. Hydrolysis Simulation of Cellulose

We performed DFTB-MD simulations for the hydrolysis cellulose in water at supercritical states. To reduce the computational resource, a minimum unit of cellulose, i.e., cellobiose (Figure. 5), was targeted. We investigated a system of a single cellobiose molecule surrounded by 100 water molecules (Figure. 6). The cell size was 17.197 × 17.197 × 17.197 A3 (0.70 g/cm3) and the periodic boundary conditions were assumed for all directions. Constant temperature MD simulations were carried out at four temperatures, T = 400, 650, 800 and 1000 K. A shorter time step of 0.25 fs was used to see the trajectories more precisely.

Figure 5. Investigated model of cellulose; cellobiose.
Figure 6. Snapshot of the system for hydrolysis simulation.

For the hydrolysis of large cellulose molecules, cleavage of the 1,4-β-glycosidic bond, shown in Figure. 5, is essential. Since this is a stochastic process, we counted the number of cleavage cases out of ten trials with different initial configurations; each trial run is 12.5 ps (50,000 steps). The results are summarized in Table 1, which indicates that the hydrolysis of cellobiose is more probable near the critical point of fluid water.

Table 1. Summary of simulation results. Each trial is a 12.5 ps MD simulation.

Temperature [K] Trial Bond cleavage
400 10 0
650 10 2
800 10 3
1000 10 1

Figure 7 shows a series of snapshots of the bond cleavage at T = 1000 K. After t = 0.25 ps, the bond breaks while a weak carbon-carbon bridge via a hydrogen atom is newly formed, temporarily maintaining the dimer. In Figure.8 comparison is made between T = 400 and 1000 K cases for the atomic distances; the glycosidic bond between C and O atoms is indicated with black lines, while the distance between the glycosidic O and the nearest hydroxyl O is shown with red, and the glycosidic O and the nearest water O with blue. At 400 K, water molecules cannot get close to C because the adjacent hydroxyl group blocks it. Under higher temperature such as 1000 K, some of the water molecules have high kinetic energy (large translational momentum) and are enabled to approach the bond, as seen in Figure.8b, leading to the bond breakup.

Figure 7. Example of the bond breakup process at T = 1000 K.

(a) T = 400 K.

(b) T = 1000 K.

Figure 8. Example of atomic distance change. C-O: Between glycosidic O and glycosidic C. O-O (hydroxy): Between glycosidic O and the nearest hydroxyl O. O-O (water): Between glycosidic O and the nearest water O.

Figures 9 and 10 show the change in the number of electrons (Mulliken charge) assigned to the C(4’) atom and the nearest neighbor water molecule to the bond. Figure 9 indicates that the Mulliken charge decreases during the bond cleavage. This decrease of electrons is roughly compensated by the increase for the nearest water molecule (Figure.10). Thus we conclude that the water molecule approaching the bond with high kinetic energy deprives the C(4’) partially of electrons and initiates the bond breakup.

Figure 9. Mulliken charge on the bonding C(4’) atom for the case of Figure.8 (b).
Figure 10. Number of electrons assigned to the nearest neighbor water molecule for the case of Figure.8 (b).

4. Conclusion

We examined a model cellulose in fluid water with quantum molecular dynamics simulations to study the mechanism of hydrolysis under supercritical conditions. In the normal states, the 1,4-β-glycosidic bond is protected by hydroxyl groups from surrounding water molecules. Under high temperature conditions, the number of water molecules with larger kinetic energy increases so that some of them can approach the glycosidic bonds of cellulose molecules. The water molecule close to the bond attracts an electron from the carbon atom, leading to the bond breakup.


[1] A. Sharma, M. Thakur, M. Bhattachara, T. Mandal and S. Goswami, “Commercial application of cellulose nano-composites – A review”, Biotechnology Reports, vol. 21, 2019. View Article

[2] D. Klemm, F. Kramer, S. Moritz, T. Lindstrom, M. Ankerfors, D. Gray and A. Dorris, “Nanocelluloses: A New Family of Nature-Based Materials”, Angewandte Chemie International Edition, vol 50, pp. 5438-5466, 2011. View Article

[3] T. Saito, Y. Nishiyama, J. L. Putaux, M. Vignon and A. Isogai, “Homogeneous Suspensions of Individualized Microfibrils from TEMPO-Catalyzed Oxidation of Native Cellulose”, Biomacromolecules, vol. 7, no. 6, pp. 1687-1691, 2006. View Article

[4] K. Abe, S. Iwamoto and H. Yano, “Obtaining Cellulose Nanofibers with a Uniform Width of 15 nm from Wood”, Biomacromolecules, vol.8, no. 10, pp. 3276-3278, 2007. View Article

[5] M. Sasaki, Z. Fang, Y. Fukushima, T. Adschiri and K. Arai, “Dissolution and Hydrolysis of Cellulose in Subcritical and Supercritical Water”, Industrial and Engineering Chemistry Research, vol. 39, no. 8, pp. 2883-2890, 2000. View Article

[6] B. Aradi, B. Hourahine and Th. Frauenheim, “DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method”, Journal of Physical Chemistry A, vol. 111, no. 26, pp. 5678-5684, 2007. View Article

[7] M. Gaus, Q. Cui and M. Elstner, “DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB)”, Journal of Chemical Theory and Computation, vol. 7, pp. 931-948, 2011. View Article

[8] M. Gaus, A. Goes and M. Elstner, “Parametrization and Benchmark of DFTB3 for Organic Molecules”, Journal of Chemical Theory and Computation, vol. 9, pp. 338-354, 2013. View Article

[9] G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, “Explicit reversible integrators for extended systems dynamics”, Molecular Physics, vol. 87, no. 5, pp. 1117-1157, 1996. View Article

[10] H. J. Monkhorst and J. D. Pack, “Special points for Brillouin-zone integrations”, Physical Review B, vol. 13, no. 12, pp. 5188-5192, 1976. View Article

[11] W. Wagner and A. Pruß, “The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use”, Journal of Physical and Chemical Reference Data, vol 31, no. 2, pp. 387-535, 2002. View Article