Wood–Moisture Relationships Studied with Molecular Simulations: Methodological Guidelines

: This paper aims at providing a methodological framework for investigating wood polymers using atomistic modeling, namely, molecular dynamics (MD) and grand canonical Monte Carlo (GCMC) simulations. Atomistic simulations are used to mimic water adsorption and desorption in amorphous polymers, make observations on swelling, mechanical softening, and on hysteresis. This hygromechanical behavior, as observed in particular from the breaking and reforming of hydrogen bonds, is related to the behavior of more complex polymeric composites. Wood is a hierarchical material, where the origin of wood-moisture relationships lies at the nanoporous material scale. As water molecules are adsorbed into the hydrophilic matrix in the cell walls, the induced ﬂuid–solid interaction forces result in swelling of these cell walls. The interaction of the composite polymeric material, that is the layer S2 of the wood cell wall, with water is known to rearrange its internal material structure, which makes it moisture sensitive, inﬂuencing its physical properties. In-depth studies of the coupled e ﬀ ects of water sorption on hygric and mechanical properties of di ﬀ erent polymeric components can be performed with atomistic modeling. The paper covers the main components of knowledge and good practice for such simulations.


Introduction and Context
Wood is a hierarchical material, where the configurations at different scales (i.e., lumber, growth ring, cellular and cell wall material) play different roles in adsorption/desorption and also in the resulting swelling/shrinkage and mechanical softening. Different modeling approaches are required at different scales and upscaling between models at different scales is required. However, the origin of wood-moisture relationships lies at the nanoporous material scale that is the sub-cellular scale. Wood polymers have different distribution and configuration in the cell wall primary (P), secondary (S1, S2 and S3) layers, and in the middle lamella ( Figure 1). The work presented here is strongly inspired by the S2 cell wall layer, the thickest layer. The nanoporous material that forms the cell wall is a nano-composite of several polymers that undergoes water sorption and swelling in co-occurrence. Simulating the behavior of water molecules interacting with polymers provides a unique probe in this coupled behavior.
The topic of this paper is the investigations of the moisture-wood relationships at the atomistic scale, using molecular dynamics and grand canonical Monte Carlo simulations. We present the different

Summary of Wood Structure, Components, and Behavior
Wood, an orthotropic cellular biomaterial, has the capacity of adsorbing water molecules from the surrounding environment into its hierarchical material structure. As water molecules attach themselves to the hydrophilic matrix in the cell walls, the induced fluid-solid interaction forces result in a swelling of the cell walls. Moisture-induced internal stresses due to the restrained swelling highly influence the hygromechanical behavior of wood and can eventually lead to cracking, as observed at the macroscale. Moisture can lead to damage and be a cause for biodegradation. Adsorption of moisture in wood, in the hygroscopic range, i.e., until around 30% moisture content mass per mass, results in swelling up to 10% volumetrically and leads to mechanical weakening, e.g., a reduction of the shear stiffness by orders of magnitude.
Wood structural hierarchical levels can be considered as lumber, growth ring, cell, and layered cell wall material. The variation of the cellular structure across the growth ring has been found responsible, at least in part, for the inhomogeneity of swelling and stiffness [1][2][3]. However, from cellular and sub-cellular investigations, it has become clear that a full understanding of the wood polymer-water interactions requires investigating at the scales of the cell wall composite material, which are the focus of this paper and is described next.
The cell wall of softwood is composed of four layers built as the exoskeleton of the cell, see Figure 1a. The primary layer is the first layer which is produced as the cell is growing. At the mature stage, the cell builds its secondary (S) layer in three phases. The thin internal and external cell wall layers (i.e., S3 and S1) act as restraining corsets against deformation due to the winding of the cellulose microfibrils around the cell. In the central and, by far, thickest cell wall layer, namely, S2, the microfibrils are almost parallel to the longitudinal axis of the cells and the presence of a small angle (called microfibril angle) results in a helicoidal organization of the fibrils. The microfibrils, laid out in quite parallel fashion, may undulate and touch each other in the deposition plane (radial direction), but apparently not among planes, as indicated from tomographic transmission electron microscopy [4] or conceptual models as presented in Reference [5]. The network of stiff cellulose microfibrils are embedded within complex matrix polymers [6] and the S2 chemical content has been in large part identified. The microfibrils (MFs) are long, thin filaments of bounded crystalline cellulose, hydrophobic or hydrophilic depending on the surfaces of the crystal, with crystal dimensions roughly 3 × 3 × 30~100 nm 3 . The crystal length depends on species, e.g., 30 nm for Norway spruce cited per Reference [5] or 100 nm for black spruce [8]. The number of chains per crystal, i.e., 18, 24, 30 or 36, is disputed although the often reported 3 × 3 nm 2 cross dimensions would tend to agree with larger numbers of chains [9,10]. The presence of amorphized or water-accessible cellulose is suspected at the crystal surface [9]. To form microfibrils, the crystals are bounded by a hydrophilic amorphous hemicellulose which, in softwood, is galactoglucomannan (GGM) [5]. The microfibrils are embedded in an amorphous matrix which displays two sub-regions: matrix A adjacent, somewhat aligned to the microfibrils, and matrix B, filling in the rest of lenticular spaces in between the microfibrils, as schematized in Figure 1a inset and in Reference [5]. The hydrophilic amorphous matrix is composed of different hemicelluloses (GGM and glucuronoarabinoxylan or xylan) and of two states of guaiacyl lignin, non-condensed and condensed, referring to the types of chemical bonds between monomers, where matrix A would be made of GGM and condensed (cross-polymerized) lignin and matrix B of xylan and uncondensed lignin [11]. The distance between cellulose microfibrils, which interspace is filled by the matrix, is estimated to be in the order of 4-10 nm in dry cell walls [12][13][14]. Overall, as summarized in Figure 1b, the cell wall material is composed of almost equal share of stiffer cellulose microfibrils, of which 50-60% is crystalline, and of a softer polymeric matrix. Of notable complexity, the chemical components of the S2 layer and their relative configuration are continuously elucidated, in order to understand the remarkable properties of wood and provide inspiration for new materials (e.g., [15][16][17]. In the thick S2 cell wall layers, cellulose nanofibrils interact via hydrogen bonds to produce microfibrils, which are embedded in an amorphous matrix of hemicellulose-lignin, thus yielding a tightly interconnected polymer network consisting of numerous non-covalent (hydrogen, hydrophobic, van der Waals) and covalent bonds [5,18]. Lignin and hemicelluloses are covalently linked through lignin-carbohydrate complexes (LCCs) involving benzyl-ether, benzyl-ester, and phenylglycoside bonds [19]. Although in low proportion [20], the covalent linkages between lignin and polysaccharides would likely favor noncovalent interactions between polymers and play a central role in determining the mechanical and physical properties of wood.
Understanding through experiments the effects of the interactions between polymers within the cell wall in dry and moist states and, subsequently, the impact of such interactions on the cell wall properties is still challenging considering that the characterization of the chemical components and the analysis of the physical properties of cell wall requires a multiscale approach ranging from nano-to macroscopic levels [21,22]. As a welcome complement, molecular simulation offers a unique lens into the S2 system. The microfibrils (MFs) are long, thin filaments of bounded crystalline cellulose, hydrophobic or hydrophilic depending on the surfaces of the crystal, with crystal dimensions roughly 3 × 3 × 30~100 nm 3 . The crystal length depends on species, e.g., 30 nm for Norway spruce cited per Reference [5] or 100 nm for black spruce [8]. The number of chains per crystal, i.e., 18, 24, 30 or 36, is disputed although the often reported 3 × 3 nm 2 cross dimensions would tend to agree with larger numbers of chains [9,10]. The presence of amorphized or water-accessible cellulose is suspected at the crystal surface [9]. To form microfibrils, the crystals are bounded by a hydrophilic amorphous hemicellulose which, in softwood, is galactoglucomannan (GGM) [5]. The microfibrils are embedded in an amorphous matrix which displays two sub-regions: matrix A adjacent, somewhat aligned to the microfibrils, and matrix B, filling in the rest of lenticular spaces in between the microfibrils, as schematized in Figure 1a inset and in Reference [5]. The hydrophilic amorphous matrix is composed of different hemicelluloses (GGM and glucuronoarabinoxylan or xylan) and of two states of guaiacyl lignin, non-condensed and condensed, referring to the types of chemical bonds between monomers, where matrix A would be made of GGM and condensed (cross-polymerized) lignin and matrix B of xylan and uncondensed lignin [11]. The distance between cellulose microfibrils, which interspace is filled by the matrix, is estimated to be in the order of 4-10 nm in dry cell walls [12][13][14]. Overall, as summarized in Figure 1b, the cell wall material is composed of almost equal share of stiffer cellulose microfibrils, of which 50%-60% is crystalline, and of a softer polymeric matrix. Of notable complexity, the chemical components of the S2 layer and their relative configuration are continuously elucidated, in order to understand the remarkable properties of wood and provide inspiration for new materials (e.g., [15][16][17]. In the thick S2 cell wall layers, cellulose nanofibrils interact via hydrogen bonds to produce microfibrils, which are embedded in an amorphous matrix of hemicellulose-lignin, thus yielding a tightly interconnected polymer network consisting of numerous non-covalent (hydrogen, hydrophobic, van der Waals) and covalent bonds [5,18]. Lignin and hemicelluloses are covalently linked through lignin-carbohydrate complexes (LCCs) involving benzyl-ether, benzyl-ester, and phenylglycoside bonds [19]. Although in low proportion [20], the covalent linkages between lignin and polysaccharides would likely favor noncovalent interactions between polymers and play a central role in determining the mechanical and physical properties of wood.
Understanding through experiments the effects of the interactions between polymers within the cell wall in dry and moist states and, subsequently, the impact of such interactions on the cell wall properties is still challenging considering that the characterization of the chemical components and the analysis of the physical properties of cell wall requires a multiscale approach ranging from nano-to macroscopic levels [21,22]. As a welcome complement, molecular simulation offers a unique lens into the S2 system.

Using Molecular Simulations for Modeling of Molecular Level Wood-Moisture Interactions
Many moisture-induced physio-chemical processes of porous materials, such as swelling and mechanical weakening, have their roots at the atomistic level. One of the main reasons is that, at various moisture contents, water molecules interfere with the polymer matrix. At atomistic level, in the material, atoms vibrate and interact with each other through various forces, which are adequately modeled in molecular simulations. From these simulations, the collective behavior of atoms can be statistically recorded. Through the theory of statistical physics, the macroscopic properties can be extracted. In the last decades, studies relying on molecular modeling are more and more performed, with a remarkable agreement with experimental results, and have achieved great successes in physical and chemical investigations. We note that, in 2013, the Nobel Prize in Chemistry was awarded to Martin Karplus, Michael Levitt, and Arieh Warshel for their ground-breaking work in molecular simulations.
Modeling of atomistic systems generally involves either of two major methods: quantum mechanics or Newton's classical mechanics. In principle, through quantum mechanics calculation, one could accurately describe any system, including chemical reactions. However, the computational costs are high and system sizes are limited to hundreds of atoms [23,24], a size too small to represent typical polymers. In essence, the processes involved in wood-moisture relationships are physical (sorption and diffusion) and do not involve chemical reactions. Given these considerations, quantum computation is not within our scope of tools for this paper. In contrast, classical mechanics simulations consider the electronic structures of atoms as fixed, meaning atoms are treated as points with fixed properties such as mass and charge, greatly reducing the number of parameters needed and expanding the size of systems to millions of atoms [25].
This atomistic approach allows to select the exact configuration of monomers and their polymerization, plus weak interactions allow to monitor the affinity with water. Molecular simulations offer several advantages. First, such an approach allows full consideration of hydrogen bonds, which play an important role in the physical properties of hydrophilic polymers [26,27]. Second it allows the consideration of dipoles, i.e., those of water molecules (e.g., [28,29]). The third advantage of molecular simulations results from their tempo-spatial resolution capacities to record exact trajectories and velocities of atoms, as well as the forces acting on each atom, at intervals of several femtoseconds. To study bulk materials, when periodic conditions are implemented, properly corrected long-range electrostatic forces can be taken into account [30,31]. We note that, for polymeric systems where hydrogen bonds account for a significant part of the non-bonded interactions between atoms and play a critical role in the hygromechanical behavior, simulations of such systems still need to be fully atomistic. Despite recent developments in coarse-grained models proposing methods to embed properly the dynamics of polymeric materials [32] and in the multiplicity of approaches and application domains [33], in our view, the nature of hydrogen bonds is best handled using atomistic molecular simulations. Lately, cellulose nanocrystals have been the object of several interesting full-atom and coarse grain MD studies, often in conjunction with different non-wood polymers (as reviewed in Reference [34] or in for example Reference [35]). Based on these considerations, all-atom Newton's classic mechanic molecular simulations are retained as the appropriate tool for the investigation of wood-moisture relationships.
Very few atomistic investigations of lignocellulosic systems include interactions with water. Apart from the innovative incursion on cellulose-xylan-lignin-xylan-cellulose systems in Reference [36] on configurations of polymers in presence of water, which included diffusion, such investigations have mainly considered small systems, such as idealized systems of cellulosic crystals lightly cross-bonded with xylan [37], surface wetting of crystalline cellulose [38], cellulose crystals solvated at two degrees of hydration [39], low-DP (degree of polymerization) lignin [40]. In terms of hygromechanical behavior, Reference [41] evaluated with MD the elastic moduli as a function of moisture content for a xylan and lignin complex and their findings were compared with nanoindentation results. Focusing on fuel-related applications, Rismiller et al. [42] used reactive force field (ReaxFF) to investigate the role of water content on the pyrolysis process of lignocellulose biomass. Their results show an enhancement in cellulose breakdown compared to dry systems. The ReaxFF approach was also employed to study the interaction between lignin depolymerization products and metal surface in water-methanol solutions [43]. Vural et al. [44] studied the effect of hydration and temperature on lignin structure and dynamics. Monitoring the temperature effect, they observed a hysteresis behavior in both structure and dynamics of lignin. In addition, Manna and Ghosh [45] studied cellulose chain separation in ionic liquid and in water mixtures and the resultant hydrogen bonding evolution in the decrystallization process. All these studies showed the potential of molecular simulations to analyze particular problems at the atomistic scale, with specific engineering applications, although a consistent approach for wood and its components is still missing.
Our group has been using MD and developing its use, as seen in [27,29,[46][47][48][49][50][51][52][53][54][55], to work on wood and its interaction with moisture. The objective of this paper is to present and summarize the methodological framework developed in our group, to study wood polymer-water interactions at the atomistic scale, and to discuss future paths for this type of work.

Introduction to Molecular Simulation
A detailed presentation of the ensemble of numerical techniques coined as molecular simulation is out of the scope of the present paper but we provide here a short non-exhaustive introduction. In the specific context of wood hydration and deformation, there are essentially two types of molecular simulation techniques that can be used: Monte Carlo and molecular dynamics. Both method families rely on statistical mechanics which can be seen as a microscopic fundamental theory of thermodynamics. This section starts with a few basic notions to guide the reader.
Comparison between the predictions from microscopic or macroscopic theories is often not straightforward as theories are usually derived for ideal situations (pores of simple geometry with ideal pore surface for instance) while experiments are carried out for real, i.e., often complex and heterogeneous, systems. With this respect, molecular simulations which will be briefly described below are very efficient at linking theory and experiment. In particular, molecular simulation allows consideration of an atomistic-level description of the processes to which experimental samples are exposed. While these molecular models cannot capture all the complexity of real samples (heterogeneities, defects, existence of different length/time scales), they are often sufficiently sophisticated to capture well many main features of the behavior observed experimentally.

Chemical Potential
Wood science relies commonly on relative humidity (RH) to document the state of water vapor in equilibrium with wood. While RH can be used for the environment to which wood is exposed, for water sitting within a porous medium, the notion of relative humidity is less convenient. Water potential has been used [56], however, it is physically most sound to use the chemical potential µ defined as the partial molar Gibbs free energy and can be expressed as: where µ, µ 0 , T, k B , and RH are the chemical potential of water, chemical potential of water at saturation, temperature, Boltzmann constant and relative humidity respectively.

Intermolecular Potentials
Properties of wood, confined water, and hydrated wood in molecular simulation studies depend on the interaction potentials used to model the water-water, water-polymer and polymer-polymer interactions. In most simulations, intermolecular energies are assumed to be the sum of pair additive potentials for all interactions. In this respect, it is important to keep in mind that these pair additive interaction potentials are "effective" potentials since they account for N-body terms only in a non-explicit fashion. The fluid-fluid interaction is usually described using potentials (often Lennard-Jones + Coulomb interactions) with parameters chosen to quantitatively capture the bulk properties of water. Similarly, the polymer-polymer interactions, like cellulose-cellulose interactions, are calculated using forcefields that allow describing the main properties of the polymer (cellulose in this example). Finally, the polymer-water interactions are usually derived from simple combination rules such as the Lorentz-Berthelot rules (to combine the like-atom Lennard-Jones parameters to obtain cross-parameters). Water-polymer, in principle, could be also derived from ab initio calculations. However, it is found to be often necessary to adjust these potentials to reproduce experimental data that depend only on the water-polymer interaction such as the isosteric heat of adsorption at very low coverage and Henry's law constant.

Statistical Ensemble
All molecular simulation methods presented in this paper are based on statistical mechanics, for which "ensemble" is one of the most important concepts. A statistical ensemble is a collection of various microstates of an equilibrium macroscopic system as determined by the constraints operating on the system. Depending on the type of constraints, the choice of the ensemble is different. Some of the well-known ensembles include [57]: Micro-Canonical Ensemble (NVE): Constraints on an isolated macroscopic system in equilibrium are constant energy (E), constant volume (V), and constant number of particles (N). This ensemble is used to address an isolated system with no exchange of matter and energy with the reservoir.
Canonical Ensemble for Closed Systems (NVT): Constraints on an isolated macroscopic system in equilibrium are constant temperature (T), constant volume (V), and constant number of particles (N). The system can exchange energy with the reservoir to keep the temperature constant.
Isobaric Isothermal Ensemble (NPT): Constraints on an isolated macroscopic system in equilibrium are constant pressure (P), constant temperature (T), and constant number of particles (N). This is one of the most widely used ensembles as most of the real experiments are carried out under controlled conditions of temperature and pressure [58].
Grand Canonical Ensemble (µVT): Constraints on an isolated macroscopic system in equilibrium are constant volume (V), constant temperature (T), and constant chemical potential (µ). Grand canonical ensemble is a natural choice for adsorption studies at constant volume, where the adsorbed gas is in equilibrium with the gas in the reservoir. The temperature, volume and chemical potential outside and inside the adsorbent are fixed but the number of gas atoms adsorbed is allowed to fluctuate.

Grand Canonical Monte Carlo (GCMC) Simulation
The most widely employed molecular simulation technique in this field regarding sorption related phenomena is grand canonical Monte Carlo (GCMC) simulation [59]. In Monte Carlo molecular simulation, one generates for a given system a set of microscopic configurations. These configurations, coined as microstates, are then accepted or rejected based on the Boltzmann factor (the exact expression for the Boltzmann factor depends on the specific ensemble considered). In general, when treating the problem of water adsorption in wood-related materials, the set of thermodynamic conditions to be considered corresponds to the Grand Canonical ensemble: constant temperature T, volume V, and chemical potential µ. This ensemble, which corresponds to the open ensemble in classical thermodynamics, is appropriate because it allows the number of molecules to fluctuate in the system. Typically, starting from a dry polymer (e.g., cellulose) matrix, adsorption at a given pressure (i.e., relative humidity) is simulated by setting the chemical potential to the corresponding value (at constant temperature, the chemical potential and pressure are related through the Gibbs-Duhem equation ρdµ = dP where ρ is the density). In this grand canonical ensemble, the thermodynamic potential that is minimum at equilibrium is the grand potential: Ω = U − TS − µN = F − µN where N, U, F, and S are the number of molecules, the internal energy, the free energy and the entropy, respectively. In this ensemble, the probability P(N,u) of a microscopic state with N molecules and an intermolecular energy u writes [60]: where Λ is de Broglie's wavelength, β = 1/k B T and Ξ is the configurational part of the grand canonical partition function. There are different Monte Carlo moves involved in GCMC governing adsorption and desorption: particle tranlations, rotations, insertions and deletions. For the insertion of a particle, a particle is inserted into a position randomly inside the simulation box. The probability of accepting the insertion is given by: where u N and u N+1 are the energies of the system before and after the insertion, respectively. The deletion is conducted by removing a randomly selected particle from the N particles in the simulation box and the probability of acceptance is given by: where u N and u N−1 are the energies of the system before and after the deletion, respectively. Many thermodynamic and structural quantities can be estimated by means of GCMC simulations such as the density/adsorption isotherm, isosteric heat, and isothermal compressibility, etc. However, GCMC simulation cannot address the deformation directly as it is performed by definition in an ensemble with constant volume. This limitation is a real obstacle for the use of GCMC in wood related materials as most of them undergo significant swelling upon sorption.

Molecular Dynamics Simulation
Molecular dynamics (MD) is a deterministic method in which one integrates Newton's classical equation of motion for all the atoms k of the system (each atom has a mass m k ): where r k (t) and f k (t) are the position of atom k and the force exerted on this atom at time t. f k (t) is related to the derivative of the interaction potential u with respect to the atom position (the treatment of intermolecular energies and corresponding forces in molecular simulations will be discussed at the end of this section). Such an interaction potential corresponds to the intermolecular energy of this atom with all the other fluid atoms as well as with all atoms of the host porous substrate. In practice, this equation is integrated with a small time step dt (typically, 0.1~1 fs) over a great number of times to reach time scales of the order of 10-100 ns, with system size typically ranging from 10 3 to 10 5 atoms. Molecular Dynamics is quite simple to implement and particularly useful/appreciated as it allows probing both the thermodynamic and dynamical properties. In addition, structural properties are also accessible in MD from a detailed analysis of the atomic configurations obtained along the MD trajectory. While the most convenient ensemble for an MD simulation is the microcanonical ensemble (constant number of molecules, volume and energy), this ensemble is often limited as it does not match the vast majority of experimental conditions. Molecular Dynamics can be modified to consider systems under various more realistic conditions: for example, simulations at constant temperature (canonical ensemble) or simulation at constant temperature and pressure (isobaric-isothermal ensemble). This is done by using algorithms which act as barostat and thermostat by controlling/imposing the convergence towards the target temperature (directly related to the particle velocity distribution) and/or pressure (external force acting on the system through the virial theorem). In practice, an MD simulation is very similar to the course of real experiments. Starting from an initial configuration, the system evolves toward equilibrium where proper statistical mechanics averages can be estimated to obtain relevant structural, thermodynamic and dynamical quantities.
In MD, the adsorption process can be mimicked by random insertion of water molecules into the polymeric system until a target moisture content is attained. Then the corresponding chemical potential can be determined by one-step perturbation methods [61]. At each moisture content level, the system should be probed. However, the major limitation of this approach is that the random insertion of water molecules into the system ignores the energy heterogeneity. In other words, the random insertion strategy is a biased sampling. As a result, this method cannot distinguish the difference between adsorption and desorption and the resulted isotherm are neither adsorption nor desorption isotherms but something lying in between.

Hybrid GCMC/MD Simulation
In the case of wood but also other compliant materials, there are other ensembles and methods that are very relevant. In particular, the osmotic ensemble is of key importance as it allows considering systems at constant temperature T, chemical potential µ, and stress σ. This ensemble can be adapted to look at the coupled effect between wood sorption (hydration) and swelling as it permits to address water sorption while allowing for the deformation of the matrix. This situation can be treated using Monte Carlo simulations in the osmotic ensemble. However, this osmotic ensemble method is challenged by the estimation of free energy differences between the solid phases of interest in the absence of adsorbed molecules, which is still an ongoing research topic [62]. Another option is to use a hybrid strategy that combines Monte Carlo simulations in the Grand Canonical ensemble with molecular dynamics simulations.
The combination of GCMC and MD is realized by simply iterating between the two methods. As shown in Figure 2, part of the trial insertion/deletion used in the GCMC simulations is replaced by an MD trajectory at constant number of water molecules N, constant pressure P and temperature T. Specifically, this means that, in addition to conventional MC steps in the GCMC algorithm (insertions and deletions), MD timesteps, i.e., at constant number of molecules, allowing for water molecule translations or rotations and for local relaxation of molecules at constant external stress, are added as the MD technique is more efficient at converging towards local equilibrium. Such hybrid strategy succeeds in capturing coupled sorption/swelling phenomena, for example, the sorption-induced phase transitions in porous materials such as Metal Organic Frameworks [63].
T. Specifically, this means that, in addition to conventional MC steps in the GCMC algorithm (insertions and deletions), MD timesteps, i.e., at constant number of molecules, allowing for water molecule translations or rotations and for local relaxation of molecules at constant external stress, are added as the MD technique is more efficient at converging towards local equilibrium. Such hybrid strategy succeeds in capturing coupled sorption/swelling phenomena, for example, the sorption-induced phase transitions in porous materials such as Metal Organic Frameworks [63].  Practically, hybrid MD/GCMC molecular simulations consist of performing a large number of blocks where one block corresponds to N GCMC insertion/deletion attempts followed by N MD timesteps. In other words, one "block" essentially corresponds to one iteration between MD and GCMC. A large number of blocks have to be performed to equilibrate the system followed by an adequate number of additional blocks to accumulate statistics. Depending on the system, the ratio between N GCMC and N MD can be adjusted for higher efficiency. For instance, in a dense system with small porosity, it can be expected that the acceptance ratio of trial insertion will be very low. That is to say most trial insertion in GCMC part will be denied and the number of molecules remains unchanged for a lot of blocks. Then it is reasonable to have a relatively high N GCMC /N MD ratio so that more GCMC trial insertion/deletion will be performed before the relaxation of MD steps. For example, for amorphous cellulose, 1000 steps GCMC are performed with 100 steps MD.
The hybrid MD/GCMC is more versatile in addressing sorption-induced deformation and the related hysteresis compared to pure MD. The cost, however, is that hybrid MD/GCMC is generally much more time-consuming. This feature mainly comes from two aspects. The first is, as mentioned, that the acceptance ratio of trial insertion/deletion of GCMC is very low because of the low porosity of wood-related polymers. The second is the poor performance of GCMC in parallel computing. An efficient GCMC in terms of parallel computing is yet to be developed.

Force Fields
All the molecular simulation methods mentioned above require a force field to describe the interactions between atoms. Based on the objectives of the molecular simulation, a number of force fields can be used. Typical choices are GROMOS [30], PCFF [64], CHARMM [65], AMBER [66], OPLS [67], etc. They are similar in terms of the forms of the interaction functions, i.e., expressions of Hamiltonian consisting of potential and kinetic energy parts, however, they differ significantly in their parameterization methodology and consequently the values of parameters. It is difficult to compare the quality of the formulation of different force fields, as they all accurately represent certain aspects of the material behavior however may have different defects. Both GROMOS and PCFF have been comprehensively applied in the investigations of cellulosic polymer, the major component of wood, making them promising choices for studying of wood-moisture relationships. In the following, we focus mainly on the force fields used for modelling cellulose, since most work has been done with respect to this wood component.

GROMOS
The Groningen Molecular Simulation (GROMOS) package and its affiliated force field were developed in the 1980s and are continuously being improved. Many researchers investigated cellulosic materials with this force field and found acceptable agreement with experimental results. Yu et al. [68] used the GROMOS force field to study the stability and solvation of cellulose fragments, where the simulation well represented the conformations of disaccharides. Kulasinski et al. [27] carried out a comparative study of crystalline, paracrystalline, and amorphous celluloses, and achieved an agreement with available experiments in multiple aspects, such as the mechanical modulus and cell parameters of cellulose crystal. Bergenstrahle et al. [69] simulated crystalline cellulose at different temperatures, and found reasonable agreement with experimental data for crystal density, corresponding packing, chain conformation, thermal expansion, etc. Charlier et al. [36] reported a model of a secondary plant cell wall where cellulose, hemicellulose, and lignin were all modeled in GROMOS force field finding a co-orientation of xylan on the cellulose surface that corresponds to experimental observations.

Polymer Consistent Force Field (PCFF)
The PCFF is a force field parameterized against a broad range of experimental observables for organic compounds. It has been applied to modeling cellulose-based materials by many researchers who have found that it captures quantitatively or semi-quantitatively most physical properties. For example, Chen et al. [70] built an amorphous cellulose model and showed that the final density obtained is consistent with experimental values (1.39 g cm −3 versus 1.48 g cm −3 ). As far as mechanical properties are concerned, Tanaka and Iwata [71] used molecular simulation with the same force field to assess Young's modulus of cellulosic crystals, finding values in the range 124-155 GPa that are in good agreement with the experimental value (~138 GPa) [72]. As for adsorption properties, in their molecular simulation of adsorption of various absorbates in cellulose, Da Silva Perez et al. [73] found that the heat of adsorption for a large variety of aromatic compounds is consistent with experimental values; 84% of the adsorbate-cellulose couples displayed differences of less than 20% between the measured and predicted heats of adsorption. Xu and Chen [74] also found that PCFF predicts formaldehyde diffusion in cellulose showing a temperature-dependent self-diffusion coefficient in good agreement with experimental data. Overall, the research presented above shows that the PCFF forcefield can provide a reasonable, at least semi-quantitative, description of cellulose in terms of its density, mechanical, and adsorption properties.

CHARMM
Chemistry at Harvard Macromolecular Mechanics (CHARMM) [75] is one of the widely used softwares and force fields for modeling biomolecules, lipids, and membranes [76] and includes united atom (CHARMM 19), all-atom (CHARMM 22), and polarizable [77] parameterization for a wide range of biomolecules. Reiling and colleagues employed the CHARMM 22 force field to obtain the physical properties of cellulose Iα, cellulose Iβ, and cellulose II. They validated their method by comparing the simulated density and Young's modulus to experimental results [78]. In a comparative study, CHARMM35, GLYCAM06, and GROMOS 45a4 force fields were employed for modeling crystalline Iβ cellulose. The CHARMM force field was successful in modeling the crystalline cellulose configuration and the results were in an acceptable agreement with experimental results [79]. Matthews et al. [80] used CHARMM to model crystals of cellulose Iβ allomorph in aqueous solution and studied the crystal-water interfaces. They reported a good agreement between their cellulose microfibrils configuration and NMR experiments. The studies of xyloglucan adsorption on cellulose microfibril surface [81] and of cellulose interface with oil-in-water emulsion [82] are other examples of cellulose microfibril molecular simulation using CHARMM force field.

Solvers
There are numerous software packages developed incorporating several molecular dynamics integration algorithms, and methods to post-process the simulation results. The LAMMPS [83], GROMACS [84], NAMD [85], Materials Studio, CHARMM [86], and AMBER [87] are some of the well-known software packages developed especially for the simulation of biomaterials. Most of these software products include a set of different force fields such as OPLS, Gromos, CHARMM, PCFF, CVFF, AMBER, and COMPASS. Almost all of these packages are capable of performing multithreading, GPU acceleration, and parallel computation with conventional integration methods such as velocity Verlet and leap-frog, and different methods of long-range calculation, energy minimization, free-energy calculation, temperature, and pressure control. Many studies have compared different MD software products using different test systems, benchmarks, and hardware [88,89]. Comparing GROMACS, LAMMPS, CHARMM, and NAMD with five protein systems from 20,000 to 3 million atoms as test systems found GROMACS to be faster for all of the tested systems sizes while GROMACS, LAMMPS, and NAMD showed similar scaling behavior for larger systems [89]. Wikipedia [90] has a detailed article comparing around 30 different programs applicable to different tasks of molecular mechanics modeling. Due to the excellent flexibility of LAMMPS and GROMACS, we mainly used these two packages for our calculations along with Material Studio and Packmol for building molecules and initial structures. While we find GROMACS faster and more compatible with the available force field parameter data banks such as SwissParam [91] and Automated Topology Builder [92], LAMMPS is more flexible and, since it is able to do hybrid MD and GCMC calculations, it was used for calculation of chemical potential for the determination of adsorption/desorption isotherms.

Single Systems and Sorption
In the last years, we developed a methodology based on atomistic simulations for the determination of the moisture-related behavior of the different components as wells as the composite behavior of the S2 cell wall layer of wood. Using the main chemical components of wood, here softwood, the nanostructures of the polymeric systems are built, their hygromechanical behavior determined and analyzed as well. The MD results are validated with available experimental data. Of particular relevance is the capacity of this method to control exactly, as input, the polymeric constituents, their polymeric configuration and the configuration of layered or composite systems. Molecular dynamics provides a wealth of data, on which statistical physics is applied and then the outcome of these analyses is compared with theoretical and continuum approaches, like thermodynamics and poromechanics, providing a unique window on the fundamental hygromechanical behavior.
In this section, we present the methodologies developed to build amorphous polymeric systems and to expose them to different moisture contents.

Building an Amorphous System
There are different methods for building an amorphous system. Here we present two common methods: a Monte Carlo method and a relaxation method.
We discuss first the Monte Carlo method, which has been implemented in some commercial software like Material Studio 8.0. The degrees of freedom, such as position, orientation and any torsion, of each repeating unit of the polymer are assigned one by one. One has to reserve a number of candidate values for the coming degree of freedom. To determine which value to pick, each candidate value will be subject to a number of geometrical checks, such as whether the associated segment will have close contacts, and/or ring spearing with the formerly placed segments. Then, the Boltzmann weight w i of each candidate that survives the geometrical check is calculated by: where k B is the Boltzmann constant, T is the absolute temperature, and E i , the energy change caused by the trial segment insertion. Then the probability of a candidate P i can be easily determined by: In this way, an amorphous molecular model with periodic boundary conditions can be constructed. One should note that, as more segments are put into the box, there will be less free space left for further insertion. As a result, there will be fewer candidates surviving the geometric check. This makes this method incapable of constructing systems with high density. In practice, one should set a relatively low density for construction, and then achieve the target density by adjusting the box dimension using other molecular simulation techniques. A model of amorphous cellulose built by this method is shown in Figure 3a. Another method of preparing amorphous dry material is through chain relaxation. In the case of, for example, amorphous cellulose, high-temperature amorphization can be used. The starting point is the initial coordinates of the cellulose crystal, which are generated by a freely available script, Cellulose Builder [102], that uses the crystallographic structures as documented by X-ray diffraction. As an example, crystalline cellulose consisting of 32 cellulose chains, 10 cellobiose units long, is subjected to temperature as high as 700 K for several nanoseconds to loosen its ordered structure. The amorphization has two stages [27]. First, the polymer is raised to a temperature above its melting point and equilibrated there in stress-free conditions until the chains lose their ordered structure, typically after 1 ns. Then, the amorphized chains are quenched to room temperature and equilibrated in stress-free conditions. The degree of crystallinity or isotropy can be measured quantitatively through the so-called Herman's orientation function, which is defined as the statistical average of the orientations of the polymer monomers: where is the angle between the axis of interest and monomer orientation. In perfectly ordered cases, i.e., crystalline cellulose, the Herman's orientation function should give values −0.5, −0.5, and 1 along the three principle directions, and for an ideal amorphous polymer, the value should be 0 in all directions showing isotropy.
It is also possible to build the amorphous polymers starting from single chains. For example, in the case of arabinoglucuronoxylan, multiple chains are inserted into a simulation box (shown in Figure 3b with different colors denoting different chains), which are then followed by relaxation Figure 3. Visualization of (a) amorphous cellulose prepared by Monte Carlo method, the box is 3 × 3 × 3 nm 3 and amorphous xylan prepared by molecular dynamic relaxation, (b) before and (c) after densification, the box is 5 × 5 × 5 nm 3 .
Another method of preparing amorphous dry material is through chain relaxation. In the case of, for example, amorphous cellulose, high-temperature amorphization can be used. The starting point is the initial coordinates of the cellulose crystal, which are generated by a freely available script, Cellulose Builder [102], that uses the crystallographic structures as documented by X-ray diffraction. As an example, crystalline cellulose consisting of 32 cellulose chains, 10 cellobiose units long, is subjected to temperature as high as 700 K for several nanoseconds to loosen its ordered structure. The amorphization has two stages [27]. First, the polymer is raised to a temperature above its melting point and equilibrated there in stress-free conditions until the chains lose their ordered structure, typically after 1 ns. Then, the amorphized chains are quenched to room temperature and equilibrated in stress-free conditions. The degree of crystallinity or isotropy can be measured quantitatively through the so-called Herman's orientation function, which is defined as the statistical average of the orientations of the polymer monomers: where θ is the angle between the axis of interest and monomer orientation. In perfectly ordered cases, i.e., crystalline cellulose, the Herman's orientation function should give values −0.5, −0.5, and 1 It is also possible to build the amorphous polymers starting from single chains. For example, in the case of arabinoglucuronoxylan, multiple chains are inserted into a simulation box (shown in Figure 3b with different colors denoting different chains), which are then followed by relaxation run of molecular dynamics at various temperatures ranging from 800 K to room temperature. The entanglement of the different chains can be visually checked (colors are mixed in Figure 3b). Density is an important quality index of the model, where simulation results are compared with experimental data. For cross-polymerized lignin, the branching density can be varied, for example between 0.1 and 0.3 in terms of the number of cross-linkages to total number of monomers. This initial structure is then equilibrated finding an equilibrium configuration. The average number of polymer atoms for such simulation has been around 3000-5000, resulting in a domain of 3 × 3 × 3 nm 3 [49]. We employ periodic boundary conditions (PBCs) in order to mimic bulk material properties.

Water Molecules
There is a number of different water models used in molecular simulations. The GROMOS force field is designed to work with single-point charge (SPC) or SPC/Extended (SPC/E). In terms of the obtained saturated vapor pressure at room temperature of 300 K, SPC (~4400 Pa) is closer to experiment (~3600 Pa) than SPC/E (~1017 Pa) [103]. However, SPC/E outperforms SPC model in terms of density and diffusion constant [104]. There is no consensus on the best water model regarding sorption. Currently both water modules are popular in the sorption community.

Varying the Moisture Content
As mentioned above, the way moisture enters the system differs between MD, GCMC, and hybrid MD/GCMC. With GCMC, a variation of moisture content is achieved by imposing a variation of chemical potential. The impact of chemical potential on moisture content can be path dependent (hysteretic). Thus, adsorption and desorption have to be addressed separately. The adsorption starts from the dry state while the desorption starts from the saturated state, similarly to what is done experimentally. We note that the insertion of moisture by controlling the chemical potential, using grand canonical Monte Carlo (e.g., [54]), is not available within the Gromacs environment.
In MD, moisture adsorption is mimicked by inserting water molecules one after another with random orientation into random locations of the polymer system. The insertions with close contacts or overlaps with existing atoms in the system are rejected, as they are singularities in terms of force field potential functions that may induce unphysical movement or collapse of the MD simulations. Every successful insertion is followed by a relaxation run which further eliminates possibly unphysical close contact. In classical MD, we obtain water adsorption curves using the One-Step Perturbation method [61,105] for the determination of chemical potential. Once the desired number of water molecules has been inserted, and after relaxation runs, the system is interrogated during production runs [29,49]. Calculated sorption curves, i.e., moisture content versus chemical potential, are compared with experimental data for validation.
Rigorously speaking, GCMC is the more proper way to introduce moisture into the system, as it evaluates the energy difference of possible moisture insertions and only the insertions with energy corresponding to the presented chemical potential are accepted. In theory, MD could reach the equivalent status, if the system would attain a sufficiently relaxed state, meaning that water molecules would have eventually diffused to favorable locations, even had they been first forcedly inserted into energetically unfavorable locations.

Probing the Systems Towards Characterization
The material can be probed to determine different hygromechanical properties at different moisture content. Determination of the swelling of the system at each water insertion allows to determine the swelling coefficient, and its dependence on moisture content. In mechanical loading (compressive or tensile, as desired), an anisotropic barostat is applied and the pressure is controlled, until the system finds its equilibrium and the corresponding strain is measured. From these data, we can determine the bulk and shear moduli, and their dependence on moisture content.

How to Determine Hygromechanical Properties: Swelling, Elastic Moduli
Hygromechanical properties such as swelling and elastic moduli can be measured through MD simulations at designed loading schemes.
The swelling strain is related to the equilibrium size of the system, which is measured after equilibration of the simulation box in NPT (p = 0 Pa, T = 300 K) ensemble for a period as long as 20 ns. The strain with respect to the dry state, i.e., the Lagrangian strain V , is: where V , m, x is the volumetric strain, moisture content and size of the system respectively. In the case of volumetric swelling strain, x refers to volume and, for uniaxial swelling strain, x refers to one length direction of the system. For wood-related polymers, uniaxial swelling strains are common as samples are frequently prepared in the form of thin films. An example of the swelling strains of three different components of wood, as well as one composite, is reported in Figure 4.  The mechanical stiffness of material is an important property of wood materials, as is wood frequently used in structural mechanics applications. It is of great interest for both academia and industry to find out the dependence of elastic constants on moisture content. Generally speaking, the most common way to measure moduli is through the stress-strain curves, where orthotropic systems are probed by loading in different directions. Methods along this path differ slightly between each other in terms of the loading protocol, the number of sampling points and the frequency of relaxation, due to the difference in magnitude of stiffness. It is also reasonable to strike a balance between accuracy and computational costs. The least computationally expensive way is to apply one level of stress/strain and collect the response after equilibration. However, to obtain reliable results, it is favorable to obtain stress-strain curves over a certain range.
The bulk, Young's and shear moduli can be measured by the different protocols. For the bulk modulus, we first simulate the system in NPT (p = 0 Pa, T = 300 K) ensemble for 20 ns to obtain the equilibrium volume at stress free state ( , 0). Then we apply an isotropic pressure, which simulates the system in NPT (p = σ, T = 300 K) ensemble for 2 ns to obtain equilibrium volume at loaded state ( , ). Then the bulk modulus is determined as:  The mechanical stiffness of material is an important property of wood materials, as is wood frequently used in structural mechanics applications. It is of great interest for both academia and industry to find out the dependence of elastic constants on moisture content. Generally speaking, the most common way to measure moduli is through the stress-strain curves, where orthotropic systems are probed by loading in different directions. Methods along this path differ slightly between each other in terms of the loading protocol, the number of sampling points and the frequency of relaxation, due to the difference in magnitude of stiffness. It is also reasonable to strike a balance between accuracy and computational costs. The least computationally expensive way is to apply one level of stress/strain and collect the response after equilibration. However, to obtain reliable results, it is favorable to obtain stress-strain curves over a certain range.
The bulk, Young's and shear moduli can be measured by the different protocols. For the bulk modulus, we first simulate the system in NPT (p = 0 Pa, T = 300 K) ensemble for 20 ns to obtain the equilibrium volume at stress free state V(m, 0). Then we apply an isotropic pressure, which simulates the system in NPT (p = σ, T = 300 K) ensemble for 2 ns to obtain equilibrium volume at loaded state V(m, σ). Then the bulk modulus K is determined as: It is worthy to note that the load value should be chosen with care. It should be large enough to minimize the fluctuation error as well as simulation time while it should be relatively small to ensure the system remains in the elastic regime. This rule also applies to Young's and shear modulus.
Young's modulus is measured by applying a directional stress. For an isotropic material, its value is independent of direction. We first simulate the system in NPT (p = 0 Pa, T = 300 K) ensemble for 20 ns to obtain the equilibrium size at stress free state X(m, 0). Then we apply an anisotropic pressure, which simulates the system in NPT (p = σ, T = 300 K) ensemble for about 100 ns to obtain equilibrium volume at loaded state X(m, σ). The Young's modulus E is determined as: Shear modulus is commonly measured through the biaxial method, which means tensioning in one direction meanwhile compressing in the two other orthogonal directions, where the absolute values of stress are σ i = 2σ j = 2σ k . The shear modulus G i is then obtained as: The shear modulus should be independent with direction for isotropic materials. As previously mentioned the magnitude of the load should be carefully chosen.
A more rigorous, however, more expensive, way to measure the mechanical stiffness is first to obtain the stress-strain curve in strain-loading, and then determine the elastic constant as the slope of the linear regime of stress-strain curve in a range of small strain. For bulk, Young's and shear modulus, volumetric, uniaxial, and shear strain are applied separately, both in tension and compression. The strain rate is chosen as low as possible considering the limitation of computational costs, which is 3.66 × 10 −5 nm/ps, 5.00 × 10 −5 nm/ps and 3.66 × 10 −5 nm/ps, for bulk, Young's and shear strain, respectively. This method is usually referred to as dynamic loading [107]. To eliminate the influence of straining rate, the strain can be applied in a stepwise manner and, after every step, the strain is maintained for a certain time, while the system is intermittently relaxed. This protocol makes the results less rate-independent, however, more computationally expensive. A typical stress-strain curve is given in Figure 5. The slope of the tangent line to the stress-strain curve over a small strain range gives the elastic modulus (simulation details can be found in Reference [47]).
Poisson's ratio ν and the moduli are related to each other for isotropic homogeneous materials: Through this relationship, the consistency of moduli measurement can be checked when isotropy is assumed. after every step, the strain is maintained for a certain time, while the system is intermittently relaxed. This protocol makes the results less rate-independent, however, more computationally expensive. A typical stress-strain curve is given in Figure 5. The slope of the tangent line to the stress-strain curve over a small strain range gives the elastic modulus (simulation details can be found in Reference [47]). Poisson's ratio and the moduli are related to each other for isotropic homogeneous materials: Through this relationship, the consistency of moduli measurement can be checked when isotropy is assumed.
Hydrogen bonds play an important role in the mechanical and thermodynamic properties of polymers. For example, it is shown that the number of hydrogen bonds correlates well with the Hydrogen bonds play an important role in the mechanical and thermodynamic properties of polymers. For example, it is shown that the number of hydrogen bonds correlates well with the mechanical stiffness [29]. In MD, the establishment of a hydrogen bond is judged by relative locations of the donor-hydrogen-acceptor triplet. The criteria are that the distance r follows r ≤ r HB = 0.35 nm and the angle α follows α ≤ α HB = 30 • where r is the distance between the donor oxygen atom and the acceptor oxygen atom, and α is the angle of acceptor oxygen atom-donor oxygen atom-donor hydrogen atom.
The material properties obtained from atomistic simulation, i.e., moisture adsorption isotherm, swelling and moisture-dependency of mechanical properties, of the different systems can be further analyzed to obtain a better insight into the hygromechanical behavior of polymers. Analysis of the nanoscopic behavior is possible given the unique level of information provided by molecular dynamics simulation. Analyses can be based on the configuration of the polymer and the distribution of water, as well as the lifetime, location, and number of hydrogen bonds in the polymeric systems, etc. For example, the analysis of the bulk and shear moduli versus the number of polymer-polymer hydrogen bonds yields information on hydrogen-bond breaking due to the water molecules' adsorption. This analysis of breakage of hydrogen bonds with increasing moisture content gives information on the mechanical weakening process in polymers. Systems properties such as monomer composition, the substitution of hemicellulose backbone, proportion, and type of linkages between monolignols and molar mass obtained from experiments can be further used to validate and analyze the MD results.

MD/GCMC Adsorption and Desorption
Although MD simulation, as presented in the last section, can show the strong coupling between sorption and deformation in wood at the molecular level at the nanoscale, the method does not allow to study hysteresis. If one wants to further include the sorption hysteresis, which is quite common in wood research, the hybrid MD/GCMC method should be applied.
In this section, we present water adsorption and desorption in amorphous cellulose as an example to illustrate the detailed method and what can be expected from this method. The hybrid molecular simulations are performed at T = 300 K using a Berendsen thermostat and an anisotropic external stress = 0 Pa (with the following relaxation times: τ T = 0.1 ps and τ s = 1.0 ps). Water molecules are described using the SPC/E water model with the SHAKE algorithm to maintain its internal structure rigid. The MD trajectory is integrated using the velocity Verlet integrator with a timestep equal to 1 fs. Hybrid MD/GCMC molecular simulations consist of performing a large number of blocks where one block corresponds to 2000 GCMC insertion/deletion attempts followed by 200 MD timesteps. In total, 1 × 10 5 blocks are first performed to equilibrate the system followed by 2 × 10 4 additional blocks to accumulate statistics. In order to check the efficiency of the phase space sampling, the pressure and volume are monitored along the hybrid GCMC/ NPT simulations. Figure 6 shows the evolution of the volume of the simulation box as a function of the number of time steps in the molecular dynamics simulation N MD . The insert in Figure 6 shows for two different chemical potentials imposed, how the volume fluctuates around its equilibrium value in the course of a block consisting of an MD segment and a GCMC segment (of course, the volume does not change during the GCMC segment as only the energy and number of molecules are allowed to change). Overall, the data in Figure 6 shows that volume sampling in such hybrid simulations is efficient and that the volume reaches its final equilibrium value provided simulations are long enough (typically, at high loadings, at least 10 7 MD timesteps are required but shorter simulations are needed for low loadings). Both water adsorption and desorption are performed on the molecular model of amorphous cellulose. For the adsorption branch, the dry amorphous cellulose samples are selected as the initial configuration and molecular simulations at different chemical potentials are conducted to obtain the adsorption amount (moisture content, m) as a function of water relative humidity RH. Then, starting from the state close to RH = 1, desorption is simulated by decreasing the water relative humidity (which is related to the water chemical potential). Figure 7a gives a comparison of the simulated water adsorption/desorption isotherm in cellulose with its experimental counterpart showing a good agreement. In Figure 7b, snapshots of molecular configurations are given for different moisture contents. Figure 6. Volume in nm 3 as a function of MD steps as obtained from hybrid GCMC/MD simulations. The data obtained for a relative humidity RH = 0.025 (black data) and RH = 1 (red data) are shown. The insert, which shows a zoom corresponding to the end of the hybrid simulation, illustrates how the volume fluctuates around its equilibrium value in the course of a block consisting of a MD segment and a GCMC segment (the volume does not change during the GCMC segment as only the energy and number of molecules are allowed to change).
Both water adsorption and desorption are performed on the molecular model of amorphous cellulose. For the adsorption branch, the dry amorphous cellulose samples are selected as the initial configuration and molecular simulations at different chemical potentials are conducted to obtain the adsorption amount (moisture content, m) as a function of water relative humidity RH. Then, starting from the state close to RH = 1, desorption is simulated by decreasing the water relative humidity (which is related to the water chemical potential). Figure 7a gives a comparison of the simulated water adsorption/desorption isotherm in cellulose with its experimental counterpart showing a good agreement. In Figure 7b, snapshots of molecular configurations are given for different moisture contents. We note this section by mentioning that atomistic investigation can also aim at documenting other aspects of material behavior and material properties such as diffusion coefficient, thermal dilation, heat capacity, friction, etc. as done in Kulasinski et al. [27,[49][50][51] and Zhang et al. [46,47,55]. A full description of all possibilities is not the aim of this paper.

Upscaling
Molecular simulations are powerful to characterize material properties and reveal the mechanism at the molecular level. However, the high resources demanded and the time-consuming nature of molecular simulation still impose limitations on both time and length scales. Therefore, if one is interested in the hygromechanical behavior at a higher scale, upscaling is generally needed. We note this section by mentioning that atomistic investigation can also aim at documenting other aspects of material behavior and material properties such as diffusion coefficient, thermal dilation, heat capacity, friction, etc. as done in Kulasinski et al. [27,[49][50][51] and Zhang et al. [46,47,55]. A full description of all possibilities is not the aim of this paper.

Upscaling
Molecular simulations are powerful to characterize material properties and reveal the mechanism at the molecular level. However, the high resources demanded and the time-consuming nature of molecular simulation still impose limitations on both time and length scales. Therefore, if one is interested in the hygromechanical behavior at a higher scale, upscaling is generally needed. All the results obtained from atomistic simulation, such as moisture sorption isotherms, swelling and moisture-dependency of mechanical properties for all the different components and composite systems form an integrated dataset which can be upscaled using a continuum approach in order to provide material properties. In this approach, the swelling coefficients and mechanical properties, and their dependence on moisture content are fed as data into classical continuum constitutive laws and the governing conservation laws can be solved using a finite-element approach and solver.
When upscaling, several special considerations have to be made for wood-related polymer. First, as polymers are micro-porous materials with pore widths smaller than 2 nm, one should generally employ an extended poromechanical approach appropriate to model microscopic materials. Traditional poromechanics addresses the coupled behavior of sorption and deformation of meso-and macroporous materials (pore size larger than 2 nm) showing multilayer adsorption (film forming) and capillary condensation. Extended poromechanics is built on the fact when the pores become filled, there is no well-defined fluid bulk phase and solid-fluid interface in micropores [109,110].
For microporous media, sorption is not described by a surface covered by a film, followed by capillary condensation, but rather by pore filling. The state of the adsorbed fluid is then fully characterized by the number of molecules adsorbed and its chemical potential. Following this line, Brochard et. al. [111] proposed a reformulation of poromechanics to account for sorption-induced stress in microporous solids in function of the number of molecules adsorbed.
The second aspect is the non-linear behavior of wood-related polymers during sorption and deformation. A clear example of this nonlinear behavior is the mechanical weakening of polymers with increasing moisture content. Thus a nonlinear poromechanics is generally needed to do a proper upscaling for wood-related polymers. The dependence of mechanical properties on sorption can be included in poromechanics by considering a higher order formulation of the free energy allowing to formulate nonlinear poroelastic constitutive models [112]. This approach can model the dependence of elastic modulus, moisture capacity, and swelling coefficient on stress and liquid pressure. As such mechanosorptive effects can be taken into account for cellulose where atomistic simulations were used to determine the material parameters of the non-linear mechanical model [53].
Third, in classical continuum mechanical models reversible adsorption without hysteresis is generally assumed, while sorption processes in wood show hysteresis, which can have an important effect on the hygromechanical behavior of wood [113]. In nanoporous materials like wood polymers, hysteresis is found to originate from the coupling between sorption and swelling [54]. Due to swelling, new sorption sites become accessible during adsorption, increasing substantially the possible sorption amount in the material. During desorption, these sites remain filled by water molecules leading to sorption hysteresis. Also hysteresis has been observed in swelling strain versus relative humidity. However, this hysteresis disappears when plotting swelling strain versus moisture content. Since a correct macroscopic hysteresis modeling on the macroscale has to be based on a coupling between sorption and swelling. The Independent Domain Model, which neglects the mechanical interaction between neighboring pores, is not adequate, although it has been used for wood [114] or another nanoporous materials [115]. Alternatively, the Dependent Domain Model (DDM) includes the mechanical influence between neighboring pores when becoming filled and swollen and is able to give a quantitative description of the coupling of sorption and deformation and related hysteresis. Briefly, DDM decomposes the porous media into a collection of porous elements following a certain pore size distribution. Deformation and sorption are connected under the framework of poromechanics. The breakage and reformation of polymer-to-polymer hydrogen bonds in the wood-related polymers can be represented by open and closed states of the pores in the continuum model. Based on the coupling between sorption and deformation, the element can switch its own states upon the influence of its neighbors, which finally leads to hysteresis in both sorption and deformation. One can refer to References [116,117] for further details. Here we show briefly an upscaling model applied to amorphous cellulose. The poromechanics extended to microporosity made in Reference [118] was employed, together with the non-linear material behavior. The DDM was adopted to account for the hysteresis. We implemented all the elements stated above into a Finite Element method [116] and built a continuum model based on molecular simulations which captures the non-linear coupling between sorption and deformation with hysteresis, with good agreement between atomistic and continuum modeling results as shown in Figure 8. Here we show briefly an upscaling model applied to amorphous cellulose. The poromechanics extended to microporosity made in Reference [118] was employed, together with the non-linear material behavior. The DDM was adopted to account for the hysteresis. We implemented all the elements stated above into a Finite Element method [116] and built a continuum model based on molecular simulations which captures the non-linear coupling between sorption and deformation with hysteresis, with good agreement between atomistic and continuum modeling results as shown in Figure 8.

Highlights of Current Molecular Simulations Research in the Group
With the aim of investigating the moisture impact in wood, the research activities in our group combine multiscale simulation and experimental methods, where molecular simulation plays an important role in extending the resolution down to the nanometer level, explaining molecular mechanisms of moisture-induced effects on wood and providing the necessary hygromechanical and poromechanical information for the upscaling into continuum models. The softwood cell wall S2 layer, the thickest wood cell wall layer, is the major modeling object. The atomistic models of individual components, including crystalline and amorphous cellulose, glalactoglucomannan, arabinoglucuronoxylan, condensed and uncondensed lignin, are built, validated and investigated. Our main findings to date are listed below.
The sorption isotherms of wood-related biopolymers are captured by molecular simulations, with hysteresis covering the entire RH range [54]. Volumetric strain stays linear with the moisture content in both adsorption and desorption branches and no hysteresis is present between them, which is explained by the fact that the swelling of biopolymers is mainly driven by the additional inter-chain space created by water molecules. The microscopic picture of different hydrogen bond network upon adsorption and desorption is found, which paves the foundation of understanding the hysteresis related to sorption and deformation. The analysis of hydrogen bonds shows that, at the same moisture content, the system accommodates the same number of water molecules but distributed according to different microscopic hydrogen bonding configurations corresponding to distinct host pore distributions. This complex coupling leads to the hysteresis observed in experimental and simulated sorption isotherms. It also leads to the hysteresis in mechanical properties such as bulk modulus against the moisture content. Further studies show that the

Highlights of Current Molecular Simulations Research in the Group
With the aim of investigating the moisture impact in wood, the research activities in our group combine multiscale simulation and experimental methods, where molecular simulation plays an important role in extending the resolution down to the nanometer level, explaining molecular mechanisms of moisture-induced effects on wood and providing the necessary hygromechanical and poromechanical information for the upscaling into continuum models. The softwood cell wall S2 layer, the thickest wood cell wall layer, is the major modeling object. The atomistic models of individual components, including crystalline and amorphous cellulose, glalactoglucomannan, arabinoglucuronoxylan, condensed and uncondensed lignin, are built, validated and investigated. Our main findings to date are listed below.
The sorption isotherms of wood-related biopolymers are captured by molecular simulations, with hysteresis covering the entire RH range [54]. Volumetric strain stays linear with the moisture content in both adsorption and desorption branches and no hysteresis is present between them, which is explained by the fact that the swelling of biopolymers is mainly driven by the additional inter-chain space created by water molecules. The microscopic picture of different hydrogen bond network upon adsorption and desorption is found, which paves the foundation of understanding the hysteresis related to sorption and deformation. The analysis of hydrogen bonds shows that, at the same moisture content, the system accommodates the same number of water molecules but distributed according to different microscopic hydrogen bonding configurations corresponding to distinct host pore distributions. This complex coupling leads to the hysteresis observed in experimental and simulated sorption isotherms. It also leads to the hysteresis in mechanical properties such as bulk modulus against the moisture content. Further studies show that the external stress and temperature may influence the coupled behavior, for instance, a compressive external stress or higher temperature may lead to less sorption (referred to a mechano-sorption) and less deformation and limited hysteresis.
Moisture-induced crossover behavior is found in a number of thermodynamic and mechanical properties of arabinoglucuronoxylan. The mechanism is explained by the formation of a double-layer adsorption film, providing a long-missing part of the general picture of polymer-moisture relationship at molecular level [47]. Through molecular simulation, weakening and swelling of an uncondensed lignin are revealed to be correlated with the reciprocal of the segmental thermal oscillations, coined as local stiffness. The phenomenological similarities between the impact of hydration and heating are examined and supported at molecular level, with the fundamental difference between heat and moisture elucidated [55]. Wood as a composite material consists of crystalline (cellulose crystal) and amorphous components (other polymers). Interesting stick-slip behavior at their interfaces is found and the microscopic details are presented. Moisture adsorption plays a significant role in weakening such interfacial frictional behavior, thus providing a possible molecular switch mechanism which might be related to recovery of wood cell wall after deformation [46].
Finally, such atomistic simulation methods can be used as part of more applied studies, for example in archeology. Experimental results show a substantial decrease in mechanical stability of water logged archaeological wood resulting in drastic distortion and collapse during drying process. As a consolidation technique used for both Swedish warship Vasa [119] and Henry VIII's warship the Mary Rose [120], polyethylene glycol (PEG) hydrogel was sprayed for decades at the surface of both shipwrecks to penetrate the wood and reinforce the wood structure [119]. Measurements show the relative recovery of mechanical strength and the elevated resistance to drying stresses. The cell walls of treated archeological wood are polymeric systems of slightly higher complexity that new wood cell walls. However, molecular dynamics simulation can be used to elucidate the interaction between PEG molecules and both new and degraded wood biopolymers, assuming the system is stable chemically. Molecular models of different components of degraded wood can be built according to available chemical characterization measurements. Methods presented in this paper can yield the microscopic mechanism by which the PEG mixture stabilizes mechanically the wood structure and decelerates adsorption-induced mechanical softening and such work is on-going. In addition to the insight towards consolidation techniques, these studies can provide a better understanding of the physics of biopolymer mixtures and their response to hydration.
For upscaling, we build a continuum model based on DDM and poromechanics to model the coupling of sorption and deformation in biopolymers. The breakage and reformation of polymer-to-polymer hydrogen bonds in the molecular systems are represented by open and closed states of the pores in the continuum model. The DDM is found to agree well with the molecular simulations and experimental results and to give a good description of sorption hysteresis. Further studies show that the sorption isotherms are greatly influenced by material properties. It was found that a small mechanical modulus of the polymer matrix and a strong adsorbent-adsorbate interaction leads to more sorption due to the higher deformation of the material.

Reflections for Future Directions
Molecular simulation brings unique perspectives on wood-moisture relations and points out future directions for the wood science community. In essence, molecular simulations extend the resolution of investigations down to the nanoscale which complements traditional experimental methods, where information at nanoscale is difficult to obtain. The molecular simulation techniques introduced in this paper, which are developed based on a number of theories, including sorption, thermodynamics and statistical mechanics, can in theory be applied to any nanoscale topic that does not involve chemical reactions, i.e., forming and breakage of covalent bonds. Many interesting physical processes, such as moisture-induced shape memory effect and cell-wall recovery after irreversible deformation, where molecular interactions play important roles in the wood-moisture relationship remain to be elucidated.
The ongoing research of deciphering the detailed chemical structure as well as material arrangements of wood is majorly driven by investigations using cutting-edge experimental instruments. However, molecular modeling can always contribute to further support the experimental results which provide, most frequently, indirect structural information. Molecular modeling offers us the freedom of investigating unlimited possibilities of different configurations and composites of wood structures, ranging from individual wood polymer materials to composite structure resembling subunits of wood cell wall, most of which still remain to be modeled and scrutinized, bequeathing a giant pool of potential modeling targets.
The mechanisms of common physical processes in wood-moisture relations, i.e., swelling, weakening, and hysteresis, are the core topics of our current research. However, the impact of better understanding natural biopolymers sorptive behavior is not limited to wood-related topics. Wood polymers, in a broader sense, can be seen as soft material with nanoporosity, for which the sorption and swelling differ greatly from macroscopic porous structures or rigid structures. The molecular interactions at the fluid-solid surface dominate the coupled deformation and sorption behavior. Integrating this knowledge is a challenge to the current theoretical frameworks, e.g., poromechanics theories, which provide no explanation of sorption-swelling coupling. The insights of these molecular studies will help to establish a proper description of other soft nanoporous materials, such as gas shale and artificial water nanochannels, to understand the mechanisms of moisture-induced phenomena related to them, and to further improve theoretical models.