PEER-REVIEW ARTICLE


 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8223 

 

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: shyfu@scut.edu.cn 

 

 
INTRODUCTION 
 

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 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8224 

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 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8225 

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. 

 

 
EXPERIMENTAL 
 

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), 

MPI=
1

𝐴
∬ |𝑉(𝑟)|𝑑𝑆
𝑆

        (1) 

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 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8226 

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 

a) b) 

c) 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8227 

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, 

Etotal=Eelec+Evdw        (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. 

 

 
RESULTS AND DISCUSSION 
 

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.  



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8228 

Table 1. Polarity Data for 16 Liquid Molecules 
 

 Atotal (Å2) pApolar 
a) 

Vmin b) 

(kcal mol-1) 
Vmax b) 

(kcal mol-1) 
MPI (kcal 
mol-1) 

Water 45.7 81.6% -36.3 42.6 22.8 

Formic acid 74.0 69.6% -31.0 52.0 16.6 

DMF 122.4 72.6% -42.8 21.9 15.4 

Acetic acid 96.5 65.8% -33.9 47.2 14.5 

Methanol 72.2 46.5% -33.4 41.5 13.2 

Propionic 
acid 

117.4 50.3% -33.8 46.6 12.1 

Ethanol 94.9 36.5% -33.7 39.8 11.0 

Pyridine 122.2 45.7% -36.2 20.7 10.4 

n-Propanol 116.4 28.9% -33.8 39.4 9.7 

H2CCl2 101.7 34.7% -12.4 27.0 9.5 

Benzene 127.1 34.7% -15.4 13.8 7.9 

HCCl3 119.3 15.5% -7.1 32.5 6.5 

CCl4 134.7 9.7% -3.2 17.5 3.8 

CS2 98.2 6.8% -1.2 14.6 2.9 

n-Hexane 171.0 0% -2.9 6.9 2.8 

Cyclohexane 145.9 0% -2.6 5.9 2.5 
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β 
 

 Cellulose Iβ a) G06 b) C36 b) 

a (Å) 7.784 7.849 (0.13) 7.894 (0.03) 

b (Å) 8.201 8.198 (0.03) 8.393 (0.02) 

c (Å) 10.380 10.696 (0.005) 10.423 (0.015) 

α (°) 90 
89.5 (0.30) c) 

89.9 (0.36) d) 

89.5 (0.17) c) 

90.2 (0.16) d) 

β (°) 90 
88.8 (0.88) c) 

89.7 (0.45) d) 

89.1 (0.26) c) 

89.6 (0.13) d) 

γ (°) 96.5 97.6 (0.66) 97.6 (0.29) 

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 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8229 

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. 

a) 

b) 

c) 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8230 

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). 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8231 

 

 
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 

a) 

c) 

b) 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8232 

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 

a) 

b) 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8233 

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 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8234 

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. 

a) 

b) 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8235 

DISCUSSION 
 

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 

 tg conformation gt conformation gg conformation Twist (°) 

Water 744 (68.9%) 191 (17.7%) 145 (13.4%) 11.23 

Formic acid 774 (71.7%) 158 (14.6%) 148 (13.7%) 9.03 

n-Hexane 832 (77.0%) 60 (5.6%) 188 (17.4%) 6.69 

Cyclohexane 857 (79.3%) 58 (5.4%) 165 (15.3%) 6.48 

 

In experiments to identify cellulose I by CP/MAS 13C NMR spectroscopy, the 

spectra of cotton staple immersed in methanol in the C-4 region showed the largest 

deviation from the spectra of cotton staple immersed in water in the spectral region typical 

of the unordered form of cellulose (80 to 86 ppm). The typical spectral region for the 

ordered form of cellulose (86 to 92 ppm) remained relatively intact. Crystalline cellulose 

was not substantially altered by solvent exchange, i.e., the crystalline form was an “in-core” 

or “in-fiber” form that was inaccessible to the surrounding solvent (Wickholm et al. 1998). 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8236 

Some experimental evidence has suggested that cellulose microfibrils are twisted. 

Narrow bands of twisted fibers were observed during the biosynthesis of bacterial cellulose 

(Benziman et al. 1980). Hanley and Gray reported the right-handed twisting of 

Micrasterias denticulata microfibrils using atomic force and transmission electron 

microscopy (Hanley et al. 1997). The authors suggested that the twisting may be an 

artificial product of sample preparation and that when microfibers are suspended in an 

aqueous solution before drying onto the substrate for study, they will exhibit a smooth, 

uniform twist along their length. The untwisting process of cellulose microfibril was more 

exciting, and although this simulation result was not necessarily correct, it should help to 

strengthen the general understanding of the structural behavior of I-HT. 

Finally, without further experimental data, it was impossible to state which force 

field simulation results were “correct definitively”. Microfibrils’ size and surface 

morphology might significantly impact the observable results obtained through molecular 

dynamics simulations, so further experimental work will be needed to understand the 

morphology of small-diameter cellulose microfibrils to get a better comparative reference 

and to be of great value for more accurate cellulose simulation. The results of this study 

should provide insight into the behavior of cellulose microfibers in different polar liquids 

and suggest how MD simulations can best be applied to cellulose studies. 

 

 

CONCLUSIONS 
 

Two force fields, G06 and C36, were used to study the behavior of cellulose Iβ 

microfibrils in non-solvent liquids, including the conformational transition and driving 

forces behind the cellulose microfibrils, the interaction with each liquid, and the fiber 

twisting. Both force field models were twisted in the equilibrium phase at the beginning of 

the simulation, and the twist direction was consistent. With a twist, the hydroxymethyl 

conformation inside the cellulose microfibril was transformed (initial conformation: tg). 

The C36 model retained more than 70% of the tg conformation in 16 liquids, and the 

retention of the tg conformation increased with decreasing liquid polarity. In fact, the 

ability of liquid polarity to affect the hydroxymethyl conformation of cellulose was only 

an appearance, and the real driving force behind it was the electrostatic interaction between 

liquid molecules and cellulose. In general, the electrostatic interaction was closely related 

to the hydrogen bonding ability of the molecule and had some connection with the 

molecular polarity. Since DMF and pyridine molecules lack hydrogen bond donors, the 

electrostatic interaction energy (absolute value) was small compared to similarly polar 

liquid molecules, but the simulation data of C36 did not show conclusive evidence in a 

short time. Therefore, another set of simulations was designed using the G06 force field 

with a 5/6 1,4 electrostatic scaling factor, with the aim of reducing the barriers to tg 

conformational rotation and obtaining an "acceleration" effect. In this case, the tg 

conformation of the model in highly polar liquids (water, formic acid, methanol, etc.) 

underwent a complete transformation, while about 60% of the tg conformation of the low-

polar liquids was transformed. The DMF and pyridine systems were amplified by 

electrostatic interactions and showed off-trend results, which further verified that the 

driving force of the conformational transformation of cellulose molecules is the 

electrostatic interaction between liquid molecules and cellulose. 

 

 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8237 

ACKNOWLEDGMENTS 
 

This work was supported by the National Key Research and Development Program 

(2021YFE0104500), the National Natural Science Foundation of China (22078114), Key 

Research and Development Program of Guangzhou Science and Technology Program 

(202103000011), and the Natural Science Foundation of Guangdong Province 

(2021A1515010360). 

 

 

REFERENCES CITED 
 

Bayly, C. I., Cieplak, P., Cornell, W., and Kollman, P. A. (1993). "A well-behaved 

electrostatic potential based method using charge restraints for deriving atomic 

charges: the RESP model," The Journal of Physical Chemistry 97(40), 10269-10280. 

DOI: 10.1021/j100142a004 

Bellesia, G., Chundawat, S. P. S., Langan, P., Dale, B. E., and Gnanakaran, S. (2011). 

"Probing the early events associated with liquid ammonia pretreatment of native 

crystalline cellulose," The Journal of Physical Chemistry B 115(32), 9782-9788. DOI: 

10.1021/jp2048844 

Benziman, M., Haigler, C. H., Brown, R. M., White, A. R., and Cooper, K. M. (1980). 

"Cellulose biogenesis: Polymerization and crystallization are coupled processes in 

Acetobacter xylinum," Proceedings of the National Academy of Sciences 77(11), 

6678-6682. DOI: 10.1073/pnas.77.11.6678 

Cosgrove, D. J. (2014). "Re-constructing our models of cellulose and primary cell wall 

assembly," Current Opinion in Plant Biology 22, 122-131. DOI: 

10.1016/j.pbi.2014.11.001 

Ding, S.-Y., Zhao, S., and Zeng, Y. (2014). "Size, shape, and arrangement of native 

cellulose fibrils in maize cell walls," Cellulose 21(2), 863-871. DOI: 10.1007/s10570-

013-0147-5 

Gardner, K. H., and Blackwell, J. (1974). "The structure of native cellulose," 

Biopolymers 13(10), 1975-2001. DOI: 10.1002/bip.1974.360131005 

Gomes, T. C. F., and Skaf, M. S. (2012). "Cellulose-Builder: A toolkit for building 

crystalline structures of cellulose," Journal of Computational Chemistry 33(14), 

1338-1346. DOI: 10.1002/jcc.22959 

Gross, A. S., Bell, A. T., and Chu, J.-W. (2011). "Thermodynamics of cellulose solvation 

in water and the ionic liquid 1-butyl-3-methylimidazolim chloride," The Journal of 

Physical Chemistry B 115(46), 13433-13440. DOI: 10.1021/jp202415v 

Gross, A. S., and Chu, J.-W. (2010). "On the molecular origins of biomass recalcitrance: 

The interaction network and solvation structures of cellulose microfibrils," The 

Journal of Physical Chemistry B 114(42), 13333-13341. DOI: 10.1021/jp106452m 

Hadden, J. A., French, A. D., and Woods, R. J. (2013). "Unraveling cellulose 

microfibrils: A twisted tale," Biopolymers 99(10), 746-756. DOI: 10.1002/bip.22279 

Hanley, S. J., Revol, J.-F., Godbout, L., and Gray, D. G. (1997). "Atomic force 

microscopy and transmission electron microscopy of cellulose from Micrasterias 

denticulata; evidence for a chiral helical microfibril twist," Cellulose 4(3), 209-220. 

DOI: 10.1023/A:1018483722417 

 

 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8238 

Horii, F., Hirai, A., and Kitamaru, R. (1987). "CP/MAS carbon-13 NMR spectra of the 

crystalline components of native celluloses," Macromolecules 20(9), 2117-2120. 

DOI: 10.1021/ma00175a012 

Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., de Groot, B. L., Grubmüller, 

H., and MacKerell, A. D. (2017). "CHARMM36: An improved force field for folded 

and intrinsically disordered proteins," Biophysical Journal 112(3, Supplement 1), 

175a-176a. DOI: 10.1016/j.bpj.2016.11.971 

Humphrey, W., Dalke, A., and Schulten, K. (1996). "VMD: Visual molecular dynamics," 

Journal of Molecular Graphics 14(1), 33-38. DOI: /10.1016/0263-7855(96)00018-5 

Imai, T., and Sugiyama, J. (1998). "Nanodomains of Iα and Iβ cellulose in algal 

microfibrils," Macromolecules 31(18), 6275-6279. DOI: 10.1021/ma980664h 

Kamide, K. (2005). “1 – Introduction,” in: Cellulose and Cellulose Derivatives, K. 

Kamide (ed.), Elsevier, Amsterdam, pp. 1-23. 

Kirschner, K. N., Yongye, A. B., Tschampel, S. M., González-Outeiriño, J., Daniels, C. 

R., Foley, B. L., and Woods, R. J. (2008). "GLYCAM06: A generalizable 

biomolecular force field. Carbohydrates," Journal of Computational Chemistry 29(4), 

622-655. DOI: /10.1002/jcc.20820 

Kutzner, C., Pall, S., Fechner, M., Eszternnann, A., de Groot, B. L., and Grubmuller, H. 

(2019). "More bang for your buck: Improved use of GPU nodes for GROMACS 

2018," Journal of Computational Chemistry 40(27), 2418-2431. DOI: 

10.1002/jcc.26011 

Lee, J., Cheng, X., Swails, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., Wei, S., 

Buckner, J., Jeong, J. C., Qi, Y., Jo, S., Pande, V. S., Case, D. A., Brooks, C. L., III, 

MacKerell, A. D., Jr., Klauda, J. B., and Im, W. (2016). "CHARMM-GUI Input 

Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM 

Simulations Using the CHARMM36 Additive Force Field," Journal of Chemical 

Theory and Computation 12(1), 405-413. DOI: 10.1021/acs.jctc.5b00935 

Liu, Z., Lu, T., and Chen, Q. (2021). "Intermolecular interaction characteristics of the all-

carboatomic ring, cyclo[18]carbon: Focusing on molecular adsorption and stacking," 

Carbon 171, 514-523. DOI: /10.1016/j.carbon.2020.09.048 

Lu, T., and Chen, F. W. (2012). "Multiwfn: A multifunctional wavefunction analyzer," 

Journal of Computational Chemistry 33(5), 580-592. 10.1002/jcc.22885 

Lu, T., and Chen, Q. (2023). "Visualization analysis of weak interactions in chemical 

systems,” in: Reference Module in Chemistry, Molecular Sciences and Chemical 

Engineering, Elsevier. DOI: 10.1016/B978-0-12-821978-2.00076-3 

MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. 

J., Fischer, S., Gao, J., Guo, H., Ha, S., Joseph-McCarthy, D., Kuchnir, L., Kuczera, 

K., Lau, F. T., Mattos, C., Michnick, S., Ngo, T., Nguyen, D. T., Prodhom, B., 

Reiher, W. E., Roux, B., Schlenkrich, M., Smith, J. C., Stote, R., Straub, J., 

Watanabe, M., Wiórkiewicz-Kuczera, J., Yin, D., and Karplus, M. (1998). "All-atom 

empirical potential for molecular modeling and dynamics studies of proteins," J Phys 

Chem B 102(18), 3586-3616. DOI: 10.1021/jp973084f 

Matthews, J. F., Beckham, G. T., Bergenstråhle-Wohlert, M., Brady, J. W., Himmel, M. 

E., and Crowley, M. F. (2012). "Comparison of cellulose Iβ simulations with three 

carbohydrate force fields," Journal of Chemical Theory and Computation 8(2), 735-

748. DOI: 10.1021/ct2007692 

Matthews, J. F., Bergenstråhle, M., Beckham, G. T., Himmel, M. E., Nimlos, M. R., 

Brady, J. W., and Crowley, M. F. (2011). "High-temperature behavior of cellulose I," 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8239 

The Journal of Physical Chemistry B 115(10), 2155-2166. DOI: 10.1021/jp1106839 

Matthews, J. F., Skopec, C. E., Mason, P. E., Zuccato, P., Torget, R. W., Sugiyama, J., 

Himmel, M. E., and Brady, J. W. (2006). "Computer simulation studies of 

microcrystalline cellulose Iβ," Carbohydrate Research 341(1), 138-152. DOI: 

10.1016/j.carres.2005.09.028 

Matthews, J. F., Himmel, M. E., and Crowley, M. F. (2012). "Conversion of cellulose Iα 

to Iβ via a high temperature intermediate (I-HT) and other cellulose phase 

transformations," Cellulose 19(1), DOI: 297-306. 10.1007/s10570-011-9608-x 

Mazeau, K., and Rivet, A. (2008). "Wetting the (110) and (100) surfaces of Iβ cellulose 

studied by molecular dynamics," Biomacromolecules 9(4), 1352-1354. DOI: 

10.1021/bm7013872 

Nishiyama, Y., Johnson, G. P., French, A. D., Forsyth, V. T., and Langan, P. (2008). 

"Neutron crystallography, molecular dynamics, and quantum mechanics studies of the 

nature of hydrogen bonding in cellulose Iβ," Biomacromolecules 9(11), 3133-3140. 

DOI: 10.1021/bm800726v 

Paavilainen, S., Róg, T., and Vattulainen, I. (2011). "Analysis of twisting of cellulose 

nanofibrils in atomistic molecular dynamics simulations," The Journal of Physical 

Chemistry B 115(14), 3747-3755. DOI: 10.1021/jp111459b 

Sarko, A., Southwick, J., and Hayashi, J. (1976). "Packing analysis of carbohydrates and 

polysaccharides. 7. Crystal structure of cellulose IIII and its relationship to other 

cellulose polymorphs," Macromolecules 9(5), 857-863. DOI: 10.1021/ma60053a028 

Shefter, E., and Trueblood, K. N. (1965). "The crystal and molecular structure of D(+)-

barium uridine5'-phosphate," Acta Crystallogr 18, 1067-1077. DOI: 

10.1107/s0365110x65002530 

Wickholm, K., Larsson, P. T., and Iversen, T. (1998). "Assignment of non-crystalline 

forms in cellulose I by CP/MAS 13C NMR spectroscopy," Carbohydrate Research 

312(3), 123-129. DOI: 10.1016/S0008-6215(98)00236-5 

Yui, T., Nishimura, S., Akiba, S., and Hayashi, S. (2006). "Swelling behavior of the 

cellulose Iβ crystal models by molecular dynamics," Carbohydrate Research 341(15), 

2521-2530. DOI: 10.1016/j.carres.2006.04.051 

Zhang, Q., Bulone, V., Ågren, H., and Tu, Y. (2011). "A molecular dynamics study of the 

thermal response of crystalline cellulose Iβ," Cellulose 18(2), 207-221. DOI: 

10.1007/s10570-010-9491-x 

Zhang, Q.-K., Zhang, X.-Q., Wan, J., Yao, N., Song, T.-L., Xie, J., Hou, L.-P., Zhou, M.-

Y., Chen, X., Li, B.-Q., Wen, R., Peng, H.-J., Zhang, Q., and Huang, J.-Q. (2023). 

"Homogeneous and mechanically stable solid–electrolyte interphase enabled by 

trioxane-modulated electrolytes for lithium metal batteries," Nature Energy 8(7), 725-

735. DOI: 10.1038/s41560-023-01275-y 

Zhang, Y., He, H., Liu, Y., Wang, Y., Huo, F., Fan, M., Adidharma, H., Li, X., and 

Zhang, S. (2019). "Recent progress in theoretical and computational studies on the 

utilization of lignocellulosic materials," Green Chemistry 21(1), 9-35. DOI: 

10.1039/C8GC02059K 

Zhao, Z., Shklyaev, O. E., Nili, A., Mohamed, M. N. A., Kubicki, J. D., Crespi, V. H., 

and Zhong, L. (2013). "Cellulose microfibril twist, mechanics, and implication for 

cellulose biosynthesis," The Journal of Physical Chemistry A 117(12), 2580-2589. 

DOI: 10.1021/jp3089929 

  



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8240 

Zhou, S., Jin, K., and Buehler, M. J. (2021). "Understanding plant biomass via 

computational modeling," Advanced Materials 33(28), 2003206. DOI: 

10.1002/adma.202003206 

 

Article submitted: August 3, 2023; Peer review completed: September 23, 2023; 

Revisions received: October 7, 2023; Revisions accepted: October 8, 2023; Published: 

October 19, 2023. 

DOI: 10.15376/biores.18.4.8223-8248 

 

 

  



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8241 

SUPPLEMENTARY 
 
APPENDIX 
 

 
Fig. S1. RMSD plots of cellulose Iβ microfibers in a) CHARMM36 and b) GLYCAM06 force fields 
  

a) 

b) 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8242 

 
Fig. S2. Snapshots of cellulose Iβ microfibers per 100 layers at 50 ns under two force fields. 
Snapshots of the CHARMM36 and GLYCAM06 force fields are shown on the left and right of 
each system, respectively. The non-reducing end is arranged at the left end. 
  



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8243 

 
Fig. S3. The plot of the number of hydroxymethyl conformations of cellulose microfibrils in 16 
solvents as a function of time. The CHARMM model is shown above the blue line, and the 
GLYCAM06 model is shown below. The gt, tg, and gg conformations are indicated by the green, 
yellow, and purple lines, respectively.  



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8244 

 
Fig. S4. Calculation of the overall twist of a) center layer and b) origin layer using C111-C111-C121-
C121 dihedral angle. The reducing end (upper end) and the non-reducing end (upper end) are 
labeled glucose #1 and #30, respectively. 
  

a) b) 

 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8245 

 
Fig. S5. Twist versus time for cellulose microfibers in water and hexane solvent system. a) and c) 
belong to the GLYCAM 06 model, b) and d) belong to the CHARMM 36 model 
 
  

a) b) 

c) d) 



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8246 

 

 
Fig. S6. Snapshots of two replicate simulations of cellulose Iβ microfibrils. The time is 50 ns, the 
solvent is water, and the force field is GLYCAM06. The differences between the two simulations 
are inside the red box line, resulting in different twist data.  



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8247 

Table S1. Lattice Parameters of Cellulose Microfibrils for All Solvent Systems 
 a (Å) b (Å) c (Å) α (°) β (°) γ (°) 

GLYC
AM06 

Water 8.000 8.189 10.701 90.0/90.6 90.1/90.4 98.4 

Formic acid 7.993 8.152 10.696 89.7/89.9 89.5/89.8 98.5 

DMF 7.942 8.212 10.693 90.1/90.2 87.9/88.4 98.0 

Acetic acid 7.979 8.174 10.690 89.5/89.9 89.3/89.7 98.3 

Methanol 7.907 8.169 10.695 89.5/90.1 89.6/90.3 98.0 

Propanoic Acid 7.820 8.195 10.696 89.7/90.3 87.4/89.6 97.5 

Ethanol 7.987 8.153 10.698 89.7/89.8 89.9/89.7 98.3 

Pyridine 7.865 8.206 10.698 89.2/90.0 87.3/89.8 96.8 

Propanol 7.989 8.158 10.697 89.7/89.8 89.8/89.7 98.2 

H2CCl2 7.896 8.189 10.688 89.0/89.1 88.6/90.2 97.9 

Benzene 7.772 8.225 10.701 89.3/89.9 87.8/89.6 97.3 

HCCl3 7.762 8.220 11.687 89.1/89.7 89.0/89.3 96.5 

CCl4 7.652 8.236 10.699 89.5/89.4 88.7/89.5 97.2 

CS2 7.642 8.226 10.703 89.4/89.7 89.4/89.9 97.4 

n-Hexane 7.734 8.234 10.694 89.3/90.1 88.4/90.0 96.8 

Cyclohexane 7.642 8.226 10.703 89.8/90.4 88.2/89.5 96.6 

CHAR
MM36 

Water 7.933 8.385 10.417 89.2/90.2 89.5/89.6 97.7 

Formic acid 7.815 8.346 10.367 89.6/89.9 89.7/90.0 97.1 

DMF 7.875 8.379 10.422 89.5/90.5 88.7/89.5 97.7 

Acetic acid 7.904 8.384 10.426 89.4/90.3 88.8/89.5 96.7 

Methanol 7.880 8.388 10.424 89.3/90.4 89.2/89.7 97.5 

Propanoic Acid 7.870 8.384 10.426 89.6/90.4 88.8/89.6 97.7 

Ethanol 7.877 8.385 10.426 89.7/90.4 89.1/89.7 97.5 

Pyridine 7.962 8.417 10.428 89.4/90.2 89.2/89.7 97.5 

Propanol 7.881 8.384 10.426 89.7/90.3 88.9/89.6 97.6 

H2CCl2 7.901 8.407 10.428 89.3/90.2 89.4/89.7 97.7 

Benzene 7.903 8.407 10.432 89.6/90.2 89.0/89.7 97.8 

HCCl3 7.945 8.407 10.429 89.1/90.0 89.2/89.4 97.8 

CCl4 7.886 8.403 10.430 89.5/90.1 89.1/89.6 97.8 

CS2 7.877 8.404 10.429 89.5/90.1 89.1/89.6 97.7 

n-Hexane 7.903 8.404 10.431 89.3/90.0 89.4/89.7 97.8 

Cyclohexane 7.895 8.400 10.429 89.6/90.1 89.0/89.7 97.8 

  



 

PEER-REVIEWED ARTICLE  bioresources.com 

 

 

Kong et al. (2023). “Cellulose I surface simulation,” BioResources 18(4), 8223-8248.  8248 

Table S2. Conformational Data from 5 Times CHARMM 36 Simulations in Formic 
Acid Solvent 

 
Gromacs 
Version 

GPU 
Acceleration 

Computer 
Operating System 

tg gt gg 

1 2022.6 Yes Win 10 35 1034 11 

2 2020.6 Yes Win 10 25 1045 10 

3 2020.6 Yes Centos 7 14 1058 8 

4 2022.6 No Linux 13 1059 8 

5 2020.6 Yes Centos 7 130 922 27