Thermal Response in Crystalline Iβ Cellulose: A Molecular Dynamics

Jul 12, 2007 - ... France, affiliated with Université Joseph Fourier, and member of the Institut de .... Malin Bergenstråhle-Wohlert , John W. Brady...
0 downloads 0 Views 462KB Size
Download PDF
9138

J. Phys. Chem. B 2007, 111, 9138-9145

Thermal Response in Crystalline Iβ Cellulose: A Molecular Dynamics Study Malin Bergenstråhle,† Lars A. Berglund,† and Karim Mazeau*,‡ KTH Department of Fibre and Polymer Technology, 10044 Stockholm, Sweden, and Centre de Recherches sur les Macromole´ cules Ve´ ge´ tales (CERMAV-CNRS), BP53, 38041 Grenoble cedex 9, France, affiliated with UniVersite´ Joseph Fourier, and member of the Institut de Chimie Mole´ culaire de Grenoble ReceiVed: March 21, 2007; In Final Form: May 25, 2007

The influence of temperature on structure and properties of the cellulose Iβ crystal was studied by molecular dynamics simulations with the GROMOS 45a4 force-field. At 300 K, the modeled crystal agreed reasonably with several sets of experimental data, including crystal density, corresponding packing and crystal unit cell dimensions, chain conformation parameters, hydrogen bonds, Young’s modulus, and thermal expansion coefficient at room temperature. At high-temperature (500 K), the cellulose chains remained in sheets, despite differences in the fine details compared to the room-temperature structure. The density decreased while the a and b cell parameters expanded by 7.4% and 6%, respectively, and the c parameter (chain axis) slightly contracted by 0.5%. Cell angles R and β divided into two populations. The hydroxymethyl groups mainly adopted the gt orientation, and the hydrogen-bonding pattern thereby changed. One intrachain hydrogen bond, O2′H2′‚‚‚O6, disappeared and consequently the Young’s modulus decreased by 25%. A transition pathway between the low- and high-temperature structures has been proposed, with an initial step being an increased intersheet separation, which allowed every second cellulose chain to rotate around its helix axis by about 30°. Second, all hydroxymethyl groups changed their orientations, from tg to gg (rotated chains) and from tg to gt (non-rotated chains). When temperature was further increased, the rotated chains returned to their original orientation and their hydroxymethyl groups again changed their conformation, from gg to gt. A transition temperature of about 450 K was suggested; however, the transition seems to be more gradual than sudden. The simulated data on temperature-induced changes in crystal unit cell dimensions and the hydrogen-bonding pattern also compared well with experimental results.

1. Introduction Cellulose provides mechanical reinforcement to plant cell walls and plant structures.1,2 During biogenesis, cellulose is synthesized in the plasma membrane of the plant cell3 and the chains associate in highly crystalline microfibrils with lateral dimensions of 3-10 nm. The high axial modulus and strength of cellulose microfibrils contribute strongly to the mechanical properties required of the plant cell wall. Cellulosic fibers in paper products also derive tensile strength and modulus from the microfibrils. Native cellulose is formed as the allomorph cellulose I, and its structure at room temperature is well documented. A major discovery was made with 13C CP/MAS solid-state NMR experiments, which revealed that native cellulose is a mixture of two distinct crystalline phases named IR and Iβ.4 The IR/Iβ ratio depends on the fiber origin, and cellulose Iβ is the more abundant form in higher plants. In these two structures, the cellulose chains are arranged in a parallel manner and are organized in sheets stabilized by interchain OH‚‚‚O hydrogen bonds, whereas the stacking of the sheets is stabilized by both Van der Waals forces and weak CH‚‚‚O hydrogen bonds.5 Recently, Yano et al. produced films and nanostructured composites based on cellulose microfibrils.6 These nanocom* Corresponding author. Phone: 33 (0)4 76 03 76 39. Fax: 33 (0)4 76 54 72 03. E-mail: [email protected]. † KTH Department of Fibre and Polymer Technology. ‡ Centre de Recherches sur les Macromole ´ cules Ve´ge´tales (CERMAVCNRS).

posites are candidate materials for flexible displays. In such application, low thermal expansion is essential and the nanocomposites meet this requirement. Low thermal expansion of microfibrils obviously derives from the structure and properties of the cellulose crystal, and it is of interest to learn more about the effect of temperature on crystal dimensions and modulus. However, molecular aspects of the thermal behavior of cellulose are not so well known. Early attempts with X-ray7-9 and solid-state NMR10 experiments had the deficiency that the samples showed low crystallinity. More recent X-ray and FTIR experimental studies on highly crystalline cellulose samples11-13 provide some insight. For all, they prove the existence of a hightemperature crystalline phase. The high-temperature phase seems to be reached at a transition temperature of about 450-500 K, but the results also indicate gradual structural changes at even lower temperatures. These experimental results induce a need of fundamental knowledge about the high-temperature structure, the temperature dependence of the structure, and its transition temperature. Molecular modeling has already been used to study crystalline cellulose. The span of modeling protocols includes both simple minimizations on frozen structures and molecular dynamics (MD) simulations at constant volume (NVT) or pressure (NPT). Results are thereby achieved either at 0 K for the frozen structures or at room temperature for the MD studies. The system structures have been either minicrystals14-20 or infinite crystals.21-26 All starting structures for the cellulose models are constituted by a limited number of cellulose chains having a low degree of polymerization, and they are systematically

10.1021/jp072258i CCC: $37.00 © 2007 American Chemical Society Published on Web 07/12/2007

Thermal Response in Crystalline Iβ Cellulose generated from duplicates of the primitive crystal structure with 2 chains of 2 residues for the Iβ allomorph. The minicrystals are either simulated isolated14-16,27 or immersed in water.17,18,20,28 In the latter case, the model is inserted in a periodic box with dimensions larger than those of the cellulose model and the empty space of the box is filled with solvent molecules. Infinite crystals are modeled within a computational box with dimensions of exact correspondence to multiples of the primitive cell, and each cellulose chain is covalently bonded to its own periodic image across the computational box. The system thereby becomes an infinite crystal because of the periodic boundary conditions.21-26 From experiments it is known that typical dimensions of a native cellulose microfibril are approximately 10-100 nm2 in cross section and the microfibril length is in the micrometer scale, corresponding to degrees of polymerization of the chains of about 10000 or more. The two different modeling approaches therefore provide different possibilities and suffer from different shortcomings. Minicrystal models made of 528 up to 48 chains17,18 are realistic in cross section. However, the lengths of the cellulose chains of 10, 20, or 40 residues correspond to an excessively short longitudinal dimension. These models have been successfully used to describe the behavior of the surface chains at the cellulose/water interface and were able to capture the fiber twist, an important and experimentally observed structural feature. The infinite crystals, on the other hand, correctly describe the almost endless longitudinal dimension of the microfibril. Therefore infinite crystals are well adapted to study bulk structures since they have infinite cross-sectional areas so that surface effects are excluded. The obvious difference from solvated minicrystals is that these models are not able to model chain end effects and twisting of the chains. Variation of crystal structure with temperature has, to the authors’ knowledge, never before been investigated by molecular modeling. The present study aims to fill this gap and provide detailed information about the thermal behavior of cellulose Iβ crystals. To this end, a cellulose model was implemented in the GROMACS 3.3 simulation software package.29 We chose an infinite crystal to focus on bulk effects. The objective was primarily to verify the ability of the force-field to correctly reproduce the crystal structure at room temperature. Second, the objective was to apply the model aiming to achieve physical insight into the correlation between changes in structure and properties of the cellulose Iβ crystal at elevated temperature. 2. Methods Throughout this study, atomistic molecular dynamics (MD) simulations were performed with the GROMACS 3.3 simulation software package.29 A large number of input parameters are allowed to be specified within the program and the most important, differing from default values, are being accounted for here. 2.1. System Setup. The cellulose Iβ bulk crystal was studied. A periodic box with initial box lengths Lxx ) 4.368 nm, Lyy ) 2.497 nm, and Lzz ) 4.347 nm, and angles xz ) yz ) 90° and xy ) 87° was filled with 8 × 4 chains, each having 8 glucose units (which corresponds to a total of 64 primitive unit cells). The zz-axis was parallel with the chain axis. An illustration of two glucose units linked together, as they are in cellulose, with corresponding atom names is presented in Figure 1. The chains were covalently bonded to their own periodic images in order to mimic chains of infinite length. Coordinates for the starting structure were taken from a structure experimentally determined by X-ray diffraction.5 The experimental

J. Phys. Chem. B, Vol. 111, No. 30, 2007 9139

Figure 1. Schematic representation of the cellobiose repeat unit with definitions of conformational torsion angles.

results show two different hydrogen-bonding patterns, and we chose hydroxyl hydrogen coordinates corresponding to the dominant one.30 A united atom force-field with explicit hydroxyl hydrogens was applied; the GROMOS 45a4 force-field for hexopyranosebased carbohydrates,31 based on GROMOS 45a3.32,33 It is the latest version available for this kind of molecules. This forcefield has been designed for simulations of small molecules in water solution but has never before been used for crystal bulk simulations. All bond lengths were subjected to constraints using the Linear Constraint Solver (LINCS) algorithm.34 Electrostatic interactions, including hydrogen bonds, were treated with a Coulomb potential and calculated with the Fast Particle Mesh Ewald method35 which, briefly, uses Ewald summation in direct space and Fast Fourier Transforms to improve performance of the reciprocal sum. Parameters used for electrostatic interactions were the electrostatic cutoff of 1.0 nm, a grid spacing of 0.1 nm, and an interpolation order 4. Weak interactions were considered with a Lennard-Jones potential within a cutoff radius of 1.0 nm, and a dispersion correction of energy and pressure was used. Both electrostatic and weak interactions were excluded for atom pairs separated by less than 3 bonds. A leapfrog algorithm was applied for integrating the equations of motion in the MD simulations36 and always with a time step of 1 fs. 2.2. MD Simulations at Different Temperatures. For the simulations on which the thermal behavior analysis is based, the crystalline experimental structure was equilibrated with NPT MD at 300 K for 1 ns. Subsequently, 10 ns NPT MD was performed at 300, 350, 400, 450, 500, and 550 K. Different temperatures were initially achieved by assigning different velocities to the atoms in the system. Velocities were generated randomly according to a Maxwell distribution. In all simulations the system was kept at constant temperature using the Nose´Hoover thermostat.37,38 The pressure was controlled with an anisotropic Berendsen barostat39 and kept constant at a reference pressure of 1 bar in the xx, yy, and zz directions and 0 bar in the off-diagonal xy/yx, xz/zx, and yz/zy directions. The compressibility in all directions was set to 1 × 10-6 bar -1. Global system properties such as energy and density and averaged values of, for example, chain conformational parameters were calculated as averages over the last 5 ns. Errors were estimated as standard deviations. Angle distributions were constructed by collecting all angle values for the last 5 ns of the simulations in bins, each with a width of 1°. Unit cell parameters were measured by locally regarding 10 different unit cells in the system. Corresponding distances and angles between O4 atoms were traced during the last 5 ns. The possible risk of artifacts, caused by the fact that the system has different x and y dimensions, was excluded by a 5 ns control simulation at 300 K with a larger system (8 × 8 chains with 8 glucose units) in which the unit cell dimensions and density were found to be the same as in the smaller system at the same temperature.

9140 J. Phys. Chem. B, Vol. 111, No. 30, 2007

Bergenstråhle et al. All MD run input parameters except temperature were identical with those for the 10 ns simulations. Statistics from the entire 200 ps trajectories were used for analysis. 3. Results and Discussion

Figure 2. Average densities as a function of temperature. Standard deviations are in the range (0.5-1.4, that is, less than 0.1%.

Hydrogen bonds were analyzed by geometric criteria which were defined by a cutoff distance between donor and acceptor oxygen of 0.35 nm and a cutoff angle donor-hydrogenacceptor larger than 150°. 2.3. Longitudinal Young’s Modulus. Crystalline Young’s modulus in chain direction was estimated for structures at 300 and 500 K. Structures well equilibrated with NPT MD at these temperatures (end structures from 10 ns MD) were placed in a box fixed at the equilibrated sizes, with Lzz ) L0. A 500 ps NVT relaxation was performed, and then the box size was modified, Lzz ) L0 ( 0.1%, 0.2%, and 0.3%. Steepest descent energy minimization was applied, followed by another 200 ps NVT MD. Young’s modulus is defined as

E ) σ/ Stress (σ) in chain direction was estimated by taking the negative average value of the internal pressure in the zz direction, and strain () was ∆Lzz/L0. The error here was estimated by dividing the trajectory into 10 parts, taking the average for each part, and applying the central limit theorem. 2.4. Thermal Expansion Coefficient at Room Temperature. For unit cell parameters a and b, the thermal expansion coefficients (TEC) were estimated at room temperature. The structure equilibrated at 300 K was used as the starting structure for 200 ps NPT MD simulations at 280, 290, 310, and 320 K.

Atomistic MD simulations were applied to investigate how the structure (and corresponding properties) of a cellulose Ιβ crystal is affected by increasing temperature, from 300 to 550 K, close to the cellulose degradation temperature. An obvious result of the simulations is that the time series of a selected characteristic (e.g., density) changes within the first nanoseconds; this occurs at all temperatures. Importantly, no further changes of the characteristic take place during the last 5 ns of the simulations, and the structures are therefore considered equilibrated. 3.1. Bulk Properties. Results presented here are estimations of density, crystal unit cell parameters, and thermal expansion coefficients at 300 K and Young’s modulus at 300 and 500 K. 3.1.1. Density. An important basic structural parameter is density. When performing MD simulations in NPT ensemble, system volume, and consequently its density, is allowed to change so that the system is kept at constant atmospheric pressure. Average density evolution with temperature is presented in Figure 2. At room temperature (300 K), the average equilibrated density of crystalline cellulose Iβ was found to be 1498 kg/m3. The polymer handbook indicates that the density of cellulose in cotton ranges between 1545 and 1585 kg/m3.40 Cotton cellulose is a mixture of IR and Iβ allomorphs and has 70% crystallinity; its density is therefore hardly comparable with that predicted in this study. Densities of the cellulose Ιβ crystal in the range of 1543-1643 kg/m3 can be estimated from its unit cell geometries measured by diffraction experiments (see Table 1). Thereby, our predicted density is 2.9-8.8% less than the experimental values. Other predicted values from molecular modeling range between 1491 and 1691 kg/m3.20,41 At higher temperatures, the results show that cellulose crystal density decreases nonlinearly to an average value of 1374 kg/ m3 at 550 K. Density is a macroscopic parameter, and it is of importance to go further into the details at the local scale of crystalline cellulose Iβ. 3.1.2. Unit Cells. A crystal structure is often described in terms of its unit cell geometry. Thermal variation of unit cell

TABLE 1: Average Unit Cell Parameters (in nm) and Angles (in degrees)a average unit cell parameters (in nm)

a

angles (in degrees)

a

b

c

R

β

γ

reference

0.8378 0.0002 0.9005 0.0055 0.778 0.801 0.778 0.785 0.794 0.786 0.771 0.790 0.803 0.806 0.847 0.760 0.760 0.750

0.8168 0.0002 0.8659 0.0069 0.820 0.817 0.820 0.814 0.839 0.817 0.834 0.835 0.813 0.851 0.811 0. 810 0.874 0.870

1.0523 0.0002 1.0470 0.0002 1.038 1.036 1.034 1.034 1.058 1.038 1.037 1.030 1.039 1.040 1.051 1.040 1.046 1.036

89.97 0.02 77.62/103.17 0.36/1.17 90 90 90 90 90 90 90 90 90 90.1 90 90 90.4 90

89.97 0.03 78.30/101.41 1.11/0.45 90 90 90 90 90 90 90 90 90 90.6 90 90 92.3 90

90.92 0.03 88.30 3.04 96.5 97.3 96.5 96.6 98.0 97.0 97.6 96.0 97.5 96.6 90 96 103.5 94.1

model 300 K SDb model 500 K SDb Nishiyama5 Sugiyama42 Woodcock43 Sarko44 Honjo45 Gardner46 Takahashi47 Meyer48 Reiling22 Neyertz23 Matthews20 Kroon-Batenburg41 Mazeau24 Vietor27

Overview of the experimental and modeled unit cells from the literature. b SD: standard deviations.

Thermal Response in Crystalline Iβ Cellulose

J. Phys. Chem. B, Vol. 111, No. 30, 2007 9141

Figure 3. The equilibrated model, view parallel to the chain axis at 300 K (top) and 500 K (bottom).

axes a, b, and c and angles R, β, and γ was evaluated; the ab projection of one of the considered unit cells in the structure equilibrated at 300 K is depicted in Figure 3 (top). At low temperature, the unit cells are well defined and easy to trace for statistically acceptable results. Values obtained at room temperature are presented in Table 1 along with the experimental ones. Parameters b and c (which corresponds to the cellulose chain axis) are within the range of reported experimental values; similarly, angles R and β are predicted close to the expected value of 90°. On the contrary, parameters a and γ are too large. The well-reproduced values for the b and c axes demonstrate that intrasheet as well as intrachain OH‚ ‚‚O hydrogen bonds are properly described. On the other hand, the discrepancy for the a axis is probably due to poorly described Van der Waals interactions in combination with too weak intersheet coulomb interactions between the united CH atoms and oxygen. This drawback is probably also responsible for the γ error. Results from other modeling of Iβ cellulose are also included in Table 1; all the tested force-fields, including CHARMM,20,22 GROMOS,41 PCFF,24 MM3,27 and the new force-field specially developed to study crystalline cellulose23 predict values of a and b systematically out of the correct ranges. Average values for all unit cell parameters at different temperatures are presented in Figure 4 I-IV. These figures clearly reveal the existence of a low-temperature structure and a high-temperature structure. Furthermore, it seems that 450 K is in the vicinity of a transition temperature where the cell parameters behave erratically (this point will be discussed in section 3.3). The graphical representation of the crystal structure equilibrated at 500 K, taken as a representative high-temperature structure, is given in Figure 3 (bottom), and its structural

characteristics are presented in Table 1. This high-temperature structure has a and b parameters that are 7.4% and 6.0% larger than in the low-temperature structure. The c parameter is only 0.5% smaller, and the γ angle also decreases by 2.9%. Concerning the R and β angles in the high-temperature structure, the model shows that the chains have moved in chain direction relative to each other. The appearance of R and β in two different unit cells at high- and low-temperature structures is shown in Figure 5, clarifying how chain sliding causes the angle distributions of R (top) and β (bottom) to separate into two populations at high temperature. These two populations appear symmetrically around 90°, and their values are about 13-14% different in the high-temperature structure compared to the lowtemperature structure. Thermal behavior of cellulose was experimentally studied on tunicate samples11 (large crystallites considered pure in the Iβ allomorph), cotton,9 and wood cellulose crystals49 (smaller crystallites with only a dominant fraction of Iβ) by X-ray diffraction. In agreement with our results, experiments suggest a crystalline phase transition where a given parameter (a or b axis) shows a sudden change in behavior. There is, however, a discrepancy about the temperature at which the transition occurs; it is measured at 453 K,49 473 K,9 or in a range between 493 and 503 K.11 Our modeling protocol precludes an exact estimate of the transition temperature; it is estimated at slightly lower values than the experimental ones, at around 450 K. However, a gradual transition could also be interpreted from the experimental results, since a starts to increase nonlinearly even before the suggested transition point. Also the FT-IR experiments by Watanabe et al.12 indicate gradual structural changes. Our study is thus

9142 J. Phys. Chem. B, Vol. 111, No. 30, 2007

Bergenstråhle et al.

Figure 4. Variation of unit cell parameters with temperature. (I) Axes a and b; (II) angles R and β; (III) chain axis c; (IV) angle γ.

Figure 6. Unit cell parameters a and b with linear fits to estimate TEC at 300 K.

TABLE 2: Predicted Thermal Expansion Coefficients (nm K-1 10-5) of the a and b Cell Parameters at Room Temperature

Figure 5. Chain sliding effect on unit cell angles R (top) and β (bottom).

reporting a thermal behavior in qualitative concordance with experimental results. The results from all experimental studies show that the a axis increases remarkably with increasing temperature. However, the behavior of the b and c axes is less clear. In the tunicate sample, b remains constant until the transition where it decreases;11 however, in the wood samples it decreases throughout the heating procedure.49 The c parameter measured for the wood samples is slightly increasing at all temperatures,49 whereas in cotton samples it decreases and increases below and above the transition point, respectively.9 It is worth noting that the unit cell parameters in those studies were calculated from d-spacings; therefore, the R and β angles, and consequently, movements in chain directions, were not considered. Our results show, in agreement with the experiments, an increasing a axis. In addition, the b axis also increases. Thermal changes in the a direction are mainly controlled by weak Van

TECa TECb a

model

exptla

7.3 1.2

6.8 0.5

Experimental values interpreted from Wada.11

der Waals interactions and CH‚‚‚O hydrogen bonds, while in the intrasheet b direction OH‚‚‚O hydrogen bonds are involved. Our results show that the cohesion between two sheets of cellulose is lower than that between cellulose chains in a single sheet; this is in qualitative agreement with recent ab initio calculations on cellulose Iβ.50 3.1.3. Thermal Expansion Coefficients at Room Temperature. Thermal expansion coefficients (TEC) were estimated only at room temperature. Figure 6 reports the estimated variations of a and b crystal cell dimensions, from 280 to 320 K. A leastsquares linear fit was performed on the results, and the TEC was calculated from the slopes. A similar fit was made on the a and b parameters reported in Figure 5 in the paper by Wada,11 and a comparison of those values of TEC and ours are presented in Table 2. In his paper, Wada presents another kind of TEC estimation based on a fit with a second-order polynomial. This might give better results over the entire temperature range from room temperature up to

Thermal Response in Crystalline Iβ Cellulose

Figure 7. Predicted stress-strain curves for the low (300 K) and high (500 K) structures.

Figure 8. Distribution (in %) of the three staggered conformations of the torsion angle ω.

473 K, but locally at 300 K a linear fit is a better choice. The anisotropic difference between a and b is apparent in both cases. 3.1.4. Longitudinal Young’s Modulus. The Young’s modulus in chain direction is an important mechanical characteristic. NVT simulations at different small strains in the chain direction were performed to estimate this property for the low-and the hightemperature structures, and the result is presented as a stressstrain graph in Figure 7. The modulus at 300 K was calculated to be 156 and 117 GPa at 500 K. The longitudinal Young’s modulus at 300 K agrees with earlier results, obtained both experimentally (120-143 GPa51-55) and theoretically (124173 GPa25,28,53,56). 3.2. Structural Aspects. In this section, the key conformational parameters of the individual cellulose chains are considered. Two important structural parameters are the conformation of the hydroxymethyl group and the conformation of the glycosidic bond. 3.2.1. Orientation of the Hydroxymethyl Group. Each glucose unit on the cellulose chains has a hydroxymethyl group. Its orientation can be described by the torsion angle ω which has three staggered positions (gt, tg, and gg) at 60°, 180°, and 300°, respectively. The simulation results are presented in Figure 8 as a histogram, which shows the percentage of torsion angles in these gt, tg, and gg states at different temperatures. At 300 K, all hydroxymethyls stay in their tg conformation during the simulation time, the number of occurrences in other conformations is negligible. This is in agreement with experimental

J. Phys. Chem. B, Vol. 111, No. 30, 2007 9143 observations from X-ray diffraction5 and solid-state NMR spectroscopy.57 At higher temperatures, the gg and gt conformations become more frequently populated, and at 550 K 83% of the angles are in gt and 16% in gg conformation; only a negligible part is still in tg. This conformational change is in accordance with experimental results from a CP/MAS 13C NMR spectroscopy study on cotton cellulose treated with steam at different temperatures.10 The experiments show that the signal corresponding to the C6 atom is affected at 553 K compared to that at room temperature, indicating a conformational change of the hydroxymethyl groups.57 3.2.2. Conformation around the Glycosidic Bonds. The glycosidic bond conformation was studied in terms of the two torsion angles, Φ and Ψ. The result indicates that the larger the temperature, the larger is the (Φ,Ψ) area explored. As a consequence, a little broader dispersion of the torsion angles at 500 K is observed when compared to the one at 300 K, which indicates a more disordered structure at high temperatures. On the other hand the center of the point collection has shifted slightly. The average (Φ, Ψ) pair is (264°, 220°) at 300 K and (271°, 216°) at 500 K. 3.2.3. Hydrogen Bonds. The OH‚‚‚O hydrogen-bonding network in cellulose Iβ is widely recognized as being of utmost importance for the stability of the crystal structure. In this study, hydrogen bonds were analyzed in the low-temperature structure (300 K) as well as in the high-temperature structure (500 K). Selected hydrogen bonds are presented in Figure 9 as the percentage number of maximum possible occurrences. Additionally, distance distributions for each hydrogen bond are reported. We were principally interested in the intrachain hydrogen bonds O2H2‚‚‚O6 and O3H3‚‚‚O5 as well as the intrasheet (but interchain) O6H6‚‚‚O3 and O6H6‚‚‚O2 ones, as they are considered responsible for single chain conformation and crystal moduli in native cellulose.25,52,56 The figure shows that increased temperature causes remarkable changes in the hydrogen-bonding pattern. At 300 K, both intrachain hydrogen bonds O2H2‚‚‚O6 and O3H3‚‚‚O5 and intrasheet O6H6‚‚‚O3 are present almost everywhere in the structure. This shows that the dominant hydrogen-bonding network determined from neutron diffraction of deuterated samples,5 chosen for the starting structure, has not changed during the simulation. On the other hand, at 500 K, the number of intrachain hydrogen bonds has decreased or completely disappeared, in the case of O2H2‚‚‚O6 (intrachain). Two new types of intrasheet (interchain) hydrogen bonds have appeared simultaneously, between O6H6‚‚‚O2 and O2H2‚‚‚O6. This shows that the hightemperature phase has a more disordered nature, since the chains are no longer showing exactly the same hydrogen-bonding pattern. Temperature-dependent changes in hydrogen bonds in cellulose Iβ have recently been studied by IR spectroscopy.12,13 A gradual disruption of the O3H3‚‚‚O5 and O2H2‚‚‚O6 intrachain hydrogen bonds in the temperature range 303-473 K was observed. In addition, around 493 K, drastic collapses of these hydrogen bonds were detected and a new structure with weak hydrogen bonds was found. Our results are in agreement with those experimental results. Distance distributions of these hydrogen bonds show that the O6H6‚‚‚O3 hydrogen bonds are slightly stronger at 300 K than at 500 K. On the other hand, the few occurrences of O6H6‚‚‚ O2 at 300 K are weaker than those at 500 K. The O2H2‚‚‚O6 (intrachain) hydrogen bond at 300 K shows a narrow distance

9144 J. Phys. Chem. B, Vol. 111, No. 30, 2007

Bergenstråhle et al.

Figure 9. Percentage of maximum possible intrachain and intrasheet hydrogen bonds and their corresponding distance distributions.

Figure 10. Structure equilibrated at 450 K.

distribution around a short distance which indicates that it is strong. At 500 K, this hydrogen bond has disappeared and has been replaced by a new intrasheet O2H2‚‚‚O6 hydrogen bond which is weaker. There is only a minor difference between O3H3‚‚‚O5 (intrachain) at 300 and 500 K. The reduced number of intrachain hydrogen bonds is one of the reasons explaining why the high-temperature structure shows a lower Young’s modulus in Figure 7. 3.3. Suggested Transition Mechanisms. The results that were discussed in the former section were those corresponding to the low-temperature structure at 300 K and the hightemperature structure at 500 K. In this section, the order in which the different structural changes were detected in the model between 300 and 500 K will be discussed. The transition between low-temperature and high-temperature structures seems to be a gradual process, according to the simulation results. At the first stage, between 300 and 400 K, the intersheet distance has increased, illustrated by an increasing unit cell a parameter. The sheets have also moved with respect to each other in the b direction and thereby slightly decreased the angle γ. At 450 K, the chains adopt two different orientations in the ab plane. In Figure 10, a graphical representation of the end structure at 450 K clarifies this issue. Because of chain sliding within the sheets, the b parameter is divided into two populations at this temperature and so is the angle R, as depicted earlier in

Figure 5. The parameters a and γ explore their highest and lowest values, respectively, at this temperature. At this temperature, every second chain has rotated about 30° with respect to the initial position in the ab plane. In addition to this rotation, the hydroxymethyl groups on the rotated chains have changed conformation from tg to gg. The hydroxymethyl groups on the remaining chains also change their conformation from tg to gt. After 450 K, the sheets slide with respect to each other in the chain direction, and thereby the third unit cell angle, β, also separates into two populations, also shown in Figure 5. The chains have then returned to their initial positions in the ab plane and the hydroxyl groups are now all in gt. 4. Conclusions The GROMOS 45a4 force-field was applied to a cellulose Iβ crystal model, and MD was used to simulate the model in the GROMACS 3.3 package. The simulation at room temperature was found to reproduce the crystal structure sufficiently well in terms of crystal density, packing and crystal unit cell dimensions, chain conformation parameters, hydrogen bonds, and Young’s modulus. Second, insight into changes at the molecular level in the cellulose Iβ crystal heated to different temperatures was achieved. As in experiments, the performed simulations suggested that the cellulose Iβ crystal will be restructured when

Thermal Response in Crystalline Iβ Cellulose heated above 450 K. The simulation results indicated an expansion of the crystal unit cell and a change in orientation of the hydroxymethyl groups from tg to gt causing a changed hydrogen-bonding pattern. Young’s modulus, which is known to be dependent on intrachain hydrogen bonds, decreased since the high-temperature structure contained less hydrogen bonds of this type. The present MD model was able to predict a complex physical phenomenon such as thermal response of the cellulose Iβ allomorph in acceptable agreement with experimental results. Comparison between predictions and experimental data will help to contribute toward an improved understanding of how crystal properties depend on specific details of molecular structure, such as the hydrogen-bonding pattern. Acknowledgment. This research project has been supported by a Marie Curie Early Stage Research Training Fellowship of the European Community’s Sixth Framework Program, under Contract No. MEST-CT-2004-503322. M.B. and L.A.B acknowledge the Swedish Agency for Innovation Systems (VINNOVA) for financial support. Olle Edholm (Royal Institute of Technology, Stockholm) is thanked for kindly granting us computer time on the PC cluster owned by the Theoretical Biological Physics group. Many thanks also to Jakob Wohlert (Royal Institute of Technology, Stockholm) for discussions and opinions on model construction and analysis. References and Notes (1) Berglund, L. Cellulose Based Nanocomposites; CRC Press: Boca Raton, FL, 2005; Vol. 26. (2) Koehler, L.; Telewski, F. W. Am. J. Bot. 2006, 93, 1433. (3) Somerville, C. Annu. ReV. Cell DeV. Biol. 2006, 22, 53. (4) VanderHart, D. L.; Atalla, R. H. Macromolecules 1984, 17, 1465. (5) Nishiyama, Y.; Langan, P.; Chanzy, H. J. Am. Chem. Soc. 2002, 124, 9074. (6) Yano, H.; Sugiyama, J.; Nakagaito, A. N.; Nogi, M.; Matsuura, T.; Hikita, M.; Handa, K. AdV. Mater. 2005, 17, 153. (7) Wakelin, J. H.; Sutherland, A.; Beck, L. R., Jr. J. Polym. Sci. 1960, 42, 278. (8) Seitsonen, S.; Mikkonen, I. Polym. J. 1973, 5, 263. (9) Takahashi, M.; Takenaka, H. Polym. J. 1982, 14, 675. (10) Horii, F.; Yamamoto, H.; Kitamaru, R.; Tanahashi, M.; Higuchi, T. Macromolecules 1987, 20, 2946. (11) Wada, M. J. Polym. Sci., Polym. Phys. 2002, 40, 1095. (12) Watanabe, A.; Morita, S.; Ozaki, Y. Biomacromolecules 2006, 7, 3164. (13) Watanabe, A.; Morita, S.; Ozaki, Y. Appl. Spectrosc. 2006, 60, 611. (14) Aabloo, A.; French, A. D.; Mikelsaar, R.-H.; Pertsin, A. J. Cellulose 1994, 1, 161. (15) Aabloo, A.; French, A. D. Macromol. Theory Simul. 1994, 3, 185. (16) French, A. D.; Miller, D. P.; Aabloo, A. Int. J. Biol. Macromol. 1993, 15, 30.

J. Phys. Chem. B, Vol. 111, No. 30, 2007 9145 (17) Yui, T.; Hayashi, S. Biomacromolecules 2007, 8, 817. (18) Yui, T.; Nishimura, S.; Akiba, S.; Hayashi, S. Carbohydr. Res. 2006, 341, 2521. (19) Tanaka, F.; Okamura, K. Cellulose 2005, 12, 243. (20) Matthews, J. F.; Skopec, C. E.; Mason, P. E.; Zuccato, P.; Torget, R. W.; Sugiyama, J.; Himmel, M. E.; Brady, J. W. Carbohydr. Res. 2005, 341, 138. (21) Hardy, B. J.; Sarko, A. Polymer 1996, 37, 1833. (22) Reiling, S.; Brickmann, J. Macromol. Theory Simul. 1995, 4, 725. (23) Neyertz, S.; Pizzi, A.; Merlin, A.; Maigret, B.; Brown, D.; Deglise, X. J. Appl. Polym. Sci. 2000, 78, 1939. (24) Mazeau, K.; Heux, L. J. Phys. Chem. B 2003, 107, 2394. (25) Kroon-Batenburg, L. M. J.; Kroon, J.; Northolt, M. G. Polym. Commun. 1986, 27, 290. (26) Heiner, A. P.; Sugiyama, J.; Teleman, O. Carbohydr. Res. 1995, 273, 207. (27) Vietor, R. J.; Mazeau, K.; Lakin, M.; Perez, S. Biopolymers 2000, 54, 342. (28) Tanaka, F.; Iwata, T. Cellulose 2006, 13, 509. (29) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701. (30) Mazeau, K. Cellulose (Dordrecht, Netherlands) 2005, 12, 339. (31) Lins, R. D.; Hunenberger, P. H. J. Comput. Chem. 2005, 26, 1400. (32) Schuler, L. D.; Van Gunsteren, W. F. Mol. Simul. 2000, 25, 301. (33) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. J. Comput. Chem. 2001, 22, 1205. (34) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463. (35) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (36) Hockney, R. W. Methods Comput. Phys. 1970, 9, 136. (37) Nose´, S. Mol. Phys. 1984, 52, 255. (38) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (39) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (40) Brandrup, J.; Immergut, E. H.; Grulke, E. A. Polymer Handbook, 4th ed.; Wiley-Interscience: New York, 1999. (41) Kroon-Batenburg, L. M. J.; Kroon, J. Glycoconjugate J. 1997, 14, 677. (42) Sugiyama, J.; Vuong, R.; Chanzy, H. Macromolecules 1991, 24, 4168. (43) Woodcock, C.; Sarko, A. Macromolecules 1980, 13, 1183. (44) Sarko, A.; Muggli, R. Macromolecules 1974, 7, 486. (45) Honjo, G.; Watanabe, M. Nature (London, United Kingdom) 1958, 181, 326. (46) Gardner, K. H.; Blackwell, J. Biopolymers 1974, 13, 1975. (47) Takahashi, Y.; Matsunaga, H. Macromolecules 1991, 24, 3968. (48) Meyer, K. H.; Misch, L. HelV. Chim. Acta 1937, 20, 232. (49) Hori, R.; Wada, M. Cellulose 2005, 12, 479. (50) Qian, X.; Ding, S.-Y.; Nimlos, M. R.; Johnson, D. K.; Himmel, M. E. Macromolecules 2005, 38, 10580. (51) Sakurada, I.; Nukushina, Y.; Ito, T. Kobunski Kagaku 1962, 19, 300. (52) Nishino, T.; Takano, K.; Nakamae, K. J. Polym. Sci., Polym. Phys. 1995, 33, 1647. (53) Sturcova, A.; Davies, G. R.; Eichhorn, S. J. Biomacromolecules 2005, 6, 1055. (54) Matsuo, M.; Sawatari, C.; Iwai, Y.; Ozaki, F. Macromolecules 1990, 23, 3266. (55) Sakurada, I.; Nukushina, Y.; Ito, T. J. Polym. Sci. 1962, 57, 651. (56) Tashiro, K.; Kobayashi, M. Polym. Bull. 1985, 14, 213. (57) Horii, F.; Hirai, A.; Kitamaru, R. Polym. Bull. 1983, 10, 357.