NC State
Kong, Y., Fu, S., Yang, X., Leu, S-Y., and Hu, C. (2023). “Cellulose Iβ behaviors in non-solvent liquid media: Molecular dynamic simulations,” BioResources 18(4), 8223-8248.


The structural changes of cellulose in non-solvent liquid media can provide insights into the high-value utilization of cellulose. This study includes molecular dynamics simulations of 36-chain cellulose Iβ microfibril model (Iβ-MF) behavior in 16 non-solvent liquids with different polarities at room temperature using two carbohydrate force fields (CHARMM36, GLYCAM06). Iβ-MF in CHARMM36 retains more than 70% of the tg conformation in 16 liquids, and the retention of the tg conformation increased with decreasing liquid polarity. Liquid polarity can affect the hydroxymethyl conformation of cellulose, which is only an appearance, and the real driving force behind is the electrostatic interaction between liquid molecules and cellulose. Furthermore, changing the 1,4 electrostatic scaling factor of GLYCAM06 can effectively affect the structural convergence of Iβ-MF. The Iβ-MF forms an alternating layer structure in the gg/gt conformation in a medium to high polarity non-solvent liquid, while the model undergoes untwisting. Model untwisting is inextricably linked to the degree of alternate layer structure formation. This paper provides a theoretical basis for the molecular study of nanocellulose structures from an energy-structure-property perspective.

Download PDF

Full Article

Cellulose Iβ Behaviors in Non-solvent Liquid Media: Molecular Dynamic Simulations

Yi Kong,a Shiyu Fu,a,b,* Xuedi Yang,b Shao-Yuan Leu,c and Chuanshuang Hu d

The structural changes of cellulose in non-solvent liquid media can provide insights into the high-value utilization of cellulose. This study includes molecular dynamics simulations of 36-chain cellulose Iβ microfibril model (Iβ-MF) behavior in 16 non-solvent liquids with different polarities at room temperature using two carbohydrate force fields (CHARMM36, GLYCAM06). Iβ-MF in CHARMM36 retains more than 70% of the tg conformation in 16 liquids, and the retention of the tg conformation increased with decreasing liquid polarity. Liquid polarity can affect the hydroxymethyl conformation of cellulose, which is only an appearance, and the real driving force behind is the electrostatic interaction between liquid molecules and cellulose. Furthermore, changing the 1,4 electrostatic scaling factor of GLYCAM06 can effectively affect the structural convergence of Iβ-MF. The Iβ-MF forms an alternating layer structure in the gg/gt conformation in a medium to high polarity non-solvent liquid, while the model undergoes untwisting. Model untwisting is inextricably linked to the degree of alternate layer structure formation. This paper provides a theoretical basis for the molecular study of nanocellulose structures from an energy-structure-property perspective.

DOI: 10.15376/biores.18.4.8223-8248

Keywords: Cellulose; Liquid polarity; Molecular dynamics; Structural transformation

Contact information: a: State Key Laboratory of Pulp and Paper Engineering, South China University of Technology, Guangzhou 510640, PR China; b: South China University of Technology-Zhuhai Institute of Modern Industrial Innovation, Zhuhai 519175, PR China; c: Department of Civil and Environmental Engineering, The Hong Kong Polytechnic University, Hong Kong SAR 999077, PR China; d: College of Materials and Energy, South China Agricultural University, Guangzhou 510640, PR China;

* Corresponding author:


Cellulose, the main component of plant cell walls, is nature’s most widely distributed and abundant homoglycan (Matthews et al. 2006; Zhao et al. 2013; Zhang et al. 2019; Zhou et al. 2021). Cellulose consists of repeating β-D-glucopyranose units linked by β (1-4) glycosidic bonds. These pyranose rings are in a chair-like conformation with the hydroxyl group in the equatorial position (Kamide 2005). In nature, wood cellulose chains’ polymerization (DP) is about 10,000 pyranose units, and cotton cellulose is about 15,000. Natural cellulose is defined as cellulose I, which is known from 13C CP/MAs NMR spectra to exist in two different forms called cellulose Iα and Iβ (Horii et al. 1987). The main difference between cellulose Iα and Iβ is the stacked arrangement of the hydrogen-bonded layers. The layers in Iα (P1 symmetry) are permanently displaced +c/4 along the c-axis, while the layers in cellulose Iβ (P21 symmetry) are alternately displaced at +c/4 and -c/4 (center and origin chains) (Gardner and Blackwell 1974). Cellulose I from naturally occurring lower plants (e.g., algae and bacteria) is rich in cellulose Iα. In contrast, cellulose Iβ, mainly from higher plants (e.g., cotton and wood), is the most studied and widely used type of cellulose (Imai and Sugiyama 1998; Nishiyama et al. 2008).

Hydrogen bonding is an essential component of the cellulose crystal structure. For cellulose I, hydrogen bonding within the cellulose chains allows a linear arrangement of cellulose, and two adjacent cellulose chains are bound together mainly by hydrogen bonding, forming cellulose sheets. Stacked sheets, conversely, are thought to have no hydrogen bonding interactions between them but rather form cellulose microfibril crystals through van der Waals interactions (Zhou et al. 2021). Notably, the orientation of the C6 hydroxymethyl group will highly affect the hydrogen bonding and chain conformation. This pyranose ring substituent (C6 hydroxymethyl) has three possible minimum energy orientations: trans–gauche (tg), gauche–trans (gt), and gauche–gauche (gg) (Shefter and Trueblood 1965). The cellulose models for each of the three orientations (tg, gt, and gg) were compared with the X-ray data to obtain the best fit, described as the model with the lowest reliability factor (R). For cellulose I, the tg direction gave the best fit with R values of 0.242, 0.292, and 0.349 for the tg, gt, and gg models, respectively (Gardner and Blackwell 1974). The majority view in the literature is that cellulose I has a tg conformation throughout the chain (Sarko et al. 1976; Nishiyama et al. 2008).

The interaction between cellulose microfibrils and solvents or the solubilization process of cellulose microfibrils is a research hotspot (Zhang et al. 2019; Zhou et al. 2021). Still, studies have yet to investigate cellulose’s behavior and conformational changes in non-solvent liquids. It is crucial to observe the effect of polarity on the internal and external structure of crystals in polar or nonpolar liquid media and to explore the mechanism behind it. Related works have provided a substantial basis for this study, including cellulose’s room- and high-temperature behavior in a single solvent (Matthews et al. 2006, 2011, 2012; Zhang et al. 2011) and cellulose phase transitions (Bellesia et al. 2011). These works have demonstrated a helpful tool – computer simulation. Molecular dynamics (MD) computational methods have gained much attention in cellulose structure studies because of their ability to obtain valuable information regarding cellulose structure and properties to have sufficient sampling to explore cellulose conformation. MD simulations were used to model the wetting of cellulose surfaces by water. The theoretical wetting limit of cellulose was solved by considering the surfaces (110) and (100) of the Iβ variant, respectively (Mazeau and Rivet 2008). Gross and Chu (2010) used MD simulations to observe cellulose microfibrils in water and concluded that inter-sheet interactions of cellulose microfibrils were the most robust component. The following year, the team used the same method to study the state of cellulose in two liquids (water and BmimCl), showing that the insolubility of cellulose in water results mainly from a reduction in liquid entropy (Gross et al. 2011). In addition, a considerable number of MD simulations focused on the distortion phenomenon of cellulose microfibrils. Matthews et al. (2006) first reported the tendency of right-handed twisting of cellulose Iβ microfibrils under Charmm force field in a short time (< 200 ps). In the same year, Yui et al. (2006) reported similar behavior of different microfibril models under Glycam force fields. Hadden et al. (2013) investigated the driving forces behind the twisting behavior, including the role of the liquid, the effect of non-bonding force field parameters, and the use of explicitly modeled oxygen lone pairs in the solute and liquid. The results showed that the twisting of microfibrils is influenced by van der Waals interactions and is counteracted by intra-chain hydrogen bonding at the microfibril surface and liquid effects.

MD calculations strongly depend on molecular force fields. The appropriate force fields can match experimental phenomena well or even be used to accurately calculate the system’s internal state to predict the experimental phenomena. MD simulations of cellulose are usually based on the molecular force fields of carbohydrates, and the most commonly used molecular force fields for cellulose are CHARMM36 and GLYCAM06 force fields (C36 and G06) (MacKerell et al. 1998; Kirschner et al. 2008; Huang et al. 2017).

In this work, two common carbohydrate force fields (C36 and G06) were used to study the behavior of cellulose Iβ microfiber in a large number of polar or nonpolar liquids, focusing on the crystal parameters, conformational changes, and fiber twisting of cellulose Iβ microfiber under the influence of liquid polarity. The effect of liquid polarity on model conformation was analyzed from an energy perspective. The reason why these liquid media cannot dissolve cellulose may be that the liquid media cannot penetrate into the crystal, cannot change the internal hydroxymethyl structure, and cannot destroy the hydrogen bond network and van der Waals forces in the crystal structure. It was noteworthy that a short time might exist to not easily discern the simulation results in liquids of similar polarity. Therefore, G06 with an electrostatic scale factor of 1,4 was chosen for the simulation (G06 is fully compatible with the Amber14 force field with a default factor of 5/6). The purpose was to “speed up” the simulation process and to observe more simulation results. The “C36 model” and “G06 model” in this paper refer to the cellulose Iβ microfiber model (named Iβ-MF) in the C36 and G06 force fields, respectively.


Cellulose Iβ Microfiber Modeling

A crystal structure file of cellulose Iβ was constructed through the Cellulose-Builder website (Gomes and Skaf 2012), and the model contained 36 glucan chains with 15 cellobiose per chain (Fig. 1). The structure was a monoclinic P21 space group with lattice parameters of a=7.784 Å, b=8.201 Å, c=10.380 Å, and y=96.5°. Current studies generally agreed that the cross-sectional shape of cellulose microfibrils was hexagonal. The number of cellulose microfibrils containing cellulose chains was controversial, and both 18-24 and 36 chains models had been reported (Cosgrove 2014; Ding et al. 2014). The 36-chain Iβ-MF was chosen because a more significant number of chains would allow more structural information and data to be obtained for better analysis of the behavior of Iβ-MF in liquids.

Liquid Polarity Calculation

The inability to disrupt the crystal structure was a prerequisite for studying the behavior of Iβ-MF in liquid media, so none of the selected liquids could dissolve cellulose. The preferred liquids were water, N,N-dimethylformamide (DMF), formic acid, acetic acid, propionic acid, methanol, ethanol, propanol, pyridine, carbon tetrachloride (CCl4), trichloromethane (HCCl3), dichloromethane (H2CCl2), benzene, carbon disulfide, n-hexane, and cyclohexane. The molecular polarity index (MPI) was an important parameter to express the polarity of a molecule and can measure the overall polarity of the molecule. MPI was defined as shown in Eq. 1 (Liu et al. 2021),


where V is the molecular electrostatic potential, A is the molecular surface area, and the function is integrated over the molecular surface S. The system’s charge distribution determines the polarity’s magnitude. Uneven charge distribution led to differences in the molecular surface’s electrostatic potential, causing a polarity change. The more inhomogeneous the charge distribution, the greater the polarity and MPI.

The computational software optimized all molecular structures at the B3LYP-D3(BJ)/def2-SVP level. The single point energy was calculated at B3LYP-D3(BJ)/def2-TZVP, followed by the quantitative analysis of the molecular surface module of Multiwfn (Lu and Chen 2012) to obtain the MPI.

Fig. 1. Iβ-MF and color scheme of hydroxymethyl group conformation. a) Iβ-MF cross-section (ab surface) and chain numbering scheme. The layers were numbered from top to bottom, left to right, e.g., 11, 12, and 13 for the first layer, 21, 22, 23, and 24 for the second layer, and so on. The blue line was a dividing line between the inner and outer cellulose chains. b) The ac surface of cellulose Iβ microfibrils, each chain contained 15 cellobiose. c) The three monomers in the single chain were shown from two orientations, TG in yellow, GG in purple, and GT in green.

Details of MD

The restrained electrostatic potential (RESP) charge (Bayly et al. 1993) was first calculated for all molecular monomers. Their RESP charges for liquid molecules could be calculated from the file in the previous step. For cellulose chains, the excessive number of atoms caused difficulties in structural optimization. To obtain the cellulose chain’s RESP charge, an oligomer was constructed containing five glucose units. The structure was geometrically optimized at the B3LYP-D3(BJ)/def2-SVP level. The single-point task was computed at the B3LYP-D3(BJ)/def2-TZVP level, and the resulting wave function file was imported into the Multiwfn program to compute the RESP charge. The RESP charge was calculated three times to constrain the total charge of the head, central, and tail units to 0, respectively. Finally, the RESP charge of the oligomer was extended to the cellulose chain.

MD simulations were performed using GROMACS 2020.6 (Kutzner et al. 2019) software and two force fields of carbohydrates (G06 and C36). Iβ-MF was placed in the middle of a 60 × 80 × 200 box and wrapped with approximately 16 molecules in each direction. When describing Iβ-MF using G06, all liquid molecules used the Gaff force field to construct molecular topology files. When Iβ-MF were defined using C36, all liquid molecules were required to generate the molecular topology files using the Glycan Reader and Modeler in CHARMM-GUI (Lee et al. 2016). Periodic boundary conditions (PBC) were applied to the three directions of the box to ensure that the number of particles inside the box was constant. The Langevin thermostat and the Nose-Hoover Langevin pressure regulator stabilized the temperature and pressure at 298.15 K and 1 atm, respectively. The TIP3P water model was used with a 2-fs time step for the dynamics. The Particle-Mesh Ewald (PME) method was used to handle the long-range electrostatic forces, and a non-bonding cutoff distance of 10 Å was applied. In the C36 and G06 simulations, the steepest descent scheme was used to minimize the system energy, requiring a maximum force of less than 100 kJ mol-1 nm-1, followed by gradual heating to the prescribed temperature. After reaching the desired temperature, the energy-minimized system did a restricted MD on the cellulose microfibers by the NVT ensemble to ensure that the cellulose microfibers would not move before liquid relaxation. Finally, a long-time MD simulation with a duration of 50 ns was performed. The trajectory data of the system were analysed using the Visual Molecular Dynamics (VMD) package (Humphrey et al. 1996).

Interaction Energies

The interaction energy is expected to be a powerful tool for studying weak interactions between molecules, as shown in Eq. 2,


where Eelec was the electrostatic interaction energy, acting as an attraction (negative) or repulsion (positive); Evdw was the dispersive interaction energy, corresponding to the long-range Coulomb correlation of electrons, acting as an attractive force. In medium-strength hydrogen bonds and dihydrogen bonds, electrostatic interactions were dominant, complemented by dispersive interactions. Both van der Waals (vdW) interactions and π-π stacking interactions are dispersive.


Liquid Polarity Data

The polarity data for all liquids are shown in Table 1, sorted by MPI from highest to lowest. Water, formic acid, and DMF molecules (MPI>15 kcal mol-1) are strongly polar, acetic acid, methanol, propionic acid, ethanol, pyridine, propanol, and H2CCl2 molecules (8<MPI<15 kcal mol-1) are moderately polar, benzene and HCCl3 molecules (4<MPI<8 kcal mol-1) are weakly polar, and CCl4, CS2, n-hexane, and cyclohexane molecules (MPI<4 kcal mol-1) could be considered as non-polar. The pApolar content, Vmin, and Vmax could be used as auxiliary indicators for the polarity of liquid molecules. Liquid molecules with a high percentage of pApolar or large absolute values of Vmin and Vmax could be considered a prerequisite for high polarity.

Table 1. Polarity Data for 16 Liquid Molecules

a) The pApolar was the molecular polar surface area percentage over the total surface area.

b) The minimum and maximum molecular surface electrostatic potentials were defined as Vmin and Vmax.

Cellulose Iβ Microfiber Crystal Structure

The RMSD of Iβ-MF (Fig. S1; see Appendix) showed that all systems reached equilibrium within 50 ns. In C36, Iβ-MF showed the most minor structural change in the formic acid liquid, considered an anomaly and discussed in detail in the subsequent subsection.

Table 2. Crystal Parameters of Crystalline Cellulose Iβ

a) Experimental data at room temperature (Nishiyama et al. 2003)

b) Lattice parameter values from the simulation were averaged over the final 10 ns. Standard deviations were shown in parentheses.

c) Lattice parameter values of the center chains

d) Lattice parameter values of the origin chains

Table 2 shows the simulated lattice parameters of Iβ-MF and the lattice parameters of Iβ determined by synchrotron radiation and neutron diffraction experiments. The lattice parameters of the Iβ-MF were generally similar in the G06 and C36. Combined with Table 2 and Fig. S1, the crystal structure of Iβ-MF did not change substantially in the liquid and could maintain the hexagonal cross-section well. Detailed data on the lattice parameters of the Iβ-MF for all liquid systems are listed in Table S1. For the G06 model, the c-lattice parameter (cellobiose length) was 0.31 Å larger than the experimental value, which could be a force field defect. Table S1 exhibited the pattern of a-lattice parameters with liquid polarity, and the Iβ-MF interlayer distance was positively correlated with polarity. There was a tendency for the Iβ-MF to swell in liquids higher than benzene polarity, while the opposite was true in liquids lower than benzene polarity. In addition, the perturbation of the α- and β-angles of the center chain by the liquid was higher than that of the origin chain. For the C36 model, the standard deviation indicated that the lattice parameters of the model were largely unaffected by the liquid.

The Behavior of Cellulose Microfibers in Liquids

Herein, 36-chain Iβ-MF (DP=30) behavior was reported in different liquids in C36 and G06. Figure S2 organizes snapshots of the cellulose microfibrils for each of the 100 layers at 50 ns for two force fields. Figure S3 depicts the variation of the number of gg, gt, and tg conformations over time for the two force fields, and the conformational data (number of gg, gt, and tg conformations) for the last moments are collated in Fig. 2.

Fig. 2. Average hydroxymethyl conformation of a) CHARMM36 and b) GLYCAM06 simulations at 50 ns. Chart c was an additional illustration of the data with large deviations in the three replicate simulations in Chart b, correlated with the twist data in Fig. 7 (explained in Fig. 7). The data sets H2CCl2_3, and CS2_1 in Chart c were not counted in chart b because of the large difference between them and the remaining two data sets in the same group.

The initial conformation of hydroxymethyl groups in all systems was tg. Several relevant observations could be derived from the simulation results for each force field. C36 simulations showed that at most 30% of the tg conformations were transformed in 50 ns. More tg conformations were retained as the liquid polarity decreased, with a negative correlation between liquid polarity and the number of tg conformations (except for the formic acid liquid system). In liquids other than formic acid, the difference in the number of gg conformations in each system was slight, ranging from 5.6% to 9.4%. The liquid media can only change the hydroxymethyl conformation on the surface and fails to change the internal conformation. It is therefore believed that these liquid media cannot dissolve cellulose because they cannot break the hydrogen bonding network and van der Waals forces of cellulose. The formic acid system is abnormal in each of the three replicate calculations. The simulation results show that the number of tg conformations in the formic acid system fluctuated very little, and there was no conversion of tg conformation to other conformations. The version number of the computer operating system and simulation software was changed to overcome this anomaly, and five calculations were carried out using the simulation parameters; the results are shown in Table S2. Most subsequent results were consistent with the earlier results, with only one set of exceptions, with 157 tg conformations transformed. According to the conclusion above, formic acid and DMF were both highly polar molecules, and the simulation results for both liquids should be similar, but they were very different. C36 might not be well described for both formic acid and cellulose microfibrils.

The anomaly of the formic acid system was not reflected in the G06 simulations. The transformation of cellulose microfibrils is evident in G06, which is in line with an expectation that the narrowing of 1,4 electrostatic scaling factors will necessarily affect the structure of cellulose substantially. The stock of tg conformation in the low-polarity liquid was less than 40%. It could be assumed that the tg conformation of cellulose microfibrils in the low polarity liquid could be transformed entirely given enough time.

The hydroxymethyl group exposed to the liquid was susceptible to conformational transitions, including odd glucose units at the upper edge, even those at the lower edge, and those at the non-reducing ends, which could be observed in the simulation results for C36. In contrast, interlayer conformational differentiation occurred in the G06 model (alternating layer pattern), with the center layers (even layers) mostly in the gg conformation and the origin layers (odd layers) mostly in the gt conformation, a phenomenon that was particularly evident in liquid systems with high polarity. In the center layer, the hydroxymethyl group exhibiting the gg conformation was essentially perpendicular to the mean plane of the pyranose ring and therefore pointed upward and downward to the original chain in the upper and lower layers.

In this alternating layer pattern, the spatial conflict between these center chain hydroxymethyl groups and the upper and lower layers (origin layers) forced the cellulose chains to tilt significantly for the plane of their layer. The center chain (gg conformation) was tilted counterclockwise, and the origin chain (gt conformation) was tilted clockwise when viewed from the non-reducing end, as shown in Fig. 3. The tilting of the chains led to changes in the hydrogen bonding network, and Figs. 4a, b, and c showed the hydrogen bonding network before and after the simulation, respectively. The extra-ring groups could form effective O6-O2 hydrogen bonds between chains in the same layer (Figs. 4b and c) or O6centerH6center-O2origin hydrogen bonds between adjacent layers. This simulated structure, in terms of chain tilt and hydroxymethyl conformation, which was very similar to the high-temperature phase structure of I-HT (Matthews et al. 2012).

Fig. 3. Schematic diagram of the alternating layer pattern of cellulose microfibers

Fig. 4. Schematic diagram of the hydrogen bonding network before and after the simulation. a) interlayer hydrogen bonding network between the center and origin layers, mainly O6centerH6center-O2origin hydrogen bonds; b) interchain hydrogen bonding network diagram inside the origin layer with the center layer in the background; c) interchain hydrogen bonding network diagram inside the center layer with the origin layer in the background

Interactions Between Iβ-MF and Liquids

The interaction between Iβ-MF and liquids was described by the interaction energy (Lu and Chen 2023), as shown in Fig. 5. The Eelec between liquids and Iβ-MF showed a positive correlation with the polarity of liquids. Due to the lack of hydrogen bond donors for DMF and pyridine liquids, the hydrogen bond formation ability was poor compared to short-chain alcohols and short-chain carboxylic acids, resulting in a weaker electrostatic interaction (attraction). The weakening of this interaction was not found to be significantly different in the simulation results for C36 for a short time (Fig. 2a). However, in the G06 model (Fig. 2b), the retention of the tg conformation confirmed that electrostatic interactions between the liquid and the Iβ-MF dominated the Iβ-MF structural changes. Dispersive interactions only played an auxiliary role, and electrostatic interactions were almost non-existent in low polarity liquids. The conformational change of Iβ-MF was guided by dispersive interactions, and about 10% of tg conformations were transformed in the C36 model. In fact, the relationship between liquid polarity and conformational change is energy-driven and dominated by electrostatic interactions between the liquid and Iβ-MF. Moreover, the Etotal of the G06 simulation was generally larger than that of the C35 simulation, which was caused by the large structural changes of the G06 model, and the Etotal difference between C35 and G06 was smaller at the beginning of the simulation.

Fig. 5. Interaction energy between liquid molecules and Iβ-MF in a) CHARMM36 and b) GLYCAM06

Large changes in atomic charge would inevitably produce large biases, and misleading simulation results would exist (Bayly et al. 1993; Zhang et al. 2023). However, small perturbations in a large system might not produce detectable biases. If the atomic charge of the molecule was near 0, then small perturbations might produce simulation results with large deviations. In the HCCl3 system, the RESP2(0.1) and RESP2(0.9) charges were calculated, when the average charges of the Cl atoms were -0.0095 and 0.004, respectively, while the RESP charges of the C atom were -0.2323 and -0.2713, respectively.

For the simulation of HCCl3 systems with different charges, the Eelec of the C36 simulation differed by 463 kJ mol-1 with no significant change in the tg conformation. The difference in Eelec for G06 was as high as 4040 kJ mol-1 and the difference in tg conformation was 30%. Atomic charge variations lead to energy changes, which result in different interactions between molecules. The C36 was less sensitive to atomic charges and could withstand some small perturbations. The G06 reduced the 1,4 electrostatic scaling factor, and small changes in the values near the atomic charge critical point produced large deviations. In addition, simulations were performed for H2CCl2, CCl4, and CS2 under the same conditions with deviation values of 50 kJ mol-1, 332 kJ mol-1, and 71 kJ mol-1, respectively.

Twist of Cellulose Microfiber

The hydrated cellulose microfibrils in both force fields were twisted at 1 ns moment, and the twisting direction is the same. The non-reduced end (distal end, ball-and-stick model) was kept in its original position. The chain orientation is indicated in Fig. 6 with blue and white spiral lines from the non-reduced end, ending at the reducing end. This finding was very similar to the work of Paavilainen et al. (2011). The distortion of the G06 model was significantly higher than that of the C36 model at 1 ns.

Fig. 6. Schematic diagram of the twisting of cellulose microfibrils at the 1ns. The non-reduced end was represented as a ball-and-stick model and placed at the distal end, and the blue and white spiral lines indicated the course of the cellulose chain, with the end being the reduced end.

Inspired by Matthews et al. (2006), the middle 5 cellobiose units were chosen to calculate the overall twist of this short oligosaccharide segment. The dihedral angle C111-C111-C121-C121 (indicated by the black line in Fig. S4) was used to define the twist, and the average twist values of the last 5 ns are summarized in Fig. 6. The difference in a twist between the center and origin layers of the model under the two force fields in a single simulation was slight and would not be discussed separately. The C36 model had little difference between the overall twist in different systems, except for the formic acid liquid system. In contrast, the overall twist in the G06 model was not related to the liquid polarity. There were data with a large degree of deviation in three replicate simulations (solid black squares in Fig. 7). The results suggest that the Iβ-MF would go through twisting and untwisting stages, and the C36 model was basically in the twisting stage.

A part of the G06 model had already been untwisted. For example, in Fig. S5a, the G06 model underwent rapid untwisting at about 20 ns in the presence of water molecules, and the tg conformation was fully transformed at the same moment (Fig. S3, Inside the black box), and it was clear that the conformational transition influences the untwisting process. Combining Fig. 7 and Fig. 2c revealed that the less tg conformation contained in the same system, the lower the twisting degree (except for the water system), which undoubtedly indicated that the degree of transformation of tg conformation was one of the essential factors affecting the twist degree. In the water liquid system, three repeated calculations showed the complete transformation of the tg conformation. Still, the overall twist of one of them was 4.79, which did not reach the state of untwisting because no alternate layer structure was formed in the 8th layer of this calculation. According to the trend, this layer should be mostly in gg conformation (e.g., Fig. 3), but in fact, it was mostly in gt conformation (Fig. S6, inside the red box), resulting in a slight increase in the overall twist. Therefore, the untwisting of cellulose microfibers requires two conditions: 1) almost complete transformation of the tg conformation and 2) the formation of the alternate layer structure.

Fig. 7. Average twist values for the last 5 ns of a) CHARMM36 and b) GLYCAM06 simulations. The anomalies in the b-plot correspond precisely to the points with large deviations in Fig. 3c.


Because the atomic charge is one of the important factors affecting electrostatic interactions, the effect of the atomic charge at the 0 points on the simulation results was explored. In the HCCl3 system, changing the Cl atomic charge (-0.0095 and 0.004) resulted in a difference of 463 kJ mol-1 in Eelec between the two simulations (C36), with no significant change in the tg conformation. However, the difference in Eelec between the two simulations of G06 was up to 4040 kJ mol-1 and the tg conformation differed by 30%. The change in atomic charge led to a change in energy and thus affected the intermolecular interactions. A charge change of about 0.05 in the G06 simulation could produce large deviations, although such significant deviations did not occur in other systems. Therefore, it was necessary to focus on the charge at the location of the critical point during the charge calculation to avoid obtaining results with large deviations.

It was clear that the reduction of the 1,4 electrostatic scaling factors weakened the torsional potential of the cellulose structure and made the whole system more chaotic. By linking this “acceleration” effect with heating, it was easy to understand why an alternating layer structure similar to the high-temperature phase intermediate I-HT appeared. This alternate layer structure had two characteristics: 1) the center layer was occupied by the gg conformation and the origin layer was occupied by the gt conformation; 2) the cellulose chains of the center and origin layers were inclined counterclockwise and clockwise, respectively. The hydroxymethyl group in the gg conformation was perpendicular to the mean plane of the pyran ring and acted like an “arm” to capture the cellulose chain of the origin layer through the O6centerH6center-O2origin hydrogen bond, forcing the cellulose chain of the origin layer to tilt clockwise. When the origin layer was tilted to a certain extent, the origin layer’s hydroxymethyl group would have a spatial conflict with the center layer, resulting in a slight counterclockwise tilt of the center layer. The formation of the alternating layer structure appears to untwist. The four groups of systems (water, formic acid, hexane, cyclohexane) were re-simulated using a 1,4 electrostatic scaling factor of 1.0, and the simulation results are summarized in Table 3. The G06 simulation results were very similar to those with C36 at the same 1,4 electrostatic scaling factor and the formic acid system could be explained. Changing the electrostatic scaling factor of GLYCAM06 could affect the thermodynamics and kinetics, both of which have a significant impact on the study of cellulose structures with MD simulations.

Table 3. G06 Simulation Data with a 1,4 Electrostatic Scale Factor of 1