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 (X-ray), 260–358 GPa (Raman), and 329 GPa (NIS) (see Table ); unfortunately, the experimental data are considerably scattered.

Table 3

Crystalline Moduli and Linear and Volume Compressibilities of PE

literature				crystalline modulus,a GPa			compressibility,b TPa–1
first author	year	method	temp, °C	Ea	Eb	Ec (E∥)	βa	βb	βc	βv
Theoretical
this study (set III)			–273	10.9	7.8	333	21.2	55.6	2.1	79
Odajima (set III)[55]	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[61]	1967		–273	13.7	11.9
Wobser[56]	1970		–273	9.5	8.6	324	42.7	53.7	2.1	99
Tashiro[58]	1978		20	6.9	8.6	316	96.8	68.3	2.4	168
Dewar[9]	1979		–273			494
Suhai[10]	1983		–273			305
Suhai[11]	1986		–273			276, 334
Horn[12]	1993		–273			340
Crist[13]	1996		–273			336
			rtc			299
Hageman[14,15]	1998		–273			334
Zhang[16]	2000		–273			360
Peeters[17]	2000		–273			340
Miao[18]	2001		–273			350–400
Cochrane[19]	2012		–273			350
Experimental
Sakurada[62]	1966	X-ray	20	3.1	3.8
Ito[59]	1971	X-ray	20				117	74.0	3.1	195
Clements[44]	1978	X-ray	rtc			150
			–165			255
Ito[60]	1982	X-ray	20				87.7	67.2	1.3	156
Matsuo[45]	1986	X-ray	20			213–229
Nakamae[5]	1991	X-ray	rtc			235
			–155			254
Mizushima[47]	1949	Raman				340
Schaufele[48]	1967	Raman	27			358
Strobl[50]	1976	Raman	rtc			290
Kobayashi[49]	1983	Raman	rtc			281
Tashiro[8]	1988	Raman	rtc			260
Pietralla[51]	1997	Raman	rtc			315
Feldkamp[52,53]	1968	neutron	20			329
Holliday[52]	1971	neutron	25	6	6	329
Twisleton[54]	1982	neutron	–196	9	8

E, E, and E 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.

β, β, and β 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] β = s11 + s12 + s13, β = s12 + s22 + s23, β = s13 + s23 + s33, and βv = s11 + s22 + s33 + 2(s12 + s23 + s13).

At room temperature.

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/d0) and the stress (σ) is estimated as F/A, where d0 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 ). 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–CH2–CH2–O part[21] 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[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, E, and E, 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 , 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).

Table 1

Lattice Constants and Crystalline Moduli of PMO

	exptl		DFT-D calcd: set
	rta	–150 °C	0b	I	II	III
lattice constant (Å)c
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 %			3.57	9.26	7.02	5.53
RMSD (−150 °C),e %			6.62	6.65	4.76	3.56
crystalline modulus (GPa)f
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.[25,26] α = β = 90° and γ = 120°.

Calculated with eq for the experimental lattice constants at room temperature.

Calculated with eq for the experimental lattice constants at −150 °C.

The experimental data were determined by vibrating reed method, X-ray diffraction, or vibrational spectroscopy.[8,34−36,39]

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 , are plotted in Figure . 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

Figure 1

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

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

Table 2

Lattice Constants and Atomic Positions of PE

	lattice constant, Å				atomic position
	a	b	c		x/a	y/b	ϕ,a deg	RMSD,b %
Calcd
set 0c	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
10d	7.120	4.842		C	0.042	0.062	45.07
				H	0.193	0.023
				H	0.016	0.284
100d	7.162	4.863		C	0.042	0.062	44.78
				H	0.192	0.024
				H	0.015	0.283
200d	7.258	4.898		C	0.041	0.062	44.35
				H	0.189	0.026
				H	0.013	0.282
300d	7.417	4.939		C	0.040	0.061	44.91
				H	0.185	0.023
				H	0.015	0.278
4e	7.121	4.851	2.548	C	0.050	0.062	49
				H	0.198	0.034
				H	0.019	0.283
90e	7.161	4.866	2.546	C	0.049	0.062	49
				H	0.197	0.030
				H	0.018	0.278
4.5f	7.128	4.852	2.553	C	0.046	0.065	45.5
rtg	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.[41]

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.

Temperature (K). By X-ray diffraction.[43]

By X-ray diffraction at room temperature.[27]

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 calculations were calculated with eq 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 , 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, E, and E) and linear compressibilities (β, β, and β) 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 . 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 X-ray 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 Ecrystapp to the crystalline modulus Ecryst may be expressed aswhere χvc is the volume fraction crystallinity and Eam is the amorphous modulus. At high temperatures where Eam ≪ Ecryst, Ecrystapp/Ecryst ratio will approach χvc. At fully low temperatures, on the other hand, the amorphous region becomes so rigid so as to enlarge Eam and hence Ecrystapp will be close to Ecryst. 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 (LLP) and volume fraction crystallinity (χvc) according to L = χvcLLP.[7] The LLP 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 in the propagation direction (x) of the phonon[52]where E 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 and E 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 . 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 = 10.9 GPa, E = 7.8 GPa, and E = 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 = 9.5 GPa, E = 8.6 GPa, and E = 324 GPa at −273 °C) of Wobser and Blasenbrey[56] and the NIS experimental data (E = 9 GPa and E = 8 GPa at −196 °C and E = 329 GPa at 20 °C)[53,54] and also very close to the three-dimensional normal coordinate calculations (E = 6.9 GPa, E = 8.6 GPa, and E = 316 GPa at 20 °C) of Tashiro et al.[58] The MO calculations (from Dewar 1979 to Cochrane 2012 in Table ),[10−19] notwithstanding at different levels of theory with various basis sets, yielded similar E (E∥) values of 300–400 GPa except 494 GPa of Dewar 1979.[9] These Young’s moduli are comparable to or consistent with our calculation of 333 GPa.

Both theoretical and experimental values of β and β are much larger than β because the PE crystal is weak (deformable) under external pressure applied in the lateral directions. However, the above-mentioned theoretical calculations yielded similar β values between the two X-ray experiments (1.3 and 3.1 GPa by Ito[59,60]): 2.1 TPa–1 (this study), 2.1 TPa–1 (Wobser and Blasenbrey[56]), and 2.4 TPa–1 (Tashiro et al.[58]).

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

Figure 2

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

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, 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 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 shows the ac projections of the optimized crystal structures, and Table 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−32,65] Accordingly, the β form is regarded as a nonequilibrium structure, from which the crystalline moduli cannot be derived with eq . For this reason, this article does not deal with the β form of PBT.

Figure 3

Optimized crystal structures (projections on the ac plane) of (a) PET, (b) PTT, and (c) PBT.

Table 4

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

	PET			PTT			PBT (α form)
	expt 1a	expt 2b	optc	expt 1d	expt 2e	optc	expt 1f	expt 2g	optc
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 %	1.7	1.2		4.4	4.0		2.3	2.8
Dihedral Angle,i 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.[20]

By Fu et al. (sample, PET4B).[66]

This study (optimized with set III parameters).

By Poulin-Dandurand et al.[28]

By Desborough et al.[29]

By Yokouchi et al.[30]

By Hall et al.[31,32]

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

See Figure . The dihedral angle is defined here according to the IUPAC recommendation (trans ∼ 180°, gauche+ ∼ −60°, and gauche– ∼ +60°).[72]

In the crystalline state, the O–CH2–CH2–O part (referred hereafter to as spacer) of PET adopts the all-trans conformation; however, the ϕ1 angle around the O–CH2 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 . 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 all-trans. 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]

Figure 4

(a) PET, (b) PTT, and (c) PBT. The dihedral angles ϕ’s are defined as illustrated.

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, E, and E) in the a, b, and c axis directions were calculated with eq and are compared with the experiment in Table . 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 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 (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%).

Table 5

Crystalline Moduli of PET, PTT, and PBT

literature		crystalline modulus,a GPa
first author	year	Ea	Eb	Ec	noteb
PET
theoretical
this study (set III)		7.2	22.3	182	–273 °C
Lyons[73]	1958			146
Treloar[74]	1960			122
Tashiro[75]	1977			95
experimentalc
Dulmage[76]	1958			137
Sakurada[62]	1966			75	χc 39.5%
Sakurada[70]	1970			108	χc 64.5%
Sakurada[77]	1975	4.6	4.1		23 °C
Thistlethwaite[78]	1988			100–129	χc 50–73%, 20 °C
				94–132	χc 50–73%, −124.4 °C
Matsuo[79]	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[80]	2000	5.9	9.2	12.21
experimentalc
Kim[81]	2003		4.2
Nishino[3]	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[82]	1980			10.1
experimentalc
Nakamae[39]	1980	2.4	2.2
Nakamae[71]	1982			13.2

E, E, and E 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. As for PTT, the (002) diffraction was used.[3]

The E 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 value of PTT is smaller than E in both DFT-D calculation and X-ray experiment. The E 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 (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 PBT (g+g+tg–g–)[30−32] show the much smaller moduli (see Figure ). Nishino et al.[3] have suggested that PTT has the smallest E 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 < E, whereas the experiments have yielded E ≥ E. In the PET crystal, for example, the benzene rings are stacked in the a axis direction (Figure ); therefore, the soft π/π attraction is formed in this direction, whereas in the b axis direction, hard steric repulsions will act. Accordingly, the relation E < E 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 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 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 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 mixing, 80% (with the modified Broyden method);[23] shrinking factor, 4. The optimized crystal structures (Figures –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 (EDFT), 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 (Edisp)The Edisp term is expressed aswhere s6 is the global scaling factor, Nat is the number of atoms included in the system, C6 is the dispersion coefficient for atom pair ij, given by C6 = (C6C6)1/2, and R is the distance between atoms i and j. Here, the damping factor is defined aswhere d represents the steepness of the damping function and Rr is the sum of van der Waals radii (RvdW) of atoms i and j. We have adopted three sets of dispersion force parameters (Table ), 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, 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.

Table 6

Parameters Adopted for Dispersion Force Correctiona

	sets
parameters	Ib	IIc	IIId
s6	1.05	1.00	1.00
RvdW (H)e	1.001	1.3013	1.3013
RvdW (C)e	1.452	1.5246	1.70
RvdW (O)e	1.342	1.4091	1.52

C6H = 0.14, C6C = 1.75, C6O = 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.[83,84]

By Civalleri, et al.[89]

By Milani, et al.[85,86]

RvdW (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,[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 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, E, and E) are given by E = s11–1, E = s22–1, and E = s33–1, respectively. Young’s modulus E(l1, l2, l3) in an arbitrary direction represented with the unit vector (l1, l2, l3) can be calculated from

The element c (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 Ecryst and V are the total energy and volume per crystal cell, respectively.

The stiffness and compliance tensors were calculated with the ELASTCON routine[90,91] 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.