Monitoring interactions through molecular dynamics simulations: effect of calcium carbonate on the mechanical properties of cellulose composites

This study describes the preparation and characterization of full atomistic models of amorphous cellulose and calcium carbonate (CaCO3) nanocomposite to assess its mechanical properties within and beyond the elastic limit via molecular dynamics simulations. The interactions by hydrogen bond and conformation of the cellulose molecules from the assessment of torsional angles were specifically monitored during the tensile stretching simulations to get deep understanding of the possible structural changes produced in the material during the deformation. On the one hand, the results showed a favorable interaction of the cellulose matrix with the calcium carbonate nanoparticle, with the electrostatic contribution being dominant over the van der Waals component. The determined mechanical elastic constants indicated that the inclusion of the CaCO3 nanoparticle provided an increase on the rigidity of the composite system of 15%, 18% and 19% in the Young, shear or bulk modulus, respectively. On the other hand, using extension and compression simulations, the recovery capacity of the material systems was also assessed in terms of plastic deformation. The elastoplastic behavior was observed for either the neat or the CaCO3 nanocomposite, with an elastic limit around 2.5%. The results also showed that the presence of the CaCO3 nanoparticle produced higher values of plastic deformation in the composite material compared to the neat cellulose system and thus decreased the flexibility of the material. A hysteresis mechanism was identified together with irreversible conformational changes on the cellulose molecules which would explain the plastic deformation observed on the cellulosic systems. It was concluded that the higher plastic deformations observed in the nanocomposite system would be a result of the disruption of the network of hydrogen bonds and the associated decrease on the number of possible interactions.

Young, shear or bulk modulus, respectively. On the other hand, using extension and compression simulations, the recovery capacity of the material systems was also assessed in terms of plastic deformation. The elastoplastic behavior was observed for either the neat or the CaCO 3 nanocomposite, with an elastic limit around 2.5%. The results also showed that the presence of the CaCO 3 nanoparticle produced higher values of plastic deformation in the composite material compared to the neat cellulose system and thus decreased the flexibility of the material. A hysteresis mechanism was identified together with irreversible conformational changes on the cellulose molecules which would explain the plastic deformation observed on the cellulosic systems. It was concluded that the higher plastic deformations observed in the nanocomposite system would be a result of the disruption of the network of hydrogen bonds and the associated decrease on the number of possible interactions.
Abstract This study describes the preparation and characterization of full atomistic models of amorphous cellulose and calcium carbonate (CaCO 3 ) nanocomposite to assess its mechanical properties within and beyond the elastic limit via molecular dynamics simulations. The interactions by hydrogen bond and conformation of the cellulose molecules from the assessment of torsional angles were specifically monitored during the tensile stretching simulations to get deep understanding of the possible structural changes produced in the material during the deformation. On the one hand, the results showed a favorable interaction of the cellulose matrix with the calcium carbonate nanoparticle, with the electrostatic contribution being dominant over the van der Waals component. The determined mechanical elastic constants indicated that the inclusion of the CaCO 3 nanoparticle provided an increase on the rigidity of the composite system of 15%, 18% and 19% in the

Introduction
Cellulose is an almost inexhaustible polymeric raw material with outstanding properties (Klemm et al. 2005) and is remarkably considered to be the most abundant organic compound mainly derived from biomass (Ummartyotin and Manuspiya 2015). With an annual production estimated to be over 7.5 × 10 10 tons (Li et al. 2018), is the most abundant biomass resource in the world. The industrial and commercial applications of cellulose and its derivatives cover a wide range of use cases. From food packaging , humidity and UV sensors (Ummartyotin and Manuspiya 2015), cancer therapy (Jian et al. 2014), wastewater treatment (Kalidhasan et al. 2012), carbon capture (Gebald et al. 2014;Sepahvand et al. 2020) or selective membranes (Hou et al. 2020), among others. Particular attention has been dedicated in the recent years to the development of nanocomposites based on cellulose and inorganic nanoparticles as fillers, such as SiO 2 (Pinto et al. 2008;Sheykhnazari et al. 2016), TiO 2 (Garusinghe et al. 2018;Hamad et al. 2020), gold (Rie and Thielemans 2017), kaolin particles (Fahmy and Mobarak 2008;Torvinen et al. 2017), among others (Alves et al. 2019). The combination of both components provides a final material with improved properties, such as mechanical (Fahmy and Mobarak 2008), antimicrobial (Rie and Thielemans 2017), optical (Hynninen et al. 2021) and thermal or electrical conductivity (Phiri et al. 2018). Among the different inorganic fillers, in the past years great effort has been paid to the synthesis and applications of cellulose nanocomposites filled with calcium carbonate (CaCO 3 ) due to its relating properties such as mechanical strength, biocompatibility, and biodegradation, and bioactivity, and potential applications including biomedical, antibacterial and water pre-treatment fields as functional materials (Ma et al. 2016). Calcium carbonate is one of the most abundant biomineral (Meldrum 2003) as formed in corals, pearls, molluscs and the exoskeleton of arthropods (Vilela et al. 2010). In the field of paper, which is the most important cellulose-based material in terms of worldwide production tons (Bajpai 2013), calcium carbonate is the dominant filler in the production of wood-free uncoated paper (Chen et al. 2011) due to its high number of benefits (economical, environmental, improves paper brightness, opacity, porosity, etc.) (Holik et al. 2013). Nevertheless, the use of calcium carbonate as filler can interfere with the fibre-fibre bonding and thus decrease the tensile strength of the filled paper (Maloney et al. 2005;Fahmy and Mobarak 2008). Some studies have reported, in contrast, the opposite effect of calcium carbonate particles on the mechanical performance of paper (Chen et al. 2011;Mohamadzadeh-Saghavaz et al. 2014). Due to this controversy, further study should be performed to clarify the real relationship between calcium carbonate and cellulose materials.
From an environmental point of view, the pollution problems associated with the use of non-degradable plastic packaging materials are demanding the finding of solutions to replace them with biodegradable alternatives. At this point, the development of eco-friendly cellulose solvents (Isobe et al. 2013;Krysztof et al. 2018) has allowed to the fabrication of various regenerated cellulose (RC) based films and hydrogels, which establish an attractive biodegradable substitute for common packaging materials (Huang and Wang 2022). This kind of cellulose derivatives, obtained from the regeneration of dissolved cellulose, have modulated chemical and physical properties and exhibits an amorphous nature compared to the natural occurring high crystalline cellulose (Ciolacu et al. 2010). Besides, their porous structure allow for the incorporation of functional fillers which can expand their performance and applications (Ardila-Diaz et al. 2020;Salama and Abdel Aziz 2020;Xu et al. 2021;Shaghaleh et al. 2021). Calcium carbonate, which was successfully incorporated in regenerated cellulose films by the precipitation method (Zhu et al. 2020), increased the tensile strength when the optimum concentration was added to the film.
Theoretical modelling of cellulose materials, which comprises ab initio calculations (Kubicki et al. 2013;Uto and Yui 2018), molecular dynamics (MD) simulations (Li et al. 2017;Zhao et al. 2019) and coarse-grained (CG) approaches (Mehandzhiyski and Zozoulenko 2021), allows for a deep understanding of the phenomena occurring at the atomistic or molecular level, which later manifests in their macroscopic properties. Molecular modelling has played an important role over the years for the understanding of molecular-scale phenomena since it offers a level of detail that surpasses what can be reached by experimental methods (Wohlert et al. 2022). Such modelling methods can be used to predict the effect of chemical modifications, such us grafting , the thermal decomposition (Paajanen et al. 2021) or the dissolution process (Bering et al. 2022) of cellulose, among others phenomena, without the need of experimental input. In the specific field of composites, MD simulations are usually used to predict the mechanical properties of the material (Pashmforoush et al. 2020;Ren et al. 2021) or the effect of different additives on its flexibility (Chen et al. 2004a(Chen et al. , b, 2007. Considering the importance of calcium carbonate as filler of cellulose-based or paper related materials and the increasing interest on the packaging industry for the regenerated cellulose based films, the present study describes the modelization performed, through MD simulations, of an amorphous cellulose (AC) composite filled with calcium carbonate nanoparticles, with the aim of studying the effect of such nanoparticles on its mechanical properties, either in the elastic region or beyond the elastic limit in the plastic regime. On the one hand, the interaction energy and radial distribution function (RDF) between the cellulose matrix and the calcium carbonate nanoparticle were determined to assess the type of interaction quantitatively and qualitatively. On the other hand, the elastic stiffness constants were determined using molecular mechanics simulations and, with the aim of assessing the recovery capacity of the materials, the modelled molecular systems were subjected to a series of extension and compression simulations from which the remanent plastic deformation could be determined. The hydrogen bond interactions and conformation of cellulose molecules through the analysis of torsion angles were also monitored during these simulations to check for possible structural changes in the material.

Computational implementation
All calculations presented in the following study were performed using the LAMMPS simulation code (Plimpton 1995). The interactions between cellulose atoms were modeled through a Class II potential as stated next: where K i , H i , V i , and F i are force constants, r is the distance between atoms, b is the bond length, θ and φ are bond and torsion angles, q i is the partial charge of a given atom, and A ij and B ij are the van der Waals parameters. The value of these parameters has been ij defined based on the COMPASS force field, specified in H. Sun's work (Sun 1998). The previous potentials describe well the interaction between atoms in organic or small inorganic molecules, and have been validated for the mechanical simulation of cellulose (Tanaka and Iwata 2006;Bazooyar et al. 2012), The molecular structure of the cellulose monomer (i.e., glucose residue) is represented and the partial atomic charge, in relative terms with respect to an electron charge, is indicated for each atom in the structure but they are not usually accurate for inorganic compounds. For ionic solids the interactions between ions are better described using other exponential potentials, such as those involved in the force field developed by , which was used in the present study to model the interactions between atoms within the calcium carbonate in the present study. The force field described by Raiteri et al. combines the sum of a Buckingham potential and coulombic interactions: where A and ρ are factors associated with the exponential term of the Buckingham potential, and C is a coefficient associated with attractive dispersion forces between oxygen atoms of carbonate ions (set to zero for all other interactions between Ca, C and O). The partial atomic charges for the calcium carbonate atoms were taken directly as defined in the Raiteri  force field. The interactions between cellulose and calcium carbonate atoms were directly modeled from van der Waals interactions using the Lennard-Jones 9-12 potential with the parameters defined in the COMPASS force field, since previous study showed good results with this combination (Ylikantola et al. 2013;Zhao et al. 2019;Pashmforoush et al. 2020). Likewise, the long-distance coulombic interactions (> 10.5 Å) were calculated by coupling an Ewald solver (Particle Mesh Ewald, PME), which uses the reciprocal images of the models. In all the simulations the Nosé-Hoover thermostat (Nosé 1984;Hoover 1985) was used to maintain the systems at 300 K of temperature, the time step was set at 0.5 fs and the Verlet algorithm was used to update the positions of the atoms.

System preparation
Initially, the molecular models of oligomers of cellulose were built by bonding 40 units of β-Dglucopyranose by O-β(1 → 4)-glycosidic bonds, resulting in a chain with a degree of polymerization with n = 40 (i.e., number of glucose residues) (French 2017). The partial atomic charges were generated considering the charge increments of the COMPASS force field (Sun 1998). Figure 1 shows the partial charges of the cellulose structure generated. In this study no ionic groups (such as carboxyl) were included in the molecules of cellulose (i.e. standard cellulose without any degree of oxidation). The model of neat amorphous cellulose matrix was initially generated by randomly adding 38 chains of the cellulose oligomers in the simulation cell under periodic boundary conditions (PBC).
To equilibrate this molecular system, first it was time integrated under the NVT canonical ensemble at 300 K during 75 ps to let the molecules of cellulose to relax and accommodate the volume of the simulation box. Further the system was subjected to 5 cycles of minimization and time integration under the isothermal-isobaric ensemble (NPT) up to a total of 4 ns at 300 K to get a fully relaxed system (i.e., with zero stress). A rendering of this model is shown in Fig. 3a. The calcium carbonate model was based on the crystalline structure of calcite, whose hexagonal unit cell together with the atomic partial charges is shown in Fig. 2. Either the cell parameters, atomic coordinates and the partial atomic charges were taken from the study described by Raiteri et al. (Raiteri et al. 2015). By replicating the unit cell of the calcite in the three directions of space it was obtained a crystal lattice that, after equilibration in the isothermal-isobaric ensemble (NPT) during 1 ns at 300 K, gives rise to the model of pure calcite (Fig. 3b). The density of the calcite model in the equilibrium state was determined as 2.611 g/cm 3 . Compared to the experimental density values determined for amorphous calcium carbonate of 1.62 g/cm 3 (Liu et al. 2010) and that for the calcite polymorph of 2.71 g/cm 3 (DeFoe and Compton 1925), the calculated density value revealed a high crystallinity degree of the developed model. On the other hand, the unit cell parameters calculated for this model were 4.99 Å, 5.16 Å and 16.95 Å for a, b and c respectively. The experimental values determined for calcite were reported as 4.99 Å and 17.062 Å for a and c respectively (Bearchell and Heyes 2002;Raiteri et al. 2015), which were very close to the values determined for the developed model and were thus in agreement with its high crystallinity.
The CaCO 3 nanoparticle model was generated by isolating a spherical region of 2 nm radius in the previously developed pure calcite lattice model. The isolation procedure was made using LAMMPS commands and checking that the stoichiometry of calcium carbonate was not altered after the isolation of atoms. The formula of the isolated nanoparticle was Ca 439 (CO 3 ) 439 . The nanoparticle was later placed on a simulation box and equilibrated under the canonical ensemble (NVT) during 1 ns at 300 K to let the atoms equilibrate under the new geometry. A picture of the stabilized CaCO 3 nanoparticle is shown in Fig. 2b. Finally, the nanocomposite model was built by placing randomly a total of 38 cellulose chains around the equilibrated CaCO 3 nanoparticle within a simulation box with PBC, which supposed a 15.1% in weight of nanofiller reinforcement. This model was relaxed by applying a total of 4 cycles of minimization and time integrated under the isothermal-isobaric ensemble (NPT) at 300 K up to a total of Interaction energy and radial distribution function (RDF) The interaction energy defines the energy of a system (in terms of its potential energy) involved in the interaction between its components, which here would be cellulose and calcium carbonate (CaCO 3 ). This interaction energy (E i ) can be calculated considering the energy of the system and the energy of the isolated components as indicated below: where E total , is the potential energy of the whole system, E cellulose is the energy of the cellulose matrix without the nanoparticle and E CaCO3 is the energy of the calcium carbonate nanoparticle without the cellulose matrix. The potential energy has contributions from bonds, angles, dihedrals, van der Waals and coulombic interactions. The radial distribution function (RDF), also called pair correlation function g(r), was calculated according to next formula: where n ij (r) is the number of pair atoms i-j within the range (r − Δr/2, r + Δr/2), and r the distance between the atoms.

Assessment of hydrogen bond interaction
According to IUPAC (International Union of Pure and Applied Chemistry) (Arunan et al. 2011), hydrogen bonds are formed between an electronegative atom (i.e. the acceptor atom) and a hydrogen atom bonded to a second electronegative atom (i.e. the donor atom) (McNaught and Wilkinson 1997). In the force field chosen in this study (i.e. a combination of COMPASS and Raiteri potentials) the hydrogen bond interaction is considered in the van der Waals and Coulomb interactions of the atoms involved in the bond (i.e. it is considered as a mean field approximation), but not in a specific energy term. For this reason, the contribution of hydrogen bonding to the system in energy terms cannot be measured directly, but otherwise it is possible to have an estimate of the number of hydrogen bonds in the system considering the distance between the acceptor atom and the hydrogen atom. In the present study, the criterion to establish the formation of a hydrogen bond interaction was based in the distance between the involved atoms. As it was established in the work from MacLeod and Rosei (MacLeod and Rosei 2011), the nature of hydrogen bond can be classified according to its energy as very strong, strong and weak (Nishiyama et al. 2008). Besides its energy, the distance between the hydrogen and donor atoms are also affected by the strength of the hydrogen bond and this fact was taken into account in the present study. Two types of hydrogen bonds were assessed from the material systems under study: strong hydrogen bonds were considered when the distance between the acceptor atom and hydrogen was lower than 2.0 Å; on the other hand, weak hydrogen bonds were considered when this distance was in the range 2.0 -3.0 Å. As was previously indicated, the acceptor and donor atoms must be of high electronegativity, which in the considered material systems could only be the oxygen (O) atoms of hydroxyl (-OH) and/or ether groups (-O-) present in the cellulose chains, and the oxygen atoms present in the carbonate ions of the CaCO 3 nanoparticle. In order to quantify the number of hydrogen bonds, the OVITO software (Stukowski 2009) was used considering a maximum cut-off of 2.0 Å or 3.0 Å between the donor and acceptor atoms, depending on the type of hydrogen bond. This software also allows to distinguish bonds between different molecules, and this fact was considered to calculate intramolecular, intermolecular and total number hydrogen bonds.

Mechanical properties
The mechanical properties in the elastic region of the different models prepared were evaluated in terms of the 6 × 6 elastic stiffness matrix, C. For any linear elastic material the stiffness matrix in contracted notation can be written written as: where C ij are the stiffness coefficients relating i stress with j strain components. In materials with symmetries expressed in their principal axis system, the stiffness matrix simplifies and some of the coefficients are equal to zero, being the rest independent (9 for an orthotropic material, 2 for an isotropic one). The calculation of the elastic stiffness matrix by means of molecular mechanics simulations was carried out considering a static or equilibrium approach. After obtaining relaxed structures by applying time integration simulations under the isothermal-isobaric assembly (NPT), the structures were subjected to a set of 6 types of deformation cycles (three uniaxial stress and compression pairs, and three pure shear pairs), controlled by the corresponding strain vectors, with the deformed component taking a small value while the other components of the vector were kept fixed at zero (i.e. changing one cell parameter while maintaining the rest of the cell parameters fixed), and applying a minimization after each deformation stage without changing the cell parameters. The deformation was applied in cycles of positive and negative deformation increments, with amplitude of 0.003 (0.3%) up to a maximum of 0.009 (0.9%). Before the elastic matrix calculation, a pre-minimization procedure was initially applied to the models using the FIRE method (Fast Inertial Relaxation Engine) (Bitzek et al. 2006) for 10,000 iterations. The response of the system at each state of deformation was sampled in terms of the stress during the equilibration time after the initial 6 ps, where the system was allowed to relax, and was defined in the form of virial stress as, where V is the total volume of the system of atoms, m k and u k denote the mass and velocity of the k-th atom, respectively, r kl refers to the distance between the k-th and l-th atoms, and f kl is the force exerted on the l-th atom by the k-th atom. The first term in the virial equation is related to the kinetic energy of the atoms, while the second is related to the potential energy (forces exerted between the atoms). In static conditions, which is the approximation considered here, the first term has no contribution to stress. Each component of the stiffness matrix was obtained from the slope of the stress versus strain plot for each σ ij component of the virial stress: The compliance matrix S was obtained from the inversion of the stiffness matrix C. Assuming an approximately orthotopic behavior, the elastic constants can be directly derived from the compliance matrix component as: On the other hand, the stress-strain curve was generated by applying small consecutive (positive) longitudinal strain increments to the structure under study. These increments were applied in two stages: initially 20 deformation increments equal to 0.003 are applied, and later 25 deformation increments equal to 0.005 were applied, which in total makes a final model deformation equal to 0.185. After each increase in deformation, the model was minimized using the FIRE (Fast Inertial Relaxation Engine) algorithm for 10,000 iterations and the model stress was calculated using the virial. Subsequently, the stress was plotted against the deformation to generate the corresponding stress-strain curve.
Plastic deformation assessment through extension and compression (recovery) simulations Molecular mechanics simulations were used to estimate the plastic deformation remaining in the systems under study after the application of a certain finite strain. The methodology was based on that described by Chen et al. (2004a). In this type of simulations, the relaxed systems are first deformed up to a given final strain value by the application of consecutive small deformation increment steps. The application of a series of small consecutive steps of deformation instead of the application a unique large deformation step allows the system to accommodate after each deformation increment, avoiding thus the generation of unrealistic overstresses that conducted to an unreliable estimation of the properties of the system. The deformation of the system is conducted by changing the box dimensions in the direction under study while keeping the rest of dimensions fixed (e.g. the box dimension X is changed in order to apply a ε 11 strain on each model). The atomic coordinates are remapped after each increment of deformation and further minimized by the application of molecular mechanics using the FIRE (Fast Inertial Relaxation Engine) algorithm (Bitzek et al. 2006) for 10,000 iterations to accommodate the atoms in an energy minimized structure. The systems were deformed up to four different final strains in separated simulations, i.e. 2.5%, 5%, 10% and 15%. Once the systems were deformed up to one of these final target deformations, then they were compressed by application of a series of deformation decrements in the opposite way back to 0% of strain. During the extension simulations, the deformation was applied initially in increments of 0.3% up to a total 4.2% and later in increments of 1.0%, as in the case of the compression simulations. The stress of the system and the number of hydrogen bond interactions were monitored at each stage of deformation to reproduce the stress-strain and hydrogen bonds versus strain curves.

Conformational analysis
The conformation of cellulose molecules was assessed from the analysis of the torsional angles involved in the glycosidic bonds, which determine the relative orientation between two adjacent glucopyranose rings. These torsional angles (i.e. dihedral) are referred as φ (phi) and ψ (psi) (Sankaranarayanan et al. 2015;Nivedha et al. 2016), and were defined as follows:  Figure 4 shows a graphic representation of the definition of φ and ψ angles, indicating the orientation for the calculus of the angle value according to Newman projection representations:

Results and discussion
Interaction between cellulose and calcium carbonate The interaction energy between cellulose and calcium carbonate was calculated from the nanocomposite model as shown in Table 1. The van der Waals and Coulombic components have also been calculated, and the results have been normalized by the area involved. The negative value of the total interaction energy indicated that the cellulose and the calcium carbonate had an attractive interaction, i.e., there is compatibility between both components. On the one hand, it is observed that the Coulombic component is dominant compared to the van der Waals component in the interaction between cellulose and calcium carbonate, which indicated that the electrostatic interactions dominated over the dispersive forces. The radial distribution function (RDF) determined between cellulose and the calcium carbonate nanoparticle (see Fig. 5) showed two specific interactions at 1.55 Å and 2.65 Å, which were due to the electrostatic interaction between the hydroxyl group in the cellulose molecules and the oxygen atom of the carbonate ion, and a third specific interaction at 2.75 Å due to the electrostatic interaction between the calcium ion and the oxygen atoms in cellulose (either in hydroxyl or ether). These results are in agreement with experimental studies described in the literature in which a good absorption and interaction between cellulose and calcium carbonate has been reported (Dalas et al. 2000;Vilela et al. 2010). On the other hand, Fimbel and Siffert (Fimbel and Siffert 1986) determined that the adsorption of carbonate on cellulose was mainly governed by the electrostatic interaction energy, which is in agreement with the results presented here and validate, at least qualitatively, the models developed in the present study. The electrostatic mechanism of interaction between calcium carbonate and cellulose would also explain the selectivity found on the precipitation of calcium carbonate towards the carboxyl group on the surface of cellulose fibers (Vilela et al. 2010), since these chemical groups are partially dissociated in water.

Evaluation of the elastic mechanical properties
The elastic stiffness matrix obtained for the pure cellulose and nanocomposite systems are shown in Fig. 6, whereas a summary of the main elastic constants obtained from the stiffness matrix are collected in Table 2. The values obtained for pure cellulose are consistent with those reported in the  literature (Chen et al. 2004a). The Young's modulus of pure cellulose on the tensile stress conditions are reported in the range 12-27 GPa (Quesada Cabrera et al. 2011), or even lower depending on the type of sample and preparation (5-9 GPa) (Hancock et al. 2000). This modulus is very dependent on the crystallinity of the cellulose, amongst other factors such as relative humidity and porosity (Hancock et al. 2000), i.e. on the relative orientation of the cellulose molecules. According to the contiguity model (Ganster et al. 1994), the volume fraction of crystals V C (i.e. crystallinity) could be calculated as where E is the mean longitudinal modulus of the material, E a is the amorphous modulus, ξ is the contiguity parameter and Ѳ is defined by the expression where E c is the crystalline modulus of cellulose. Taking the values of 128GPa (Sakurada et al. 1962) and 5 GPa (Eichhorn and Young 2001) for the crystalline and amorphous moduli and a value of 5 for the contiguity parameter (Eichhorn and Young 2001), the volume fraction of crystals V c for the pure cellulose model developed in the present study would be of 0.19, i.e. a crystallinity of 19%, considering an average elastic modulus E of 10.3 GPa. In other words, the pure cellulose model developed in this study is essentially amorphous with low degree of orientation or alignment between the molecules, which is the source of crystallinity in cellulose materials. In the same direction, the Poisson's ratio is also very dependent on the crystallinity and phase as stated by Nakamura et al. (Nakamura et al. 2004) who determined values of 0.38, 0.30 and 0.46 for the Poisson's ratio for different cellulose crystalline forms. In the present study, an average value of 0.28 was determined for the Poisson's ratio, which would be in agreement with the reported experimental values considering also the low crystallinity of 19% estimated for this model. The inclusion of the calcium carbonate nanoparticle on the cellulose matrix clearly had an impact on its elastic mechanical properties. Either the Young's, shear or bulk modulus were increased to 11.9 GPa, 4.8 GPa and 9.4 GPa in average, which supposed an increment of 15%, 18% and 19% respectively compared to the neat cellulose. The Poisson's ratio was not changed and the nanocomposite system showed a value of 0.28 as the neat cellulose model. In the papermaking industry, it is usually recognized that the strength properties of paper, whose structure is a network of cellulose fibers, are reduced by increasing the paper ash content (  where RBA m is the relative bonded area between the fibers of the material system. This factor establishes that the strength in a fibrous material is a combination of the area in contact between the fibers and the number of interactions per unit of area in the material. The MD simulations showed that there exists a favorable interaction between the cellulose and the calcium carbonate, so the limiting parameter in the strength factor would be the RBA m (i.e. the area in contact). Precisely, Chen et al. (2011) indicated that the form of the calcium carbonate affected the strength properties of paper depending on how the calcium carbonate particles adapted and embedded between the fibers. In other words, the form of calcium carbonate (e.g. whiskers, spheres) affected thus to the effective area in contact between fibers which would explain, depending on the global balance of the area in contact, why calcium carbonate could improve or deteriorate the strength properties of cellulose-based fibrous materials. For this reason, several experimental studies also found an improvement on the mechanical properties of paper due to the reinforcing of calcium carbonate particles (Chen

Strain recovery assessment through extension and compression simulations
The elastic mechanical properties shown previously were obtained at low deformations of the material. The next type of simulations were obtained by extending the simulations to higher deformations and then compressed back to the initial state. In Fig. 7 the stress-strain curves generated during these extension and compression simulations for the pure cellulose and nanocomposite systems are shown. As it was previously described, the systems were initially deformed (strained) in separated extension simulations up to four different final deformations (2.5%, 5.0%, 10% and 15%) and after that the systems were compressed back to the initial dimensions. During the extension simulations at low levels of deformation (below 2.5%) either the neat cellulose or the nanocomposite systems exhibited a linear stress-strain behavior (i.e., elasticity), which was also observed in the compression curve, since it was very close to the extension curve (i.e., no hysteresis mechanism). When the applied strain in the extension simulation was higher than 2.5% the compression curve did not follow the extension curve and tended to be further away (i.e., not superposed), which indicated that there existed a hysteresis mechanism. This fact also implied that the state at which the stress of the system became zero during the compression simulations had a different length than the original, exhibiting thus a permanent or plastic deformation. In Fig. 8 the results of plastic deformation versus applied strain are plotted. Below applied strains of 5.0% the remanent plastic deformation of both material systems was below 0.5%, which indicated an elastic behavior at this range of deformations. Nevertheless, when the applied strain increased the plastic deformation increased as well over 1.0%, which indicated that the elastic limit had been overcome and the materials experienced an elastoplastic behavior. Besides, the presence of the calcium carbonate nanoparticle produced an additional increase on the remanent plastic deformation of the material after the extension and compression simulations respect to the neat cellulose. In other words, the presence of calcium carbonate nanoparticle Fig. 8 Plastic deformation after extension and compression simulations. In this graph the plastic deformation that remain after the extension and compression simulations of either the neat cellulose or calcium carbonate nanocomposite model is represented for each level of applied strain decreased the flexibility of the cellulosic material system. This type of simulations were also carried out by Chen et al. (2004aChen et al. ( , b, 2007 to assess the flexibility of cellulose models and the effect of different crosslinking moieties. These authors obtained a similar curve for the neat cellulose system either in terms of calculated stress and strain values, and even predicted similar plastic deformation of the material as has been obtained in the present study after the extension and compression simulations using a different force field (Chen et al. 2004a). Together with the stress, the number of hydrogen bonds was also monitored during the extension and compression simulations. As was pointed out previously, two types of hydrogen bonds (i.e. weak and strong) were considered depending on the distance between the acceptor atom and hydrogen. In Fig. 9 a graphic representation of the evolution of hydrogen bonds, considering the total number obtained by the sum of strong and weak interactions, at different strain states during the extension simulations is shown for the neat and nanocomposite systems. In this representation it is possible to observe the formation of failure regions in form of cavities in which there were no hydrogen bonding when increasing the deformation of the material systems, and which were responsible of the decrease of the stress in the system. From Fig. 10, where the total number of hydrogen bonds is summarized, it can be observed that the majority of the hydrogen bonds corresponded to the intramolecular interactions, which agreed with other molecular dynamic simulation studies (Ren et al. 2021). On the other hand, it was observed a decrease of around 3% in the total number of hydrogen bonds, which comes from the sum of either strong and weak bonds, due to the presence of the CaCO 3 nanoparticle in comparison with the neat cellulose, indicating thus that the calcium carbonate was less active towards the hydrogen bond formation than the cellulose itself. This effect of calcium carbonate on the hydrogen bonding of cellulose has been applied experimentally in other studies found on the scientific literature to increase the adsorption capacity of fluff pulp . The evolution of the number of hydrogen bonds during the recovery simulations is plotted in Fig. 11 and Fig. 12 for the neat and nanocomposite systems respectively. In these figures up to six plots can be found, since the assessment of hydrogen bonds was split into strong and weak, regarding to its length distance, and then divided into total, intramolecular and intermolecular attending to the position of the hydrogen bond in the molecules. The tendency observed in the evolution of the number of hydrogen bonds was similar in either the neat cellulose or calcium carbonate composite, except for the fact that the number of bonds in the composite system was below the neat cellulose system as indicated previously. The important points to check in these curves were the tendency on the number of bonds respect to the deformation (i.e., negative or positive), and how the extension and compression curves superpose each other (i.e., presence or not of a hysteresis mechanism). In the case of the strong hydrogen bonds (Fig. 11a and Fig. 12a) it was observed that the intramolecular bonds tended to increase while the deformation increased, while the opposite tendency was observed for the intermolecular bonds. When the systems were compressed back to the relaxed state, a hysteresis effect was observed in the intramolecular bonds, while the evolution of intermolecular bonds was parallel to the original curve, although not totally superposed. The positive tendency and hysteresis mechanism found in the intramolecular bonds indicated that the deformation of the material produced a reorganization of the molecules into a new conformation which was maintained even after the systems were compressed back (i.e., the original geometry was no recovered). This fact was especially important at deformations higher than 5%. The evolution of the total number of bonds was dominated by the intermolecular ones, since they showed a similar tendency. When considering the weak bonds ( Fig. 11b and Fig. 12b), the situation was simpler than that found with the strong bonds. The number of hydrogen bonds, either intramolecular or intermolecular, always decreased while the extension deformation increased, and no hysteresis effect was found later since the compression curve was superposed with the extension curve. Given that the hysteresis mechanism is related to a non-reversible process, the fact that the hysteresis was mainly found only in the strong intramolecular hydrogen bonds indicated in this way that the deformation process would have produced irreversible conformational changes within the cellulose molecules, since these type of hydrogen bonds are mainly located between hydroxyl and/or ether groups located at close positions in the molecule (e.g., O2-H···O3 or O2-H···O5) (Nishiyama et al. 2008). The conformational analysis was assessed trough the phi (φ) and psi (ψ) angles as it is described later in the text.
Comparing the initial and final number of hydrogen bonds it was possible to calculate the recovery of hydrogen bonds as: where N is the number of hydrogen bonds of the system where the stress became zero during the compression simulation and N 0 is the initial number of hydrogen bonds of the system. A value of recovery over 100% indicates that new bonds were created after the deformation process, while a value below 100% reflects that there was produced a loss of hydrogen bonds respect to the original number. A summary of the recovery of hydrogen bonds in terms of the total number of bonds, obtained from the sum of weak and strong interactions, is plotted in Fig. 13. At deformations below 5% the curves of both materials (i.e. neat cellulose and nanocomposite) were statistically similar, but at deformations higher than 5% the nanocomposite system exhibited lower recovery values, coherent with the results of plastic deformation described previously A deeper analysis of hydrogen bonds recovery can be found in the supplementary material (Fig. S1). In summary, the calcium carbonate nanoparticle did not affected significatively the behavior of the hydrogen bonds towards the deformation as has been described previously, so the higher plastic deformation and lower hydrogen bond recovery values observed in the nanocomposite system in comparison with the neat cellulose would be explain due to a disruption of the network of hydrogen bonds and the associated decrease on the number of possible interactions with the calcium carbonate nanoparticle. Fig. 13 Summary of recovery of total (weak + strong) number of hydrogen bonds after extension and compression simulations. In this graph the recovery of hydrogen bonds at each strain is plotted for the neat cellulose an nanocomposite systems

Conformational analysis
The conformation analysis was made in the basis of the values of φ (phi) and ψ (psi) torsional angles, as were defined in Fig. 4. Initially, the phi-psi map (i.e. the Ramachandran plots) were obtained from each of the relaxed structures, as it is shown in Fig. S2 (see supplementary material). A total of 31,121 points make up each plot and each data point represents the combination of φ and ψ angles occurring in a glycosidic bond of the cellulose molecules. The plots were similar in both systems, indicating that the calcium carbonate nanoparticle had no significant effect on the conformation of cellulose molecules in the composite system. In both systems the torsional angles were mainly found at two consecutive oval regions in the ranges of φ/ψ of approximately (− 225, − 120)/ (− 150, − 15) and (− 120, − 15)/(− 150, − 15). These values were similar to those found for cellobiose using quantum mechanical calculations at B3LYP/6-31+G(d) level (French et al. 2021). During the extension and compression simulations, the values of φ and ψ were also monitored at the different deformation states to check for possible conformation changes of the cellulose molecules. These changes were calculated as delta values (ΔѲ i ) respect to the initial values of the φ and ψ angles as stated next, where Ѳ i refers to φ or ψ angles, ε is the deformation, abs indicate the absolute value, Ѳ 0 the angle value at the relaxed state (i.e. ε = 0) and Ѳ ε is the angle value at the ε deformation. These delta values were averaged over all the monitored φ/ψ angles and the average was plotted versus the applied deformation as shown in Fig. 14. As can be observed, during the extension simulation the delta value ΔѲ i increased while the strain increased in both systems and both torsional angles, indicating that there was produced a finite change in the conformation of the cellulose molecules due to this deformation. The variations found in φ angles were higher than those produced in the ψ angles. When the systems were compressed back, the delta value ΔѲ i decreased, i.e. the torsional angles tended to their initial values. Nevertheless, as occurred with the hydrogen bonds, the final values of ΔѲ i (i.e. the recovery of the torsional angles) depended on the applied deformation. The final values of ΔѲ i after the extension and compression simulations have been summarized in Fig. 15. As can be seen, the final ΔѲ i increased with the deformation or, in other words, the recovery of the conformation of the molecules (i.e. their original torsional angle values) was lower when the applied deformation increased, with an exponential tendency for both types of torsional angles. These results agreed with the evolution observed in the strong intramolecular hydrogen bonds (i.e., those mainly located at close positions within a molecule) and explained the hysteresis mechanism associated with this type of bonds at deformations higher than 5%. No significative differences were observed between the neat and nanocomposite systems, which indicated that the calcium carbonate nanoparticle not affected to the conformation of cellulose molecules itself. In brief, the results of the assessment of torsional angles revealed that the tensile deformation of the cellulosic materials, either neat and nanocomposite, produced irreversible conformational changes within the cellulose molecules which contributed to the plastic deformation in the systems.

Concluding remarks
The present study addressed the preparation and characterization of amorphous cellulose composites molecular models to analyze the effect of the presence of a calcium carbonate (CaCO 3 ) nanoparticle on its mechanical properties, either in the elastic region or beyond the elastic limit in the plastic regime, for which the hydrogen bond interactions and torsional angles were monitored to understand in a molecular level the deformation process. It was observed that the interaction of the cellulose matrix with the calcium carbonate nanoparticle was favorable (i.e., attractive), with the electrostatic contribution being dominant Fig. 15 Final values of ΔѲ i for φ and ψ angles after extension and compression simulations. In this graph the final values of ΔѲ i for φ and ψ torsional angles after the deformation process are plotted over the van der Waals component as was obtained either from the analysis of the interaction energy or RDF. The crystallinity of 19% for the neat cellulose model was calculated from the values of the mechanical elastic constants using the contiguity model equations, which was in agreement with the calculated Poisson's ratio value. In addition, the elastic constants determined for the neat calcium carbonate model also agreed with the reported values in the literature. The inclusion of the CaCO 3 nanoparticle gave place to an increase on the rigidity of the composite system, which indicated that the calcium carbonate would potentially contribute to an improvement of the mechanical properties of the cellulose composite. Extension and compression simulations were in the same way conducted to study the effect of tensile deformation on the recovery capacity of the material systems in terms of remanent plastic deformation. At strains below 5% either the neat cellulose or CaCO 3 nanocomposite exhibited an elastic behavior (i.e., lineal stress-strain relationship), with the remanent plastic deformation taking nearly-zero values. Nevertheless, when the applied strain was higher than 5%, an elastoplastic behavior was observed in both material systems, and the plastic deformation increased. The presence of the CaCO 3 nanoparticle produced higher values of plastic deformation in the composite material compared to the neat cellulose system and decreased thus the flexibility of the material. Through the analysis of the evolution of hydrogen bond interactions and torsional angle conformation during the extension and compression simulations, it was possible to identify the presence of a hysteresis mechanism on the strong intramolecular bonds (i.e., short hydrogen bonds located at close positions within the same molecule), and thus the presence of irreversible conformational changes on the cellulose molecules which would explain the plastic deformation observed after the deformation process. The CaCO 3 nanoparticle affected mainly to the number of overall hydrogen bond interactions due to its lower hydrogen bond formation capacity, but not to the evolution of hydrogen bond interactions or to the conformation of cellulose molecules during the deformation process, so the higher plastic deformations observed in the nanocomposite material would be a result of the disruption of the network of hydrogen bonds and the associated decrease on the number of possible interactions.
Author contributions All authors contributed to the study conception and design. Modeling, data collection and analysis were performed by CSE. The first draft of the manuscript was written by CSE and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work has been granted by Spanish Ministry of Science, Innovation and Universities-State Research Agency-PREFIJO project RTC-2017-6258-5.

Declarations
Competing interest The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.