Article pubs.acs.org/Macromolecules
Atomistic Simulation of the Structure and Mechanics of a Semicrystalline Polyether Nikolaos Lempesis,† Pieter J. in ‘t Veld,‡ and Gregory C. Rutledge*,† †
Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139, United States ‡ Polymer Physics, BASF SE, 38 Carl-Bosch Street, Ludwigshafen D-67056, Germany S Supporting Information *
ABSTRACT: We report the use of atomistic simulation to study semicrystalline poly(tetramethylene oxide) (PTMO), which is one of the major components of thermoplastic polyurethanes. This work reports the first application of an Interphase Monte Carlo model previously developed for polyethylene to a more complex chemistry involving heteroatoms, about which much less is known experimentally. The interface between the crystalline and amorphous domains of PTMO has been modeled in detail, complete with the equilibrium distributions of tails, loops and bridges. In doing so, a criterion has been established for selecting the relevant interface between domains, and a methodology developed that identifies the energetically most favorable interface in a heterogeneous material. A representative sample of configurations was then simulated by molecular dynamics, and analysis of deformation to small strains at different strain rates is described. Estimation of the full stiffness matrix of semicrystalline PTMO is reported for the first time. neutron scattering,12 and calorimetry,13 indicates that the interphase between ordered and disorder regions is of nanometer dimensions and may contain a variety of topological structures; nevertheless, none of these techniques has revealed how such structures are organized within the interphase or which combinations are responsible for the material’s response to deformation. Theoretical, multiscale simulation approaches can be a good starting point to examine the behavior of these materials across different length and time scales. Molecular simulation offers a complementary approach to gain detailed, mechanistic insight into the structure and response of complex, polymeric materials. Of particular relevance to the current work, we have previously reported the development of a thermodynamically rigorous interphase Monte Carlo (IMC) method to describe the topology and packing of chain segments within the noncrystalline region sandwiched between two semi-infinite lamellar crystallites, representative of a one-dimensional lamellar stack motif, and its application to the study of relatively simple, semicrystalline hydrocarbons such as polyethylene (PE)14−24 and isotactic polypropylene (iPP).25,26 Structural,14−16,21 thermomechanical,20 and interfacial properties18,19,27 for PE have been evaluated by this method and validated against the extensive experimental record for this prototypical polymer. More recently, the method has been shown to provide a good
1. INTRODUCTION Crystallizable polymer melts rarely crystallize completely. Instead, such materials are generally found to be semicrystalline, characterized by small crystallites interspersed with noncrystalline, amorphous regions. Consequently, the thermodynamic, dynamic, and mechanical properties of a semicrystalline polymeric material are sensitive to the amount and properties of the different domains it comprises, as well as to the interphase1 transition zones through which they interact. The crystal domains in these materials are characterized by long-range crystalline order, whereas the adjacent amorphous regions are more typically characterized as entangled, elastic networks. Hence, semicrystalline polymers are structurally and mechanically heterogeneous. The coexistence and interaction of the two structurally different domains in response to external stimuli is responsible for their useful properties such as toughness and strength. Such semicrystalline materials are of interest both scientifically and technologically, and represent the majority of plastics in general use today. Ever since the seminal studies of Keller, Fisher, and Dill2−4 that revealed the chain-folded nature of the semicrystalline morphology in linear polymers,5 various experimental techniques have been used to shed light on the structure−property relationships that govern these materials. Despite these decades of effort, important details regarding the structural organization and mechanistic understanding of the properties of these materials remain unresolved. In particular, experimental evidence from various sources, including X-ray scattering,6,7 nuclear magnetic resonance,8,9 Raman spectroscopy,6,10,11 © XXXX American Chemical Society
Received: March 18, 2016 Revised: July 8, 2016
A
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules Table 1. TraPPE−UA Force Field Parameters and Functional Forms Used in This Worka
a The bond stretching potential parameters were taken from ref 43. Where equivalent, CH3 and CH2 groups are denoted by CHn. The Lorentz− Berthelot mixing rules were used for Lennard-Jones interactions between dissimilar sites.
This report is organized as follows. Section 2 describes the model employed to simulate each system and is divided into three subsections, each of which describes one of the three systems: PTMO melt (section 2.1), PTMO crystal (section 2.2) and finally semicrystalline PTMO (section 2.3). In section 3, nonequilibrium deformation simulations are considered for both crystalline and semicrystalline PTMO. The calculations are described in section 3.1, and the corresponding findings are reported, compared and discussed in detail in section 3.2. Finally, section 4 summarizes the basic conclusions of this work and compares them to analogous findings for chemically related crystalline and semicrystalline systems.
starting point for molecular dynamics simulations of the semicrystalline lamellar stack, permitting detailed examination of the response mechanisms underlying large strain deformations in extension, compression, and shear.22−24 The onedimensional lamellar stack representation, consisting of alternating crystalline and noncrystalline domains, is the simplest motif to represent the structure of a semicrystalline polymer; more complex morphologies such as spherulites are observed in semicrystalline polymers, but these may be envisioned to be built up by organizing such one-dimensional alternating motifs in three-dimensional space. In this work, we extend the methods developed on these simpler, chemically homogeneous semicrystalline polymers to a more complicated, chemically heterogeneous system for which there is significantly less information available in the experimental record. As a result, more attention must be given to construction of the model from basic principles. We focus our attention on semicrystalline poly(tetramethylene oxide) (PTMO), also known as poly(tetrahydrofuran), poly(butylene oxide), poly(oxytetramethylene), or sometimes poly(tetrahydrofuran). Our approach to understand better the structure−property relationships of semicrystalline PTMO entails the development and application of theoretical and computational modeling at the atomistic length scale, where the chemistry, topological structure, crystal lattice geometry, and excluded volume interactions are all explicitly taken into account. In preparation for simulation of the semicrystalline PTMO system, we first examine a pair of auxiliary systems: the melt and the perfect crystal. Construction of the semicrystalline system is then based on the information collected from the auxiliary crystal and melt simulations.
2. SYSTEMS AND METHODS 2.1. PTMO Melt. The PTMO melt was studied first, for the purpose of validating a suitable force field for this material. Validation requires that there exist ample, unambiguous thermodynamic data for the material in question. Fortunately, there is a relatively large amount of experimentally measured properties for molten PTMO in the literature. For this reason, simulation of the melt was a good starting point for the selection of a suitable force field. An ensemble of ten independent PTMO melt configurations was constructed at T = 453 K and P = 1 atm using the Enhanced Monte Carlo (EMC) software developed by In ’t Veld and Rutledge.21,28 These conditions corresponded to a temperature substantially higher than the melting temperature of PTMO at atmospheric pressure (Tm = 316 K).29−31 Each one of the ten configurations contained five PTMO chains with 555 repeat units apiece. A repeat unit consisted of the sequence − (CH2CH2OCH2CH2)−, where every oxygen and methylene B
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules group was treated as an atom or united atom, respectively. Thus, every polymer chain had a total of 2775 sites, corresponding to a number-average molecular weight of Mn = 39960 g/mol, and the total number of sites in the system was 13875. The specifications of the system (number of atoms, number of chains) were chosen such that the resulting numberaverage molecular weight Mn resembled as much as possible that of a system investigated experimentally by Tsujita al.,32 in order to eliminate molecular weight dependencies of properties.33−35 The simulation box was cubic, with dimensions a = b = c = 72.50 Å and α = β = γ = 90°, in order to reproduce the experimentally measured density at this temperature (ρ = 0.87 g/cm3).32 Each of the 10 initial configurations of the PTMO melt was simulated by molecular dynamics (MD) using the LAMMPS MD simulator36 in the NPT ensemble with P = 1 atm and temperatures ranging from 453 K down to 313 K in increments of 20 K. The equations of motion were integrated using the velocity−Verlet method,37 with a 1 fs time step. The deterministic Nosé−Hoover thermostat and barostat, with constants of 0.1 and 1 ps, respectively, were used to maintain isothermal and isobaric conditions.38,39 Each simulation ran for a total of 30 ns, the last 25 ns of which were sampled every 500 fs to compute ensemble averages. A slight modification of the TraPPE−united atom (UA)40,41 force field was tested for the quantification of the inter- and intramolecular interactions. Recently, this force field has been successfully implemented for the study of polyols and polyethers.42,43 Whereas the original TraPPE−UA force field considers a stiff bond stretching potential, in the current work a harmonic bond stretch potential was chosen for simulations of melt PTMO. The parameters for the harmonic bond stretching potential term listed in Table 1 were taken from a slight alteration of the OPLS-UA force field previously tested by Chen and co-workers for polyethers.44 The force field parameters used in this work, along with the corresponding functional forms, are listed in Table 1. Figure 1 compares the volumetric properties of the PTMO melt obtained by simulation to the available experimental values.32 Figure 1a shows the specific volume as a function of temperature. The agreement between experiment and simulation results is gratifying. Figure 1b shows a comparison between the calculated thermal expansion coefficient and the corresponding experimental measurements at different temperatures.32 The simulation points of Figure 1b were calculated from the local slopes of Figure 1a, using a 5-point estimator of the derivative at each temperature. The simulated isothermal compressibility was calculated from volume fluctuations45 using the following equation: βT =
1 ⎛ ∂⟨V ⟩ ⎞ 1 ⎜ ⎟ = (⟨V 2⟩ − ⟨V ⟩2 ) ⟨V ⟩ ⎝ ∂P ⎠T ⟨V ⟩kBT
Figure 1. (a) Comparison between simulated volumetric behavior (diamonds) and experimental measurements (squares), with data taken from ref 32. (b) Comparison of thermal expansion coefficient (diamonds) to experimental measurements (squares) for melt PTMO at P = 1 atm.
Figure 2. Temperature dependence of isothermal compressibility calculated via volume fluctuations45 for PTMO melt at P = 1 atm. The experimental points were taken from ref 32.
Figure 3, along with corresponding experimental values available in the literature.46 The agreement between experiment
(1)
The computed values were plotted as a function of temperature for comparison to the experimental values, as shown in Figure 2. Finally, heat capacities at constant pressure were calculated from energy fluctuations, using Cp =
1 (⟨E2⟩ − ⟨E⟩2 ), where E = U + K + PV kBT 2
Figure 3. Comparison between calculated heat capacities under constant pressure (rhombuses) to experimental data (squares) with data from ref 46.
and simulation results is very good for high temperatures, but begins to deviate in the close vicinity of the melting temperature, where the heat capacity diverges. Taken all together, these comparisons for PVT behavior and isobaric heat
(2)
In eq 2, E represents the total, U the potential, and K the kinetic energy of the system. The obtained results are plotted in C
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules capacity serve to validate the TraPPE−UA force field for simulation of the equilibrium state of molten PTMO. 2.2. PTMO Crystal. In this section, we examine the suitability of the TraPPE−UA force field to simulate the structure and properties of the PTMO perfect crystal. This step is necessary since the semicrystalline PTMO contains both crystalline and noncrystalline domains, and it is important to confirm that the force field correctly describes the thermodynamics of each phase separately, before proceeding to simulations of the heterogeneous semicrystalline state. The perfect PTMO crystal was constructed using the polymer builder of the EMC code21,28 based on the size and symmetry of the crystallographic unit cell and the coordinates of the asymmetric unit. In the PTMO crystal, the chains adopt a linear, planar all-trans conformation.47,48 From investigation of oriented PTMO samples, Imada and co-workers49,50 proposed a monoclinic unit cell. The monoclinic crystal system was subsequently confirmed by IR and Raman spectroscopies.49−51 Two other research groups have also reaffirmed the monoclinic symmetry of the PTMO unit cell, albeit with slightly different lattice parameters.29,52 Vainshtein and co-workers51 subsequently proposed a rhombic unit cell with the oxygen atoms all in register but this alternate unit cell has not been confirmed independently. In this work, we consider the crystal structure of PTMO having monoclinic symmetry with lattice parameters reported by Imada et al.49 The lattice parameters are a = 0.559 nm, b = 0.89 nm, c = 1.207 nm, α = γ = 90° and β = 134.2°. There are two planar zigzag chains running parallel to the c crystallographic axis. The zigzag planes of the molecular chains are parallel to the bc plane, and the second chain is displaced by (a + b)/2 relative to the first chain. Figure 4 shows the monoclinic unit cell in three different perspectives. The fractional coordinates of the atoms in the
conformer and the bc-crystallographic plane is zero for both chains. This indicates that both chains have the same alignment and have molecular planes parallel to each other. This observation is in accordance with experimental findings.48,49 Interestingly, the angle β of the monoclinic unit cell is such that the electronegative oxygen atoms of one chain are nearly equidistant (in projection along the c-axis) between those of its neighbor in the a-direction. To create the initial configuration of the perfect PTMO crystal for simulation, the unit cell of Figure 4 was replicated in the a, b, and c directions, for a total of (8 × 8 × 15) unit cells. This resulted in a simulation box of total size 4.472 × 7.12 × 18.105 nm, containing 128 chains of 150 sites each, for a total of 19200 sites in the simulation box. Periodic boundaries were applied in all three directions, and all of the chains were bonded with their images in the c-direction. Three molecular dynamics simulations were run at P = 1 atm and T = 300 K, below the melting point but above the glass transition temperature of PTMO (Tg = 187 K),31,53−57 for a total of 5 ns, since the crystals equilibrated relatively quickly. All other specifications for the MD run of the crystal were identical to those of the melt. After equilibration of the three crystal configurations, the average density and lattice parameters of the simulated crystal phase were calculated and compared to the experimental values for the PTMO crystal.33,48,51 The results are shown in Table 2. Table 2. Comparison between Experimental and Equilibrated Simulated Crystal Density and Lattice Parameters for Crystal PTMO at T = 300 K and P = 1 atma 3
crystal density (g/cm ) a (nm) b (nm) c (nm) α (deg) β (deg) γ (deg) a
simulation
experiment
ref
1.17 ± 0.03 0.551 ± 0.13 0.879 ± 0.16 1.215 ± 0.2 88.9 ± 1.6 133.1 ± 2.1 90.3 ± 1.3
1.15 ± 0.03 0.559 0.89 1.207 90 134.2 90
32, 47, and 48 49 49 49 49 49 49
Uncertainty in all simulated results is about 2%.
Agreement is generally quite good, confirming the utility of the modified TraPPE−UA force field for the crystal phase as well as for the melt of PTMO. The uncertainty in the experimental density reflects the deviation in values reported by different groups, using different evaluation methods. To make sure that the considered force field captures well not only the structure, but also the energetics of the perfect PTMO crystal, heat capacities under constant pressure were calculated using eq 2 for a set of temperatures in the vicinity of Tm for the three equilibrated configurations of crystal PTMO. The calculated average heat capacities are compared with the corresponding experimental values46 in Figure 5. Here also, the agreement is quite good. We note that the agreement observed here is justified at least in part by the use of a united atom model for PTMO, and by the temperature range considered. The exclusion of hydrogen atoms from explicit treatment in the UA model avoids the errors that would arise from the treatment of such high frequency vibrational modes as classical oscillators by MD. Additionally, at temperatures close to Tm, anharmonicities become important;58 such anharmonicities are captured better when treated classically by MD than by the
Figure 4. Projections along the (a) a-, (b) b-, and (c) ccrystallographic axis of the unit cell structure of crystalline PTMO. In all parts, united atoms are represented with spheres. Cyan spheres correspond to methylene groups, and red spheres describe oxygen atoms.
PTMO unit cell are given in Appendix A of the Supporting Information. Figure 4a shows a projection of the unit cell on the bc-plane, whereas parts b and c of Figure 4 show projections of the unit cell on the ac- and ab-planes, respectively. The setting angle between the plane of the zigzag D
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
−(C βH2 C αH2 OCα H2 Cβ H2 )−) of PTMO was preserved throughout the simulation. This was accomplished by cutting chains only between the β-methylenes when excising sites from the simulation, so that every tail terminates with the −(CβH2CαH2OCαH2CβH3) sequence, and then allowing only those end-bridging moves that attack a neighboring segment (loop or bridge) at the bond between the β-methylenes. Unlike previous simulations of polyethylene and polypropylene, the reptation move was not used in this work because the success rate was very low, due to the size of the repeat unit insertion required to maintain the correct chemistry. Upon equilibration in the IMC method, Boltzmann distributions of loops, bridges, and tails of various lengths were obtained in the noncrsytalline region of the simulation, and a nonuniform density profile was obtained that includes the transitional, or interphase, region between crystalline and fully amorphized domains. To estimate the density of the amorphous region, and thereby determine the number of sites to excise prior to amorphization, we extrapolated the simulated values for the melt in Figure 1a to the desired temperature T < Tm. The rationale for doing so is the experimental observation that conformers in the melt and in the amorphous domain in a semicrsytalline polymer are similar.60 The extrapolated density of amorphous PTMO at 300 K was ρamorphous = 0.98 g/cm3. Thus, this value was the target density for the amorphous domain, not including the transition zones. The details of the semicrystalline PTMO simulation, then, were as follows. The initial simulation box was identical to that used to simulate the perfect PTMO crystal, as described in the preceding section. The two fixed crystal regions each had the dimensions 4.472 × 7.12 × 3.5 nm. Thus, the domain subjected to the MC moves, located in the middle of the box, had the dimensions 4.472 × 7.12 × 11.105 nm, wherein the gradients of the transition zone (interface) between the two neighboring domains lay in the z-dimension. Next, a total of 1000 sites were removed from the central portion of the polymer crystal in increments of 100 sites, in a series of NNeVT Monte Carlo simulation runs. Amorphization of the middle region after removal of each increment was achieved using the IMC algorithm and allowing the MC moves described previously to populate the system with tails, loops and bridges. An ensemble of 5 independent semicrystalline configurations was equilibrated with the help of MD simulations in the NNePT ensemble for 20 ns, using a time step of 2 fs. During the course of the MD simulations, the density of the system was monitored. The average density in the amorphous region was then measured, where caution was taken to exclude the transition regions. The result is shown in Figure 6. In this way, we determined that removal of about 900 united atoms obtains the desired value of ρamorphous. 2.3.2. Selection of Interface Orientation. A particularly complicated issue in the construction of a semicrystalline polymer model is the selection of which crystallographic facet best approximates the plane parallel to the transition zone between the crystalline and amorphous domains. Different choices of this facet correspond to different angles of tilt between the chain axes in the crystal and the plane of transition, which in turn has been shown to result in different distributions of loop and tail length.27 Unlike semicrystalline polyethylene, where this interface selection has been examined both experimentally and computationally,16,18,19,25,59 there is no information related to it in the case of PTMO, to our
Figure 5. Comparison between calculated heat capacities under constant pressure (rhombuses) to experimental data (squares) with data from ref 46.
quasiharmonic approximation, in which thermal motion is treated as a collection of essentially independent oscillators. 2.3. Semicrystalline PTMO. 2.3.1. Construction of the Simulation. Atomistic configurations of semicrystalline PTMO were created using the interphase Monte Carlo (IMC) method, as previously reported in considerable detail for polyethylene and polypropylene.14−17,20,21,27,59 In this method, one starts with a perfect crystalline configuration, in some instances rotated so that a desired crystallographic plane is placed parallel to the xy-plane of the simulation box; that plane ultimately becomes the plane through which the transition from crystalline to amorphous domains is realized. As a first guess, the crystal phase in this work was constructed such that the chain axis, which is parallel to the c-axis of the monoclinic lattice, was also parallel to the z-axis of the simulation box. Layers of sites were “frozen” in their crystallographic positions at the top and bottom (i.e., in the z-direction) of the simulation, and then a number of sites were removed from the layer of material in between, so that the density of the middle region approached that of the amorphous melt; 80 chain end sites (Ne) were introduced to cap the segments where sites were excised. Next, an NNeVT Monte Carlo simulation with simulated annealing was used to amorphize and equilibrate the noncrystalline material in the middle region of the simulation. For this purpose the following Monte Carlo moves were used: (a) single-site displacement, wherein a single site is displaced from its current position, without altering any bonded connections, (b) end-rotation, wherein rotation of one to three sites at a chain end is attempted, (c) rebridging move (or concerted-rotation), wherein three consecutive sites in a chain are excised and reinserted as a different conformer within the same chain, without changing bond lengths and angles, and finally (d) end-bridging move, where the end of a chain “attacks” a set of three sites in a nearby loop or bridge segment, resulting in a new loop or bridge and tail. The first three moves change the orientation and position of sites, without any change of connectivity or topology. In contrast to that, the fourth move enables changes in the connectivity of the chains. As a result of the end-bridging move, upon amorphization the middle region becomes populated with three types of chain segments: “bridges”, which span the noncrystalline region and connect two different lamellae, “loops”, which exit and enter the same crystal boundary, and “tails”, which originate from one crystal lamella and terminate somewhere within the noncrystalline domain. Short loops may be interpreted as “folds”, whereas short bridges are analogous to “tie chains”. Additional restrictions on the end-bridging move were implemented to ensure that the chemistry (i.e., E
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
Figure 6. Density of amorphous region as a function of the number of removed atoms (diamonds). The extrapolated (cf. Figure 1a) target value of the amorphous density at 300 K is shown as the dashed line.
knowledge. For purpose of making a reasonable determination of the interface orientation, we assume that the preferred facet is that which minimizes the interfacial energy between the crystalline and amorphous domains. One simple proxy for interfacial energy associated with a certain plane is the areal chain density at that plane, i.e. the number of polymer chain segments that intersect a unit area of the plane. In principle, an energetically favorable interface would incur the least change in chain density between the two neighboring regions of a semicrystalline polymer. It has been argued elsewhere that such a change in areal chain density (sometimes referred to as “chain flux”61) over a thin transition zone imposes certain topological constraints on the organization of the interphase,59,62−64 with concurrent entropic and/ or enthalpic costs.54,1,65,66 The advantage of such a proxy is that it can be computed up front for any crystallographic plane, without resort to simulation. It thus serves as a useful screening tool in the absence of experimental data. At this point a clarification regarding our counting algorithm is in order. Even in a system containing only one polymer chain, the backbone of that chain may cross and recross any given plane numerous times. Thus, any algorithm that simply counts backbone bonds that intersect the plane may overcount, due to corrugation of the backbone in the vicinity of the plane. In our analysis, the algorithm was formulated so that such “trivial” crossings are discounted. A crossing event is deemed “trivial” when the distance between crossing points is less than a user-defined threshold, which we took to be 0.5 nm, about 30% larger than the van der Waals length σ. The chain density was calculated for a total of 2196 crystallographic planes as specified by their Miller indices (h k l) ranging from −6 to 6, for the perfect crystal of PTMO. The number of crossing events through any defined (h k l) plane was evaluated and then divided by the area of that plane. In general, the total number of planes belonging to a family of planes {h k l} is called the multiplicity factor. The multiplicity factor may vary between 2 and 48, depending upon whether the Miller indices are different from one another or zero, and the symmetry of the crystal lattice. For a monoclinic crystal system, in the most general case, where all three Miller indices are different from one another, the multiplicity factor is equal to four.67 For a more detailed description of the chain density calculation methodology the reader is referred to Appendix B of the Supporting Information. Figure 7a shows a “thermal plot” of the areal chain density as a function of the three Miller indices for each crystallographic plane. This plot illustrates the highly anisotropic “chain flux” of the PTMO crystal, as indicated by the appearance of “warm”
Figure 7. (a) Areal chain density analysis for the case of PTMO perfect crystal at T = 300 K and P = 1 atm. The three axes of the cubic heat map describe the three Miller indices h, k ,and l ranging between −6 and 6, whereas the color bar describes the intensity of the calculated areal chain density, in crossings per Å2. The cubic area in part a corresponding to planes for which all three Miller indices are negative is rendered in dark blue. Use of the origin as a reference point in our algorithm results in zero intersections with planes of this type; these planes are ignored in the current work. Part b shows the probability distribution of sampled chain densities. dw and dun are the weighted and unweighted averages of areal chain density, respectively (see text for details).
and “cold” zones in the chain density heat map. The maximum calculated areal chain density, about 0.083 Å−2, can be attributed to the dense, highly packed chain arrangement in the crystal structure. The histogram of observed areal chain densities for the calculated sets of (h k l) planes is shown in Figure 7(b). The weighted (dw) and unweighted (dun) averages of areal chain density are also shown. The weighted average d w = ∑i Pi × cfi is weighted by the probability of appearance Pi of every encountered chain density cf i. By contrast to that, the unweighted average is simply the average of all chain densities, wherein every certain chain density value is counted only once, and thus it is located in the middle of the spanned areal chain density interval. The distance between those two averages can be deemed a measure of anisotropy; the farther apart these two quantities are from each other, the more anisotropic the system is. As Figure 7b shows, the weighted and unweighted areal chain density averages are far apart from one another in PTMO, a fact that confirms further the anisotropy of the crystal, in terms of areal chain density. The same analysis for areal chain density was implemented for the case of PTMO melt, assigning Miller indices based on the cubic simulation box. Figure 8a confirms, as expected, that the areal chain density is more isotropic in the amorphous melt. Although the maximum and minimum areal chain densities F
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
Table 3. Low Index Crystal Planes Exhibiting Areal Chain Densities Similar to the Average Chain Density of the Melt (0.0453 Å−2) and Corresponding Interfacial Energiesa crystal plane (1̅ (1̅ (1̅ (1̅ (0 (2 (2 (2
2̅ 1̅ 0 1 0 1̅ 0 1
1) 2) 1) 1) 1) 2) 1) 1)̅
chain density [Å‑2] 0.0459 0.0468 0.0464 0.0459 0.0465 0.0468 0.0463 0.0450
interfacial energy [mJ/m2] 116 126 131 109 91 115 110 107
± ± ± ± ± ± ± ±
8 8 6 5 5 7 7 5
insertion angle [deg] 73.6 40.8 63.9 67.8 44.2 76.1 81.9 100.8
a
The insertion angles of the tilted polymer chains with respect to the xy-plane are also displayed.
Table 3. Among the eight considered facets, the (001) facet exhibited the lowest interfacial energy, within the uncertainties of the simulation. The insertion angles formed by the polymer chain backbone and the plane normal to the lamellar stack direction are also shown in Table 3. On the basis of this analysis, semicrystalline PTMO systems with the interfacial region oriented parallel to the (001) crystallographic facet were used for all subsequent studies. 2.3.3. Morphological Analysis. Next, a set of 10 uncorrelated initial configurations of semicrystalline PTMO with the interfacial region parallel to the (001) plane were created using the IMC algorithm. After simulated annealing and topological equilibration by IMC, each configuration was relaxed by MD in the NNePT ensemble for 20 ns, using a time step of 2 fs. Figure 9 shows a representative snapshot of the simulation cell, in which the polymer chains in the crystalline domains that are tilted with respect to the xy-plane, as well as the different topological entities that make up the amorphous domain (bridges, tails and loops), can be distinguished. During the course of the MD simulations, the density of the system was monitored. Each simulation ran for a total of 20 ns, the last 15 ns of which were sampled every 500 fs to compute ensemble averages. The average density profile over the 10 configurations as a function of the distance from the centerplane of the crystal is shown in Figure 10. The difference in density between the crystalline and noncrystalline domains, as well as the transition zone in between them, is readily apparent. The average density of the crystalline region was ρcryst = 1.16 ± 0.01 g/cm3, similar to the value obtained for the perfect PTMO crystal. The average density in the amorphous region spanning z = [4.8, 13.31] nm, was 0.99 ± 0.01 g/cm3, similar to the target melt density of ρamorphous = 0.98 g/cm3. The probability distribution of the orientational order parameter P2,i = (3⟨cos2 θi⟩ − 1)/2 is shown in Figure 11. Angle θi is the angle between the vector from the (i − 1)th united atom to the (i + 1)th united atom and the z direction, and the angle brackets indicate averages taken over all j ≠ i pairs within a cutoff radius rij ≤ rP2 of the ith united atom. Figure 11 exhibits two peaks, one corresponding to a highly ordered region with P2,i > 0.9 (crystal domain) and the other region with P2,i ∼ 0 (amorphous domain). For intermediate values (0.6 < P2,i < 0.9), there seems to be a broad plateau, which we attribute to united atoms lying in the interfacial region where some degree of local order is maintained. It has
Figure 8. (a) Areal chain density analysis for the case of PTMO melt at T = 300 K and P = 1 atm. The three axes of the cubic heat map describe the three Miller indices h, k, and l ranging between −6 and 6, whereas the color bar describes the intensity of the calculated areal chain density, in crossings per Å2. Part b shows the probability distribution of sampled areal chain densities.
observed in the melt were comparable to those observed in the crystal, the histogram of areal chain density with direction in the melt was more sharply peaked than that in the crystal. The narrower distribution of the sampled melt chain densities indicates that unusually high and low values of areal chain density are relatively rare in the melt. A quantitative indication of isotropy in the melt is given by the close proximity of weighted and unweighted areal chain density averages, shown in Figure 8b. On the basis of the principle that equal chain flux through a plane in the crystal and an (arbitrary) plane in the amorphous melt is a good proxy for identifying the crystal plane having low interfacial energy, we searched for crystal planes exhibiting areal chain density comparable to the weighted average areal chain density (henceforth “average”) of the melt, dw = 0.0453 Å−2. Eight low index planes were thus identified and are listed in Table 3: In an effort to refine the selection of crystal orientation for the transition zone, three semicrystalline configurations were generated by the IMC method as described previously, but with the crystal rotated in each case so that one of the desired planes (cf. Table 3) aligns with the xy-plane of the simulation box. After simulated annealing and equilibration, the potential energy contribution to interfacial energy was calculated for each of these eight tabulated configurations using the method of Gautam et al.27 The results are shown in the third column of G
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
Figure 10. Average density profile of equilibrated semicrystalline PTMO as a function of the distance measured from the center plane of the crystal. The simulation results have been mirrored across the center plane of the crystal for purposes of computing the average density profile.
Figure 11. Probability distribution of the order parameter P2, calculated using as cutoff radius rP2 = 2.5σ.
Figure 12. Population distributions of bridges, loops, and tails in the amorphous region of semicrystalline PTMO. Figure 9. Semicrystalline PTMO atomistic model used in this work, viewed down the laboratory axis. United atoms within the simulation cell are colored according to the type of segments to which they belong: stems (cyan) in the crystal domain, or bridges (green), tails (blue), and loops (red) in the noncrystalline domain. Periodic boundary conditions were applied in all directions. The outline of a PTMO unit cell (red box) viewed along the b-axis, as well as the tilt angle αtilt, is shown in the upper crystalline domain.
potential μ*(Nb > Nb,min) = 0 and μideal(Nb) = Nbμbead, where μbead is constant. Furthermore, for chain segments longer than 100 beads, the populations of loops and bridges overlap, indicative of equilibration of these two populations, for which μbead is the same. From the slopes of the distributions in semilog format of Figure 12, incremental chemical potentials μbead for loops, bridges, and tails were estimated, as described by Kuppa al.,25 to be −0.0297 ± 0.0016, −0.0315 ± 0.0024, and −0.0337 ± 0.0023 kJ/mol, respectively. While there are no experimental measurements available in the literature for the chemical potential required to insert or remove a bead of PTMO, nevertheless, the reported values herein are very similar to values previously reported for iPP. This observation is important, since it confirms both the credibility of the force field, and the competence of the modified IMC algorithm to sample the PTMO interphase ensemble effectively. The
been shown previously,68 that the general shape of Figure 11 does not depend upon the cutoff radius rP2. Figure 12 shows the probability distribution of tails, loops and bridges in the amorphous region as functions of their numbers of constituent united atoms, Nb. For Nb > 100, these distributions decay exponentially, consistent with a most probable distribution of segment lengths, as defined by Krishna-Pant and Theordorou,69 having reduced chemical H
DOI: 10.1021/acs.macromol.6b00555 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
Figure 13. Stress−strain behavior of semicrystalline PTMO under simple straining (uniaxial and shear) at T = 300 K and P = 1 atm. The slopes of these plots for strains
