Introduction {#sec1}
============

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).^[@ref1],[@ref2]^ For example, the *E*~∥~ value of polyethylene (PE) has been reported as 150--254 GPa (X-ray), 260--358 GPa (Raman), and 329 GPa (NIS) (see [Table [3](#tbl3){ref-type="other"}](#tbl3){ref-type="other"}); unfortunately, the experimental data are considerably scattered.

Of the experimental methods, X-ray diffraction is the most common.^[@ref3]^ 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,^[@ref4]−[@ref6]^ whereas not all researchers have accepted the assumption.^[@ref2],[@ref7],[@ref8]^ 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](#tbl3){ref-type="other"}](#tbl3){ref-type="other"}). 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,^[@ref9]−[@ref19]^ 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,^[@ref20]^ 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^[@ref21]^ but nevertheless crystallizes in the all-trans 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^[@ref22],[@ref23]^ for representative polymers with different crystal structures and chain conformations, namely, poly(methylene oxide) (PMO, 9/5 or 29/16 helix),^[@ref24]−[@ref26]^ polyethylene (PE, all-trans),^[@ref27]^ poly(ethylene terephthalate) (PET, all-trans),^[@ref20]^ poly(trimethylene terephthalate) (PTT, tggt),^[@ref28],[@ref29]^ and poly(butylene terephthalate) (PBT, g^+^g^+^tg^--^g^--^),^[@ref30]−[@ref32]^ 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 {#sec2}
======================

Poly(methylene oxide) {#sec2.1}
---------------------

It is well known that poly(methylene oxide) whisker, exceptionally, includes no amorphous phase.^[@ref33]^ 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.^[@ref34]^ 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.^[@ref35]^ 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.^[@ref8],[@ref36]^

The *E*~∥~ value of highly oriented PMO sample was calculated from its vibrational spectroscopic data to be 109 GPa.^[@ref36]^ From molecular force-field calculations using the COMPASS parameters,^[@ref37]^ the *E*~∥~ and *E*~⊥~ values were evaluated to be 122 and 28 GPa, respectively.^[@ref38]^ 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](#tbl1){ref-type="other"}](#tbl1){ref-type="other"}, 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.^[@ref26]^ 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](http://pubs.acs.org/doi/suppl/10.1021/acsomega.8b00506/suppl_file/ao8b00506_si_001.pdf) (Appendix B, Supporting Information).

###### Lattice Constants and Crystalline Moduli of PMO

                                                              exptl    DFT-D calcd: set                              
  ----------------------------------------------------------- -------- ------------------ -------- -------- -------- --------
  lattice constant (Å)[c](#t1fn3){ref-type="table-fn"}                                                               
  *a*, *b*                                                    4.464    4.373              4.567    4.173    4.252    4.302
  *c*                                                         17.389   17.274             17.640   17.541   17.742   17.745
  RMSD (rt),[d](#t1fn4){ref-type="table-fn"} %                                            3.57     9.26     7.02     5.53
  RMSD (−150 °C),[e](#t1fn5){ref-type="table-fn"} %                                       6.62     6.65     4.76     3.56
  crystalline modulus (GPa)[f](#t1fn6){ref-type="table-fn"}                                                          
  *E*~∥~                                                      44, 75   100--105           105      132      125      114
  *E*~⊥~                                                      8.0                         5.0      32       20       16

rt: room temperature.

Without the dispersion force correction.

Either (9/5) or (29/16) helical structure.^[@ref25],[@ref26]^ α = β = 90° and γ = 120°.

Calculated with [eq [11](#eq11){ref-type="disp-formula"}](#eq11){ref-type="disp-formula"} for the experimental lattice constants at room temperature.

Calculated with [eq [11](#eq11){ref-type="disp-formula"}](#eq11){ref-type="disp-formula"} for the experimental lattice constants at −150 °C.

The experimental data were determined by vibrating reed method, X-ray diffraction, or vibrational spectroscopy.^[@ref8],[@ref34]−[@ref36],[@ref39]^

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^[@ref34]^ and that (105 GPa) measured at −150 °C.^[@ref35]^ 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;^[@ref39]^ therefore, Young's moduli at cryogenic temperatures must be larger than the value. The *E*~⊥~ distributions on the *ab* plane, calculated from [eq [16](#eq16){ref-type="disp-formula"}](#eq16){ref-type="disp-formula"}, are plotted in [Figure [1](#fig1){ref-type="fig"}](#fig1){ref-type="fig"}. 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^)^[@ref40]^ and hence forms a uniform helical structure. The stiffness (*C*) and compliance (*S*) tensors derived from the set III calculation are as follows

![(Left) *E*~⊥~ distributions on the *ab* plane of the PMO crystal (right): purple, set 0 (without the dispersion force correction); yellow, set I; green, set II; blue, set III; red, the experiment at room temperature. The grid spacing corresponds to 10 GPa.](ao-2018-005067_0005){#fig1}

Polyethylene {#sec2.2}
------------

The crystal structure of PE at room temperature was earlier determined by Bunn:^[@ref27]^ 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](#tbl2){ref-type="other"}](#tbl2){ref-type="other"} 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.

###### Lattice Constants and Atomic Positions of PE

                                          lattice constant, Å           atomic position                                   
  --------------------------------------- --------------------- ------- ----------------- --- ---------- -------- ------- ------
  Calcd                                                                                                                   
  set 0[c](#t2fn3){ref-type="table-fn"}   8.592                 5.795   2.569             C   0.0297     0.0574   37.2    28.5
                                                                                          H   0.1577     0.0585           28.4
                                                                                          H   --0.0041   0.2404           28.3
  set I                                   6.481                 4.448   2.551             C   0.0480     0.0636   47.5    12.1
                                                                                          H   0.2143     0.0187           12.2
                                                                                          H   0.0316     0.3091           12.3
  set II                                  6.941                 4.773   2.569             C   0.0428     0.0618   44.9    2.9
                                                                                          H   0.1994     0.0305           3.1
                                                                                          H   0.0209     0.2898           3.2
  set III                                 7.004                 4.849   2.567             C   0.0424     0.0605   45.5    1.6
                                                                                          H   0.1967     0.0284           1.8
                                                                                          H   0.0217     0.2851           1.8
  Exptl                                                                                                                   
  10[d](#t2fn4){ref-type="table-fn"}      7.120                 4.842                     C   0.042      0.062    45.07    
                                                                                          H   0.193      0.023             
                                                                                          H   0.016      0.284             
  100[d](#t2fn4){ref-type="table-fn"}     7.162                 4.863                     C   0.042      0.062    44.78    
                                                                                          H   0.192      0.024             
                                                                                          H   0.015      0.283             
  200[d](#t2fn4){ref-type="table-fn"}     7.258                 4.898                     C   0.041      0.062    44.35    
                                                                                          H   0.189      0.026             
                                                                                          H   0.013      0.282             
  300[d](#t2fn4){ref-type="table-fn"}     7.417                 4.939                     C   0.040      0.061    44.91    
                                                                                          H   0.185      0.023             
                                                                                          H   0.015      0.278             
  4[e](#t2fn5){ref-type="table-fn"}       7.121                 4.851   2.548             C   0.050      0.062    49       
                                                                                          H   0.198      0.034             
                                                                                          H   0.019      0.283             
  90[e](#t2fn5){ref-type="table-fn"}      7.161                 4.866   2.546             C   0.049      0.062    49       
                                                                                          H   0.197      0.030             
                                                                                          H   0.018      0.278             
  4.5[f](#t2fn6){ref-type="table-fn"}     7.128                 4.852   2.553             C   0.046      0.065    45.5     
  rt[g](#t2fn7){ref-type="table-fn"}      7.40                  4.93    2.534             C   0.038      0.065    41.2     

The setting angle, namely, the angle between the molecular plane and the *b* axis.

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).

Without the dispersion force correction.

Temperature (K). By neutron diffraction for deuterated PE.^[@ref41]^

Temperature (K). By neutron diffraction for deuterated PE.^[@ref42]^ The original atomic coordinates are so shifted as to be consistent with those of the other studies.

Temperature (K). By X-ray diffraction.^[@ref43]^

By X-ray diffraction at room temperature.^[@ref27]^

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.^[@ref41]−[@ref43]^ The RMSD values of the set 0 (without the dispersion force correction), I, II, and III calculations were calculated with [eq [11](#eq11){ref-type="disp-formula"}](#eq11){ref-type="disp-formula"} 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](#tbl2){ref-type="other"}](#tbl2){ref-type="other"}, 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°)^[@ref41]^ and 4.5 K (45.5°).^[@ref43]^ 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 crystal at 0 K were derived as shown in [Table [3](#tbl3){ref-type="other"}](#tbl3){ref-type="other"}. For the mechanical properties obtained from all the parameter sets, see [Table S1](http://pubs.acs.org/doi/suppl/10.1021/acsomega.8b00506/suppl_file/ao8b00506_si_001.pdf) (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 X-ray method was already outlined in the [Introduction](#sec1){ref-type="other"}. By this method, for example, Clements et al.^[@ref44]^ 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.^[@ref35]^ The ratio of the observed modulus *E*~cryst~^app^ to the crystalline modulus *E*~cryst~ may be expressed aswhere χ~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.^[@ref5]^ 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.^[@ref44]^ On this basis, they have concluded that the assumption of homogeneous stress distribution should be well founded. Matsuo and Sawatari^[@ref45]^ 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^[@ref1]^ estimated Young's modulus (as 2.8 GPa) of the isotropic amorphous region of PE from the shear modulus and Poisson's ratio.^[@ref46]^ 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.

###### Crystalline Moduli and Linear and Volume Compressibilities of PE

  literature                                     crystalline modulus,[a](#t3fn1){ref-type="table-fn"} GPa   compressibility,[b](#t3fn2){ref-type="table-fn"} TPa^--1^                                          
  ----------------------------- ------ --------- ---------------------------------------------------------- ----------------------------------------------------------- ------ ---------- ------- ------ ----- -----
  Theoretical                                                                                                                                                                                                  
  this study (set III)                           --273                                                      10.9                                                        7.8    333        21.2    55.6   2.1   79
  Odajima (set III)^[@ref55]^   1966             20                                                         5.9                                                         9.1    256        133.0   75.6   0.6   209
                                                 --196                                                      8.0                                                         9.4    253        82.4    63.1   1.1   147
  Williams^[@ref61]^            1967             --273                                                      13.7                                                        11.9                                    
  Wobser^[@ref56]^              1970             --273                                                      9.5                                                         8.6    324        42.7    53.7   2.1   99
  Tashiro^[@ref58]^             1978             20                                                         6.9                                                         8.6    316        96.8    68.3   2.4   168
  Dewar^[@ref9]^                1979             --273                                                                                                                         494                              
  Suhai^[@ref10]^               1983             --273                                                                                                                         305                              
  Suhai^[@ref11]^               1986             --273                                                                                                                         276, 334                         
  Horn^[@ref12]^                1993             --273                                                                                                                         340                              
  Crist^[@ref13]^               1996             --273                                                                                                                         336                              
                                                 rt[c](#t3fn3){ref-type="table-fn"}                                                                                            299                              
  Hageman^[@ref14],[@ref15]^    1998             --273                                                                                                                         334                              
  Zhang^[@ref16]^               2000             --273                                                                                                                         360                              
  Peeters^[@ref17]^             2000             --273                                                                                                                         340                              
  Miao^[@ref18]^                2001             --273                                                                                                                         350--400                         
  Cochrane^[@ref19]^            2012             --273                                                                                                                         350                              
  Experimental                                                                                                                                                                                                 
  Sakurada^[@ref62]^            1966   X-ray     20                                                         3.1                                                         3.8                                     
  Ito^[@ref59]^                 1971   X-ray     20                                                                                                                                       117     74.0   3.1   195
  Clements^[@ref44]^            1978   X-ray     rt[c](#t3fn3){ref-type="table-fn"}                                                                                            150                              
                                                 --165                                                                                                                         255                              
  Ito^[@ref60]^                 1982   X-ray     20                                                                                                                                       87.7    67.2   1.3   156
  Matsuo^[@ref45]^              1986   X-ray     20                                                                                                                            213--229                         
  Nakamae^[@ref5]^              1991   X-ray     rt[c](#t3fn3){ref-type="table-fn"}                                                                                            235                              
                                                 --155                                                                                                                         254                              
  Mizushima^[@ref47]^           1949   Raman                                                                                                                                   340                              
  Schaufele^[@ref48]^           1967   Raman     27                                                                                                                            358                              
  Strobl^[@ref50]^              1976   Raman     rt[c](#t3fn3){ref-type="table-fn"}                                                                                            290                              
  Kobayashi^[@ref49]^           1983   Raman     rt[c](#t3fn3){ref-type="table-fn"}                                                                                            281                              
  Tashiro^[@ref8]^              1988   Raman     rt[c](#t3fn3){ref-type="table-fn"}                                                                                            260                              
  Pietralla^[@ref51]^           1997   Raman     rt[c](#t3fn3){ref-type="table-fn"}                                                                                            315                              
  Feldkamp^[@ref52],[@ref53]^   1968   neutron   20                                                                                                                            329                              
  Holliday^[@ref52]^            1971   neutron   25                                                         6                                                           6      329                              
  Twisleton^[@ref54]^           1982   neutron   --196                                                      9                                                           8                                       

*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.

β*~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:^[@ref63]^ β*~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~).

At room temperature.

The Raman spectroscopic method is in principle based on the following relation^[@ref2]^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~.^[@ref7]^ 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^[@ref47],[@ref48]^ and afterward, smaller *E*~∥~ values of 260--315 GPa were evaluated for PE crystals.^[@ref8],[@ref49]−[@ref51]^

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^[@ref52]^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,^[@ref52],[@ref53]^ the *E*~∥~ and *E*~⊥~ values at room temperature were determined to be 329 and 6 GPa, respectively. In a later article,^[@ref54]^ 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^[@ref7]^ 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^[@ref55]^ conducted lattice dynamics calculations for PE and derived the crystalline moduli, as shown in [Table [3](#tbl3){ref-type="other"}](#tbl3){ref-type="other"}. Wobser and Blasenbrey^[@ref56]^ 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^[@ref57]^) extrapolated to 0 K. Therefore, McCullough^[@ref1]^ 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^[@ref56]^ 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)^[@ref53],[@ref54]^ 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.^[@ref58]^ The MO calculations (from Dewar 1979 to Cochrane 2012 in [Table [3](#tbl3){ref-type="other"}](#tbl3){ref-type="other"}),^[@ref10]−[@ref19]^ notwithstanding at different levels of theory with various basis sets, yielded similar *E*~*c*~ (*E*~∥~) values of 300--400 GPa except 494 GPa of Dewar 1979.^[@ref9]^ These Young's moduli are comparable to or consistent with our calculation of 333 GPa.

Both theoretical and experimental values of β~*a*~ and β~*b*~ are much larger than β~*c*~ because the PE crystal is weak (deformable) under external pressure applied in the lateral directions. However, the above-mentioned theoretical calculations yielded similar β~*c*~ values between the two X-ray experiments (1.3 and 3.1 GPa by Ito^[@ref59],[@ref60]^): 2.1 TPa^--1^ (this study), 2.1 TPa^--1^ (Wobser and Blasenbrey^[@ref56]^), and 2.4 TPa^--1^ (Tashiro et al.^[@ref58]^).

[Figure [2](#fig2){ref-type="fig"}](#fig2){ref-type="fig"} 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.

![(Left) *E*~⊥~ distribution on the *ab* plane of the PE crystal (right), obtained from the set III parameters. The grid spacing corresponds to 10 GPa.](ao-2018-005067_0004){#fig2}

Aromatic Polyesters {#sec2.3}
-------------------

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](http://pubs.acs.org/doi/suppl/10.1021/acsomega.8b00506/suppl_file/ao8b00506_si_001.pdf), 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,^[@ref64]^ which is a list vdW radii determined from X-ray diffraction data so as to reproduce crystal densities at 0 K, gas kinetic collision cross sections, critical densities, and liquid-state 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](#fig3){ref-type="fig"}](#fig3){ref-type="fig"} shows the *ac* projections of the optimized crystal structures, and [Table [4](#tbl4){ref-type="other"}](#tbl4){ref-type="other"} compares the calculated and observed lattice constants. For the atomic coordinates, see [Tables S5--S7](http://pubs.acs.org/doi/suppl/10.1021/acsomega.8b00506/suppl_file/ao8b00506_si_001.pdf) (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.^[@ref30]−[@ref32],[@ref65]^ Accordingly, the β form is regarded as a nonequilibrium structure, from which the crystalline moduli cannot be derived with [eq [17](#eq17){ref-type="disp-formula"}](#eq17){ref-type="disp-formula"}. 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°).^[@ref20],[@ref66]^ For the ϕ angles, see [Figure [4](#fig4){ref-type="fig"}](#fig4){ref-type="fig"}. The most stable conformation of the PET spacer was proved to be trans--gauche--trans.^[@ref21]^ However, the interchain π/π attraction between the benzene rings renders the spacer all-trans. In a previous study, Galimberti and Milani^[@ref67]^ 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.^[@ref68],[@ref69]^

![Optimized crystal structures (projections on the *ac* plane) of (a) PET, (b) PTT, and (c) PBT.](ao-2018-005067_0003){#fig3}

![(a) PET, (b) PTT, and (c) PBT. The dihedral angles ϕ's are defined as illustrated.](ao-2018-005067_0001){#fig4}

###### Lattice Constants and Chain Conformations of PET, PTT, and PBT

                                                        PET     PTT       PBT (α form)                                                 
  ----------------------------------------------------- ------- --------- -------------- -------- --------- --------- ------- -------- --------
  Lattice Constant, Å, deg                                                                                                             
  *a*                                                   4.56    4.51      4.45           4.64     4.59      4.38      4.83    4.89     4.59
  *b*                                                   5.94    5.88      5.93           6.27     6.21      6.00      5.94    5.95     5.88
  *c*                                                   10.75   10.79     10.76          18.64    18.31     18.35     11.59   11.67    11.81
  α                                                     98.5    100.0     101.6          98.4     98.0      101.9     99.7    98.9     100.9
  β                                                     118     118.4     119.2          93.0     90.0      90.1      115.2   116.6    114.3
  γ                                                     112     110.6     112.7          111.1    111.7     103.8     110.8   110.9    111.2
  RMSD/√*N*,[h](#t4fn8){ref-type="table-fn"} %          1.7     1.2                      4.4      4.0                 2.3     2.8       
  Dihedral Angle,[i](#t4fn9){ref-type="table-fn"} deg                                                                                  
  ϕ~1~                                                  159     126.0     141.9          180.0    --152.8   --166.0   --88    --92.5   --94.4
  ϕ~2~                                                  180     180.0     180.0          --67.0   --61.9    --58.0    --68    --76.2   --66.7
  ϕ~3~                                                  --159   --126.0   --141.9        --67.0   --73.7    --67.0    180     180      180.0
  ϕ~4~                                                                                   180.0    167.4     --169.9   68      76.2     66.7
  ϕ~5~                                                                                                                88      92.5     94.4

By Daubeny et al.^[@ref20]^

By Fu et al. (sample, PET4B).^[@ref66]^

This study (optimized with set III parameters).

By Poulin-Dandurand et al.^[@ref28]^

By Desborough et al.^[@ref29]^

By Yokouchi et al.^[@ref30]^

By Hall et al.^[@ref31],[@ref32]^

The standard deviation of the lattice constants for each experiment. *N* is the number of lattice constants (*N* = 6).

See [Figure [4](#fig4){ref-type="fig"}](#fig4){ref-type="fig"}. The dihedral angle is defined here according to the IUPAC recommendation (trans ∼ 180°, gauche^+^ ∼ −60°, and gauche^--^ ∼ +60°).^[@ref72]^

The crystalline moduli were calculated with the sets 0--III parameters, and the results are listed in [Table S3](http://pubs.acs.org/doi/suppl/10.1021/acsomega.8b00506/suppl_file/ao8b00506_si_001.pdf) (Appendix A, Supporting Information). Only the stiffness and compliance tensors from set III are presented in Appendix C ([Supporting Information](http://pubs.acs.org/doi/suppl/10.1021/acsomega.8b00506/suppl_file/ao8b00506_si_001.pdf)). 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](#eq16){ref-type="disp-formula"}](#eq16){ref-type="disp-formula"} and are compared with the experiment in [Table [5](#tbl5){ref-type="other"}](#tbl5){ref-type="other"}. 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](#tbl5){ref-type="other"}](#tbl5){ref-type="other"} 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.^[@ref62],[@ref70]^ reported two different data: 75 GPa (crystallinity χ~c~, 39.5%) and 108 GPa (64.5%).

###### Crystalline Moduli of PET, PTT, and PBT

  literature                                     crystalline modulus,[a](#t5fn1){ref-type="table-fn"} GPa                           
  ---------------------------------------------- ---------------------------------------------------------- ----- ------ ---------- -------------------------
  PET                                                                                                                               
  theoretical                                                                                                                       
  this study (set III)                                                                                      7.2   22.3   182        --273 °C
  Lyons^[@ref73]^                                1958                                                                    146         
  Treloar^[@ref74]^                              1960                                                                    122         
  Tashiro^[@ref75]^                              1977                                                                    95          
  experimental[c](#t5fn3){ref-type="table-fn"}                                                                                      
  Dulmage^[@ref76]^                              1958                                                                    137         
  Sakurada^[@ref62]^                             1966                                                                    75         χ~c~ 39.5%
  Sakurada^[@ref70]^                             1970                                                                    108        χ~c~ 64.5%
  Sakurada^[@ref77]^                             1975                                                       4.6   4.1               23 °C
  Thistlethwaite^[@ref78]^                       1988                                                                    100--129   χ~c~ 50--73%, 20 °C
                                                                                                                         94--132    χ~c~ 50--73%, −124.4 °C
  Matsuo^[@ref79]^                               1990                                                                    123        \[1̅ 0 5\] direction
                                                                                                                         118        *c* axis direction
  PTT                                                                                                                               
  theoretical                                                                                                                       
  this study (set III)                                                                                      6.9   18.4   7.1        --273 °C
  Jang^[@ref80]^                                 2000                                                       5.9   9.2    12.21       
  experimental[c](#t5fn3){ref-type="table-fn"}                                                                                      
  Kim^[@ref81]^                                  2003                                                             4.2                
  Nishino^[@ref3]^                               2011                                                             3.7    2.59       27 °C
                                                                                                                  7.4    5.39       --255 °C
  PBT (α Form)                                                                                                                      
  theoretical                                                                                                                       
  this study (set III)                                                                                      4.8   11.6   20.8       --273 °C
  Tashiro^[@ref82]^                              1980                                                                    10.1        
  experimental[c](#t5fn3){ref-type="table-fn"}                                                                                      
  Nakamae^[@ref39]^                              1980                                                       2.4   2.2                
  Nakamae^[@ref71]^                              1982                                                                    13.2        

*E*~*a*~, *E*~*b*~, and *E*~*c*~ are Young's moduli parallel to the *a*, *b*, and *c* axes, respectively.

Crystallinity (χ~c~), measurement temperature, and direction.

In most X-ray diffraction experiments on PET and PBT, the modulus in the \[1̅0 5\] (PET) or \[1̅ 0 4\] (PBT) direction has been substituted for *E*~*c*~. As for PTT, the (002) diffraction was used.^[@ref3]^

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,^[@ref3]^ 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.^[@ref71]^ 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^[@ref20]^ exhibits the large *E*~∥~ value, whereas the distorted ones of PTT (tggt)^[@ref28],[@ref29]^ and PBT (g^+^g^+^tg^--^g^--^)^[@ref30]−[@ref32]^ show the much smaller moduli (see [Figure [3](#fig3){ref-type="fig"}](#fig3){ref-type="fig"}). Nishino et al.^[@ref3]^ 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](#fig3){ref-type="fig"}](#fig3){ref-type="fig"}); 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 {#sec3}
===========

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^[@ref83],[@ref84]^ was investigated with different parameter sets. As a result, the parameters proposed by Milani et al.^[@ref85],[@ref86]^ 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^[@ref34]^ and that (105 GPa) observed at −150 °C from a highly oriented PMO film.^[@ref35]^ The calculated *E~c~* value (333 GPa) of PE falls within the range of those (260--358 GPa) derived from Raman spectroscopy^[@ref8],[@ref47]−[@ref50]^ and exactly agrees with that (329 GPa) determined by neutron inelastic scattering^[@ref52],[@ref53]^ but is much larger than those (235--255 GPa) estimated by X-ray diffraction.^[@ref5],[@ref44],[@ref45]^ 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,^[@ref3],[@ref62],[@ref70],[@ref71],[@ref76],[@ref78],[@ref79]^ 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^[@ref62],[@ref70]^ and temperature,^[@ref3],[@ref35],[@ref44]^ the X-ray data must be significantly perturbed by the amorphous phase. Therefore, it seems reasonable to conclude that the homogeneous stress model^[@ref4],[@ref5]^ 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 {#sec4}
=======

Optimization of Crystal Structure {#sec4.1}
---------------------------------

Density functional theory calculations under periodic boundary conditions were conducted with the CRYSTAL14 code.^[@ref22],[@ref23]^ 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 mixing, 80% (with the modified Broyden method);^[@ref23]^ shrinking factor, 4. The optimized crystal structures ([Figures [1](#fig1){ref-type="fig"}](#fig1){ref-type="fig"}--[3](#fig3){ref-type="fig"}) were depicted by Mercury CDC 3.9 program.^[@ref87]^

Dispersion Force Correction {#sec4.2}
---------------------------

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.^[@ref88]^ For instance, Grimme^[@ref83],[@ref84]^ has proposed the DFT-D calculation by introducing a correction term (*E*~disp~)The *E*~disp~ term is expressed aswhere *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 aswhere *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](#tbl6){ref-type="other"}](#tbl6){ref-type="other"}), applied them to the structural optimization, and compared the resultant lattice constants with the experimental data. The agreement with experiment was checked by the root-mean-square deviation (RMSD)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](#sec1){ref-type="other"}, 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.

###### Parameters Adopted for Dispersion Force Correction[a](#t6fn1){ref-type="table-fn"}

                                                 sets             
  ---------------------------------------------- ------- -------- --------
  *s*~6~                                         1.05    1.00     1.00
  *R*~vdW~ (H)[e](#t6fn5){ref-type="table-fn"}   1.001   1.3013   1.3013
  *R*~vdW~ (C)[e](#t6fn5){ref-type="table-fn"}   1.452   1.5246   1.70
  *R*~vdW~ (O)[e](#t6fn5){ref-type="table-fn"}   1.342   1.4091   1.52

*C*~6~^H^ = 0.14, *C*~6~^C^ = 1.75, *C*~6~^O^ = 0.70, and *d* = 20.0 are common to the three sets. The cutoff distance to truncate direct lattice summation was set equal to 25.0 Å.

By Grimme.^[@ref83],[@ref84]^

By Civalleri, et al.^[@ref89]^

By Milani, et al.^[@ref85],[@ref86]^

*R*~vdW~ (X): van der Waals radius of element X (X = H, C, or O).

Recently, a modification of the above DFT-D method, the so-called D3 approach,^[@ref88]^ 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.^[@ref67],[@ref85],[@ref86]^ Thus, this study has employed the B3LYP-D/6-31G(d,p) calculations exclusively.

Crystal Elasticity^[@ref63]^ {#sec4.3}
----------------------------

The stress tensor is defined asInasmuch 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* bywhereThe 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 asYoung'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

The element *c*~*uv*~ (*u* and *v* = 1--6, Voigt's notation) of the stiffness tensor can be calculated fromwhere ϵ is the rank-2 symmetric tensor of pure strain and *E*~cryst~ and *V* are the total energy and volume per crystal cell, respectively.

The stiffness and compliance tensors were calculated with the ELASTCON routine^[@ref90],[@ref91]^ included in CRYSTAL14 code. The typical conditions were as follows. The SCF convergence criterion was 10^--8^, the number of points for calculation of the numerical second derivative was 5, and the strain-step size along a given deformation was 0.005.

The Supporting Information is available free of charge on the [ACS Publications website](http://pubs.acs.org) at DOI: [10.1021/acsomega.8b00506](http://pubs.acs.org/doi/abs/10.1021/acsomega.8b00506).(Appendix A, Tables S1--S3) crystalline moduli and linear and volume compressibilities of PE and lattice constants and crystalline moduli of PET, PTT, and α form of PBT; (Appendix B, Tables S4--S7) fractional coordinates of PMO, PET, PTT, and α form of PBT; (Appendix C) stiffness and compliance tensors of PET, PTT, and α form of PBT ([PDF](http://pubs.acs.org/doi/suppl/10.1021/acsomega.8b00506/suppl_file/ao8b00506_si_001.pdf))

Supplementary Material
======================

###### 

ao8b00506_si_001.pdf

^‡^ Iwatani Co. Ltd., Osaka 541-0053, Japan (M.T.).

^†^ Unimatec Co. Ltd., Kitaibaraki-shi 319-1593, Japan (T.K.).

The authors declare no competing financial interest.

We thank Professor Tashiro of Toyota Technological Institute and Professor Nishino of Kobe University for giving us helpful comments at an early stage of this study and the late Professor Roberto Orlando of University of Torino for helpful advice on computations with CRYSTAL14 program. This study was partially supported by the Grants-in-Aid for Scientific Research (C) (16K05906) from the Japan Society for the Promotion of Science.
