All-atomistic molecular dynamics study of the glass transition of amorphous polymers
Glass transition is an important phenomenon of polymer materials and it has been intensively studied over the past a few decades. However, the influencing factors arising from the chemical structures of the polymers are often ignored due to a continuous or coarsegrained description of the polymer. Here, we approached this phenomenon using all-atomistic molecular dynamics (MD) simulations and two conventionally used polymer materials, polycarbonate (PC) and poly-(methyl methacrylate) (PMMA). We reproduced the glass transition temperatures (Tg) of the two materials reasonably well. Then we characterized and investigated the glass transition process by looking at the changes of potential energy, dihedral transition, and thermal fluctuation of the individual degrees of freedom in the systems, over the entire temperature range of glass transition. As previously reported, the dihedral angles stop their conformational changes gradually at the Tg, especially for the main chain dihedrals, and sidechain rotations immediately rooting from the main chain. The volumetric change during the temperature decrease is confirmed to be because of conformational adjustment, probably due to the tendency of chain stretching for the maintenance of the radius of gyration, and the loss of thermal energy. The strength of motions of single degrees of freedom and polymer chains, and overall slow motions obtained by normal mode analysis (NMA) shows that different motions at different spatial scale may gradually stop at distinct temperature in the MD simulation temporal and spatial scales. Presumably, the small spatial scale do not contribute to the glass transition at the experimental scale since the timescale is much longer than their relaxation time.
Introduction<!>Potential energy<!>Radius of gyration<!>Dihedral rotation<!>Structural fluctuations<!>Correlation coefficient between individual degrees of freedom<!>Normal modes<!>Chain terminal effect<!>Results and Discussion<!>Potential energy<!>Radius of gyration<!>Correlation coefficients<!>Normal modes<!>End effect<!>Conclusions
<p>Glass transition is a fundamentally important phenomenon of amorphous materials and intensively studied for a long time, but still open topic. It is a widely-observed event that inorganic, organic, and polymeric materials undergo a second-order transition accompanied by various property changes such as the thermal expansion coefficient, heat capacity, bulk modulus, dielectric constant, and viscosity, without obvious conformational or structural reorganization [1]. In a cooling process, the structural relaxation time grows for orders at this temperature [2]. Usually this phenomenon is characterized by techniques such as the differential scanning calorimetry (DSC) [3], dynamic mechanical thermal analysis (DMTA) [4], and nuclear magnetic resonance (NMR) [5]. Some of the interesting but puzzling experimental observations, for example, the kinetic slowdown without obvious structural changes and the spatially heterogeneous distribution of dynamic properties that characterize a sudden loss of ergodicity during the cooling process, stimulated the theoretical research on the glass transition. Most of them are the studies of super-cooled liquids, and resulted in two main categories of theories from the thermodynamic and kinetic perspectives, respectively [6,7].</p><p>In the thermodynamic approach, the random first-order transition (RFOT) theory depicts that the entropy drives the generation and annihilation of randomly isolated glassy clusters in the supercooled phase. The glassy clusters merge with each other as the entropic driving force vanishes when the temperature decreases and approaches the ideal transition temperature [8]. The mode-coupling theory (MCT), on the other hand, starts from the structural factor and predicts the dynamics of the intermediate scattering function using the assumption of the couplings between density modes of different wavevectors, and successfully predicts the glass transition, though it overestimated the caging effect [9,10]. A locally correlated lattice model (LCL), which is based on the concept of free volume [11] that can be measured by positron annihilation [12], also found good correlation with the experimental data [13].</p><p>The kinetic approach is based on the dynamic heterogeneity and the dynamic facilitation, and assumes that the kinetic constraints are the reason for the loss of ergodicity [7]. Following this assumption, the glass transition is claimed to be a first order phase transition in space-and-time space, or space time for short [14,15]. This theory predicts that the discontinuous transition broadens into a crossover with the system size and observation time decrease and become finite. [15] and provides a plausible explanation for the widely observed broadening of glass transition in a typical molecular dynamics (MD) simulation.</p><p>With the development of computational techniques, simulation methods have also been devoted to the research of glass transition. As mentioned above, various efforts have been devoted to the theoretical investigations of the glass transition using the super-cooled liquid as model systems [15,16]. In addition to these, various properties of polymers have been examined and investigated across the glass transition temperature (Tg). The most common way to characterize glass transition used in a molecular dynamic simulation is the dilatometric measurement. Tg measured by MD was shown to deviate from the experimental values because of the fast cooling rate following the Williams-Landel-Ferry equation or equivalently the Vogel-Fulcher-Tamman equation [17]. Various physical properties, such as dihedral transition [18], correlations between motions [19], bulk modulus [20], were calculated and compared with experiments. Collective motions of the polycarbonate (PC) near Tg was found [19], where the segmental motions and the local motions may have correlations with the αand β-relaxations respectively, as suggested by an IR study [21]. A delicate study about dihedral relaxations also suggested that the α-and β-relaxations may correspond to the segmental motion which requires the collective relaxation of all dihedrals, and individual relaxations of the dihedrals, respectively [22]. The relatively easy-to-simulate β-relaxations was also shown to be the local cooperative motions of neighboring particles [23]. The co-existence of small molecules was found to decreases Tg by employing molecular dynamics [24]. The α-relaxation process was examined by simulations and compared with MCT [25]. Although with all these efforts, the glass transition has been relatively thoroughly examined by using molecular dynamics, it still lacks the investigation into the contributions from various degrees of freedom embedded in realistic polymer materials.</p><p>In this work, we tried to approach the glass transition phenomenon using all-atomistic molecular dynamics simulations, especially aiming at clarifying individual degrees of freedom that contribute to the glass transition of polycarbonate (PC) and poly-(methyl methacrylate) (PMMA). The previously published [26] and newly performed cooling simulations of these two materials were analyzed. Firstly, some basic quantifications of the cooling simulation were performed, including the variations of internal energy, radii of gyration, and dihedral transitions along the cooling process. Then the fluctuations of all degrees of freedom (DoFs), i.e. the translation and the rotation of the polymer chains, and the bond stretching, angle bending, and dihedral rotation, were obtained by groups (defined in the method section), as well as the correlations between the DoFs within the same group. From these analyses, we will show that the various degrees of freedom begin to have larger displacement, and their correlations reach separate peaks at temperatures slightly different from the glass transition temperature. Collective motions were examined by using normal mode analysis (NMA). These results imply that MD simulation captures both the segmental and local motions that are believed to contribute the αand β-relaxations, and at the time scale of an MD simulation the β-relaxation might contribute to the glass transition and broadens the glass transition temperature range because the various local degrees of freedom experience an increase of movement at different temperatures. The chain end effect was investigated, where high flexibility and large local volumes of the terminals were confirmed in both PC and PMMA. We investigated the glass transition phenomenon by examining two amorphous polymer materials, polycarbonate (PC) and poly-(methyl methacrylate) (PMMA) (Fig. 1). We utilized the trajectories previously published for the 1.05 million-atom PC model [26] which were generated using the standard cooling protocol to equilibrate the system from 640K to 200K at a rate of 20ns/20K. New MD simulations have been performed for the 1.05 million-atom PMMA model using the same protocol from 800K instead of previous 600K to 200K at the same rate of 20ns/20K to cover a larger temperature range. All trajectories during the cooling process were saved at every 100ps. Details of the MD calculations can be found in the citation [26]. The molecular weight distribution and other details of the two systems are provided in Table S1. Additional simulations performed for specific analysis purposes will be described later accordingly.</p><!><p>The internal energy change of PC and PMMA during cooling process was characterized using the potential energy, and was compared with the potential energies after energy optimization and after energy optimization preceded by scaling the volume at higher temperatures to that at the lowest 200K. The energy optimizations were performed using a limitedmemory Broyden-Fletcher-Goldfarb-Shanno (l-bfgs) minimizer [27] until the potential energy converges or reach the machine precision. The volume scaling was performed by scaling the atom positions uniformly according to the unit cell size ratio.</p><!><p>We characterized the decrease of the radii of gyration in the cooling process by comparing it to the extent of volumetric shrinking. The radii of gyration of the simulated polymer chains were set according to the experimental values measured in the melt state. In this analysis, they were computed for the two model systems at each temperature. Because the model systems were prepared according to the experimental molecular weight distributions and, thus, contain various chains with different degrees of polymerization, we calculated the average <RgAi> of the coefficient RgAi in the empirical formula Rg,i = RgAi • Mi RgB for the radius of gyration, where Mi is the molecular weight of the i th chain and RgB equals to 0.5 and 0.49 for PC and PMMA, respectively [28,29]. Furthermore, we compared <RgAi> with the unit cell dimension by compute the ratio <RgAi>/L, where L is the cubic root of the unit cell volume, for the purpose of comparison to the volumetric shrinking.</p><!><p>We computed the conformational transition rate and the autocorrelation of the pseudo dihedral rotation C1-O2-O2-C1 of PC, and the dihedral rotations C1-C2-C1-C2 (or equivalently C2-C1-C2-C1), C1-C2-C4-O2, C2-C4-O2-C5, C4-O2-C5-H5 of PMMA. For this purpose, we computed distributions of dihedral angles at all temperatures. We defined each peak in the calculated distributions as one distinct conformation. The data of dihedral distributions are provided in Fig. S1. The conformational transition rate was calculated from the percentage of dihedral angles that travel to different conformations in 100 ps. The dihedral angle autocorrelation was computed using Eq. 1 [30]. These analyses were performed for the two model systems at all available temperatures. The Kohlrausch-Williams-Watts (KWW) function [31], Eq. 2, was used to fit the autocorrelation function and the relaxation time τ was evaluated for certain temperatures in which the timescale is long enough for the fitting to be meaningful. Additional 2 nslong MD simulations saved at a frequency of 1ps from the final conformations of the trajectories at each temperature of PC and PMMA was performed for the evaluation of the dihedral transition rate at 1 ps frequency and the results are provided in SI.</p><p>Eq. 2</p><!><p>The 3N degrees of freedom of the N atoms in the model systems were grouped into the 3M translational and 3M rotational degrees of freedom of the M polymer chains, and the 3N-6M intrachain degrees of freedom including the bond stretching, angle bending, and dihedral rotations in the molecules. The intrachain degrees of freedom were selected using a Z-matrix approach. The variances of these 3N degrees of freedom were then computed separately. Note that the rotational and the dihedral degrees of freedom were represented by two variables sin(θ) and cos(θ) to account for the periodicity [32]. Finally, the group sums of the variances were obtained, where the groups were defined as follows: the chains were divided into nine groups by their molecular weight evenly from low molecular weight to high molecular weight, and the nine groups were also used for the analyses of translational and rotational degrees of freedom; the bond, angle, and dihedral degrees of freedom were grouped by atom types.</p><!><p>We computed the correlation coefficient between degrees of freedom i and j within the groups defined in 2.4. For example, the correlation coefficients between translations of chains in group 1 were calculated and averaged as the correlation coefficient of translational DoF of group 1. Similarly, the correlation coefficients between any two of the C1-C2-C1-C2 back bone dihedrals in PMMA were calculated and the average of these correlations were taken as the correlation coefficient of this DoF. This correlation coefficient indicates the strength of the collective motion of the corresponding DoF. The correlation coefficient between degrees of freedom i and j was computed using Eq. 3 and Eq. 4, where T is the total time at each temperature, xi and xj denote the values of the degrees of freedom i and j, x i ̅ and x j ̅ their averages, and σ(xi) denotes the standard deviation of xi. The group averages were then obtained by ignoring the self-corrections in which i = j.</p><p>T Eq. 4</p><!><p>We evaluated the number of slow modes found at the lowest temperature, 200K, that persisted at higher temperatures. For this purpose, a small model system, which is referred to as PC(S), with 5.5k atoms was prepared for PC. It was equilibrated from 640K to 200K at a rate of 40ns/40K using the same settings as above. The normal mode analysis (NMA) was performed by using a translation and rotation of rigid block approach (RTB) [33], where the monomer units were treated as the rigid blocks and the three translational and three rotational degrees of freedom were included in the reduced hessian matrix. The energy-optimized final conformations (the inherent conformation) at each temperature of 200-640K were used for the NMA and then the eigenvectors of 200K were multiplied by the eigenvectors at each temperature to evaluate their similarities with the 200K modes. The first 100 slow modes at 200K were considered to exist at temperature T if their product with any modes of the first 100 or 500 slow modes, or any modes at temperature T was greater than a threshold chosen to be 0.8.</p><!><p>The atomic averages of the Voronoi volumes and the root mean square fluctuation (RMSF), and the dihedral conformational change rates, were calculated and averaged by the monomer index i which was defined as the number of its preceding monomers from the nearer one of the polymer's two chain ends, with i = 0 represents the two terminal monomers. The atomic average of the Voronoi volume was obtained using LAMMPS [34]. We didn't distinguish the elements because the volumes of the hydrogen atoms, to which the interchain spaces were partitioned, turned out to be larger than the heavy atoms by doing so and this is meaningless. The dihedral transition rates were obtained as in section 2.1. Then the ratio of the monomers' dihedral change rate to the average dihedral change rate of the corresponding chain was obtained. The RMSF of the monomers were also calculated after removing the periodic boundary condition of the polymer chains and superimposing the chains' conformations along the trajectory to the initial ones.</p><!><p>We performed analyses on the cooling simulations of three all-atomistic polymer model systems for PC and PMMA and investigated the glass transition of the polymers from a molecular perspective. From the density measurement, we determined the dilatometric Tg of these model systems, although the glass transition occurs over a span of more than 100K (Fig. 2). PC shows a glass transition regime from 320K to 520K with a Tg at 417K, which is slightly lower than the experimental value 422K; PMMA shows a glass transition regime from 360K to 580K with a Tg of 481K compared to the experimental value of 383K. As discussed previously, the large discrepancy between calculated and experimental Tg of PMMA possibly originates from the lower flexibility of the simulated PMMA at the spatial scale of ~20 nm. The analyses in this work also confirm that the simulated PMMA shows lower flexibility than PC. This broad glass transition regime compared to the relatively narrow range observed by the experiments [35][36][37] has been reported by many other simulation studies [38,39]. This can be attributed to the difference in the spatial [40,41] and temporal [42] scales between the calculations and the experiments: the dimension of the model systems are ~10 7  times smaller than the real experimental samples and the quench rates are ~10 7-8 times faster than the experimental conditions. The small scales are known to raise Tg for ~25K [43] and broaden the glass transition [15]. The density variation rates of these systems, which are similar but not identical with the thermal expansion rates, all show two plateau regions at the low temperatures and high temperatures regions. Interestingly, the low temperature regions are less flat than the high temperature regions, suggesting that some degrees of freedom may still be gradually freezing when the temperature approaches 200K from above. Because of this reason and the extra data points in the 200-300K range, the currently reported Tg of the same PC system is slightly lower than the previously reported value, 440K [26].</p><p>Figure 2. Densities (lines with squares) and density variation rates (broken lines) at each temperature.</p><!><p>We characterized the internal energy change of the PC and PMMA systems in the cooling process using the potential energy, and compared it with the potential energies after energy-optimization and after energy-optimization with the volume rescaled to the volume at 200K. We confirmed that the cooling process indeed explores the polymer's configurational space and drives the polymer to more stable states corresponding to each temperature even in the short 20 ns-time scale in our MD calculations (Fig. 3). The potential energy of PC from the MD trajectories and its decomposition by force field energy terms show mostly straight lines over the entire temperature range except the vdW term (Fig. 3 a). The vdW term, which is sensitive to the interatomic distances which is closely related to the specific volume, shows a bend at temperature near the glass transition. Because the potential energies contain a substantial portion from the thermal fluctuation of the structure which also strongly depends on the temperature and behaves as a noise, we performed energy optimization to remove it until the force converged to the machine precision and plotted the minimized energy (Fig. 3 b). By doing so, all the energy terms decrease to zero except for the vdW term. The vdW term decreases by ~7000 kcal/mol and becomes the only term that contributes to the potential energy drop along the cooling process. This confirms that the vdW interaction is optimized during the cooling process. However, this vdW interaction optimization could be a result of purely the specific volume change. To test this, we simply scaled the volume of all configurations to the volume at 200K, and performed energy minimization again (Fig. 3 c). It is interesting to find that, even in this process, the vdW energy still decreases by ~1800 kcal/mol and still contribute to a large portion of the total energy change. The behavior of PMMA is identical to that of PC. Therefore, the potential energy optimization occurring spontaneously in the cooling process consists of two aspects, 1) the system is trapped in the best possible metastable state visited within the given quench time and this metastable state becomes the glassy state below Tg and 2) the enhancement of non-bonded interaction due to the shrink of the specific volume, which may in turn alter the conformations of the metastable states at different temperatures. The efficiency of the exploration of the metastable states apparently depends on the temperature and decreases significantly at temperatures around and below Tg. For example, in glassy states at low temperature, this contributes largely to the aging process which occurs extremely slowly [44]. Therefore, the fast cooling rate in MD tends to make the calculated Tg higher [42]. By examining the relaxation of potential energy and density at each temperature in the cooling simulations, it is found that the relaxation is likely to be fast enough for the system to find its best metastable state at 640K for PC and 800K for PMMA which are higher compared to their Tg. However, the configurational exploration is not complete at temperatures near Tg within the current simulation time. Thus, further extension of the MD simulation at temperatures around Tg might provide better description of the glass transition, in particular, for PMMA.</p><!><p>We examined the radius of gyration to confirm whether the conformational changes accompanying the volumetric shrinking due to temperature decrease is affine or not. If it is affine, the radius of gyration should decrease at the same rate of the lateral dimension of the unit cell. The result shows that the radius of gyration does shrink during the cooling process, but the <RgAi>/L ratio increases except at the higher temperatures where the random conformational changes are severe (Fig. 4). This suggests that in the simulation the radius of gyration has a tendency to relax towards some equilibrium value, which is determined by the experimental value together with some bias from the force field in use. This tendency of the radius of gyration relaxation is probably one of the main driving forces of relaxation of the stress generated due to the volumetric shrinking when the temperature is lowered, which may include other factors such as the local environment optimization of the bonded terms as indicated in Fig. 3C.  The 100ps frequency conformational transitions of the major dihedral angles in the systems slowly cease around Tg. The main chain dihedral angle rotation C1-O2-O2-C1 of the PC system starts to increase at around 320K and reaches the asymptotic increase rate at around 500K, with an intermediate value at around 410K, which is close to the dilatometric Tg. The behavior of PMMA is similar to that of PC: the main chain dihedral C1-C2-C1-C2 and the side chain dihedral C1-C2-C4-O2 both show the increase at around 530K, except that the side chain dihedral transition does not vanish at lower temperatures. These behaviors agree with those found in the previous studies, and it could be shown that the activation energy of the dihedral rotations obtained from the Arrhenius diagram varies at different temperatures and peaks around the Tg [18,38]. These indicate that the dihedral rotation gradually ceases at temperatures around the Tg, and the transition cessation does not occur suddenly at a specific temperature. This smooth transition might be exaggerated by the fast quench in MD simulation but it does mean that no phase transition should occur at this temperature range [6,45,46]. The transition rates measured every 1ps (Fig. S2) show two major differences from those found in the 100 ps measurement: the percentage of transitions is small in the 1ps measurement, and the temperatures at which the transition rate decreases shifts to higher temperatures, especially for the methyl rotation. Because the aim of the present analysis is not to quantify Tg using the dihedral rotation, the identical trend found in the two measurements is satisfactory. This gradual and smooth change was also observed in the autocorrelation functions of the dihedral angles (Fig. 6). The autocorrelation functions of all dihedral rotations show that the relaxation becomes slower with decreasing temperature suggesting that free energy barriers of the dihedral rotations are overcome gradually with the constantly increasing thermal fluctuation. This also agrees with the previous study finding that the dihedral rotation barrier changes along with the temperature [18,38]. The lengths of the simulations are not long enough for the fitting using Kohlrausch-Williams-Watts (KWW) function [31] at all temperatures. So the evaluations were only performed at certain temperatures. By interpolation of the fitting, the temperatures at which the relaxation time, τ, equals to 20 ns, are roughly 444K, 541K, 498K, for the PC C1-O2-O2-C1, PMMA C1-C2-C1-C2 (or equivalently C2-C1-C2-C1), and C1-C2-C4-O2 dihedrals, respectively. These values suggest that the PC main chain rotation is more flexible than that of PMMA, which agrees with the comparison between the calculated Tgs of PC and PMMA and the larger discrepancy between the calculated and experimental Tg of PMMA. At these temperatures, the present polymer systems are all in the melt state, and therefore the 20 ns-long simulations cover both the αand β-relaxation processes [22]. Below these temperatures, the systems begin to pass through the glass transition regime, where the exponentially increasing relaxation time becomes longer than the simulation time. The fluctuations, or the variances, of individual degrees of freedom indicate that these degrees of freedom experience the increase in the relevant molecular motions at different temperature ranges in a heating process. In the case of PC, the fluctuation of all the polymer translational degrees of freedom of the nine groups of molecules begin to increase around 450K. The motion of centers of mass of the molecules with smaller molecular weights increases more than those with the larger molecular weights. Similar behaviors were observed in the rotational degrees of freedom while the increase starts at around 435K. Only two angle bending terms, the C4-C5-C6 angle and C3-C4-C5 angle, experience an abrupt increase at the temperature range 400-480K and 340K-440K, respectively, while the remaining angles show a linear increasing trend. Discrepancy found between the temperature ranges is more apparent in the case of dihedral angles. The three different main chain dihedral angles, C4-C5 rotation, C1-O2 rotation, and O2-C7 rotations, experience abrupt increase in the fluctuations roughly at the temperature ranges, 300-500K, 320-480K, and 340-540K, respectively. Similarly in the case of PMMA, the translational and rotational degrees of freedom start to increase at 540K and 470K, respectively. The main chain angle bending and two other angle bending terms involving the main chain/side chain angles experience the abrupt increase at the range of 470-560K, and the main chain and side chain dihedral rotations roughly at the range of 450-600K. The fluctuations indicate the extent of their motions and are the origins of the volumetric expansion in Fig. 2, and therefore should correlate with the dilatometric Tg. This is true in the translational, rotational, bending, and dihedral degrees of freedom, despite that the translational degrees of freedom start to experience a large boost in the motion only 10-20K higher than the corresponding dilatometric Tg but still within the glass transition regime. The fact that different DoFs thaw at different temperature raises the possibility that these degrees of freedom all experience their own individual glass transition at slightly different temperatures, at the timescale of MD simulations. To examine such possibility, we continue to investigate the correlation coefficients of these degrees of freedom.</p><!><p>Fig. 8 presents the calculated correlation coefficients between the DoFs within predefined groups which were defined as follow: the polymer chains were divided evenly into nine groups from low molecular weight to high molecular weight for the translational and rotational DoFs; the bond, angle, and dihedral DoFs were grouped by atom types. The correlation coefficient is basically the normalized cross correlation with a time lag of zero, and it indicates the strength of the concerted motion between the two degrees of freedom. Certain correlations show peaks at different temperatures. In the case of PC, well-defined local maxima appear in the molecular rotations, angles, and dihedral rotations. In contrast, with respect to the translational degrees of freedom, the correlation between polymer chains become significantly strong at low temperatures [47]. Such strong correlation that equals to nearly 1.0 suggests a perfectly concerted motion between these polymer chains in our 20 ns simulations. The rotational degrees of freedom reach the maxima at 500K to 520K depending on the molecular weights, with the shorter chains maximize at lower temperatures. Two prominent angles, C3-C4-C5 and C4-C5-C6, peak at 415K and 465K, respectively. These two angles describe the relative motion between the phenylene groups and the methyl groups and represent the extent of the bending of PC's main chain. Two other angles, C4-C5-C4 and C2-C1-O2, which also govern the bending of the main chain, show almost negligible humps at 395K and 435K, respectively. The scenario of dihedral rotation is more complicated: three types of main chain rotations, C1-O2, C4-C5, and O2-C7, peak at 410K, 425K, and 465K, respectively; the methyl group rotation reaches its maximum at a low temperature, 240K; the C1-C2 rotation on the phenylene group shows a minor protuberance at 400K. Similarly, the correlation coefficients show peaks at various temperatures in the case of PMMA: the translational degrees of freedom show a broad and, thus, unclear bump at 680K, despite of the noise at low temperatures; the rotational ones show peaks from 560K to 620K, again, depending on the molecular weight; the maxima of the angles, C1-C2-C1 (C2-C1-C2), C1-C2-C3, and C1-C2-C4, are at 540K, 550K, and 540K, respectively; the dihedral rotations, the main chain C1-C2, and the sidechain C2-C4 show two almost overlapping peaks at 540K and 550K, respectively; the methyl rotation C2-C4 maximizes at 340K; the main chain rotation C1-C2 involving the hydrogen atom peaks at 275K. Note that the bond stretch of both PC and PMMA does not show any observable protrusions at the entire temperature range in the current calculation and thus the data are not shown. The peak temperatures of these degrees of freedom in both PC and PMMA lies on a broad range covering the entire glass transition temperature regime observed in the measurement by MD simulations. The above results show that the concerted motion between different ones of the same type of degrees of freedom reaches a maximum along the increase of temperature around the glass transition regime. This is because the molecular motion is too slow to show any concerted motion at lower temperatures and too fast such that the thermal fluctuation overwhelms the converted motion at higher temperatures. Therefore, the peak is an indication of the defreezing temperature of the degree of freedom under consideration. Interesting distinct types of degrees of freedom starts to experience the glass transition at slightly different temperatures. This suggests that these local degrees of freedom might have their own intrinsic temperature, which is equivalent to the macroscopic glass transition temperature.</p><p>With these computational observations, it is natural to reach the consideration as follows. The broadening of glass transition observed in the MD simulations is a result of two origins. The first is the broadening of glass transition of each individual degrees of freedom to the short observation time as discussed in the citation [14]. The second is the wide range of Tgs of the translational, rotational, various angles, and various dihedral rotations, which contribute to the over glass transition at different temperatures, as discussed above. In the timescale of hundreds of nanoseconds, the translational and rotational degrees of freedom of the polymer chains indeed play a role in the overall glass transition, but the internal degrees of freedom, such as main chain and some of the side chain angles and dihedrals, which can be naturally categorized as local motions corresponding to the β-relaxation, also contribute to a substantial portion of the observed glass transition via volumetric changes.</p><p>As the timescale approaches the experimental conditions, presumably the contributions from multiple degrees of freedom become weak or diverge from the Tg in common sense. The concerted motion of local DoFs, for example, the dihedral rotation in Fig S3, occurs on a short timescale of pico-or nanoseconds, even at low temperature in the glassy state. It is expected that such local DoFs, which correspond to the βrelaxation, are fully relaxed at experimental timescales. Under those conditions, β-relaxation is no longer a temporal obstacle for finding the best meta-state for the local DoFs, and therefore has little influence on the overall glass transition. For the above reasons, the experimentally measured glass transition, for example, via the differential scanning calorimetry (DSC) usually only exhibit a much narrower broadening of about 15K-20K [35][36][37]. In sum, an MD simulation differs from an experiment mainly in the sense that the short timescale available in MD amplifies the contribution to Tg from the βrelaxation of local DoFs, which requires additional attention when interpreting the data. However, this does not undermine MD as a useful tool for the study of glass transition of polymers because the transition due to collective motion is included, although at slightly higher temperature.</p><!><p>We further investigated the normal modes because the individual degrees of freedom lack the knowledge of the collective motions. Since normal mode analysis (NMA) is intractable for large systems with more than one million atoms, we performed the cooling simulation for a PC model with 54887 atoms referred to as PC(S) and performed NMA using a RTB approach. The PC(S) models shows a glass transition regime from 340K to 480K with a Tg at 393K (Fig. S4). We chose the first 100 slow modes, and quantified the number of these modes that survive at higher temperatures (Fig. 9). Note that these global slow modes are different from the local modes that may correspond to the spatial heterogeneity [48]. These 100 modes already start to be gradually replaced by higher temperature modes at the next temperature, 220K. This process is accelerated at 360K, where the density increment becomes larger (Fig. 2). Finally, many of these modes disappear at around 500K, where the glass transition regime ends. It means that 80-90% of the 100 slow normal modes found in the 200K glassy state persist at room temperature and are gradually replaced by the high temperature modes during the glass transition regime. At the temperatures of fully melt state, none of these modes survives. This is an indication that the global slow modes could survive through the timescale of the current MD simulations over a large temperature span.</p><!><p>The flexibility of the polymer chain ends was found to be higher than the middle part of the chain by investigating the average atomic RMSF, average atomic Voronoi volume, and main chain dihedral transition rate ratio of monomers as a function of the monomer index defined as the number of the preceding monomers from the nearer chain ends (Fig. 10). Note that the x-axis values of PC and PMMA are chosen so that the monomer indices represent roughly similar stretched length where the chain end effect completely disappears in Fig. 10. The average atomic RMSFs of the polymer end monomers are systematically higher than those of the monomers in the middle of the polymer chains. The number of monomers at the chain ends with larger RMSFs than the average decreases with the decreasing temperature. For example, in PMMA this number is ~18 at 800K and only ~4 at 600K. Such effects exist even at temperatures below Tg. The effect of chain ends is stronger in PMMA because the PMMA's monomer is smaller than that of PC and is more coiled and compacted in space. The dihedral rotations at the ends also tend to undergo more frequent transitions. In the case of PMMA, the main chain dihedral angles in the end monomers have systematically larger transition rates than the monomers in the middles of the chains, and the end effect is more pronounced at lower temperatures approaching the glass transition. This is because the flexibility of the ends is less reduced than the middle parts when the temperature decreases. The data below Tg is too noisy due to the low frequency of transition and not shown in Fig. 10. In the case of PC, only the first end monomers have higher dihedral transition rates than the middles, again because the large size of the PC's monomers. This effect is not amplified by the lower temperature, probably because of the internal rigidity of the PC monomer together with its large size. The middle monomers show higher dihedral transition rates because the long PC chains tend to undergo more conformational changes than the shorter chains (Fig. S5), and the middle monomers are only averaged on these longer chains. Such enhanced dihedral rotation by chain ends is consistent with previous observations [18].</p><p>The average atomic Voronoi volume also shows similar trend. In the case of PMMA, the volume of monomers at the chain end is also larger than that found in the middle of the chain. In the case of PC, such effect exists in the melt state but much smaller than PMMA. The effect vanishes in the glassy state, again, because of the large size of the PC monomer. The large volume at the chain ends correlates with the higher mobility of the chain end monomers and is consistent with previous united atom study [49]. These results suggest that the chain ends have higher flexibility enhanced by the larger local volume than the middles, especially in the melt state. However, this may not be interpreted as the evidence of the chain end effects on Tg, because the correlation was shown to be weak or negative [50,51].</p><!><p>Microscopic picture of the glass transition was investigated by performing all-atomistic MD simulations for two amorphous glassy polymer materials, PC and PMMA. The volumetric change that has usually been used to characterize the glass transition temperature in MD simulations was found to consist of two key factors, the short-time simple thermal expansion and the long-time conformational optimization. As has been previously reported, the conformations cease to change below the glass transition temperature, and this change is smooth and gradual. Different degrees of freedoms, namely the translational, rotational, angle bending, and dihedral rotation, start to be activated at distinct temperatures, and experience their own characteristic glass transition at these temperatures. In the submicrosecond level MD simulation timescale, all the abovementioned degrees of freedom play significant role in the overall glass transition. Therefore the glass transition regime is broadened by the individual types of degrees of freedom that begin to thaw at distinct temperatures and the broadening effect from the short observation time. The slow normal modes found in the glassy state gradually disappear during the glass transition and vanishes completely in the melt state. These results give an overall picture of the glass transition in MD simulation that the motions ranging from the scale of a couple of nanometers to a single degree of freedom contributes to the glass transition and the component from the local motions which might correspond to the β-relaxation broaden the overall glass transition; this contrasts that at the experimental scale, where the contribution from the local degrees of freedom is absent due to a longer timescale. The chain ends have high mobility enhanced by larger local volumes.</p>