Crystalline Moduli of Polymers, Evaluated from Density Functional Theory Calculations under Periodic Boundary Conditions

A theoretical methodology based on quantum chemistry to calculate mechanical properties of polymer crystals has been developed and applied to representative polymers. By density functional theory calculations including a dispersion force correction under three-dimensional periodic boundary conditions, crystal structures of poly(methylene oxide) (PMO), polyethylene (PE), poly(ethylene terephthalate) (PET), poly(trimethylene terephthalate) (PTT), and poly(butylene terephthalate) (PBT) were optimized and their mechanical properties, such as crystalline moduli and linear and volume compressibilities, were calculated. The optimized crystal structures were proved to be fully consistent with those determined by X-ray and neutron diffraction. The crystalline moduli (E∥) parallel to the chain axis were calculated to be 114 GPa (PMO), 333 GPa (PE), 182 GPa (PET), 7.1 GPa (PTT), and 20.8 GPa (PBT) and compared with those determined from X-ray diffraction, Raman spectroscopy, and neutron inelastic scattering experiments. Herein, the E∥ values thus determined are interpreted in terms of conformational characteristics of the polymeric chains and the validity of the homogeneous stress hypothesis adopted in the X-ray diffraction method is also discussed.


■ INTRODUCTION
Mechanical properties are crucial characteristics of polymeric materials, especially from a practical viewpoint. Solid polymers are composed of crystalline and amorphous regions, and the former is more rigid than the latter. In general, polymer crystals exhibit the largest modulus in the chain axis direction; therefore, the crystalline modulus (E ∥ ) parallel to the chain axis at 0 K may correspond to the ultimate hardness expected from the polymer. By reference to the ultimate modulus, one can determine whether the polymer is suitable for one's practical application and estimate how much the mechanical properties will be improved by processings such as stretching and annealing. Therefore, it is of particular significance to determine the crystalline moduli of polymers as precisely as possible. The E ∥ values of polymers have been estimated by experiments such as X-ray diffraction, Raman spectroscopy, and neutron inelastic scattering (NIS). 1,2 For example, the E ∥ value of polyethylene (PE) has been reported as 150−254 GPa (Xray), 260−358 GPa (Raman), and 329 GPa (NIS) (see Table  3); unfortunately, the experimental data are considerably scattered.
Of the experimental methods, X-ray diffraction is the most common. 3 The change (Δd) in lattice spacing due to a load (F) is observed by X-ray diffraction to evaluate the strain (ϵ = Δd/ d 0 ) and the stress (σ) is estimated as F/A, where d 0 is the lattice spacing without load and A is the cross-sectional area of the sample. The slope of the stress−strain curve yields Young's modulus of the crystal. On the assumption that both crystalline and amorphous regions are subjected to a homogeneous stress, the microscopic strain ϵ is directly related to the macroscopic stress σ. The X-ray experimentalists have argued that the homogeneous stress hypothesis was already demonstrated, 4−6 whereas not all researchers have accepted the assumption. 2,7,8 In fact, as referred to PE above, the E ∥ values from X-ray diffraction tend to be smaller than the other experimental and theoretical data (see Table 3). Therefore, the validity of the homogeneous stress model has occasionally been disputed for these several decades in the community of polymer science. For the above reasons, the principle of natural science, i.e., theory must be proved by experiment, would not be simply applied to this case; therefore, a universal methodology to yield the ultimate crystalline moduli of polymers has been desired, regardless of whether it is experimental or theoretical.
From molecular orbital (MO) and density function theory (DFT) calculations, 9−19 Young's moduli (E ∥ ) of PE have been estimated and fall within the range of 300−400 GPa. The models have mostly been a single chain of finite length, namely, n-alkanes; therefore, weak interchain interactions of van der Waals (vdW) forces were not considered and hence the lateral modulus, E ⊥ , cannot be derived from the simplified model. In the poly(ethylene terephthalate) (PET) crystal, 20 interchain π/ π interactions between the benzene rings fulfill a significant role to form and stabilize the crystal structure because the free PET chain strongly prefers a trans−gauche−trans conformation (by approximately −1.2 kcal mol −1 ) to trans−trans−trans in the O−CH 2 −CH 2 −O part 21 but nevertheless crystallizes in the alltrans structure at the expense of the conformational stability so as to form the effective π/π stacking, which must significantly influence both E ∥ and E ⊥ magnitudes. As for such polymer crystals, we are required to accurately evaluate the interchain interactions of vdW forces.
For the above reasons, we have conducted density functional theory calculations with a dispersion force correction under three-dimensional periodic boundary conditions 22,23 for representative polymers with different crystal structures and chain conformations, namely, poly(methylene oxide) (PMO, 9/5 or 29/16 helix), 24−26 polyethylene (PE, all-trans), 27 poly(ethylene terephthalate) (PET, all-trans), 20 poly-(trimethylene terephthalate) (PTT, tggt), 28,29 and poly-(butylene terephthalate) (PBT, g + g + tg − g − ), 30−32 where the abbreviation and crystal conformation are written in the parenthesis. The crystalline moduli in the directions parallel (E ∥ ) and perpendicular (E ⊥ ) to the chain axis or in the crystal a, b, and c axis directions (E a , E b , and E c , respectively) of these polymers have been derived from the calculated compliance tensors and compared with the experimental and theoretical data reported so far. Herein, the theoretical methodology and results are explained in detail and the validity of the experimental data is discussed.

■ RESULTS AND DISCUSSION
Poly(methylene oxide). It is well known that poly-(methylene oxide) whisker, exceptionally, includes no amorphous phase. 33 Therefore, the PMO whisker is expected to give the true crystalline moduli directly; however, the whisker crystal itself is so small that it has to be embedded in resin matrices and undergo the mechanical measurements and Young's modulus (E ∥ ) of the whisker at room temperature was indirectly evaluated to be about 100 GPa. 34 Young's modulus of a highly oriented PMO film was estimated by X-ray diffraction to be 44 GPa at room temperature and 105 GPa at −150°C. 35 From Raman and infrared frequency shifts induced by tensile stress, Young's moduli of highly oriented PMO samples were determined to be 75 GPa at room temperature and approximately 100 GPa at −150°C. 8,36 The E ∥ value of highly oriented PMO sample was calculated from its vibrational spectroscopic data to be 109 GPa. 36 From molecular force-field calculations using the COMPASS parameters, 37 the E ∥ and E ⊥ values were evaluated to be 122 and 28 GPa, respectively. 38 In summary, Young's modulus in the chain axis direction may depend on temperature. However, at as low temperatures as −150°C, the experimental E ∥ data stay within a narrow range of 100−105 GPa, being comparable to the theoretical values (109−122 GPa).
In Table 1, the results of the geometrical optimization at the B3LYP-D/6-31G(d,p) level are summarized. The DFT-D calculations represent the structure and properties at 0 K, thus being preferably compared with the crystal structure at −150°C. 26 Of the four kinds of computations, the set III parameters yielded the best agreement with the cryogenic structure. For the atomic coordinates, see Table S4 (Appendix B, Supporting Information).
The DFT-D calculated crystalline moduli in the chain axis direction range from 105 to 132 GPa, thus being comparable to that (100 GPa) of the PMO whisker 34 and that (105 GPa) measured at −150°C. 35 In contrast, the E ⊥ value depends significantly on the dispersion force correction (5.0−32 GPa). This is because the E ⊥ value reflects interchain forces due mainly to van der Waals (vdW) interactions. The experimental E ⊥ value at room temperature was reported as 8 GPa; 39 therefore, Young's moduli at cryogenic temperatures must be larger than the value. The E ⊥ distributions on the ab plane, calculated from eq 16, are plotted in Figure 1. All of the distributions are perfect circles, whose radii correspond to the magnitudes of E ⊥ . This is because the PMO chain strongly prefers the gauche conformation (the energy difference from the trans state is −1.5 to −2.5 kcal mol −1 ) 40 and hence forms a uniform helical structure. The stiffness (C) and compliance (S) tensors derived from the set III calculation are as follows     (2) Polyethylene. The crystal structure of PE at room temperature was earlier determined by Bunn: 27 an orthorhombic lattice of space group Pnam with a = 7.40, b = 4.93, and c = 2.534 Å. Afterward, its crystal structures at cryogenic temperatures have also been determined by X-ray and neutron diffraction. The neutron diffraction may precisely define atomic coordinates of hydrogen as well as carbon atoms. Table 2 shows the lattice constants and fractional coordinates determined at different temperatures. The a and b dimensions are seen to depend largely on temperature and show a tendency to increase with temperature, whereas the c value is insensitive to temperature. This is because the lateral dimensions (a and b) are influenced by interchain interactions of weak vdW interactions. The variation in the c dimension must be accompanied by changes in bond length and bond angle of the extended all-trans PE chain and hence needs a much larger energy.

ACS Omega
Inasmuch as the DFT-D calculations here reflect the electronic structure at 0 K, it is preferable that the optimized crystal structures should be compared with those at extremely low temperatures. Therefore, the reliability of the DFT-D computations was checked by the RMSD value for the lattice constants at 10, 4, and 4.5 K. 41−43 The RMSD values of the set 0 (without the dispersion force correction), I, II, and III a The setting angle, namely, the angle between the molecular plane and the b axis. b The root-mean-square deviation (RMSD) with respect to the lattice constants observed at 10 K (first), 4 K (second), and 4.5 K (third). c Without the dispersion force correction. d Temperature (K). By neutron diffraction for deuterated PE. 41 e Temperature (K). By neutron diffraction for deuterated PE. 42 The original atomic coordinates are so shifted as to be consistent with those of the other studies. f Temperature (K). By X-ray diffraction. 43 g By X-ray diffraction at room temperature. 27

ACS Omega
Article calculations were calculated with eq 11 to be 28, 12, 3, and below 2%, respectively. Compared with the c value, the a and b lengths are largely dependent on the employed parameter set. This is because, as mentioned above, the dispersion force correction influences the lateral dimensions via the interchain interactions. As shown in Table 2, the set III parameters gave the minimum RMSD and, furthermore, satisfactorily reproduced the atomic coordinates determined by neutron diffraction at 10 and 4 K. The setting angle (45.5°), namely, the angle formed between the molecular plane and the b axis was also in good agreement with those at 10 K (45.07°) 41 and 4.5 K (45.5°). 43 On this basis, the results from the set III calculation will be exclusively referred to in the following discussion. The stiffness and compliance tensors are as follows From the compliance tensor, the crystalline moduli (E a , E b , and E c ) and linear compressibilities (β a , β b , and β c ) in the a, b, and c axis directions and volume compressibility (β v ) of the PE a E a , E b , and E c are Young's moduli parallel to the a, b, and c axes, respectively. When the data are not given explicitly in the literature, Young's moduli were calculated from the compliance tensor. b β a , β b , and β c are, respectively, the linear compressibilities in the a, b, and c directions, and β v is the volume compressibility. When the data are not given explicitly in the literature, the compressibilities were calculated from the compliance tensor: 63 β a = s 11 + s 12 + s 13 , β b = s 12 + s 22 + s 23 , β c = s 13 + s 23 + s 33 , and β v = s 11 + s 22 + s 33 + 2(s 12 + s 23 + s 13 ). c At room temperature.

ACS Omega
Article crystal at 0 K were derived as shown in Table 3. For the mechanical properties obtained from all the parameter sets, see Table S1 (Appendix A, Supporting Information). Fortunately, we can find a number of the experimental and theoretical data on PE crystals. The experiment methods to obtain the mechanical properties are classified into X-ray diffraction, Raman spectroscopy, and neutron inelastic scattering. The Xray method was already outlined in the Introduction. By this method, for example, Clements et al. 44 reported the E ∥ values of PE at room temperature and at −165°C as 150 and 255 GPa, respectively. The sizable temperature dependence was explained as follows. 35 The ratio of the observed modulus E cryst app to the crystalline modulus E cryst may be expressed as  (5) where χ vc is the volume fraction crystallinity and E am is the amorphous modulus. At high temperatures where E am ≪ E cryst , E cryst app /E cryst ratio will approach χ vc . At fully low temperatures, on the other hand, the amorphous region becomes so rigid so as to enlarge E am and hence E cryst app will be close to E cryst . Nakamae et al. 5 measured the E ∥ values of various PE samples such as ultramodulus, linear low-density, high-pressure-crystallized, and heat-drawn PEs and, consequently, obtained the same E ∥ values of 235 GPa at room temperature and 254 GPa at −155°C. The latter E ∥ value exactly agrees with that (255 GPa at −165°C) of Clements et al. 44 On this basis, they have concluded that the assumption of homogeneous stress distribution should be well founded. Matsuo and Sawatari 45 have also argued that the homogeneous stress hypothesis should be valid because the E ∥ values that they observed from PE were 213−229 GPa at 20°C and virtually independent of the elongation ratio. McCullough 1 estimated Young's modulus (as 2.8 GPa) of the isotropic amorphous region of PE from the shear modulus and Poisson's ratio. 46 On the basis of these data, the homogeneous stress model implies that the strain of the amorphous component would be about 90 (=254/2.8) times as much as that of the crystalline region.
The Raman spectroscopic method is in principle based on the following relation 2 where ν is the vibrational frequency of the longitudinal acoustic mode (abbreviated as LAM) and L is the molecular length that can be estimated from the long period (L LP ) and volume fraction crystallinity (χ vc ) according to L = χ vc L LP . 7 The L LP value is determined from small-angle X-ray scattering. In early studies on n-paraffins with definite molecular lengths, comparatively large E ∥ values of 340 and 358 GPa were reported 47,48 and afterward, smaller E ∥ values of 260−315 GPa were evaluated for PE crystals. 8,49−51 Neutron inelastic scattering (NIS) measurements provide phonon frequencies and wave vectors simultaneously. The slope of the frequency vs wave vector plot yields the sound velocity v x in the propagation direction (x) of the phonon 52 where E x is Young's modulus in the x direction. This method has the advantage of informing us of both phonon frequency and wave vector in a specific direction. Therefore, the crystalline moduli along a, b, and c axes may be evaluated separately, whereas the Raman scattering yields the unidirectional LAM, from which only the E ∥ value can be derived. In an early study on NIS of PE, 52,53 the E ∥ and E ⊥ values at room temperature were determined to be 329 and 6 GPa, respectively. In a later article, 54 the E a and E b values at −196°C were reported as 9 and 8 GPa, respectively. These data are larger than those estimated by X-ray diffraction (3.1 and 3.8 GPa). On the grounds that the LAM will not be perturbed by the noncrystalline region, Fanconi and Rabolt 7 have concluded that the crystalline moduli derived from NIS should be the ultimate values and furthermore that the E ∥ data obtained from X-ray diffraction are too low. At least, the Raman and NIS methods have the decisive advantage of selective observations only for the crystallites.
As a theoretical approach, Odajima and Maeda 55 conducted lattice dynamics calculations for PE and derived the crystalline moduli, as shown in Table 3. Wobser and Blasenbrey 56 also calculated the crystalline moduli in a similar manner. Odajima and Maeda used a Urey−Bradley force field for the intramolecular interactions and intermolecular neighboring atom− atom interactions without adapting the potential functions for the crystal structure, in other words, under nonequilibrium conditions. On the other hand, Wobser and Blasenbrey minimized the total energy using their own potential functions and successfully reproduced the experimental crystal structure (by Swan 57 ) extrapolated to 0 K. Therefore, McCullough 1 suggested that the elastic constants of Wobser and Blasenbrey should be more reliable.
The present DFT-D calculations (E a = 10.9 GPa, E b = 7.8 GPa, and E c = 333 GPa) represent the ultimate mechanical properties of the perfect PE crystal of infinite size at 0 K, thus being, naturally, superior to any experimental observation. The crystalline moduli thus obtained are in almost exact agreement with the theoretical values (E a = 9.5 GPa, E b = 8.6 GPa, and E c = 324 GPa at −273°C) of Wobser and Blasenbrey 56 and the NIS experimental data (E a = 9 GPa and E b = 8 GPa at −196°C and E c = 329 GPa at 20°C) 53,54 and also very close to the three-dimensional normal coordinate calculations (E a = 6.9 GPa, E b = 8.6 GPa, and E c = 316 GPa at 20°C) of Tashiro et al. 58 The MO calculations (from Dewar 1979 to Cochrane 2012 in Table 3 58 ). Figure 2 shows the E ⊥ distributions on the ab plane of the PE crystal. The maxima can be found in the directions parallel to the molecular plane.
Aromatic Polyesters. This section deals with three aromatic polyesters, PET, PTT, and PBT (α form). For all of the polyesters, the structural optimizations were carried out with the set 0 (without the dispersion force correction), I, II, and III parameters (the results are compared in Table S2,

ACS Omega
Article Appendix A, Supporting Information). The set III parameters have most satisfactorily represented all of the crystal structures treated here. The principal difference between sets I−III consists in the van der Waals (vdW) radii used. In set III, the vdW radii of heavy atoms except hydrogen were selected from the Bondi table, 64 which is a list vdW radii determined from Xray diffraction data so as to reproduce crystal densities at 0 K, gas kinetic collision cross sections, critical densities, and liquidstate properties. As far as the polymer crystals are concerned, the empirical vdW radii seem to match the Grimme expression.
Therefore, the set III results are mainly discussed here. Figure 3 shows the ac projections of the optimized crystal structures, and Table 4 compares the calculated and observed lattice constants. For the atomic coordinates, see Tables S5−S7 (Appendix B, Supporting Information). The agreement between theory and experiment seems to be satisfactory in every case. It is known that the α form of PBT reversibly changes to the so-called β form under elongation. [30][31][32]65 Accordingly, the β form is regarded as a nonequilibrium structure, from which the crystalline moduli cannot be derived with eq 17. For this reason, this article does not deal with the β form of PBT.
In the crystalline state, the O−CH 2 −CH 2 −O part (referred hereafter to as spacer) of PET adopts the all-trans conformation; however, the ϕ 1 angle around the O−CH 2 bond is somewhat shifted from the normal value (∼180°) to 141.9°, which is intermediate between the two experimental values (159 and 126.0°). 20,66 For the ϕ angles, see Figure 4. The most stable conformation of the PET spacer was proved to be trans−gauche−trans. 21 However, the interchain π/π attraction between the benzene rings renders the spacer alltrans. In a previous study, Galimberti and Milani 67 determined the crystal structure of PTT at the B3LYP-D/6-31G(d,p) level: a = 4.39 Å, b = 6.54 Å, c = 18.38 Å, α = 101.0°, β = 90.0°, γ = 116.5°, ϕ 1 = −164°, ϕ 2 = −58°, ϕ 3 = −67°, and ϕ 4 = −172°. Slight differences between their data and ours would be allowances of the structural optimization. The ϕ angles show that the spacers of PTT and PBT lie in a trans−gauche− gauche−trans and a gauche + −gauche + −trans−gauche − − gauche − conformations, respectively. Both conformations are not the most stable but likely to exist in the solutions and melts. 68, 69 The crystalline moduli were calculated with the sets 0−III parameters, and the results are listed in Table S3 (Appendix A, Supporting Information). Only the stiffness and compliance tensors from set III are presented in Appendix C (Supporting Information). From the compliance tensors, the crystalline moduli (E a , E b , and E c ) in the a, b, and c axis directions were calculated with eq 16 and are compared with the experiment in Table 5. It should be emphasized again that the calculated moduli correspond to the ultimate values at 0 K and would be larger than any experimental data. All of the experimental data in Table 5 were obtained by the X-ray diffraction method. Especially for PET, a number of the data have been reported but they are considerably scattered. The E c (E ∥ ) modulus can also be found to depend largely on crystallinity; for example, Sakurada et al. 62,70 reported two different data: 75 GPa (crystallinity χ c , 39.5%) and 108 GPa (64.5%).
The E c value of PTT was calculated here to be as small as 7.1 GPa and comparable to an experimental value of 5.39 GPa at −255°C, 3 and the X-ray data depend on temperature: 2.59 GPa at 27°C. Interestingly, the E c value of PTT is smaller than E b in both DFT-D calculation and X-ray experiment. The E c magnitude of PBT was calculated to be 20.8 GPa, whereas the experimental E ∥ value was reported as 13.2 GPa. 71 Nevertheless, both DFT-D and X-ray data on the three polyesters give a similar relation between E c (E ∥ ) and spacer conformation; the all-trans spacer of PET 20 exhibits the large E ∥ value, whereas the distorted ones of PTT (tggt) 28,29 and

ACS Omega
Article PBT (g + g + tg − g − ) 30−32 show the much smaller moduli (see Figure 3). Nishino et al. 3 have suggested that PTT has the smallest E c of all polymers investigated so far.
For the lateral crystalline moduli of the three polyesters, the DFT-D calculations show a common tendency of E a < E b , whereas the experiments have yielded E a ≥ E b . In the PET crystal, for example, the benzene rings are stacked in the a axis direction ( Figure 3); therefore, the soft π/π attraction is formed in this direction, whereas in the b axis direction, hard steric repulsions will act. Accordingly, the relation E a < E b seems to be more acceptable.

■ CONCLUSIONS
Crystalline moduli of representative polymers, PMO, PE, PET, PTT, and PBT, were calculated by the density functional theory at the B3LYP-D/6-31G(d,p) level under periodic boundary conditions. The dispersion force correction formulated by Grimme 83,84 was investigated with different parameter sets. As a result, the parameters proposed by Milani et al. 85,86 most satisfactorily reproduced all of the experimental structures. The calculated E ∥ (114 GPa) value of PMO is comparable to that (100 GPa) of the PMO whisker 34 and that (105 GPa) observed at −150°C from a highly oriented PMO film. 35 The calculated E c value (333 GPa) of PE falls within the range of those (260− 358 GPa) derived from Raman spectroscopy 8,47−50 and exactly agrees with that (329 GPa) determined by neutron inelastic scattering 52,53 but is much larger than those (235−255 GPa) estimated by X-ray diffraction. 5,44,45 The DFT-D calculations on the aromatic polyesters yielded the following E c values: 182 GPa (PET), 7.1 GPa (PTT), and 20.8 GPa (PBT). These magnitudes are approximately 1.3−1.7 times as large as X-ray diffraction data, 3,62,70,71,76,78,79 whereas the tendency of PET ≫ PBT > PTT in E c is fully consistent to that found for the X-ray data. Inasmuch as the crystalline moduli obtained from X-ray diffraction depend greatly on crystallinity 62,70 and temperature, 3,35,44 the X-ray data must be significantly perturbed by the amorphous phase. Therefore, it seems reasonable to conclude that the homogeneous stress model 4,5 presupposed in the X-ray experiments would be a coarse approximation except when highly crystalline samples are measured at cryogenic temperatures. Nevertheless, the X-ray method should be still a convenient and useful tool to estimate mechanical properties of polymer crystals. This study has also shown that the E ∥ value depends largely on the conformation of the polymeric chain: all-trans (PE, 333 GPa and PET, 182 GPa) > helical (PMO, 114 GPa) > g + g + tg − g − (PBT, 20.8 GPa) > tggt (PTT, 7.1 GPa). Interestingly, even in the crystalline state, both structures and properties of polymers are strongly influenced by the conformational characteristics.

■ METHODS
Optimization of Crystal Structure. Density functional theory calculations under periodic boundary conditions were conducted with the CRYSTAL14 code. 22,23 The lattice constants and atomic coordinates determined mainly by X-ray diffraction were set as the initial structure and fully optimized at the B3LYP/6-31G(d,p) level with the dispersion force correction to be described below, and the computational conditions were typically as follows: self-consistent field (SCF) convergence threshold, 10 −7 ; truncation criteria for bielectronic integrals, 10 −7 , 10 −7 , 10 −7 , 10 −7 , and 10 −14 ; integration grid, 55 radial and 434 angular points; Fock/Kohn−Sham matrix

ACS Omega
Article mixing, 80% (with the modified Broyden method); 23 shrinking factor, 4. The optimized crystal structures (Figures 1−3) were depicted by Mercury CDC 3.9 program. 87 Dispersion Force Correction. Mean-field electronic structure methods such as the Hartree−Fock and Kohn− Sham DFT theories include no long-range electronic correlation effects, thus being unable to represent the so-called London dispersion interactions. To improve the DFT calculations (E DFT ), a variety of tools for dispersion force corrections have been developed. 88 For instance, Grimme 83,84 has proposed the DFT-D calculation by introducing a correction term (E disp ) The E disp term is expressed as  (9) where s 6 is the global scaling factor, N at is the number of atoms included in the system, C 6 ij is the dispersion coefficient for atom pair ij, given by C 6 ij = (C 6 i C 6 j ) 1/2 , and R ij is the distance between atoms i and j. Here, the damping factor is defined as where d represents the steepness of the damping function and R r is the sum of van der Waals radii (R vdW ) of atoms i and j. We have adopted three sets of dispersion force parameters (Table  6), applied them to the structural optimization, and compared the resultant lattice constants with the experimental data. The

ACS Omega
Article agreement with experiment was checked by the root-meansquare deviation (RMSD) x calcd obsd obsd 2 1/2 (11) or the standard deviation, σ, defined as σ = RMSD/√N (%), where x represents lattice constant (a, b, c, α, β, or γ, dependent on the crystal system) and N is the number of lattice constants. As stated in the Introduction, inasmuch as the experimental crystalline moduli are scattered and hence the comparison with the experimental moduli does not always prove the reliability of the calculation, instead, the validity of the parameter sets has been checked by comparison with experiment in lattice constants throughout this study.
Recently, a modification of the above DFT-D method, the so-called D3 approach, 88 has been developed. As far as polymer crystals are concerned, however, the B3LYP-D/6-31G(d,p) calculations have been extensively used and yielded favorable results if combined with the parameter set suitable for the dispersion force correction. 67,85,86 Thus, this study has employed the B3LYP-D/6-31G(d,p) calculations exclusively.
Crystal Elasticity. 63 The stress tensor is defined as  (12) Inasmuch as the stress tensor is a symmetric matrix, only six elements are irreducible and numerically designated according to Voigt's notation: 1 = xx, 2 = yy, 3 = zz, 4 = yz, 5= xz, and 6 = xy. Then, the stress σ and strain ϵ are expressed as 6 × 1 column vectors and σ is related to ϵ through the stiffness tensor C by The stiffness tensor is a 6 × 6 matrix with 21 independent elements. The crystal symmetry further simplifies the C tensor. The compliance tensor S is related to the ϵ and σ vectors by ϵ = Sσ; therefore, S is the inverse matrix of C, being represented as Young's moduli in the x, y, and z directions (E x , E y , and E z ) are given by E x = s 11 −1 , E y = s 22 −1 , and E z = s 33 −1 , respectively. Young's modulus E(l 1 , l 2 , l 3 ) in an arbitrary direction represented with the unit vector (l 1 , l 2 , l 3 ) can be calculated from