DFT Optimization of Isolated Molecular Chain Sheet Models Constituting Native Cellulose Crystal Structures.

Because of high crystallinity and natural abundance, the crystal structures of the native cellulose allomorphs have been theoretically investigated to elucidate the cellulose chain packing schemes. Here, we report systematic structure optimization of cellulose chain sheet models isolated from the cellulose Iα and Iβ crystals by density functional theory (DFT). For each allomorph, the three-dimensional chain packing structure was partitioned along each of the three main crystal planes to construct either a flat chain sheet model or two stacked chain sheet models, each consisting of four cello-octamers. Various combinations of the basis set and DFT functional were investigated. The flat chain sheet models constituting the cellulose Iα (110) and Iβ (100) planes, where the cellulose chains are mainly linked by intermolecular hydrogen bonds, exhibit a right-handed twist. More uniform and symmetrical sheet twists are observed when the flat chain sheet models are optimized using a basis set with diffuse functions (6-31+G(d,p)). The intermolecular interactions are more stable when the chain sheet models are optimized with the two hybrid functionals CAM-B3LYP and M06-2X. Optimization of the two stacked chain sheet models, where van der Waals interactions predominated between adjacent chains, gave differing results; those retaining the initial structures and those losing the sheet appearance, corresponding to the cellulose Iα/Iβ (010)/(11̅0) and (100)/(110) chain sheet models, respectively. The cellulose Iβ (11̅0) chain sheet model is more stable using the M06-2X functional than using the CAM-B3LYP functional.

Introduction

Cellulose ((1 → 4)-β-d-glucan), the most abundant biopolymer, is the main component of plant cell walls. Cellulose does not exist in a single molecular chain in organisms, but several cellulose chains assemble into high-crystallinity microfibrils.[1−3] In addition to the conventional uses of cellulose, it is expected that cellulose can potentially act as an advanced material.[4−7] In particular, cellulose nanofibers prepared from wood by refinement techniques have attracted attention as a next-generation material because of their high strength, high elasticity, light weight, and low environmental impact.[8−13] Such increasing attention has evoked interest in the atomistic details of the higher order structure of cellulose.

Cellulose forms various crystalline allomorphs derived from the two native crystalline phases: cellulose Iα and Iβ.[14−16] One important physical aspect of these allomorphs is that the Iα phase is metastable and can be converted into the Iβ phase by hydrothermal treatment[17] or high-temperature treatment in organic solvents or even helium gas.[18] It has been suggested that cellulose Iβ is the more stable form of the two native allomorphs. In the crystal structures of the native allomorphs, all of the molecular chains are aligned with the same polarity, which is known as parallel chain packing.[19−21] Cellulose IIII can be obtained by treatment of cellulose Iβ with liquid ammonia or amines,[22] and it readily converts back to cellulose Iβ by hot-water treatment[23] while retaining the parallel chain polarity.[24] In contrast, the unit cell of cellulose II has chains oriented in opposite directions.[25−28] The native crystal forms are irreversibly converted into cellulose II by either mercerization or regeneration processes, which is accompanied by a complete change in the cellulose chain polarity from parallel to antiparallel.[29]

In the crystalline form, the cellulose molecular chains have a flat ribbonlike conformation with twofold helical symmetry and amphiphilic nature, where the polar functional groups mostly align along the side of the flat ribbon and the ribbon face consists of pyranose faces with hydrophobic nature. The amphiphilic nature two-dimensionally expands as the cellulose chains assemble into the flat chain sheet that characterizes the chain packing manner of native cellulose crystals.[30,31] The cellulose chains are bound together by O3···O6 intermolecular hydrogen bonds, the representative intermolecular hydrogen bond observed in both cellulose Iα and Iβ allomorphs in common, in the chain sheets, and the O2 and O3 hydroxyl groups participate in intramolecular hydrogen bonds to stabilize the twofold helical molecular structure (see Figure ).[24,32−34] Each allomorph involves the two alternative hydrogen bond networks that are correlated to the diverse rotational positions of O2 and O6 hydroxyl groups.[32,33] The polar functional groups are neutralized by hydrogen bond formation, whereas van der Waals interactions may dominate at the interface of the chain sheets,[35] as can be seen from the ab projection of the native crystal model in Figure . The multiple weak hydrogen bonds were also suggested between adjacent chain sheets. Cellulose Iα and Iβ have similar chain sheet structures, but the chain packing structures differ in the relative displacement of the chain sheets along the chain axis.[16] Adjacent chain sheets are displaced with a quarter fiber repeat length, which is called one-quarter chain staggering.[19,20] In the cellulose Iα crystal,[33] one-quarter chain staggering progressively occurs in one direction to give a triclinic unit cell, whereas every other chain sheet is displaced with quarter staggering in the monoclinic unit cell of the cellulose Iβ crystal[32] (see Figure ). Jarvis suggested the unique features of these hydrogen bonding networks in the native cellulose crystals in comparison with the other biopolymers.[34] The zigzag, repeating O–H···O–H··· networks from the stringed hydrogen bonds running continuously along the fiber axis, which, coupled with the interconversion between the two alternative networks, may serve a proton conducting path in the crystal.

Figure 1

Chain sheet constituting the cellulose Iβ (100) plane. The blue dotted lines indicate hydrogen bonds.

Figure 2

ab projections of the native crystal models. The labels indicate the relative translational shift of the flat chain sheet along the fiber axis (c is the fiber repeat distance). The blue rhombic frames show the ab projections of the unit cells.

Density functional theory (DFT) calculations of the native cellulose crystal structures have been performed to analyze the hydrogen bond networks,[36,37] intermolecular forces,[38] vibrational spectra,[39,40] NMR spectra,[41,42] and electronic properties,[43,44] where supercell models consisting of one to several unit cells were subjected to DFT optimization. Recently, we proposed an alternative strategy to conventional supercell models for DFT modeling of the cellulose crystals, in which the structures of cellulose chain sheet models isolated from the parent crystal structures were optimized until the calculations reached convergence.[45,46] In contrast to supercell calculations to reproduce the original crystal structures and physical properties, we intended to observe the sheet structure deformation resulting from DFT optimization and assess the structural stability of the target chain sheet model from its deformation behavior. Generally, if one wants to evaluate intermolecular interactions with respect to each dimension, the interactions will be decomposed along a lower dimensional fragment, whereas the energy decomposition may have less physical legitimacy in the strict sense. It is likely that the chain sheet in isolation, free from the crystal packing force along the other dimension, is on the slope of the potential energy landscape. The structure would be deformed toward lower potential surface by DFT optimizations, revealing the intrinsic stress involved in the chain sheet. In addition to the cellulose Iα and Iβ chain sheet models, two chain sheet models from cellulose IIII and three from cellulose II were investigated.[46] The most striking result was that both of the native chain sheet models developed a right-handed twist while retaining the sheet form. A similar sheet twist was observed for the two sheet models derived from the cellulose II crystal with parallel chain polarity. The remaining cellulose II sheet model with antiparallel chain polarity deformed from the initial sheet form. In the two cellulose IIII sheet models, one essentially retained and the other lost its initial structure. Deringer et al. reported a DFT-D study of isolated building blocks corresponding to the molecular chain sheets constituting α-chitin, and they proposed a building block bundled with strong intermolecular hydrogen bonds to elucidate the nature of the hydrogen bonding network.[47] In their study, the concept of chain sheets was extended to a stacked molecular chain sheet in which the chains stacked through the pyranose ring faces.

The plans of the present study are as follows: (1) DFT optimization of the stacked chain sheet models along the (100) and (010) planes of the Iα crystal and (110) and (11̅0) planes of the Iβ crystal and re-investigation of the two flat chain sheet models along the Iα (110) and Iβ (100) planes (see Figure ) and (2) benchmarking of the various levels of theory defined by different DFT functional/basis set combinations for the present supermolecular systems. As mentioned above, different interactions dominate at the molecular chain interfaces for flat and stacked chain sheets, where polar interactions (represented by hydrogen bonds) mainly contribute to chain sheet formation and van der Waals interactions are the most significant at the interfaces between adjacent pyranose faces, respectively. A different DFT functional/basis set combination could be suitable for the two types of chain sheet models.

Computational Methods

The six types of molecular chain sheet models, the (110), (100), and (010) chain sheets of cellulose Iα and the (100), (110), and (11̅0) chain sheets of cellulose Iβ, were constructed based on the set of atomic coordinates reported by X-ray and neutron diffraction studies.[32,33] Each molecular chain sheet model consists of four cello-oligomers with the degree of polymerization of 8. As mentioned above, we classified the Iα (110) and Iβ (100) chain sheet models as the flat chain sheet model and the Iα (100), Iα (010), Iβ (110), and Iβ (11̅0) chain sheet models as the stacked chain sheet model.

The positions of the hydrogen atoms were initially optimized by classical mechanics calculations using the Glycam 06 parameters,[48] followed by quantum mechanics optimization of all of the atom positions. The level of theory for full optimization of the flat chain sheet models was either Hartree–Fock (HF) or DFT with either the B3LYP,[49,50] CAM-B3LYP,[51] LC-ωPBE,[52] or M06-2X[53,54] functional coupled with either the 6-31G(d), 6-31G(d,p), 6-311G(d,p), or 6-31+G(d,p) basis set. The stacked chain sheet models were optimized using either the CAM-B3LYP or M06-2X functional coupled with the 6-311G(d,p) basis set. The single-point energy calculations for all of the optimized models were performed with the 6-311+G(d,p) basis set.

The binding energy (ΔEbind) between adjacent oligomers is defined as the difference between the DFT energy of the total model (Etotal) and its separate parts (E1 and E2)where E1 and E2 are the energies of an interior oligomer and the remainder of the chain sheet, respectively. Because the four-chain sheet models contain two interior oligomers, the two ΔEbind values were averaged to give the final ΔEbind value. The basis set superposition error (BSSE) for the binding energy was estimated by counterpoise correction[55,56] with either the 6-311+G(d,p) basis set or the same basis set as those used for the optimizations.

The classical and quantum mechanics calculations were performed with the AMBER14 program[57] and Gaussian 09 program package,[58] respectively. The images were generated with PyMOL software (version 1.7.0.1, Schrödinger, LLC).[59]

In the flat chain sheet models, the degree of twisting deformation was estimated. As shown in Figure , the virtual bonds connecting the gravity centers of the residues define the torsion angle θ, which represents the twisting angle of the chain sheet along the fiber axis.[46,60,61] The θ value is the average of the θ values between the positive and negative residue positions. Obviously, θ depends on the residue position along the fiber axis. The largest θ value is expected at either of the terminal ends (θ±4), and θ decreases to nearly 0° on approaching the middle position, or θ±1 if an overall uniform twist occurs. When the chain sheet twists in the clockwise direction to give a right-handed twist, the θ value is defined as positive.

Figure 3

Molecular chain sheet and residue position labeling. The sheet-twist torsion angles are defined by the virtual bonds connecting the centers of gravity (G) of the residue(s): θ+4 = G(a,+4) – G[(b,+4), (c,+4)] – G[(b,–4), (c,–4)] – G(d,–4), θ–4 = G(a,–4) – G[(b,–4), (c,–4)] – G[(b,+4), (c,+4)] – G(d,+4), and so forth.

Results and Discussion

Optimized Structures of the Chain Sheet Models

Figure shows a comparison of the superimposed structures of the initial and CAM-B3LYP/6-311G(d,p)-optimized four cello-octamer chain sheet models along with the root mean square deviation (RMSd) values between the initial and DFT-optimized structures. In the case of the stacked chain sheet models, the neighboring chains slide by one-quarter staggering along the fiber axis (1/4c). Coupled with the twofold helical symmetry of the finite length cellulose oligomer, this gives the four edge patterns of the chain sheet models. These variations of the edge structures were also examined for the Iβ (110) and (11̅0) chain sheet models, and their CAM-B3LYP/6-311G(d,p) DFT-optimized structures are shown in Figure S1.

Figure 4

Superimposed structures of the crystal (blue) and CAM-B3LYP/6-311G(d,p) DFT-optimized (red) four cello-octamer chain sheet models of cellulose. The RMSd value between the initial and DFT-optimized structure is shown below each chain sheet model.

In Figure , the chain sheet models can be classified into three types with respect to their optimized structures: (1) the Iα (110) and Iβ (100) flat chain sheet models with a right-handed twist, (2) the Iα (100) and Iβ (110) stacked chain sheets with disordered deformation, and (3) the Iα (010) and Iβ (11̅0) stacked chain sheet models that essentially retain their initial structures. The differences in the edge structures do not have a significant effect on the optimized structures in terms of the above three classifications (see Figure S1). Movies S1–S3 show the DFT-optimization processes of chain sheet models (1)–(3), respectively. The right-handed twist of the flat chain sheet models was observed in our previous DFT studies, and we discussed the relation to the similar right-handed twist of the parent crystal models.[46,60,61] A right-handed fiber twist is a significant issue in real cellulose microfibrils,[62] and it has been observed by atomic force and transmission electron microscopy observations of native cellulose microfibrils.[63] Twisted cellulose ribbons have been observed in vivo during biosynthesis of bacterial cellulose.[64−66] In a wide-angle X-ray scattering study, Fernandes et al. suggested that the twisted cellulose microfibrils may decrease the distance along the fiber axis with their coalescence at the crystal faces of spruce wood.[67] The present results thus confirm that the fiber twist of the native cellulose crystal is driven by the inherent twisting stress in the flat chain sheet.

The deformed stacked chain sheet models can also be distinguished by their large RMSd values (>2 Å). When comparing the ab projection structures of the stacked chain sheet models, the molecular chains overlap more in the retained stacked chain sheets than in the deformed stacked chain sheets. We have reported similar DFT optimization behavior in the cellulose IIII (100) and (11̅0) chain sheet models, where all chains are connected by hydrogen bonds.[45,46,68] The latter chain sheet model has a more overlapped chain arrangement and retains the initial structure, whereas the former loses the sheet structure. From both X-ray diffraction measurements[23] and our previous molecular dynamics (MD) calculations,[45,69] the lines of cellulose chains comprising the cellulose IIII (11̅0) chain sheet are conserved in the cellulose Iβ crystal during crystalline conversion. In this respect, conservation of the cellulose Iβ (11̅0) chain sheet structure in the present study seems to be reasonable. In crystalline conversion to the cellulose I–ammonia complex, ammonia molecules penetrate between the (11̅0) chain sheets, leaving the cellulose chain arrays in the complex phase.[70−72] An MD study of the cellulose Iβ crystal models dissolved in imidazolium-based ionic liquids found similar diffusion behavior of the solvent molecules, where the (11̅0) chain sheet resisted penetration of the ionic liquid from the dispersion in the early stage of dissolution.[73−75] It should be also noted that, when comparing the retained stacked chain sheet models between the two allomorphs, the Iα (010) chain sheet gives larger RMSd values than those of the Iβ (11̅0) chain sheet, which may have partly resulted from less stable interactions between the adjacent chains in the former chain sheet.

Effect of the DFT Functional/Basis Set Combination on the Flat Chain Sheet Twist

Figures and 6 show the variation of the sheet twist angle (θ) with respect to the residue position for the Iα (110) and Iβ (100) chain sheet models optimized with different DFT functional/basis set combinations, respectively (see Tables S1 and S2 for the θ values). For comparison, the results of HF combined with different basis set optimizations are shown. Increasing the basis set size leads to a smoother and more symmetrical sheet twist angle curve, suggesting that a symmetrical chain sheet twist is indicative of a good description of the intermolecular interactions. The twist angle curves calculated with different functionals mostly coincide when using the basis set with diffuse functions (6-31+G(d,p)). Without diffuse functions, the DFT optimizations with the M06-2X functional tend to show an asymmetric sheet twist. The RMSd and difference in the sheet twist angle per cellobiose unit (Δθ) values for the two optimized flat chain sheet models are listed in Table S3. The results show that the larger the basis set, the more twisting there is of the flat chain sheet model. The twist is increased (from the least to most) in the order of HF, B3LYP, CAM-B3LYP, LC-ωPBE, and M06-2X. The cellulose chain conformations essentially retained the initial twofold helical symmetry (see Tables S4 and S5 for the glycosidic rotation angles of the optimized chain sheet structures), indicating that the apparent sheet twist has been mostly resulted from rearrangement of the cellulose chains. On the basis of the representative value of Δθ, 7°/cellobiose unit, the length of the complete fiber twist is estimated to be about 50 nm. On the other hand, the twist lengths of real cellulose crystalline fibers were observed to be the order of several hundreds of nanometers[63,67,76] and even a few micrometers for tunicin whiskers.[77] Possible interpretations for this difference may be twofold: (1) the chain sheet structures were optimized in vacuum where electrostatic interactions, especially hydrogen bonding interactions, were more emphasized than in the crystalline environment and (2) the interactions between the chain sheets should have restrained excess sheet twisting in the crystal packing. We have proposed in the earlier MD studies that the crystal fiber twist was depressed with increasing dimensions of the crystal models.[60,61]

Figure 5

Variation of the sheet twist angle (θ) with respect to the residue position for the DFT-optimized structures of the Iα (110) chain sheet model. Crystal; the sheet twist angle of the central chain sheet of the crystal model obtained from our previous MD calculations.[44]

Figure 6

Variation of the sheet twist angle (θ) with respect to the residue position for the DFT-optimized structures of the Iβ (110) chain sheet model. Crystal; the sheet twist angle of the central chain sheet of the crystal model obtained from our previous MD calculations.[44]

The ΔEbind and BSSE values between adjacent oligomers of the optimized Iα (110) and Iβ (100) flat chain sheet models are given in Table . These values were evaluated by single-point calculations of the optimized structures with the 6-311+G(d,p) basis set. The ΔEbind values obtained using the M06-2X and CAM-B3LYP functionals are the lowest for all of the basis sets, followed by those calculated with the B3LYP functional and HF. All of the HF energies are weaker than the DFT energies because of the poor estimation of the electron correlation energy in HF. Marom et al. and Kruse et al. suggested that the B3LYP functional alone cannot accurately estimate weak dispersion interactions.[78,79] Functionals with long-range correction, such as CAM-B3LYP and LC-ωPBE, and dispersion correction, such as the M06 family, have thus been developed,[51−54] and they enable relatively good estimation of the noncovalent interactions for weakly bonded molecular systems.[80−84] In the present cellulose chain sheet systems, the weak intermolecular interactions between the cellulose chains are well-estimated by the CAM-B3LYP and M06-2X functionals, giving the lowest ΔEbind values, followed by those obtained with the LC-ωPBE functional. Unlike the ΔEbind values, the choice of the basis set seems to have less effect on the BSSE irrespective of the functional. The ΔEbind and BSSE values derived from single-point calculations using the same basis sets as those used for the optimizations are given in Table . The BSSE values account for 24.3–34.4% of the uncorrected ΔEbind values without diffuse functions and 9.8–14.5% for the 6-31+G(d,p) optimizations. Allinger et al. reported that that hydrogen bonding energies between a water dimer and monosaccharides are overestimated using a basis set without diffuse functions.[85] In DFT calculations of carbohydrate molecules, although diffuse functions on the hydrogen atoms appear to be necessary, the 6-31+G(d) basis set is sufficient for geometry calculations and the 6-31+G(d,p) basis set gives good energy values when combined with pure or hybrid density functionals.[86−89] The intermolecular interactions of the present flat chain sheet models, which are dominated by hydrogen bonds, are more accurately estimated using basis sets with diffuse functions. It should be noted that the 6-31+G(d,p)-optimized flat chain sheet models also exhibit more symmetric twisting, as shown in Figures and 6.

Table 1

ΔEbind (kcal/mol per Residue) and BSSE (kcal/mol per Residue) Values for the HF- and DFT-Optimized Structures of the Flat Sheet Models Derived from Iα (110) and Iβ (100) Obtained by Single-Point Calculations with the 6-311+G(d,p) Basis Set

	6-311+G(d,p)//6-31G(d)		6-311+G(d,p)//6-31G(d,p)		6-311+G(d,p)//6-311G(d,p)		6-311+G(d,p)//6-31+G(d,p)
	ΔEbind	BSSE	ΔEbind	BSSE	ΔEbind	BSSE	ΔEbind	BSSE
Iα (110)
HF	–7.56	1.27	–7.51	1.36	–7.50	1.23	–7.65	1.22
B3LYP	–9.73	1.22	–9.68	1.43	–9.62	1.20	–9.88	1.17
CAM-B3LYP	–11.27	1.30	–11.25	1.30	–11.19	1.27	–11.46	1.25
LC-ωPBE	–9.86	1.43	–9.82	1.22	–9.79	1.39	–10.00	1.36
M06-2X	–11.39	1.43	–11.29	1.25	–11.26	1.34	–11.52	1.31
Iβ (100)
HF	–7.49	1.27	–7.43	1.22	–7.42	1.23	–7.58	1.22
B3LYP	–9.62	1.22	–9.57	1.23	–9.51	1.20	–9.78	1.17
CAM-B3LYP	–11.15	1.30	–11.13	1.30	–11.07	1.28	–11.35	1.25
LC-ωPBE	–9.74	1.43	–9.71	1.43	–9.74	1.39	n.d.a	n.d.a
M06-2X	–11.28	1.39	–11.21	1.38	–11.18	1.35	–11.44	1.31

Not determined. The optimization failed to reach convergence.

Table 2

ΔEbind and BSSE Values (kcal/mol per Residue) for the HF- and DFT-Optimized Structures of Flat Sheet Models Derived from Iα (110) and Iβ (100) Obtained by Single-Point Calculations at the Same Level as That of Optimization

	6-31G(d)		6-31G(d,p)		6-311G(d,p)		6-31+G(d,p)
	ΔEbind	BSSE	ΔEbind	BSSE	ΔEbind	BSSE	ΔEbind	BSSE
Iα (110)
HF	–7.93	2.60	–7.82	2.51	–7.25	2.60	–7.74	1.42
B3LYP	–9.67	4.75	–9.43	4.88	–8.91	4.32	–9.91	1.18
CAM-B3LYP	–11.16	5.02	–10.91	5.17	–10.40	4.57	–11.49	1.22
LC-ωPBE	–9.99	4.46	–9.79	4.59	–9.21	4.11	–10.14	1.24
M06-2X	–11.36	4.12	–11.04	4.07	–10.65	3.79	–11.49	1.37
Iβ (100)
HF	–7.84	2.60	–7.74	2.51	–7.17	2.60	–7.66	1.43
B3LYP	–9.56	4.75	–9.31	4.88	–8.79	2.34	–9.80	1.18
CAM-B3LYP	–11.02	5.03	–10.78	5.17	–10.27	4.60	–11.38	1.22
LC-ωPBE	–9.85	4.47	–9.67	4.59	–9.13	4.13	n.d.a	n.d.a
M06-2X	–11.30	4.02	–10.95	4.10	–10.57	3.86	–11.40	1.37

Not determined. The optimization failed to reach convergence.

Effect of the DFT Functional/Basis Set Combination on Optimization of the Iβ (11̅0) Stacked Chain Sheet Model

In addition to the CAM-B3LYP/6-311+G(d,p)//CAM-B3LYP/6-311G(d,p) calculations, DFT optimization of the four Iβ (11̅0) stacked chain sheet models, which differ in their edge structures, was performed with the M06-2X functional. There is no significant difference between the optimized structures obtained with CAM-B3LYP/6-311G(d,p) and M06-2X/6-311G(d,p). The ΔEsheet and BSSE values derived from single-point calculations using the same basis set (6-311G(d,p)) and those with diffuse functions (6-311+G(d,p)) are listed in Table . As observed in the previous DFT-optimized flat chain sheet models, introduction of diffuse functions lowers the ΔEsheet values and effectively reduces the BSSE from 49–73 to 28–39% of the uncorrected ΔEsheet values. This indicates a certain contribution of the hydrogen bonding interactions between the stacked cellulose chains. In fact, O6–H···O5 intermolecular hydrogen bonds are present in both the B3LYP-optimized Iα (010) and Iβ (11̅0) stacked chain sheet models. The same hydrogen bond is present in the MD structures of cellulose Iβ, but it is less frequent.[90] The ΔEsheet values are considerably stronger using the M06-2X functional compared with the corresponding CAM-B3LYP calculations. Dispersion interactions (i.e., van der Waals interactions) are more significant between the cellulose chains in the stacked chain sheets, so the M06-2X functional effectively corrects the dispersion interactions in the present DFT calculations.

Table 3

ΔEbind (kcal/mol per Residue) and BSSE (kcal/mol per Residue) Values of the DFT-Optimized Iβ (11̅0) Stacked Chain Sheet Modela

	6-311G(d,p)		6-311+G(d,p)//6-311G(d,p)
	ΔEbind	BSSE	ΔEbind	BSSE
CAM-B3LYP	–2.93	7.83	–4.41	2.83
M06-2X	–7.65	7.45	–8.67	3.31

The average values from the calculations of the four Iβ (11̅0) stacked chain sheet models with different edge structures.

Conclusions

In previous theoretical calculations of cellulose crystals, their structural stability has been mostly assessed on the basis of the thermochemistry, such as the binding and lattice energies, and their changes during crystalline transformations.[36−44] In these calculations, structure optimization was performed to refine the initial structures, resulting in slight adjustments of the atomic positions, such as the exocyclic group conformations. In the present study, the native cellulose crystal structures were dimensionally decomposed into three chain sheet models of finite sizes and each of the models was optimized by DFT with various DFT functional/basis set combinations. Being isolated from the original crystalline environment, the chain sheet models often exhibit large structural deviations from the initial structures and a large amount of computational resources is required for the optimizations to reach convergence. However, this has been found to be worthwhile because the present results shed light on the unusual aspects of native cellulose crystal packing features. The twist of the flat chain sheet models, which was first reported in our previous DFT studies,[45,46] is reproduced for both the cellulose Iα (110) and Iβ (100) chain sheet models by DFT calculations using various functionals and basis sets. Conley et al. attributed the fiber twist of cellulose crystals to the intrinsic propensity of a cellulose chain to deviate from the flat ribbon conformation of exact twofold helical symmetry.[91] In contrast, the present DFT calculations reveal that the cellulose chains essentially retain the twofold helical symmetry. The higher the level used for description of the hydrogen bonding and van der Waals interactions, the more symmetrical and larger the sheet twist of the chain sheets and the lower the binding energy, suggesting that the intermolecular interactions may cooperatively contribute to the chain sheet twist. We should emphasize that the inherent twisting of the chain sheet is likely to be a main cause of the parent crystalline fiber twist. In the stacked chain sheet models, a slight difference in chain overlap results in the optimized structures either retaining the initial chain structures or losing the sheet appearances. The cellulose Iα (010) and Iβ (11̅0) chain sheet models with more overlap retain the initial chain sheet structures, whereas the initial sheet appearance is lost in the DFT-optimized Iα (100) and Iβ (110) chain sheet models. It should be noted that the cellulose chain array constituting the cellulose Iβ (11̅0) chain sheet has been found to be conserved during crystalline conversions of the real crystalline phases.[70−72] Similarly, the cellulose IIII (11̅0) chain sheet has also been suggested to be “a conserved chain sheet” in crystalline conversion from the cellulose IIII phase to the Iβ phase.[23,45,46] We thus conclude that the chain sheet twist and the self-support of the conserved chain sheet are inherited from the parent crystals. As a future DFT study, the chain sheet models can be reunited to build a three-dimensional supercell model. For example, the two-layered flat chain sheet models are expected to twist with smaller extent owing to the intersheet interactions. Similarly, the two retained stacked chain sheets are combined to estimate collective interactions of the intermolecular hydrogen bonds.

In the present study, we used the thermochemical properties (herein, the binding energy) to perform a benchmark test of different levels of theory rather than for assessing the stability of the crystal packing. Inclusion of diffuse functions in the basis set appropriately describes the hydrogen bonding interactions in the flat chain sheet models. The results of 6-31+G(d,p) optimization are independent of the choice of the DFT functional and mostly give similar results in terms of the chain sheet geometry, whereas CAM-B3LYP and M06-2X calculations give lower values of the binding energy. In contrast, the use of diffuse functions is likely to decrease the effective convergence of the self-consistent field calculation and geometry optimization cycles (the 6-31+G(d,p) optimizations of the stacked chain sheets required a large number of minimization cycles only to result in convergence failure). It should be noted that M06-2X calculations coupled with basis sets without diffuse functions tend to result in less symmetrical twists in the flat chain sheet models. Single-point calculations with the 6-311+G(d,p) basis set slightly improve the binding energy values of the optimized structure over those obtained using the basis set without diffuse functions, but they are accompanied by relatively large BSSE correction errors. The above factors are more obvious for DFT optimization of the cellulose Iβ (11̅0) stacked chain sheet model, where van der Waals interactions (dispersion forces) are dominant between the cellulose chains. Single-point calculations with the 6-311+G(d,p) basis set greatly improve the final binding energy values. DFT calculations with the M06-2X functional, which has been suggested to well describe the dispersion forces and hydrogen bonds in the cellulose crystal models,[37] give much lower binding energies than those with the CAM-B3LYP functional, whereas the binding energies obtained with both functionals are essentially the same for DFT calculations of the flat chain sheet models.

It seems likely that the energy landscapes for the present chain sheet models consist of a shallow potential well including multiple local minima. Especially, optimizations of the flat chain sheet models, while displaying exclusively a right-handed twist, may have been partly trapped in one of the local minima. Therefore, we were unwilling to describe the further details of the optimized structures such as the hydrogen bonding geometries and the side group orientations. To discuss structural details of the native cellulose crystal of a fiber twist, the DFT results of the chain sheet and three-dimensional models should be integrated with those of thorough and systematic calculations of the related small molecules along with the experimental results.[92]