Thermal Response in Cellulose Iβ Based on Molecular Dynamics

The structural details of cellulose Iβ were discussed according to molecular dynamics simulations with the GLYCAM-06 force eld. The simulation outcomes were in agreement with previous experimental data, including structural parameters and hydrogen bond pattern at 298 K. We found a new conformation of cellulose Iβ existed at the intermediate temperature that is between the low and high temperatures. Partial chain rotations along the backbone direction were found and conformations of hydroxymethyl groups that alternated from tg to either gt or gg were observed when the temperature increased from 298 K to 400 K. In addition, the gg conformation is preferred than gt. For the structure adopted at high temperature of 500 K, major chains were twisted and two chains detached from each plain. In contrast to the observation under intermediate temperature, the population of hydroxymethyl groups in gt exceeded that in gg conformation at high temperature. In addition, three patterns of hydrogen bonding were identi ed at low, intermediate and high temperatures in the simulations. The provided structural information indicated the transitions occurred around 350 K and 450 K, considered as the transitional temperatures of cellulose Iβ in this work.


Introduction
The research tendency in cellulose and its applications has increased extensively during the past few years As an important natural polymer to industrialization, cellulose can be extracted from various plants and has been widely used to produce materials such as textiles, pulp, biofuels, etc. (Liu et al. 2011;Mazeau 2005;Wyman 2007). In addition, cellulosic ber and regenerated cellulosic ber are principal raw materials for textile industry and medical application (Mazeau 2005). Therefore, research in the e ect of temperature on the structure and mechanical properties of the cellulose crystal can be signi cant not only to understanding the fundamental mechanisms of the complex cellulosic self-assembly formation, but also for redesigning the sustainable regenerated cellulosic materials with desired performance (Bergenstrahle et al. 2007; Zhang et al. 2011).
Structure wise, cellulose is linear-chain polymer that consists of D-glucose repeat units. The hydroxyl groups between adjacent chains contribute to the hydrogen bond patterns in crystalline structure. Divergent morphologies of cellulose crystals were adopted in nature (Atalla and Vanderhart 1984), and cellulose with crystalline form Iβ constitutes the secondary wall of cotton ber, which is one of the most abundant nature bers on earth (Horii et al. 1987; Sugiyama et al. 1991). Cellulose crystal in Iβ form has a typical structure in which cellulosic chains are arranged to form parallel sheets (Gardner and Blackwell 1974). The sheets are stabilized by hydrogen bonds while the stacking of the sheets are stabilized by both Van der Waals forces and weak intersheet hydrogen bond interactions (Nishiyama et al. 2008; Nishiyama et al. 2002). In a polymersolvent system, the hydrogen bonds are critical to the dissolution and stability of the polymer molecules (Cai et al. 2012;Chen et al. 2011;Makhatadze and Privalov 1992;Pace 1986;Swatloski et al. 2002;Zou et al. 2002). Therefore, the computational and experimental works on cellulose currently concentrated on the analysis of the hydrogen bonds in order to study the mechanical properties of cellulosic materials (Cai et  Among all the physical properties of cellulose, thermal behavior is one of the most critical properties for the material processing and research has been done in studying cellulose thermal stability with a temperature swing (Zhang et al. 2011). Cellulosic materials were subjected to dynamic thermogravimetry analysis in both nitrogen and air by Huang et al. (Huang and Li 1998), giving that an approximated range of initial degradation temperatures of cellulose is between 546 and 612 K. The measured values of initial thermal degradation temperature and the temperature where the maximum degradation rate occurs were used to describe the thermal stability of the cellulose and its esters. The thermal behavior of cellulose Iβ has been studied experimentally by Wada et al. (Wada 2002;Wada et al. 2010) and Watanabe et al. (Watanabe et al. 2006a;Watanabe et al. 2006b) using X-ray di raction and temperature-dependent infrared spectroscopy (IR) measurements. According to the temperature dependent IR measurements, thermolysis of cellulose in Iβ crystalline form has not been detected at the temperature from 303K to 533K, although signi cant changes of hydrogen bonds were observed (Watanabe et al. 2006a). These divergent experimental results suggest a necessary investigation of the thermal response of cellulose Iβ crystals, the existence of transition temperature, the temperature-dependent hydrogen bond formation and the hydrogen bond patterns.
However, experimental studies have limitations in tracking the thermal response of cellulose structure along a temperature change at molecular level. Instead, molecular dynamics (MD) simulations can act as powerful tools in monitoring the structural alternation of cellulose under heating. The previous works have studied the exibility and interaction strength of cellulose molecules at ambient conditions using MD simulations (Nishiyama et al. 2008;Zhang et al. 2011). In MD simulation, the accuracy of outcomes is highly relied on the selection of force eld. Stortz and his co-workers (Stortz et al. 2009) compared empirical force elds and the semi-empirical quantum methods when studying carbohydrates, such as β-cellobiose, α-maltose, and α-galabiose. They concluded that the force elds GLYCAM-06, GROMOS, and MM3 are preferred with above carbohydrate molecules. In addition, the GLYCAM-06 force eld is consistent and transferable for modeling carbohydrates (Kirschner et al. 2008 (Spiwok et al. 2010) found that the results of the surface of free energy of β-D-glucopyranose using that force eld were agree with the results of quantum metadynamics studies, providing a higher accuracy of the simulation process.
Here in the present work, we applied the GLYCAM-06 force eld to study the thermal behavior of cellulose Iβ with a temperature swing from 298 K to 550K. We carried out MD simulations with 36 cellulose chains under the NPT ensemble. Beyond the room temperature phase and high temperature phase of cellulose Iβ that were discussed in literature (Wada et al. 2010), our simulation results suggested the existence of an intermediate structure of cellulose Iβ. At this intermediate state, the hydroxymethyl groups adopt either gt or gg conformation and twisting was observed within some polymer chains. To fully understand the thermal behavior of cellulose, further analysis of the structure-thermal property correlation of cellulose Iβ was carried out.

. Initial structure and conformational parameters
Cellulose elementary bril with 36 cellulose chains ( Fig. 1a) proposed by Ding and Himmel (Ding and Himmel 2006) was used to build the model of cellulose Iβ here. The raw atomic coordinates and the positions of hydroxyl groups are assigned according to the crystallographic unit cell data and dominant hydrogen bonding pattern, which was determined by X-ray and neutron ber di raction re nement experimental methods (Nishiyama et al. 2002). The initial structure was generated by a user-friendly program named cellulosebuilder (Gomes and Skaf 2012). The unit cell belongs to the monoclinic space group P2 and conformational parameters are a = 7.784Å, b = 8.201Å, c = 10.38Å and γ = 96.5 • , where c is the chain direction and unique axis. Cellulose chains are packed parallel and each chain contains eight glucoses units.
The conformational parameters of the cellulose chains are shown in Fig. 1(b). The unit cell a, b, c, α, β, and γ are commonly used to describe the crystal structure. The ab projection of one of considered unit cells is shown in Fig. 1a. The unit cell c represents the length of one cellobiose unit alone chain direction. The unit α is the angle between directions of b and c, whereas β is the angle between directions of a and c. Thermal characters are evaluated according to thermal variations of these unit cell parameters. Torsion angle O5-C5-C6-O6 at the C6 primary alcohol group is de ned as ω (Fig. 1b). The dihedral angle ω has three staggered positions at 60 • , 180 • , and 300 • and the conformations of the hydroxymethyl groups are labeled gt, tg and gg respectively.

. System optimization
The solvent box was set as a cube with a given distance of 12 Å and the TIP3P water model was used (Jorgensen et al. 1983). Non-bonded interactions were calculated with a cuto radius of 12 Å. GLYCAM-06 force eld was applied in this study and it has been extensively used in simulations of carbohydrate (Kirschner et al. 2008). Some previous works also shown that the force eld GLYCAM-06 was suitable for the MD simulation at high temperature (Bergenstrahle et al. 2007;Zhang et al. 2011).
The energy minimization procedure consists of two stages, rst for the solvation only then for the whole system. In the energy minimization of solvation, 2000 steps of steepest descent minimization were performed followed by 3000 steps of conjugate gradient minimization. Prior to the minimization of energy of the solvated system, cellulose molecules were xed with positional restraints. To fully optimize the structure, the energy of entire system is minimized with 2000 steps of steepest descent minimization and then 3000 steps of conjugate gradient minimization.

. Molecular dynamics simulation
To investigate the thermal response of cellulose in crystalline Iβ form, we carried out MD simulations at different temperatures (298 K, 350 K, 400 K, 450 K, 500 K and 550 K). All simulations were conducted by using Amber software (D.A. Case et al. 2012). The equilibration process was performed prior to the production therefore the system can be heated up from 0 K to the target temperature. In order to avoid wild uctuations and hot solvent cold solute etc., the Langevin temperature equilibration scheme was used for the equilibrium of the system temperature (Sindhikara et al. 2009). Due to the inaccuracy of the pressure calculation at rst few picoseconds (ps) under low temperatures, one nanosecond (ns) of simulations at the NVT ensemble with restraints were performed initially for all conditions. All bonding involved hydrogen atoms were constrained by the SHAKE algorithm simultaneously (Miyamoto and Kollman 1992). Then following simulations were performed under NPT ensemble with the pre-determined temperature. The Berendsen barostat was used to maintain the pressure around 1 atm and the relaxation time was set as 2 ps (Berendsen et al. 1984). Same to the equilibrium steps, SHAKE algorithm was applied to constrain binding involved hydrogen atoms and the Langevin dynamics was used for the temperature controlling (Miyamoto and Kollman 1992; Sindhikara et al. 2009). The time step was set as 2 femtoseconds (fs) and the snapshots were saved every 10 ps at desired temperature. 100 ns MD simulations were performed at 298 K, 350 K, 400 K and 450 K. In order to ensure the equilibrium was achieved with su cient time period, 200 ns simulations were run for 500 K and 550K.
Additional 10 sets of 5 ns simulations with di erent starting con gurations based on 100 ns simulation trajectories were done under each temperature to enhance the validity of our simulations. The random number generators were used in each MD simulation. From these MD simulation sets, errors can be estimated for di erent parameters (density, unit cell, hydrogen bonding and torsion angle etc.) calculated in simulations and the average values were used for further analysis.

Results and Discussions . Equilibration of system
The root means square deviation (RMSD) of the atom C1, C2, C3, C4, C5 and O5 and total energy were used to evident the equilibration of system at di erent temperatures. A stable RMSD and total energy value indicates the equilibration was achieved in the system. Di erent time lengths were required to achieve stable values under di erent temperatures (Fig. 2).
The curve of simulation at 298 K in Fig. 2(a) demonstrates that the mean of RMSD is around 0.49 Å and the total energy is around − kcal/mol ( Fig. 2(b)). The system reached equilibrium immediately, indicating that the backbone structure changes negligibly to initial structure at the room temperature. The mean The mean of RMSD is nearly 1.46 Å at temperature 500 K and the total energy is nearly 50 kcal/mol. This indicates that the system achieves equilibrium after 50 ns (Fig.2, a). According to the nonlinear increasing of RMSD from 450 K to 500 K, an alteration of cellulose Iβ molecular structure occurred under higher temperature. At 550 K, the system reached equilibrium after 75 ns with a mean of RMSD about 3.72 Å and the total energy about 20000 kcal/mol. This further increased RMSD indicates that the molecular structure of cellulose Iβ experienced signi cant changes under this temperature with a slow pace in achieving equilibrium. In general, the cellulose Iβ system reached a temporary equilibrium in rst 20 ns regardless the simulation temperature, which is consistent with the previous study (Zhang et al. 2011). However, the structure of system changes signi cantly within rst 75 ns at high temperature.

. Structural properties
The density, unit cell parameters and torsion angles were used to represent the structural proprieties of cellulose Iβ at di erent temperatures in the simulation. The MD simulation trajectories were utilized to investigate in uence of temperature on the structure and correlated properties from 298 K to 550 K.
The crystal density is a macroscopic parameter of cellulose in Iβ crystalline form and is a critical indicator of its structural characters. The density is convertible so that the constant pressure is retained while performing MD simulation in NPT ensemble. The average density evolution with temperature was shown in Fig. 3. At the room temperature (298K), the average equilibrated density was found about 1587 kg/m with a standard deviation of 2.92 kg/m . This value is less than the experimental density of 1636 kg/m estimated from the crystal cell parameters (Nishiyama et al. 2002). The di erence between the density obtained from simulation and the experimental value is about 3.0 %. In addition to evaluate the accuracy of our simulation by comparing with experimental data at 298 K, we monitored the density change of cellulose Iβ when raising of the temperature. The general trend shows that the density of cellulose molecules decreased when increasing the temperature. The density decreases about 3.4% as the temperature changes from 298 K to 350 K. This is greater than that when the temperature raised from 350K to 450K, which is about 1.3%. However, the density decreases nonlinearly by about 3.2% as the temperature increases to 500 K from 450 K and a sudden drop occurs by about 22.5% after heating the system to 550K. The sudden drop of density of cellulose was reported at 453K-523K in experimental studies for various samples (Swatloski et al. 2002;Zou et al. 2002) and a range of 450-500K was reported in MD simulation of cellulose Iβ (Wyman 2007). According to that reported in previous studies by both experiments and molecular dynamics simulations, a temperature where a sudden drop of cellulose density happened is considered as the transition temperature (Swatloski et al. 2002;Wyman 2007;Zou et al. 2002). At this temperature, the structural properties of cellulose have signi cant changes. Our simulation results (Fig. 3) shown a similar trend when comparing with these literatures, indicating that the structure of cellulose is becoming slightly looser with the temperature increasing and discussion ensues.
Besides the density estimation above, the unit cell parameters of each glucose units were calculated in this work under various temperatures. Unit cell parameters were measured by locally regarding 60 di erent unit cells for di erent temperature. Corresponding distances and angles between O4 atoms were traced during the last 4 ns temperature (Bergenstrahle et al. 2007;Zhang et al. 2011). The unit cell parameters obtained from simulation are listed in Table 1  The structure alternation within temperature range is unveiled by the variation of the average values for all unit cell parameters at di erent temperatures as shown in Fig. 4. Signi cant increases of unit cell a, b and β were observed when temperature raised from 298 K to 500 K. Therefore, structure at 500 K can be deemed as a representative of cellulose Iβ structure at a high temperature. In addition, change of the unit cell angles (α, β, and γ) at high temperature are smaller than 0.7% when compared with that at 298 K.
The unit cell parameter a is 8.1% longer than that at the low temperature, which is comparable with the 5.3% detected by experiment (Wada 2002) and 7.5% by another MD simulation (Bergenstrahle et al. 2007). The parameter b increases merely 1.8% from 298K to 500 K in our simulation, but it was reported with an about 6% increase from 298K in previous work (   the b axis in unit cell demonstrated di erence in values, indeed anisotropic expansion thermal expansion was con rmed in our simulation results (Wada 2002). A sudden change in unit cell parameter is speculated on a transition in structure along the temperature change. Therefore, with the observation of the increased axis lengths a, b and angle β, it can be concluded that the transition temperature of cellulose Iβ lies in the temperature around 450 K at which the signi cant structural alternation occurred (Fig. 4). An intermediate structure evidenced by the unit cell parameters was observed in between low and high temperature. The equilibrated structure at temperature 400 K was considered as the intermediate structure in our study and the structural parameters are listed in Table 1. The trend of unit cell parameters changes from 298K to 550K is plotted in Fig. 4. At 400 K, the parameter a was increase about 4% while both parameters b and γ were increased about 0.9%, compared to that at 298 K. When system was heated to 450 K, only parameters a and b demonstrated a slight growth compared to that at 400 K. Therefore, the changes of unit cell parameters from 298 K to 400 K are less than that from 400 K to 500 K. In addition, a rearrangement of structure was observed from temperature 298 K to 400 K as shown in Fig. 5a and Fig. 5b. The minimal of unit cell parameter α was observed at 350 K (Fig. 4), suggesting that a temperature around 350 K can be considered as another transition temperature in addition to the 400 K. Because the changes of the surface chains at the cellulose/water interface are caused by the water environment and are able to capture the ber twist. This phenomena was considered as surface e ects (Bergenstrahle et al. 2007;Yui et al. 2006). Therefore, the structure transition at 350K is caused by the water environment and surface e ects.
The graphical structures of the crystal structure equilibrated at temperature 298 K, 400 K and 500 K are shown in Fig. 5 as the representatives of structures at low (Fig. 5a), intermediate (Fig. 5b) and high temperature ( Fig. 5c and Fig. 5d). It appears that the intermediate structure at 400 K tends to twist clockwise. The number of twisting chains is mounting progressively with the temperature increasing from 350 K to 450 K. This phenomena was also observed in experiments (HANLEY et al. 1997) and the previous MD simulations (Paavilainen et al. 2011). The twisting in structure results in an increment of the unit cell angle β after 450 K as shown in Fig. 4. At 500 K, two cellulose chains in adjacent sheets were further away from each other compared to that at lower temperature, suggesting that chains were separated from each other along the c axis (Fig. 5d). This separation contributes to about 1.1% increasing of unit cell angle β compared to that at 400 K ( Table 1). The change in a parameter is consist of two factors: the expanded distance between di erent sheets and the slippage along the molecule chains. If the distance was not enlarged between di erent sheets, a slight change in unit cell angle β should cause only 0.04% increase of unit cell parameter a mathematically with the sine function. However, about 4% increase of parameter a was found in our study. Therefore, the increase of a parameter found in our study is mainly derived from the increase of expanded distance between di erent sheets.
We also analyzed the conformation of molecule chains in order to characterize the thermal response of the cellulose Iβ. The population of gt, tg and gg conformations at di erent temperature can be estimated by the distribution of the torsion angle ω as shown in Fig. 6. Without translational symmetry along the chain direction in the initial cellulose Iβ structure, the end residues of every chain were unable to be used in estimating the distribution of hydroxymethyl groups. At 298 K, most of hydroxymethyls are in tg conformation, and the population of gg and gt conformations are minority. This is in agreement with the X-ray di raction data (Nishiyama et al. 2002).
A drastic change emerged at 350 K where about half of the hydroxymethyl groups convert from tg to the gt or gg conformation. In addition, the population of gg conformation surpass that of gt conformation. The gt and gg conformations represent rotated and non-rotated backbone chains respectively, according to the previous work (Zhang et al. 2011). Therefore, the observed twisting of cellulose chains can be attributed to bond rotations in order to adopt gt conformation. The gt conformation became more populated when the temperature ascending from 400 K to 550 K, while the population of gg conformation slightly decline. This result agrees with the work of Bergenstråhle (Bergenstrahle et al. 2007) that the population of gt conformation exceeds that of gg conformation at 550 K, whereas original tg conformation become negligible. It was reported by former works that bond rotations often occur near the transition temperature (Bergenstrahle et al. 2007;Zhang et al. 2011). Although larger quantities of twisted chain along the c axis were observed in our study compared to others results (Bergenstrahle et al. 2007;Zhang et al. 2011); however, this character agrees with experimental data that the C6 atom is more impressionable at high temperature compared to that at room temperature (Watanabe et al. 2006b).

. Hydrogen Bonds
The change in the orientation of the hydroxymethyl groups and glycosidic bond may lead to a di erent hydrogen bond network between cellulose molecule chains. Fig. 7 demonstrates the characteristics of hydrogen bonds in interior area that still maintain crystalline structure. The hydrogen bond cuto was set as a bond length of 3.0 Å and a bond angle of 150 • . According to the studies of Nishiyama and others (Nishiyama et al. 2008;Nishiyama et al. 2002;Zhang et al. 2011), the hydrogen bonds in cellulose crystal structure can be classi ed into three types: intrachain, interchain and intersheet hydrogen bonds. The hydrogen bonds O3-HO3. . . O5 and O2-HO2. . . O6 are the major types of intrachain and the hydrogen bond O6-HO6. . . O3 is the major type of interchain in the native cellulose Iβ (Fig. 7a). The intrachain hydrogen bond O3-HO3. . . O5 is constructed by hydroxyl group 3OH at the C3 position with the O5 atom in the adjacent ring. The intrachain hydrogen bond O2-HO2. . . O6 is formed by group 2OH at the C2 with the O6 atom in the adjacent ring, while the hydroxyl group 6OH at the C6 position and the O3 atom in the neighboring chain induce interchain hydrogen bond formation O6-HO6. . . O3 .
Above results indicate that the hydrogen bonds are critical to stability of intrachain and interchain but not to intersheet.
We also analysis the alternation of hydrogen bond network pattern under various temperatures in this study. In general, the populations of the intrachain hydrogen bonds O3-HO3. . . O5, O2-HO2. . . O6 and interchain hydrogen bond O6-HO6. . . O3' decrease with the temperature rising (Fig. 7b). Combing with the structural alternation of cellulose above, it is obvious that hydrogen bond network contributed to the chain rearrangement under di erent temperatures. When the system was heated from 298 K to 350 K, the population of the hydrogen bond O3-HO3. . . O5 decreased by 28% and the hydrogen bond O2-H2. . . O6 decreased by 36%. Although about half of the hydroxymethyl groups convert from tg to the gt or gg conformations at 350 K, most of them occurred at the cellulose/water interface. In the addition, the hydrogen bond was estimated by crystal internal 16 di erent unit cells in the system. This leads to about 60% of the intramolecular hydrogen bond of C2OH-O6 still remained. These results are partially in agreement with the experimental data provided by studies using IR spectroscopy, in which the alternation of intrachain hydrogen bond O3-HO3. . . O5 and O2-HO2. . . O6 was observed within the temperature range of 303-473 K (Watanabe et al. 2006a;Watanabe et al. 2006b). This is also comparable with previous MD simulations, in which a 40% decreasing of the interchain hydrogen bond O6HO6. . . O3 was observed when the temperature was raised from 298 K to 350 K (Zhang et al. 2011). In addition, a new type of interchain hydrogen bond O2-HO2. . . O6 was formed with a population of 72% at 350 K. This can be attributed to the conformational changes of the hydroxymethyl groups in interior area. The change of populations of intrachain and interchain hydrogen bonds indicates that hydrogen bonding pattern at 298 K translated to another one when the system was heated to 350 K.
On the other hand, at 500 K, the population of the hydrogen bond O3-HO3. . . O5 decreased to about 36% and the hydrogen bond O2-HO2. . . O6 has only about 3% left. The interchain hydrogen bond O6HO6. . . O3 and O2-HO2. . . O6 decreased to 13% and 61% respectively, which are in line with the previous MD simulations (Bergenstrahle et al. 2007;Zhang et al. 2011). The population of O2-HO2. . . O6 hydrogen bond is 13% less at 500 K compare that at 400 K. These results imply the changing of structural pattern of cellulose crystals. Therefore, it can be concluded that the structure of cellulose Iβ at high temperature introduced a less ordered nature.
Minor di erence in hydrogen bonds network in response to the temperature change is demonstrated as bond distance distribution of the hydrogen bonds as shown in Fig. 8. The most probable hydrogen acceptor distance of O3-HO3. . . O5 is 2.76 Å at 298 K, which is slightly shorter than that under other temperatures. The population of O3-HO3. . . O5 hydrogen bond is 16.1% at 298 K (Fig. 8a), which is slightly larger than that under other temperatures. These phenomenon of O3-HO3. . . O5 hydrogen bond are in agreement with ndings of previous MD simulations ( Bergenstrahle et al. 2007). The ranges of distance from 2.5 Å to 3.1 Å (Fig. 8b) showing that the intrachain hydrogen bond O2-HO2. . . O6 at lower temperature range (298 K and 350 K) are stronger than that under higher temperature (400 K to 550 K). Likewise, the intrachain hydrogen bonds O3-HO3. . . O5 (Fig. 8a) and O2-HO2. . . O6 (Fig. 8b) are stronger than the interchain hydrogen bonds O6HO6. . . O3 (Fig. 8c) and O2-HO2. . . O6 (Fig. 8d) at lower temperature range. These results are consistent with the previous works (Bergenstrahle et al. 2007;Zhang et al. 2011). The most probable hydrogen acceptor distance of the newly formed O2-HO2. . . O6 (Fig. 8d) hydrogen bond is about 2.76 Å at 400 K, shorter than that (2.84 Å) at 500 K and the corresponding population of O2-HO2. . . O6 hydrogen bond decreases from 400 K to 500 K. These can be explained by the fact that a few hydroxymethyl groups transform to the gt conformation from gg conformation, contributing to the fact of weakening the hydrogen bond strength.

Conclusions
In a summary, thermal response of cellulose Iβ was reported by molecular dynamics simulations with the GLYCAM06 force eld in the present paper. All simulations reached equilibrium with the RMSD and total energy used as the indicator of equilibrated system. With less than 3.9% deviations of density and unit cell parameters between the data measured in previous experiments and that in our simulation performed at 298 K, the computational methods and selection of force eld in our study was proved with validity for further analysis of the structure. In addition, the hydrogen bonding pattern calculated at temperature 298 K was in su cient accordance with that measured experimentally, herein as another support of this computational simulation.
A determined structure transition temperature (about 450 K) was consistent with the experimental data, meanwhile, another transition temperature around 350 K was detected in our simulation. Three structures at low, intermediate and high temperature was identi ed and demonstrated individually in this work. In the intermediate temperature structure, twisting emerged in a few chains along the chain direction. We also observed the sudden changing of density and unit cell parameters a, b, α and γ at this intermediate temperature.
In addition, the conformational changes of hydroxymethyl groups from tg to gt and gg were observed when raising the system temperature. Note that the gg conformation was dominant into conformational alternation. This also contributes to the transitions of the hydrogen bonding pattern under various temperature. When the temperature was further increased to 400 K, majority of the hydrogen bonds O2-HO2. . . O6 and O6-HO6. . . O3 were disrupted and substituted by interchain hydrogen bonds O2-HO2. . . O6 . The structural responses of cellulose Iβ at higher temperature were illustrated as consequences of chain sliding along c axis, bond rotations, drastic changing of density and unit cell parameters a, b and β, then the population of hydroxymethyl groups in gt and gg exchanged, therefore a conversion of hydrogen bonding pattern.
The structures under low and high temperature in our work were in reasonable agreement with the results obtained by the previous MD simulations and experiments. The discovery and illustration of three structures of cellulose Iβ that existed simultaneously under di erent temperature by molecular dynamic simulation will contribute to understanding the details of structural properties at molecular level. The simulated trajectories demonstrate details of thermal response of cellulose Iβ by providing the structures from atoms to hydrogen bond patterns, explaining the physical phenomenon that observed in experimental studies at a molecular level.