Carbohydrate Research 388 (2014) 61–66
Contents lists available at ScienceDirect
Carbohydrate Research journal homepage: www.elsevier.com/locate/carres
Investigation of the structure and interaction of cellulose triacetate I crystal using ab initio calculations Takanori Kobayashi, Daichi Hayakawa, Tegshjargal Khishigjargal, Kazuyoshi Ueda ⇑ Department of Chemistry, Graduate School of Engineering, Yokohama National University, 79-5 Tokiwadai, Hodogaya-Ku, Yokohama 240-8501, Japan
a r t i c l e
i n f o
Article history: Received 26 December 2013 Received in revised form 10 February 2014 Accepted 11 February 2014 Available online 21 February 2014 Keywords: Cellulose triacetate CTA I crystal Dissolution CH/O interaction AIM analysis Density functional calculation
a b s t r a c t The crystal structure of cellulose triacetate I (CTA I) was investigated using first-principles density functional theory (DFT) calculations. The results are in good agreement with the experimental structure obtained by Sikorski et al. when performing the calculation with inclusion of the dispersion correction. However, the cell parameters calculated with inclusion of the dispersion correction are slightly smaller than those experimentally obtained, especially along the a-axis. This smaller cell parameter could be reasonably explained by considering thermal expansion effects, since optimization with the density functional calculation gives the structure without inclusion of thermal effects. The atoms-in-molecules (AIM) theory is also employed to identify and characterize interatomic interactions in the CTA I crystal. CH/O interactions sites are shown to exist in the crystal structure of CTA I. Moreover, CH/O interactions are considered the main interactions in operation to maintain the crystal structure of CTA I. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Cellulose triacetate (CTA) is one of the more common industrially manufactured cellulose derivatives in use for many years in the fiber and textile industries and currently used in a variety of commercial products such as films and fibers. The applications of CTA have seen significant expansion not only to these mass production fields, but also to high-technology fields. For example, CTA has been recognized as a powerful chiral polymeric sorbent for chromatographic separation of enantiomers, and its industrial application is widespread in chemical and pharmaceutical fields. As for its structure, CTA is known to adopt two main types of crystalline polymorphs, CTA I and CTA II, which are derived from the corresponding cellulose I and II, respectively.1,2 According to Sprague et al.,1 CTA I can only be obtained by heterogeneous acetylation of native cellulose I fibers. Numerous structural studies of CTA I have been conducted over the years; the crystal structure of CTA I has been investigated by X-ray diffraction, NMR and IR spectroscopic techniques, and computational theoretical calculation methods, and several unit cell models have been proposed.3–5 Although the crystal structure of CTA I has been refined according to these models, its exact determination is still a matter of controversy. Recently, Sikorski et al.6 proposed a further refined crystal structure of CTA I with cell dimensions of a = 0.5939 nm, ⇑ Corresponding author. Fax: +81 45 339 3945. E-mail address: [email protected] (K. Ueda). http://dx.doi.org/10.1016/j.carres.2014.02.015 0008-6215/Ó 2014 Elsevier Ltd. All rights reserved.
b = 1.1431 nm, c (chain axis) = 1.046 nm, and c = 95.4 from X-ray diffraction experiments using a highly oriented fiber sample of CTA I. This model was in good agreement with the observed X-ray diffraction intensities, with a crystallographic reliability index of R = 0.224. Moreover, cross polarization–magic angle spinning (CP/MAS) 13C NMR spectroscopic studies demonstrated that each ring carbon atom of CTA I only shows a singlet resonance;7,8 this indicates that the asymmetric unit of CTA I should be composed of only a single CTA chemical unit. The crystal structure determined by Sikorski et al.6 can provide an explanation for these NMR results. Recent reports have described crystal structure determinations of cellulose using ab initio quantum chemical calculations.9,10 These reports indicated that the polymer crystal structure could be calculated with a reliable accuracy. In our previous report,11 we applied a similar approach to investigate the crystal structures of CTA I using the coordinates obtained by Sikorski et al.6 as an initial geometry. Sikorski et al.6 used a molecular mechanics method to optimize their model with the condition of cell size fixed to the experimental value in the refinement process of their CTA I crystal structure. In our calculation, we performed a full optimization of the crystal structure, which included relaxation of the cell parameters using density functional theory (DFT) calculations; the results of our study showed that the ab initio calculation supports their crystal structure. Here, we present the full work for the optimization of the CTA I crystal structure. The effects of crystal thermal expansion are discussed in detail along with a comparison of the
62
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
experimental data. The interactions involved in the formation of crystal structure of CTA I are also evaluated. 2. Computational details A schematic illustration and the nomenclature of the CTA chain are shown in Figure 1. Initial parameters of the unit cell and the atomic coordinates of CTA I were used as determined by Sikorski et al.6 The variable-cell relaxation (abbreviated later as vc-relax), which optimizes both the cell parameters and internal geometry optimization with monoclinic P1 space group, was performed with the Quantum ESPRESSO program package.12 We used the ultra-soft pseudopotentials (US-PP) plane wave basis sets, which were downloaded from the Quantum ESPRESSO Homepage.13 The convergences of energies and the structure were tested with different cutoff energy values of 50, 70, 100, and 150 Ry. Monkhorst–Pack k-point grids14 were tested with 2 2 2, 3 3 3, and 4 4 4. The Perdew–Burke–Ernzerhof (PBE) gradient-corrected density functional15 was used. The long-range van der Waals (vdW)-type correction term was included by using the DFT-D approach (PBE-D).16 NMR chemical shift calculations were performed using the GIPAW method17 as a single-point calculation on both the crystal structure obtained experimentally by Sikorski et al. and the theoretically optimized structures using the PBE and PBE-D methods. After the test calculations with different k-point spacings and cutoff energies, a cutoff energy (Ecut) of 70 Ry and the 3 3 3 Monkhorst–Pack k-point grid was chosen for all chemical shift calculations. The GIPAW module is available in the Quantum ESPRESSO program package. Data visualization was done using VMD1.8.7.18 Input files for atoms-in-molecules (AIM) analysis19 were generated by Gaussian 0920 and analyzed using the AIM 2000 package.21 The wave functions generated from the quantum chemical calculations were calculated using the M062X/631+G(d,p) basis set and used in the AIM calculation. AIM analysis was used to calculate the electron density, q, and the Laplacian of the electron density, Dq, at the bond critical points. 3. Results and discussion 3.1. Optimized conformation of CTA I Before discussing the crystal structure in detail, it would be important to check the convergence of optimized cell parameters
according to the selection of the parameters in the calculation such as the functional, kinetic energy cutoff, k-points, etc. The optimized cell parameters calculated using different kinetic energy cutoff values with the PBE functional are shown in Table 1. As can be seen in Table 1, the cell parameters of b, c, and c were very similar, even with variations in the Ecut value. However, cell parameter a calculated using Ecut = 50 Ry showed a slightly smaller value compared to others obtained with larger Ecut values. Above an Ecut value of 70 Ry, cell parameter a showed almost same value irrespective of the Ecut value increments, that is, this cell parameter can be considered to be converged at Ecut = 70 Ry for this crystal. A similar calculation was performed on the cell parameters with the PBE-D functional, the results of which are shown in Table 2. A similar dependence of the cell parameter on the Ecut value was observed. That is, cell parameters of b, c, and c were similar to each other even if the Ecut value varied, and cell parameter a calculated using Ecut = 50 Ry showed a slightly larger value compared to others obtained with larger Ecut values. These results again suggest that Ecut = 70 Ry is a maximum economical choice for the calculation of cell parameters of this crystal. The Monkhorst–Pack k-point grids11 selected for testing were 2 2 2, 3 3 3, and 4 4 4. The results calculated using PBE and PBE-D functional with the Ecut value of 70 Ry are shown in Tables 3 and 4, respectively. All cell parameters (a, b, c, and c) did not show any dependence on the Monkhorst–Pack k-point grids selected. As the model system used in this calculation is large, the selection of the k-point grid does not affect the calculation results. According to this result, the 2 2 2 Monkhorst–Pack k-point grid was used in subsequent crystal structure optimization.
Table 1 Optimized cell parameters of CTA I obtained from DFT calculations with PBE functional at several values of kinetic energy cutoff, Ecut Ecut/Ry
a (Å) b (Å) c (Å) c (°)
50
70
100
150
6.19 11.70 10.61 96.34
6.45 11.78 10.64 95.57
6.52 11.72 10.65 94.07
6.45 11.77 10.64 95.58
The values of k-points are 2 2 2.
Table 2 Optimized cell parameters of CTA I obtained from DFT calculations with PBE-D functional at several values of kinetic energy cutoff, Ecut Ecut/Ry
a (Å) b (Å) c (Å) c (°)
50
70
100
150
5.65 11.44 10.53 93.45
5.61 11.42 10.55 94.77
5.60 11.40 11.55 94.77
5.61 11.40 10.55 94.79
The values of k-points are 2 2 2.
Table 3 Optimized cell parameters of CTA I obtained from DFT calculations with PBE functional at several k-points with kinetic energy cutoff, Ecut = 70 Ry k-Points
Figure 1. Nomenclature of cellulose triacetate.
a (Å) b (Å) c (Å) c (°)
222
333
444
6.45 11.78 10.64 95.57
6.44 11.77 10.63 95.52
6.40 11.77 10.63 95.52
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66 Table 4 Optimized cell parameters of CTA I obtained from DFT calculations with PBE-D functional at several k-points with kinetic energy cutoff, Ecut = 70 Ry k-Points
a (Å) b (Å) c (Å) c (°)
222
333
444
5.61 11.42 10.55 94.77
5.61 11.40 10.55 94.74
5.61 11.40 10.55 94.74
The cell parameters obtained by using Ecut = 70 Ry with PBE and PBE-D methods are listed again and compared with the experimental values in Table 5, in addition to some torsion angles obtained from their structures. The calculated cell parameters of b, c, and c were almost the same compared to the experimental values. On the other hand, the calculated value of a (6.45 Å) with the PBE functional was considerably larger than the experimental value of 5.939 Å. However, the calculated result with the inclusion of a dispersion correction (PBE-D) gave a value of cell parameter a considerably shrunk from 6.45 Å to 5.61 Å; the value was slightly smaller but much closer to the experimentally obtained value. As can be seen in Figure 2, the direction of cell parameter a is in the inter-sheet direction where two CTA chains are stacked with a face-to-face orientation of their glucopyranose planes. In such a
Table 5 Comparison between experimentally and theoretically calculated cell parameters and some structure values
a (Å) b (Å) c (Å) c (°) / (°) w (°) s (°) v5 (°) v6 (°)
PBE
PBE-D
Exp.6
6.45 11.78 10.64 95.57 99.3 145.4 119.2 59.6 162.7
5.61 11.42 10.55 94.77 99.6 145.4 118.0 61.3 158.2
5.939 11.431 10.460 95.20 100.0 144.7 118.8 57.9 165.6
The calculated values by PBE and PBE-D with kinetic energy cutoff at 70 Ry were shown. The values of the main chain torsion angles (/, w) are defined as O5-C1-O40 C40 and C1-O40 -C40 -C50 . k-Points 2 2 2. s is an angle of C1-O40 -C40 , and v5, v6 are torsion angles at C4-C5-C6-O6 and C5-C6-O6-CA6, respectively.
63
chain arrangement, weak CH/O and/or vdW-type interactions between sheets would be expected;22 that is, the dispersion effect would be effective in the direction of cell parameter a. As a result, the calculation with the PBE functional underestimates such a dispersion effect, which would be properly corrected with inclusion of the empirical dispersion correction of the PBE-D functional. A similar effect in the cellulose crystal was previously observed, and the weak interaction has been discussed by Li et al.10 Compared to the cell parameter a, the extent of the shrinkage in cell parameters b and c were small when comparing results with and without the dispersion term. As is seen in Figure 2, cell parameter b corresponds to the direction of the CTA I orientation where CH/O and/or vdW-type interactions would be weak between the neighboring cellulose chains. Similarly, cell parameter c should not have space for variance, since it corresponds to the direction of the cellulose chains. Therefore, it is a reasonable result that the dispersion does not affect the cell parameters in both cases. Some CTA I torsion angles after the calculation were compared with experimental values (Table 5), and were all found to be in good agreement with each other; the values of v5 and v6 indicate that the side chain conformation at the 6-position adopts the gg (gauche–gauche) conformation after the minimization. 3.2. Thermal effects on cell size As was discussed in the above sections, the results obtained with the empirical dispersion correction suggested that all cell parameters are slightly smaller than the experimentally obtained values. However, it would be reasonable to consider that the thermal fluctuation of the crystal would expand its cell structure under experimental conditions, since the vc-relax calculation corresponds to the optimization of the crystal structure at 0 K. Wada and Hori experimentally observed the expansion of cell parameters of the CTA I crystal accompanied by an increment in temperature from 20 °C to 250 °C.23 The thermal expansion coefficient of the a-axis can be estimated as 12.5 105 nm/°C from their data. As the calculated value of 5.61 Å is that found at 0 K, we should add 0.366 Å to estimate the cell size at room temperature (20 °C); the estimated value of 5.976 Å at 20 °C is very close to the experimental value of 5.939 Å. Similarly, the thermal expansion coefficient of the b-axis was given as 0.3 105 nm/°C. Correction of the calculated value of 11.42 Å by addition of 0.00879 Å gave a value of 11.429 Å. This value again is in good accordance with the experimental value of 11.431 Å. Moreover, the thermal expansion
Figure 2. Crystal structure of CTA I after optimization using DFT calculations. Views are orthogonal to the (a) ab- and (b) bc- planes, respectively.
64
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
coefficient of the c-axis was given as 1.9 105 nm/°C. A correction of 0.056 Å was added to the calculated value of 10.55 Å to give the corrected value of 10.494 Å, which again is in good accordance with the experimental value of 10.460 Å. Although Wada and Hori obtained the thermal expansion coefficient of CTA I in the temperature range from 20 °C to 250 °C, this can be applied to the low-temperature range by linear extrapolation to a temperature of absolute zero. 3.3. Calculation of
13
C NMR chemical shifts
The 13C NMR chemical shifts of the CTA I crystal structure were next calculated, the results of which are summarized in Table 6. The chemical shifts of CTA I reported by Kono et al.8 are also shown in the same Table. In this Table, all carbon chemical shifts are normalized, with carbon C6 as a reference using the experimental value of 60.80 ppm. Although the calculated values are a little higher than the experimental values,8 the agreement between them can be considered to be acceptable. A comparison between the experimental and calculated values is shown in Figure 3. A linear relationship can be observed between the experimental and the calculated values with the slope close to one. This again indicates that the calculated results are in good agreement with the experimental values. Experimental results of CP/MAS 13C NMR spectroscopic analysis demonstrated that each ring carbon atom of CTA I only shows a singlet resonance,8 indicating that the asymmetric unit of CTA I should consist of only a single CTA chemical unit. The crystal structure determined by Sikorski et al.6 provides an explanation of these NMR results. However, NMR experiments showed that the chemical shifts of the carbonyl carbons of acetyl groups were split into doublet signals with more than 2 ppm difference. The calculated chemical shifts reproduced the doublet signals of the acetyl groups; specifically, carbonyl carbons at the 2, 3, and 6 positions were split into two groups with values around 175 and 178 ppm, although the initial structure of CTA I shows almost the same chemical shift values at around 173 ppm. In contrast, the experimental chemical shift values of the methyl carbons of acetyl groups were also split into doublet signals, but with a difference in signals of only ca. 1 ppm. The calculated values of methyl carbons at the 2, 3, and 6 positions were also split in the range between 18.4 and 19.7 ppm, indicating that the calculated structure using PBE-D can adequately explain the experimental NMR values for both acetyl carbonyl and methyl carbons.
Figure 3. Relationship between the calculated 13C NMR chemical shift and experimental results.8 Horizontal line shows the experimental chemical shift values and the vertical shows the calculated values. Calculated chemical shift values were evaluated based on the conformations experimentally obtained from X-ray structure analysis by Sikorski et al.6 (s), and PBE-D optimized structure (h). Calculated values were normalized with C6 carbon as a reference (experimental value, 60.80 ppm).
Table 7 AIM analysis of CTA I crystal structure, which was obtained from DFT calculations with PDE-D functional Type
Atom—atom
Distance [Å]
q [e/Å3]
Dq [e/Å5]
I
C1-H1—OA2 C3-H3—OA2 C5-H5—OA2
2.18 2.52 2.71
0.1084 0.0522 0.0449
1.3290 0.7210 0.5515
II III
CM3-HM3—OA3 CM3-HM3—O4
2.66 2.69
0.0421 0.0417
0.5674 0.5637
IV
CM6-HM6—OA6
2.35
0.0723
0.9380
V
CM2-HM2—OA6
2.47
0.0729
0.9443
VI
CM6-HM6—OA3
2.52
0.0575
0.8155
Interacting atom pairs and their distance, electron density q and laplacian of electron density Dq, at the critical points for intermolecular interactions are listed.
3.4. Analysis of the interaction between the CTA chains in the crystal Table 6 Calculated
13
C NMR chemical shifts (ppm) of the asymmetric unit of CTA crystal 8
Initial
PBE-D
Experimental
C1 C2 C3 C4 C5 C6
105.29 76.07 75.80 79.79 71.21 60.80
107.41 76.62 76.03 83.68 72.35 60.80
101.60 74.50 74.50 78.60 71.00 60.80
CA2 CA3 CA6
173.70 172.31 172.78
177.70 175.15 178.29
171.4, 169.10
CM2 CM3 CM6
20.32 18.97 17.76
19.69 18.38 17.93
21.2, 20.1
‘PBE-D’ shows the chemical shifts calculated for the optimized structure using PBE-D method and ‘Initial’ is the single point calculation using the X-ray structure obtained by Sikorski et al. k-Points are 3 3 3. Experimental data8 of methyl group (CA) and carbonyl group (CM) shows two values of the doublet chemical shifts.
A quantum chemical analysis of the interaction between the sheets of CTA chains in the crystal structure of CTA I has not yet been reported. We applied AIM (atoms-in-molecules) analysis to the optimized structure from the DFT-D calculation and elucidated the mechanism of the interactions among the CTA chains. The positions of possible interacting sites in the CTA crystal were roughly estimated by evaluating the distances between the atoms belonging to different CTA chains. The possible interacting atom pairs were selected having interatomic distances shorter than 2.7 Å; these pairs are listed in Table 7. The CTA chains can interact with one another in three types of orientations in the crystal structure of CTA I, as is shown in Figure 4. In these orientation groups, the six possible interacting atom pairs are located and labeled from I to VI. Their positions in the CTA I crystal structures are shown schematically in the same Figure. All of the identified pairs were found to form between hydrogen atoms of CH groups and oxygen atoms, strongly implying the existence of the so-called CH/O-type weak hydrogen-bond interactions between these atom pairs;
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
65
Figure 4. Three types of interacting CTA chain pairs in the crystal structure of CTA I (surrounded by dotted lines in the upper part of the figure). Six interacting atom pairs (I to VI), which included the three CTA chain pairs, are shown in the lower part of the picture with an exaggerated representation of the interacting atom pairs.
however, if these interacting forces are actually in effect among these CTA chains, any evidence of the existence of these interactions has not yet been reported. Proof of the existence of such interactions could be obtained using AIM analysis, which can evaluate the through-space bond formation between the atoms using the electron-density distribution obtained by quantum chemical calculations using the Gaussian program.21 However, the Gaussian calculation is difficult to apply to as large a molecular system as the CTA I crystal system with infinite repeating units. Therefore, three pairs of glucose units, which are shown in Figure 4, were cut out from the optimized structure of CTA I, and electron density calculations were performed on these dimer pairs. Subsequently, the AIM analysis was performed on these systems using the electron density obtained in the above calculation, the results of which are shown in Table 7. AIM analysis can calculate the critical bonding points from the electron density landscape between the atom pairs. Electron densities and the Laplacian of electron densities at the critical points between the atom pairs are listed in Table 7. All the atom pairs showed both positive electron densities and positive Laplacian of electron densities at the critical points. The electron density at the critical point indicates the strength of the bond and the sign of Laplacian indicates the types of the interactions;24,25 specifically, negative values of the Laplacian indicate a covalent bond nature and positive values indicate closed-shell
interactions as found in ionic bonds, hydrogen bonds, and vdW-type interactions. Koch et al.24 and Parthasarathi et al.26 suggested that the standard value of the CH/O interaction should be around 103 au. The positive values obtained for CTA I in the present study show the same order with their criteria, indicating the existence of the weak hydrogen-bond-type interaction between the CH and O atom pairs listed in Table 7. It should be emphasized that no other atom pairs, whose values are larger than those listed in Table 7, were found in the AIM analysis. Among the atom pairs in the Table 7, the type I atom pair showed the largest value of electron density at the critical point, which indicates the importance of the interaction between these atoms on the formation of the CTA I crystal structure; the carbonyl oxygen interacts with the ring CH groups of the glucose ring, which would be the largest interaction site in the crystal structure of CTA I. 4. Conclusions We have investigated the crystal structure of CTA I using density functional theory calculations. Without the dispersion correction, the cell parameters are overestimated; however, with inclusion of the dispersion term, all cell lengths become smaller than the experimental values. Especially, cell parameter a is considerably shrunk. As CTA chains orient in parallel with face-to-face
66
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
interactions of their hydrophobic sides in the direction of cell parameter a, there exists a weak interaction such as CH/O (hydrogen bonding). In such a case, inclusion of the dispersion term is important in the calculation. After intensive investigations of the choice of computational parameter, optimized geometry and the interaction mechanisms present in the crystal structure were investigated. AIM analysis was performed on the optimized structure from the DFT-D calculation the existence of the interactions among the CTA chains was elucidated. Positive electron densities between some atom pairs were observed, which indicate the existence of non-bonded interactions between the atom pairs. The analysis showed that the carbonyl oxygen interacts with ring CH groups of the glucose ring, which would be the largest interaction site in the crystal structure of CTA I. This work clearly shows that the first-principle density functional calculation is a powerful method to investigate the structure and the interactions of the CTA crystal. We are now trying to investigate the crystal structure of CTA II whose crystal structure is still under controversy. Acknowledgement The authors wish to thank the Research Center for Computational Science, Okazaki, Japan for the use of their computer to perform part of the calculation. We greatly acknowledge Prof. M. Wada at the Tokyo University for sending us their structural data of CTA I. References 1. Sprague, B. S.; Riley, J. L.; Noether, H. D. Text. Res. J. 1958, 28, 275–287. 2. Roche, E.; Chanzy, H.; Boudeulle, M.; Marchessault, R. H.; Sundararajan, P. Macromolecules 1978, 11, 86–94. 3. Stipanovic, A. J.; Sarko, A. Polymer 1978, 19, 3–8. 4. VanderHart, D. L.; Hyatt, J. A.; Atalla, R. H.; Tirumalai, V. C. Macromolecules 1996, 29, 730–739. 5. Wolf, R. M.; Francotte, E.; Glasser, L.; Simon, I.; Scheraga, H. A. Macromolecules 1992, 25, 709–720. 6. Sikorski, P.; Wada, M.; Heux, L.; Shintani, H.; Stokke, B. T. Macromolecules 2004, 37, 4547–4553.
7. Kono, H.; Erata, T.; Takai, M. J. Am. Chem. Soc. 2002, 124, 7512–7518. 8. Kono, H.; Numata, Y.; Erata, T.; Takai, M. Polymer 2004, 45, 2843–2852. 9. Bucˇko, T.; Tunega, D.; Ángyán, J. G.; Hafner, J. J. Phys. Chem. A 2011, 115, 10097– 10105. 10. Li, Y.; Lin, M.; Davenport, J. W. J. Phys. Chem. C 2011, 115, 11533–11539. 11. Kobayashi, T.; Hayakawa, D.; Khishigjargal, T.; Ueda, K. MRS Online Proc. Lib. 2012, 1470. mrss12-1470-xx03-05. 12. Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. J. Phys.: Condens. Matter 2009, 21, 395502. 13. http://www.quantum-espresso.org/index.php. 14. Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13, 5188–5192. 15. Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. 16. Grimme, S. J. J. Comput. Chem. 2006, 27, 1787–1799. 17. Pickard, C. J.; Mauri, F. Phys. Rev. B 2001, 63, 245101. 18. Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38. 19. Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, 1990. 20. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T., ; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision, A.02 ed.; Gaussian: Wallingford CT, 2009. 21. AIM2000, Designed by F. Biegler-Konig, University of Applied Sciences, Bielefeld, Germany, http://www.AIM2000.de. 22. Hayakawa, D.; Ueda, K.; Yamane, C.; Miyamoto, H.; Horii, F. Carbohydr. Res. 2011, 346, 2940–2947. 23. Wada, M.; Hori, R. J. Polym. Sci., Part B: Polym. Phys. 2009, 47, 517–523. 24. Koch, U.; Popelier, P. L. A. J. Phys. Chem. 1995, 99, 9747–9754. 25. Kobayashi, T.; Kohno, Y.; Takayanagi, T.; Seki, K.; Ueda, K. Comput. Theor. Chem. 2012, 991, 48–55. 26. Parthasarathi, R.; Bellesia, G.; Chundawat, S. P. S.; Dale, B. E.; Langan, P.; Gnanakaran, S. J. Phys. Chem. A 2011, 115, 14191–14202.
Contents lists available at ScienceDirect
Carbohydrate Research journal homepage: www.elsevier.com/locate/carres
Investigation of the structure and interaction of cellulose triacetate I crystal using ab initio calculations Takanori Kobayashi, Daichi Hayakawa, Tegshjargal Khishigjargal, Kazuyoshi Ueda ⇑ Department of Chemistry, Graduate School of Engineering, Yokohama National University, 79-5 Tokiwadai, Hodogaya-Ku, Yokohama 240-8501, Japan
a r t i c l e
i n f o
Article history: Received 26 December 2013 Received in revised form 10 February 2014 Accepted 11 February 2014 Available online 21 February 2014 Keywords: Cellulose triacetate CTA I crystal Dissolution CH/O interaction AIM analysis Density functional calculation
a b s t r a c t The crystal structure of cellulose triacetate I (CTA I) was investigated using first-principles density functional theory (DFT) calculations. The results are in good agreement with the experimental structure obtained by Sikorski et al. when performing the calculation with inclusion of the dispersion correction. However, the cell parameters calculated with inclusion of the dispersion correction are slightly smaller than those experimentally obtained, especially along the a-axis. This smaller cell parameter could be reasonably explained by considering thermal expansion effects, since optimization with the density functional calculation gives the structure without inclusion of thermal effects. The atoms-in-molecules (AIM) theory is also employed to identify and characterize interatomic interactions in the CTA I crystal. CH/O interactions sites are shown to exist in the crystal structure of CTA I. Moreover, CH/O interactions are considered the main interactions in operation to maintain the crystal structure of CTA I. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Cellulose triacetate (CTA) is one of the more common industrially manufactured cellulose derivatives in use for many years in the fiber and textile industries and currently used in a variety of commercial products such as films and fibers. The applications of CTA have seen significant expansion not only to these mass production fields, but also to high-technology fields. For example, CTA has been recognized as a powerful chiral polymeric sorbent for chromatographic separation of enantiomers, and its industrial application is widespread in chemical and pharmaceutical fields. As for its structure, CTA is known to adopt two main types of crystalline polymorphs, CTA I and CTA II, which are derived from the corresponding cellulose I and II, respectively.1,2 According to Sprague et al.,1 CTA I can only be obtained by heterogeneous acetylation of native cellulose I fibers. Numerous structural studies of CTA I have been conducted over the years; the crystal structure of CTA I has been investigated by X-ray diffraction, NMR and IR spectroscopic techniques, and computational theoretical calculation methods, and several unit cell models have been proposed.3–5 Although the crystal structure of CTA I has been refined according to these models, its exact determination is still a matter of controversy. Recently, Sikorski et al.6 proposed a further refined crystal structure of CTA I with cell dimensions of a = 0.5939 nm, ⇑ Corresponding author. Fax: +81 45 339 3945. E-mail address: [email protected] (K. Ueda). http://dx.doi.org/10.1016/j.carres.2014.02.015 0008-6215/Ó 2014 Elsevier Ltd. All rights reserved.
b = 1.1431 nm, c (chain axis) = 1.046 nm, and c = 95.4 from X-ray diffraction experiments using a highly oriented fiber sample of CTA I. This model was in good agreement with the observed X-ray diffraction intensities, with a crystallographic reliability index of R = 0.224. Moreover, cross polarization–magic angle spinning (CP/MAS) 13C NMR spectroscopic studies demonstrated that each ring carbon atom of CTA I only shows a singlet resonance;7,8 this indicates that the asymmetric unit of CTA I should be composed of only a single CTA chemical unit. The crystal structure determined by Sikorski et al.6 can provide an explanation for these NMR results. Recent reports have described crystal structure determinations of cellulose using ab initio quantum chemical calculations.9,10 These reports indicated that the polymer crystal structure could be calculated with a reliable accuracy. In our previous report,11 we applied a similar approach to investigate the crystal structures of CTA I using the coordinates obtained by Sikorski et al.6 as an initial geometry. Sikorski et al.6 used a molecular mechanics method to optimize their model with the condition of cell size fixed to the experimental value in the refinement process of their CTA I crystal structure. In our calculation, we performed a full optimization of the crystal structure, which included relaxation of the cell parameters using density functional theory (DFT) calculations; the results of our study showed that the ab initio calculation supports their crystal structure. Here, we present the full work for the optimization of the CTA I crystal structure. The effects of crystal thermal expansion are discussed in detail along with a comparison of the
62
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
experimental data. The interactions involved in the formation of crystal structure of CTA I are also evaluated. 2. Computational details A schematic illustration and the nomenclature of the CTA chain are shown in Figure 1. Initial parameters of the unit cell and the atomic coordinates of CTA I were used as determined by Sikorski et al.6 The variable-cell relaxation (abbreviated later as vc-relax), which optimizes both the cell parameters and internal geometry optimization with monoclinic P1 space group, was performed with the Quantum ESPRESSO program package.12 We used the ultra-soft pseudopotentials (US-PP) plane wave basis sets, which were downloaded from the Quantum ESPRESSO Homepage.13 The convergences of energies and the structure were tested with different cutoff energy values of 50, 70, 100, and 150 Ry. Monkhorst–Pack k-point grids14 were tested with 2 2 2, 3 3 3, and 4 4 4. The Perdew–Burke–Ernzerhof (PBE) gradient-corrected density functional15 was used. The long-range van der Waals (vdW)-type correction term was included by using the DFT-D approach (PBE-D).16 NMR chemical shift calculations were performed using the GIPAW method17 as a single-point calculation on both the crystal structure obtained experimentally by Sikorski et al. and the theoretically optimized structures using the PBE and PBE-D methods. After the test calculations with different k-point spacings and cutoff energies, a cutoff energy (Ecut) of 70 Ry and the 3 3 3 Monkhorst–Pack k-point grid was chosen for all chemical shift calculations. The GIPAW module is available in the Quantum ESPRESSO program package. Data visualization was done using VMD1.8.7.18 Input files for atoms-in-molecules (AIM) analysis19 were generated by Gaussian 0920 and analyzed using the AIM 2000 package.21 The wave functions generated from the quantum chemical calculations were calculated using the M062X/631+G(d,p) basis set and used in the AIM calculation. AIM analysis was used to calculate the electron density, q, and the Laplacian of the electron density, Dq, at the bond critical points. 3. Results and discussion 3.1. Optimized conformation of CTA I Before discussing the crystal structure in detail, it would be important to check the convergence of optimized cell parameters
according to the selection of the parameters in the calculation such as the functional, kinetic energy cutoff, k-points, etc. The optimized cell parameters calculated using different kinetic energy cutoff values with the PBE functional are shown in Table 1. As can be seen in Table 1, the cell parameters of b, c, and c were very similar, even with variations in the Ecut value. However, cell parameter a calculated using Ecut = 50 Ry showed a slightly smaller value compared to others obtained with larger Ecut values. Above an Ecut value of 70 Ry, cell parameter a showed almost same value irrespective of the Ecut value increments, that is, this cell parameter can be considered to be converged at Ecut = 70 Ry for this crystal. A similar calculation was performed on the cell parameters with the PBE-D functional, the results of which are shown in Table 2. A similar dependence of the cell parameter on the Ecut value was observed. That is, cell parameters of b, c, and c were similar to each other even if the Ecut value varied, and cell parameter a calculated using Ecut = 50 Ry showed a slightly larger value compared to others obtained with larger Ecut values. These results again suggest that Ecut = 70 Ry is a maximum economical choice for the calculation of cell parameters of this crystal. The Monkhorst–Pack k-point grids11 selected for testing were 2 2 2, 3 3 3, and 4 4 4. The results calculated using PBE and PBE-D functional with the Ecut value of 70 Ry are shown in Tables 3 and 4, respectively. All cell parameters (a, b, c, and c) did not show any dependence on the Monkhorst–Pack k-point grids selected. As the model system used in this calculation is large, the selection of the k-point grid does not affect the calculation results. According to this result, the 2 2 2 Monkhorst–Pack k-point grid was used in subsequent crystal structure optimization.
Table 1 Optimized cell parameters of CTA I obtained from DFT calculations with PBE functional at several values of kinetic energy cutoff, Ecut Ecut/Ry
a (Å) b (Å) c (Å) c (°)
50
70
100
150
6.19 11.70 10.61 96.34
6.45 11.78 10.64 95.57
6.52 11.72 10.65 94.07
6.45 11.77 10.64 95.58
The values of k-points are 2 2 2.
Table 2 Optimized cell parameters of CTA I obtained from DFT calculations with PBE-D functional at several values of kinetic energy cutoff, Ecut Ecut/Ry
a (Å) b (Å) c (Å) c (°)
50
70
100
150
5.65 11.44 10.53 93.45
5.61 11.42 10.55 94.77
5.60 11.40 11.55 94.77
5.61 11.40 10.55 94.79
The values of k-points are 2 2 2.
Table 3 Optimized cell parameters of CTA I obtained from DFT calculations with PBE functional at several k-points with kinetic energy cutoff, Ecut = 70 Ry k-Points
Figure 1. Nomenclature of cellulose triacetate.
a (Å) b (Å) c (Å) c (°)
222
333
444
6.45 11.78 10.64 95.57
6.44 11.77 10.63 95.52
6.40 11.77 10.63 95.52
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66 Table 4 Optimized cell parameters of CTA I obtained from DFT calculations with PBE-D functional at several k-points with kinetic energy cutoff, Ecut = 70 Ry k-Points
a (Å) b (Å) c (Å) c (°)
222
333
444
5.61 11.42 10.55 94.77
5.61 11.40 10.55 94.74
5.61 11.40 10.55 94.74
The cell parameters obtained by using Ecut = 70 Ry with PBE and PBE-D methods are listed again and compared with the experimental values in Table 5, in addition to some torsion angles obtained from their structures. The calculated cell parameters of b, c, and c were almost the same compared to the experimental values. On the other hand, the calculated value of a (6.45 Å) with the PBE functional was considerably larger than the experimental value of 5.939 Å. However, the calculated result with the inclusion of a dispersion correction (PBE-D) gave a value of cell parameter a considerably shrunk from 6.45 Å to 5.61 Å; the value was slightly smaller but much closer to the experimentally obtained value. As can be seen in Figure 2, the direction of cell parameter a is in the inter-sheet direction where two CTA chains are stacked with a face-to-face orientation of their glucopyranose planes. In such a
Table 5 Comparison between experimentally and theoretically calculated cell parameters and some structure values
a (Å) b (Å) c (Å) c (°) / (°) w (°) s (°) v5 (°) v6 (°)
PBE
PBE-D
Exp.6
6.45 11.78 10.64 95.57 99.3 145.4 119.2 59.6 162.7
5.61 11.42 10.55 94.77 99.6 145.4 118.0 61.3 158.2
5.939 11.431 10.460 95.20 100.0 144.7 118.8 57.9 165.6
The calculated values by PBE and PBE-D with kinetic energy cutoff at 70 Ry were shown. The values of the main chain torsion angles (/, w) are defined as O5-C1-O40 C40 and C1-O40 -C40 -C50 . k-Points 2 2 2. s is an angle of C1-O40 -C40 , and v5, v6 are torsion angles at C4-C5-C6-O6 and C5-C6-O6-CA6, respectively.
63
chain arrangement, weak CH/O and/or vdW-type interactions between sheets would be expected;22 that is, the dispersion effect would be effective in the direction of cell parameter a. As a result, the calculation with the PBE functional underestimates such a dispersion effect, which would be properly corrected with inclusion of the empirical dispersion correction of the PBE-D functional. A similar effect in the cellulose crystal was previously observed, and the weak interaction has been discussed by Li et al.10 Compared to the cell parameter a, the extent of the shrinkage in cell parameters b and c were small when comparing results with and without the dispersion term. As is seen in Figure 2, cell parameter b corresponds to the direction of the CTA I orientation where CH/O and/or vdW-type interactions would be weak between the neighboring cellulose chains. Similarly, cell parameter c should not have space for variance, since it corresponds to the direction of the cellulose chains. Therefore, it is a reasonable result that the dispersion does not affect the cell parameters in both cases. Some CTA I torsion angles after the calculation were compared with experimental values (Table 5), and were all found to be in good agreement with each other; the values of v5 and v6 indicate that the side chain conformation at the 6-position adopts the gg (gauche–gauche) conformation after the minimization. 3.2. Thermal effects on cell size As was discussed in the above sections, the results obtained with the empirical dispersion correction suggested that all cell parameters are slightly smaller than the experimentally obtained values. However, it would be reasonable to consider that the thermal fluctuation of the crystal would expand its cell structure under experimental conditions, since the vc-relax calculation corresponds to the optimization of the crystal structure at 0 K. Wada and Hori experimentally observed the expansion of cell parameters of the CTA I crystal accompanied by an increment in temperature from 20 °C to 250 °C.23 The thermal expansion coefficient of the a-axis can be estimated as 12.5 105 nm/°C from their data. As the calculated value of 5.61 Å is that found at 0 K, we should add 0.366 Å to estimate the cell size at room temperature (20 °C); the estimated value of 5.976 Å at 20 °C is very close to the experimental value of 5.939 Å. Similarly, the thermal expansion coefficient of the b-axis was given as 0.3 105 nm/°C. Correction of the calculated value of 11.42 Å by addition of 0.00879 Å gave a value of 11.429 Å. This value again is in good accordance with the experimental value of 11.431 Å. Moreover, the thermal expansion
Figure 2. Crystal structure of CTA I after optimization using DFT calculations. Views are orthogonal to the (a) ab- and (b) bc- planes, respectively.
64
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
coefficient of the c-axis was given as 1.9 105 nm/°C. A correction of 0.056 Å was added to the calculated value of 10.55 Å to give the corrected value of 10.494 Å, which again is in good accordance with the experimental value of 10.460 Å. Although Wada and Hori obtained the thermal expansion coefficient of CTA I in the temperature range from 20 °C to 250 °C, this can be applied to the low-temperature range by linear extrapolation to a temperature of absolute zero. 3.3. Calculation of
13
C NMR chemical shifts
The 13C NMR chemical shifts of the CTA I crystal structure were next calculated, the results of which are summarized in Table 6. The chemical shifts of CTA I reported by Kono et al.8 are also shown in the same Table. In this Table, all carbon chemical shifts are normalized, with carbon C6 as a reference using the experimental value of 60.80 ppm. Although the calculated values are a little higher than the experimental values,8 the agreement between them can be considered to be acceptable. A comparison between the experimental and calculated values is shown in Figure 3. A linear relationship can be observed between the experimental and the calculated values with the slope close to one. This again indicates that the calculated results are in good agreement with the experimental values. Experimental results of CP/MAS 13C NMR spectroscopic analysis demonstrated that each ring carbon atom of CTA I only shows a singlet resonance,8 indicating that the asymmetric unit of CTA I should consist of only a single CTA chemical unit. The crystal structure determined by Sikorski et al.6 provides an explanation of these NMR results. However, NMR experiments showed that the chemical shifts of the carbonyl carbons of acetyl groups were split into doublet signals with more than 2 ppm difference. The calculated chemical shifts reproduced the doublet signals of the acetyl groups; specifically, carbonyl carbons at the 2, 3, and 6 positions were split into two groups with values around 175 and 178 ppm, although the initial structure of CTA I shows almost the same chemical shift values at around 173 ppm. In contrast, the experimental chemical shift values of the methyl carbons of acetyl groups were also split into doublet signals, but with a difference in signals of only ca. 1 ppm. The calculated values of methyl carbons at the 2, 3, and 6 positions were also split in the range between 18.4 and 19.7 ppm, indicating that the calculated structure using PBE-D can adequately explain the experimental NMR values for both acetyl carbonyl and methyl carbons.
Figure 3. Relationship between the calculated 13C NMR chemical shift and experimental results.8 Horizontal line shows the experimental chemical shift values and the vertical shows the calculated values. Calculated chemical shift values were evaluated based on the conformations experimentally obtained from X-ray structure analysis by Sikorski et al.6 (s), and PBE-D optimized structure (h). Calculated values were normalized with C6 carbon as a reference (experimental value, 60.80 ppm).
Table 7 AIM analysis of CTA I crystal structure, which was obtained from DFT calculations with PDE-D functional Type
Atom—atom
Distance [Å]
q [e/Å3]
Dq [e/Å5]
I
C1-H1—OA2 C3-H3—OA2 C5-H5—OA2
2.18 2.52 2.71
0.1084 0.0522 0.0449
1.3290 0.7210 0.5515
II III
CM3-HM3—OA3 CM3-HM3—O4
2.66 2.69
0.0421 0.0417
0.5674 0.5637
IV
CM6-HM6—OA6
2.35
0.0723
0.9380
V
CM2-HM2—OA6
2.47
0.0729
0.9443
VI
CM6-HM6—OA3
2.52
0.0575
0.8155
Interacting atom pairs and their distance, electron density q and laplacian of electron density Dq, at the critical points for intermolecular interactions are listed.
3.4. Analysis of the interaction between the CTA chains in the crystal Table 6 Calculated
13
C NMR chemical shifts (ppm) of the asymmetric unit of CTA crystal 8
Initial
PBE-D
Experimental
C1 C2 C3 C4 C5 C6
105.29 76.07 75.80 79.79 71.21 60.80
107.41 76.62 76.03 83.68 72.35 60.80
101.60 74.50 74.50 78.60 71.00 60.80
CA2 CA3 CA6
173.70 172.31 172.78
177.70 175.15 178.29
171.4, 169.10
CM2 CM3 CM6
20.32 18.97 17.76
19.69 18.38 17.93
21.2, 20.1
‘PBE-D’ shows the chemical shifts calculated for the optimized structure using PBE-D method and ‘Initial’ is the single point calculation using the X-ray structure obtained by Sikorski et al. k-Points are 3 3 3. Experimental data8 of methyl group (CA) and carbonyl group (CM) shows two values of the doublet chemical shifts.
A quantum chemical analysis of the interaction between the sheets of CTA chains in the crystal structure of CTA I has not yet been reported. We applied AIM (atoms-in-molecules) analysis to the optimized structure from the DFT-D calculation and elucidated the mechanism of the interactions among the CTA chains. The positions of possible interacting sites in the CTA crystal were roughly estimated by evaluating the distances between the atoms belonging to different CTA chains. The possible interacting atom pairs were selected having interatomic distances shorter than 2.7 Å; these pairs are listed in Table 7. The CTA chains can interact with one another in three types of orientations in the crystal structure of CTA I, as is shown in Figure 4. In these orientation groups, the six possible interacting atom pairs are located and labeled from I to VI. Their positions in the CTA I crystal structures are shown schematically in the same Figure. All of the identified pairs were found to form between hydrogen atoms of CH groups and oxygen atoms, strongly implying the existence of the so-called CH/O-type weak hydrogen-bond interactions between these atom pairs;
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
65
Figure 4. Three types of interacting CTA chain pairs in the crystal structure of CTA I (surrounded by dotted lines in the upper part of the figure). Six interacting atom pairs (I to VI), which included the three CTA chain pairs, are shown in the lower part of the picture with an exaggerated representation of the interacting atom pairs.
however, if these interacting forces are actually in effect among these CTA chains, any evidence of the existence of these interactions has not yet been reported. Proof of the existence of such interactions could be obtained using AIM analysis, which can evaluate the through-space bond formation between the atoms using the electron-density distribution obtained by quantum chemical calculations using the Gaussian program.21 However, the Gaussian calculation is difficult to apply to as large a molecular system as the CTA I crystal system with infinite repeating units. Therefore, three pairs of glucose units, which are shown in Figure 4, were cut out from the optimized structure of CTA I, and electron density calculations were performed on these dimer pairs. Subsequently, the AIM analysis was performed on these systems using the electron density obtained in the above calculation, the results of which are shown in Table 7. AIM analysis can calculate the critical bonding points from the electron density landscape between the atom pairs. Electron densities and the Laplacian of electron densities at the critical points between the atom pairs are listed in Table 7. All the atom pairs showed both positive electron densities and positive Laplacian of electron densities at the critical points. The electron density at the critical point indicates the strength of the bond and the sign of Laplacian indicates the types of the interactions;24,25 specifically, negative values of the Laplacian indicate a covalent bond nature and positive values indicate closed-shell
interactions as found in ionic bonds, hydrogen bonds, and vdW-type interactions. Koch et al.24 and Parthasarathi et al.26 suggested that the standard value of the CH/O interaction should be around 103 au. The positive values obtained for CTA I in the present study show the same order with their criteria, indicating the existence of the weak hydrogen-bond-type interaction between the CH and O atom pairs listed in Table 7. It should be emphasized that no other atom pairs, whose values are larger than those listed in Table 7, were found in the AIM analysis. Among the atom pairs in the Table 7, the type I atom pair showed the largest value of electron density at the critical point, which indicates the importance of the interaction between these atoms on the formation of the CTA I crystal structure; the carbonyl oxygen interacts with the ring CH groups of the glucose ring, which would be the largest interaction site in the crystal structure of CTA I. 4. Conclusions We have investigated the crystal structure of CTA I using density functional theory calculations. Without the dispersion correction, the cell parameters are overestimated; however, with inclusion of the dispersion term, all cell lengths become smaller than the experimental values. Especially, cell parameter a is considerably shrunk. As CTA chains orient in parallel with face-to-face
66
T. Kobayashi et al. / Carbohydrate Research 388 (2014) 61–66
interactions of their hydrophobic sides in the direction of cell parameter a, there exists a weak interaction such as CH/O (hydrogen bonding). In such a case, inclusion of the dispersion term is important in the calculation. After intensive investigations of the choice of computational parameter, optimized geometry and the interaction mechanisms present in the crystal structure were investigated. AIM analysis was performed on the optimized structure from the DFT-D calculation the existence of the interactions among the CTA chains was elucidated. Positive electron densities between some atom pairs were observed, which indicate the existence of non-bonded interactions between the atom pairs. The analysis showed that the carbonyl oxygen interacts with ring CH groups of the glucose ring, which would be the largest interaction site in the crystal structure of CTA I. This work clearly shows that the first-principle density functional calculation is a powerful method to investigate the structure and the interactions of the CTA crystal. We are now trying to investigate the crystal structure of CTA II whose crystal structure is still under controversy. Acknowledgement The authors wish to thank the Research Center for Computational Science, Okazaki, Japan for the use of their computer to perform part of the calculation. We greatly acknowledge Prof. M. Wada at the Tokyo University for sending us their structural data of CTA I. References 1. Sprague, B. S.; Riley, J. L.; Noether, H. D. Text. Res. J. 1958, 28, 275–287. 2. Roche, E.; Chanzy, H.; Boudeulle, M.; Marchessault, R. H.; Sundararajan, P. Macromolecules 1978, 11, 86–94. 3. Stipanovic, A. J.; Sarko, A. Polymer 1978, 19, 3–8. 4. VanderHart, D. L.; Hyatt, J. A.; Atalla, R. H.; Tirumalai, V. C. Macromolecules 1996, 29, 730–739. 5. Wolf, R. M.; Francotte, E.; Glasser, L.; Simon, I.; Scheraga, H. A. Macromolecules 1992, 25, 709–720. 6. Sikorski, P.; Wada, M.; Heux, L.; Shintani, H.; Stokke, B. T. Macromolecules 2004, 37, 4547–4553.
7. Kono, H.; Erata, T.; Takai, M. J. Am. Chem. Soc. 2002, 124, 7512–7518. 8. Kono, H.; Numata, Y.; Erata, T.; Takai, M. Polymer 2004, 45, 2843–2852. 9. Bucˇko, T.; Tunega, D.; Ángyán, J. G.; Hafner, J. J. Phys. Chem. A 2011, 115, 10097– 10105. 10. Li, Y.; Lin, M.; Davenport, J. W. J. Phys. Chem. C 2011, 115, 11533–11539. 11. Kobayashi, T.; Hayakawa, D.; Khishigjargal, T.; Ueda, K. MRS Online Proc. Lib. 2012, 1470. mrss12-1470-xx03-05. 12. Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. J. Phys.: Condens. Matter 2009, 21, 395502. 13. http://www.quantum-espresso.org/index.php. 14. Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13, 5188–5192. 15. Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. 16. Grimme, S. J. J. Comput. Chem. 2006, 27, 1787–1799. 17. Pickard, C. J.; Mauri, F. Phys. Rev. B 2001, 63, 245101. 18. Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38. 19. Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, 1990. 20. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T., ; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision, A.02 ed.; Gaussian: Wallingford CT, 2009. 21. AIM2000, Designed by F. Biegler-Konig, University of Applied Sciences, Bielefeld, Germany, http://www.AIM2000.de. 22. Hayakawa, D.; Ueda, K.; Yamane, C.; Miyamoto, H.; Horii, F. Carbohydr. Res. 2011, 346, 2940–2947. 23. Wada, M.; Hori, R. J. Polym. Sci., Part B: Polym. Phys. 2009, 47, 517–523. 24. Koch, U.; Popelier, P. L. A. J. Phys. Chem. 1995, 99, 9747–9754. 25. Kobayashi, T.; Kohno, Y.; Takayanagi, T.; Seki, K.; Ueda, K. Comput. Theor. Chem. 2012, 991, 48–55. 26. Parthasarathi, R.; Bellesia, G.; Chundawat, S. P. S.; Dale, B. E.; Langan, P.; Gnanakaran, S. J. Phys. Chem. A 2011, 115, 14191–14202.
