Effect of moisture content on the microscopic properties of amorphous cellulose: a molecular dynamics simulations

Performance degradation of cellulose and cellulose-based materials caused by water is an inevitable problem in application processes. In most studies, this was attributed to the fracture and rebuilding of the hydrogen bond network in the system; however, limited attention was paid to the movement, aggregation state, and specific property evolution of cellulose and water during this process. In this study, based on molecular dynamics simulations, the effects of moisture content on the microscopic properties of cellulose are investigated, including the mechanical properties, diffusion coefficient, glass transition temperature, microscopic motion of water molecules, and preferred hydration sites on cellulose. The results show that the mechanical properties of the system increase and then decrease as the water content in the system increases. When the moisture content is 4%, the mechanical properties of cellulose are the best, and the elastic modulus and shear modulus increase by 7.6% and 9.4%, respectively, compared with those of dry cellulose. The glass transition temperature of the system with 4% moisture content increases by 72 K compared with that of dry cellulose. The mean square displacement and diffusion coefficient of water in the system is affected by the water molecules’ polymerisation state and the free water content. In the entire range of water contents studied, hydroxyl groups O2, O3, and O6 of cellulose dominate the reaction with water compared with acetal oxygens O4 and O5. In the system with 4% moisture content, the number of water molecules around the glycosidic bonds O4 are the most minor and cause the least damage to the cellulose structure. A critical water content point of 4% is recommended, and this result is expected to provide a reference for maintaining the excellent and stable properties of cellulose and cellulose-based materials.


Introduction
The energy demand is increasing daily due to the advent of the information age and the rapid development of the industry. It has been reported that the total energy consumption in the 20th century was almost equal to half of the energy consumed in the previous 19 centuries, and the amount of energy consumed daily by humans is unimaginable. Moreover, the massive use of fossil fuels has caused severe environmental damage, resulting in acid rain and global warming, which have seriously damaged human habitats. An increasing number of people now realise that only by developing green, sustainable energy and materials we can avoid a series of problems caused by resource depletion and environmental degradation and enable human beings to survive and develop for a long time. Cellulose is the most abundant renewable biomass resource in the world and has contributed significantly to energy. Cellulosic biofuels are considered the most promising alternative to traditional fossil fuels. In addition, research on cellulose has continued in the past decades because cellulose exhibits excellent usage performance and can be used as a substrate to synthesise composite materials with superior properties [1][2][3]. However, in the practical application of cellulosic materials and composites built on cellulose substrates, Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. material stability and performance degradation due to environmental humidity is inevitable [4]. This is due to the unique structure of the cellulose chain, which easily forms hydrogen bonds with the exposed hydroxyl (O2 and O3) and hydroxymethyl groups (O6) on the cellulose chain, causing water molecules to cluster around the cellulose chains (the oxygen atoms labelled in figure 1(a)) [5], which destroys the original structure of cellulose chains, resulting in reduced reliability and performance in the process of material used. Existing research is devoted to studying cellulose degradation under water-containing conditions from the perspective of hydrogen bonds. Ebrahimzadeh et al [6] analysed the dynamic mechanical response of a sample paper during the adsorption and desorption processes. Moisture adsorption leads to sample swelling, and a mechanical loss peak is detected during this process. This is due to the disruption of load-bearing hydrogen bonds by water molecules in the sample. Kulasinski et al [5] reported the relationship between volumetric strain, water content, and porosity of cellulose during the hygroscopic process. The study results revealed that swelling is directly related to the space created by adsorbed water molecules, and the interaction of water with the cellulose structure leads to the degradation of the mechanical properties of the system. When the water content was increased from 0% to 50%, the amorphous cellulose's bulk modulus and shear modulus decreased by 49% and 96%, respectively. This is attributed to weakening the mechanical properties due to breaking hydrogen bonds within the system. Wang et al [7] conducted molecular dynamics simulations on pure amorphous cellulose and amorphous cellulosewater, and calculated both the glass transition temperature (Tg) and mechanical properties. The glass transition temperature of the aqueous cellulose decreased by 30 K. The presence of water reduced the mechanical strength of the cellulose, lowering the elastic modulus from 13.64 GPa to 13.06 GPa while slightly increasing its flexibility and plasticity. The Poisson's ratio and C 12 -C 44 ratio increased by 1% and 19.5%, respectively. This change is attributed to the breakdown of hydrogen bonds between cellulose chains. The presence of water reduces the number of hydrogen bonds between the cellulose chains from 140 to 120 approximately. Most existing studies attribute the problem of moisture-induced cellulose degradation to the disruption and reformation of the hydrogen bond lattice in the system. Still, Malin et al [8] argued that the role of hydrogen bonds in cellulose research had been exaggerated. In many studies, hydrogen bonds have been frequently used to explain various phenomena and properties of cellulose and cellulose-based materials. However, in Malin's view, hydrogen bonds are only one of the weak forces among several interactions, and the intrinsic causes of changes in the properties of cellulosic materials cannot be entirely attributed to differences in hydrogen bonds within the system and should partly separate the focus from hydrogen bonds.
Mazeau [9] constructed a composite model of amorphous cellulose with different water combinations and studied the state of aggregation of water molecules, the mobility of water molecules, and the effect of water on the glass transition temperature of cellulose. They found that the state of water in the cellulose system ranged from single to aggregated into clusters and then formed continuous channels. The agglomeration phenomenon of water molecules caused the diffusion coefficient of water to decrease from 1.08 × 10 -11 m 2 s −1 to 0.13 × 10 −11 m 2 s −1 during the increase of water content of cellulose from 0.83% to 18.14%. Dan et al [10] investigated the effect of the water state on the mechanical properties of cellulose using molecular dynamics simulations. The hydrophilic sites of cellulose have been completely occupied as water content increases, and the water molecules exist more in the form of free water. The powerful diffusion ability of free water destroys the stable structure of the cellulose chains to a greater extent, forming a plasticization effect that causes a decrease of 31% and 26% in the elastic modulus and shear modulus of the cellulose system, respectively. The study of the effect of water content is of interest not only for the further development of cellulosic materials but also for the development of new composite materials. Bulut et al [11] focused on water adsorption tests in their research to develop carbon/ aramid fibre-reinforced hybrid composites to ensure the stability of materials for use in hygrothermal environments. It is necessary to report the motion of cellulose chains and water and the state of water aggregation at the molecular level to understand better the mechanical degradation behaviour of cellulose in the presence of moisture.
In this study, molecular dynamics simulations were carried out for amorphous cellulose systems with different moisture contents. Materials Studio (MS) software was used to simulate the molecular dynamics processes. The microscopic properties of the different moisture content models were obtained, including mechanical properties, diffusion coefficients, and glass transition temperature. Next, the form of water in the composite system and the preferred hydration sites of cellulose under different moisture contents were analysed. Finally, the coordination number of water molecules around the glycosidic bonds O4 was calculated to evaluate the mechanical degradation mechanism of cellulose further. Combined with the research results, the moisture content point (4%) for cellulose materials to maintain stable and superior performance in practical applications was derived.

Details of the simulation 2.1. Model
Cellulose is a linear chain of glucose molecules linked together through a β 1-4 glucosidic bond. The repeat unit (figure 1(a)) is comprised of two anhydroglucose rings ((C 6 H 10 O 5 )n, n = 10 000 to 15 000, where n is determined by the cellulose source material) [12]. Natural cellulose has a two-phase structure that includes crystalline and amorphous regions. The cellulose chain is arranged tightly and neatly in the crystalline region, and the inter-and intra-molecular interactions are significant. Therefore, water molecules can only be adsorbed on its surface, which does not effect the properties of cellulose [9]. Water molecules can easily act on the amorphous areas of cellulose chains and destroy the original molecular structure of cellulose, resulting in uncontrollable deformation and degradation of the mechanical properties. Therefore, the amorphous region of cellulose was selected as the model material in this study. Mazeau et al [13] performed molecular dynamics simulations for cellulose with degrees of polymerisation of 10, 20, and 40. The results showed no significant differences in the physicochemical properties exhibited by the cellulose models with different degrees of polymerisation. Wang et al [14] reported a model polymerisation degree greater than 10, and the physicochemical properties obtained from simulations were generally consistent with those of natural cellulose. Considering the simulation time and the actual state of cellulose, a cellulose chain with a degree of polymerisation of 20 was first constructed (figure 1(b)), and simulated boxes containing two cellulose chains were constructed using the amorphous cell module in MS software. Periodic boundary conditions were used for the system. The cellulose content in cotton can achieve more than 90%, making it the best material source for natural cellulose and an essential raw material in industrial and production activities. In this study, the moisture content range was determined according to cotton's standard moisture regain rate under atmospheric conditions below 8.5% [15]. Five composite models with different moisture contents (0%, 2%, 4%, 6%, and 8%) were established to explore the effect of moisture content on cellulose. The moisture content of the model was regulated by subjectively varying the number of water molecules. The amorphous cellulose-water model was constructed using the method proposed by Theodorou et al [16]. The initial target density of the model was set to be 1.5 g cm −3 [17].

Simulation details
In this study, the COMPASS force field was used to describe the interaction in the cellulose water model, which has been proven to be suitable for describing the mechanical properties of cellulose systems [14,18,19]. In the COMPASS force field, the energy of the simulated system follows the following equation [20]: Before molecular dynamics simulations, the internal energy of the preliminarily established amorphous molecular model was too high and in a high-energy state. A smart algorithm was used to minimise the model's energy, and the energy convergence accuracy was set to 0.0001 kcal mol −1 . A molecular dynamics annealing simulation was performed to make the simulation system closer to the natural state. The initial temperature was set to 300 K with a step temperature of 100 K. The target temperature was set to 900 K, much higher than the glassy transition temperature of amorphous cellulose. The system runs 10 ps of molecular dynamics equilibrium at each temperature under the NVT ensemble, and the cycle is continued 10 times. During annealing, the molecular system can be fully relaxed, and the configuration of the system is reasonable. The lowest energy configuration after 10 annealing cycles was considered the configuration at the beginning of the kinetic simulation. On this basis, kinetic equilibrium simulations were performed for each model under the NVT and NPT ensembles for 500 ps (temperature, 300 K). During the entire simulation process, the timestep was 1 fs, and the kinetic trajectories of all atoms in the system were saved every 1 ps for data analysis. Temperature and pressure were controlled using the Andersen [21] and Berendsen [22] function methods, respectively. The initial velocities of the individual atoms in the system were determined by the Boltzmann distribution, calculated by van der Waals using an atom-based approach. In addition, the velocity Verlet algorithm was used for integration, and the electrostatic interactions were calculated by the Ewald method.

Elastic modulus and shear modulus
The stress-strain behaviour of solid linear-elastic materials obeys Hooke's law [23]: Where: C ij denotes the Voigt matrix representation of the elastic stiffness tensor. σ i and ε j denote the stress and strain vectors, respectively. For static models in which the elasticity coefficients are calculated at the microscopic level of the material, the stress tensor can be calculated using the virial equation, which is expressed as: Where: m i , v i , r i, and F i denote the mass, velocity, position, and force of the particle i, respectively. V denotes the volume. N denotes the number of particles. For a parallel hexahedron simulation system with side lengths a, b, and c, the initial column vectors were a0, b0, and c0, respectively. When the system is stressed, the vectors become a, b, and c, and the strain tensor is expressed as: Where: h 0 and h denote the matrices composed of column vectors a0, b0, c0, and a, b, and c, respectively.
The elastic coefficient matrix C ij of the system can be obtained according to the above formula. The lame coefficients λ and μ can be calculated as follows:  The elastic modulus and shear modulus can be calculated as follows:

Diffusion coefficients
According to Einstein's diffusion law, the diffusion coefficient can be determined by the mean square displacement (MSD) of the corresponding molecule, which is expressed as: Where: r(0) and r(t) represent the initial and final positions of the particles in the simulation process, respectively, and brackets < > represent the systematic average of the mean square displacement. When the time is sufficiently large, equations (10) and (11) can be abbreviated as follows: Where: k is the slope of the line obtained by least-squares fitting of the curve obtained from the molecular simulation of mean square displacement at time.

Radial distribution function (RDF)
The expression of the radial distribution function is given as the following equation: Where N B is the number of type B atoms, V is the total volume, and ρ is the numerical density. g(r) indicates the local probability density of finding type B atoms at a distance from type A atoms. The coordination number can be calculated using g(r).

Effect of moisture content on the stress-strain behaviour
The stress-strain curves of the systems with different moisture contents were calculated using the Souza-Martins method [24] ( figure 2(a)). Cellulose had the best tensile properties and the highest elastic modulus at a water content of 4%. It is worth mentioning that the effect of moisture content on the cellulose system was not negative every time. When the moisture content was between 0% and 4%, the effect of moisture on the mechanical properties of the system was positive. As the moisture content increased, the tensile properties of the system decreased to different degrees; i.e., when the moisture contents were 6% and 8%. The calculated stress-strain curves alone can only qualitatively analyse the effect of moisture content on the system properties. For further quantitative analysis, the elastic modulus and shear modulus were calculated for the simulated system using the mechanical properties module in the MS ( figure 2(b)). When the water content is 0%, the elastic modulus and shear modulus of the system are 10.96 GPa and 4.34 GPa, respectively. The calculation of elastic modulus and shear modulus for the dry cellulose system was consistent with the results reported in the literature [18,25]. This shows that the simulation method and force field selected for the study are suitable, and the difference is due to the difference between the chosen force field and the size of the model. According to the data in the figure 2(b), elastic modulus and shear modulus of the system first increased and then decreased. Compared to dry cellulose, the systems with 2% and 4% moisture content had 3.6% and 7.6% higher elastic modulus and 8.3% and 9.4% higher shear modulus, respectively. Notably, with a further increase in water content, the system's elastic modulus and shear modulus began to decrease. When the water content was 8%, the elastic modulus and shear modulus of the system were 8.35 GPa and 2.89 GPa, respectively, and decreased by 23.81% and 33.4%, respectively, when compared with dry cellulose, and decreased by 29.17% and 35.15%, respectively, when compared to the 4% water content. The effect of water molecules on the mechanical properties of the simulated system is significant. To further investigate the reasons for the degradation of the mechanical properties of the cellulose system, the mean square displacement and diffusion coefficient of the cellulose chains in the system at different moisture contents were calculated ( figure 3 and table 2). The mean square displacement and diffusion coefficient of cellulose in the model system decreased when the water content increased from 0% to 4%, and the corresponding cellulose chains became more rigid, increasing the elastic modulus and shear modulus of the system. However, the ratio of increase in elastic modulus and shear modulus and the ratio of decrease in the diffusion coefficient is only correlated but not strictly corresponding. This shows that the elastic modulus and shear modulus of cellulose are influenced by the coupling effect of the cellulose chain stiffness and the water-induced cracks and voids in the system rather than by a single factor. The mean square displacement of the cellulose chain first decreases and then increases, which is related to the form of the water molecules in the system. Dan et al [10] proposed that water molecules in an amorphous cellulose simulation system can be divided into three forms according to the number of hydrogen bonds formed: type I and type II water (one or more than two hydrogen bonds), and free water (no hydrogen bonds). Compared to free water, type I and type II water makes the cellulose chains more stable, reduces flexibility, and increases the system's stiffness. When the water content increased from 0% to 4%, type I or type II water molecules adsorb to the cellulose chains and restrict their movement, leading to a decrease in the mean square displacement and diffusion coefficient of cellulose. The mean square displacement of the cellulose chains in the system at 6% and 8% water content was higher than that of dry cellulose. With an increase  in water content, the available hydrophilic sites on the cellulose chains were occupied, and the water molecules were more in the form of free water, which increased the mean square displacement of the cellulose chains to some extent. As shown in figure 3(b), when the water content is 4%, the diffusion coefficient of cellulose is the smallest, and the mechanical properties of the corresponding system are the best. This is due to the low water content, which results in fewer pores and cracks caused by water molecules. In addition, the reduction in the diffusion coefficient leads to the enhancement of cellulose chain rigidity.

Effect of moisture content on the glass transition temperature
The glass transition temperature corresponds to the interconversion of amorphous polymer material between glassy and highly elastic states. The properties of cellulose materials are very stable below the glass transition temperature, similar to glass, and undergo only minimal deformations under external forces. When the external temperature exceeds the glass transition temperature, the properties of the molecules inside the material, such as mobility, freedom, and strength, change. The effect of the water content on the glass transition temperature of cellulose should be studied to ensure the stability of its properties. In this study, the specific volume method was used to measure the glass transition temperature, which has been proven suitable for measuring the glass transition temperature of polymers [26]. The specific simulation process was as follows. The simulated system after kinetic equilibrium was followed by 300 ps of kinetic equilibrium at a quantitative temperature of 650 K to relax the model under the NPT ensemble completely. Then, the relaxed model was cooled from 650 K to 250 K in temperature steps of 50 K, and the kinetic equilibrium was run for 300 ps at each temperature under the NPT ensemble. The calculated structure at the latter temperature was based on the final configuration of the structure calculated at the previous temperature. During data analysis, to prevent errors caused by system fluctuations, the last 150 ps of the data were extracted for analysis. The scatter points in the figure were fitted using the regional linear least-square method. It should be noted that two linear fitting lines with different slopes can be obtained for each model, and the intersection of the two fitting lines is the glass transition temperature of the model.
Based on the simulated data, when the water content was 0%, 2%, 4%, 6%, and 8%, the glass transition temperature values of the simulated system were 424, 432, 496, 422, and 381 K, respectively (figure 4). The changing trend of the glass transition temperature was the same as that of the mechanical properties, which first increased and then decreased. When the water content was 4%, the glass transition temperature of the system reached a maximum of 496 K, which was 17% higher than that of dry cellulose (424 K). A higher glass transition temperature suggests superior mechanical properties of the cellulose material. Figure 5(a) shows that when the water content is 2% and 4%, the mean square displacement of water molecules is about the same. Then, the mean square displacement of water molecules increases with increasing water content. In theory, the mean square displacement of water in the 4% water content model is more significant than that in the 2% water content model. Still, the result obtained was the opposite, which may be related to the aggregation state of water molecules in the model. The water molecules present in the constructed model can be in the form of individuals or together, and the two forms are distinguished by analysing the trajectory of the molecular dynamics process. If the distance between two water molecules exceeds 4 Å, the two water molecules are deemed independent of each other; in contrast, if the distance is less than 4 Å, they are considered clustered together [9]. The equilibrium process of molecular dynamics is dynamic, and it is not easy to quantitatively analyse the state of water molecules. To investigate the influence of the aggregation state on the diffusion coefficient of water molecules, the radial distribution function between oxygen atoms, a characteristic atom in water molecules, was calculated to characterize the distance between water molecules ( figure 5(b)). Wang et al [27] proposed that the distance between atoms corresponds to different interaction types: hydrogen bonding (0.26 nm-0.31 nm), strong van der Waals (0.31 nm-0.50 nm), and weak van der Waals (above 0.50 nm). The peak value at 2.76 Å in figure 5(b) needs particular attention, which is a typical hydrogen bond interaction distance and provides a basis for illustrating the aggregation state of water molecules. When the water content is 2%, the peak at 2.76 Å is the highest. When the content of the water molecules in the system is small (due to hydrogen bonds), it results in a small value for the diffusion coefficient of water molecules, with a value of 4.183 × 10 −11 m 2 s −1 (table 3). The peaks formed at 2.76 Å in the 4% water content system and the 2% water content system are about the same height, which indicates that the possibility of water clusters formed by hydrogen bonds is the same in both systems. The water molecules in the 4% water content system form two peaks between 3.1 Å and 4 Å due to van der Waals forces, while the 2% water content system forms only one peak, and the height of this peak is lower than that of the 4% water content system. This means that the possibility of forming water clusters in the 4% water content system due to van der Waals forces is higher than that of the 2% water content system. The presence of water clusters inhibits the movement of water molecules, which is reflected in the lower diffusion coefficient of water molecules in the system. When the water content of the system is 6% and 8%, although the peak value of 6% at 2.76 Å is lower than that of 8%, it is evident that the peak width of 6% at 2.76 Å is larger than that of 8%, indicating that the molecular water agglomeration formed by hydrogen bonding in the system is more than that of 8%. With increasing water content, the number of free water molecules in the 8% system increased. These two reasons lead to the mean square displacement and diffusion coefficient of the 8% system being higher than those of the 6% water content.

Preferred hydration sites on the cellulose
In addition, the exposed hydroxyl groups O2, O3, and O6 on the cellulose chain easily form hydrogen bonds with water molecules and destroy the original structure of the cellulose chain. The ether groups O4 and O5 can also adsorb water molecules in the system. Figures 6(a)-(d) shows the radial distribution function of the oxygen atoms and water molecules on the cellulose chain. These radial distribution functions illustrate the priority of the reaction order of hydration sites on the cellulose chain under different water contents. As depicted in figures 6(a)-(d), regardless of the water content of the system, there are two obvious peaks: one is the peak formed by the hydroxyl groups O2, O3, and O6 with water molecules in the system under hydrogen bonds, and the other is the peak formed by the acetal groups O4 and O5 with water molecules in the system through van der Waals forces. This peak is very low compared to the peak formed by hydroxyl groups O2, O3, and O6 with water molecules. This indicates that all five hydration sites of the cellulose monomer adsorb water molecules in the system. Still, the adsorption of the hydroxyl groups O2, O3, and O6 with water is dominant over the O4 and O5 acetal groups. As shown in figure 6(a), when the water content is 2%, the number of water molecules in the system is minimal, and the hydroxy O2 adsorption of water is the strongest, followed by O6 and O3. In other words, if there is only one water molecule in the system, it must be adsorbed in the vicinity of hydroxy O2. When the water content was 4%, 6%, and 8%, the number of water molecules increased, and the peaks formed by hydroxyl groups O2, O3, O6, and water had different degrees of increase or decrease. However, the hydration order did not change as O6 > O2 > O3. At 2% and 4% water content, the acetal groups O4 and O5 adsorbed water in the order O4 > O5. As the water content increased, the system was at 6% and 8%, and the acetal groups O4 and O5 adsorbed water in the priority order of O5 > O4.
Usually, several water molecules gather near the glycosidic bonds O4 destroys the oxygen bridge, forms free hydroxyl groups, breaks the cellulose chain, and causes the degradation of cellulose material properties. The coordination numbers of water molecules around O4 in the range of hydrogen bonds (0.26 nm-0.31 nm) were calculated, as shown in figure 7. The number of water molecules coordinated to the glycosidic bond O4  increased approximately 3.26 times with increasing water content. The cellulose properties decreased as the water content increased. It is worth mentioning that in the model with 4% water content, the coordination number of water molecules near the glycosidic bond O4 is the smallest, which again confirms that 4% water content is the crucial point for cellulose to maintain its superior performance.