Predicting Raman Spectra of Condensed Polymer Phases from MD

Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

Predicting Raman Spectra of Condensed Polymer Phases from MD Simulations Qin Chen and Scott T. Milner* Department of Chemical Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States ABSTRACT: Raman spectroscopy is useful for probing conformational disorder and detecting rotator phases in nalkanes and polyethylene. Raman spectra do not give direct information about molecular conformations or phase structure, which are inferred from spectral signatures of previously identified phases. Here, we use molecular dynamics simulations to model polymer systems and predict their Raman spectra, so spectral and structural information can be obtained at once. Central to our method is the bond polarizability model, which represents molecular polarizability as a sum of individual bond contributions. We use two approaches for generating the Raman spectrum. First is normal-mode analysis, useful for identifying Raman modes in ordered phases. More generally, we compute the time-dependent polarizability and Raman spectra directly from simulation trajectories. Our results capture the main features of experimental spectra for crystalline, rotator, and melt phases of polyethylene.

■

INTRODUCTION Raman spectroscopy is a useful tool for probing molecular conformations in condensed polymer phases. In particular, it can be applied to the investigation of mesomorphic phase transitions during polymer crystallization. Strobl proposed that polymers crystallize from the melt via a mesomorphic phase before transforming to the final lowest energy crystalline state.1,2 In polyethylene (PE) and n-alkanes, this mesomorphic phase has increased rotational disorder in the chains compared to the orthorhombic crystalline structure and is fittingly named the rotator phase.3,4 Time-dependent Raman spectroscopy has been used to examine melting and crystallization transitions in PE and nalkanes, at both high pressure5 and atmospheric conditions.6−9 At high pressure, in combination with X-ray scattering to determine phase structure, characteristic Raman signatures have been identified for orthorhombic, hexagonal (identified with the rotator phase), monoclinic, and melt phases of PE.5 At atmospheric pressure, time-dependent spectra collected during heating of n-alkanes resembled spectra obtained for mesomorphic phases at high pressures. But because PE rotator phases at atmospheric pressure are only metastable, simultaneous structural characterizations by X-ray scattering are difficult to achieve. As a result, transient mesomorphic structures in PE are only tentatively identified, inferred by comparison to high pressure data. In Figure 1, we present experimental Raman spectra of Migler et al.,8 for alkanes in the orthorhombic, rotator, and melt phases at atmospheric pressure. Peaks associated with CC stretching, CH twisting, and CH bending vibrations are used as spectral signatures to identify alkane and PE phase structures. The main distinguishing feature of the orthorhombic phase is © XXXX American Chemical Society

Figure 1. Experimental Raman spectra of alkane phases at atmospheric pressure. Reproduced from Migler et al.8 Red, orthorhombic; blue, rotator, and black, melt. Dashed lines divide the x axis into regions corresponding to each type of vibrational motion.

the peak splitting observed at 1417 cm−1. Heating the crystal produces a Raman spectrum very similar to the rotator phase spectrum at high pressure.5 By this comparison, one can infer that the sample is largely populated by “consecutive trans” conformations, while the absence of peak splitting at 1417 cm−1 suggests that orthorhombic packing is no longer present. Further heating produces the melt phase, characterized by Received: June 7, 2017 Revised: November 27, 2017

A

DOI: 10.1021/acs.macromol.7b01202 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

results with published experimental data. Some previous efforts have been made to simulate Raman spectra of PE crystals21,22 and alkane liquids.33 However, there has been no comprehensive simulation study of the effect of phase structure on Raman vibrational signatures. More broadly, our approach can predict Raman spectra of proposed polymer phase structures or link an unknown spectral feature to a corresponding molecular structure. In particular, we are interested in the Raman spectral signatures of iPP rotator phases, which presently are not wellcharacterized.

broad peaks throughout the Raman spectrum, and the emergence of the amorphous CC stretching peak at 1080 cm−1. Compared to PE and n-alkanes, disordered phase structures of isotactic polypropylene (iPP) are not well-characterized. Rapid quenching of iPP produces a conformationally disordered phase, which may be a rotator phase.10−13 Other studies looked at optical14 and Raman properties of iPP spherulite cross sections15 to characterize rotator phase formation at the crystallite growth front. However, we still do not know the structure of iPP rotator phases or how that structure influences the Raman spectrum. A difficulty in using Raman spectroscopy to probe polymer phases is that Raman spectra do not give direct information about molecular conformations or phase structure. Instead, structural information is inferred from spectral signatures previously observed in samples of known structure (e.g., by Xray scattering).5,9 However, sometimes structural characterization is impossible because the samples are not stable or can only be achieved under extreme conditions. Much effort has gone into developing analytical and simulation tools to predict vibrational signatures of molecular conformations then using these predictions to interpret experimental spectra. There has been extensive work on how the fraction and distribution of trans and gauche bonds on alkanes and iPP molecules affects the shape, intensity, and frequency of Raman peaks.16−19 In contrast, less has been done to predict how bulk phase structure can affect Raman signatures. The biggest challenge to make such predictions is the computational effort required for large systems such as polymer bulk phases. In general, the Raman spectrum of a polymer can be simulated by two approaches. The first is to calculate all vibrational normal modes of the system and selecting those that are Raman-active based on molecular symmetry.20 The Raman intensity can be computed from the amplitude of molecular polarizability oscillations as the mode vibrates. The second approach is to obtain the time-varying molecular polarizability directly from a simulation trajectory and compute the power spectrum; this method has been applied only a few times, to crystalline systems.21−24 The necessary polarizability can either be obtained either from ab initio calculations25−27 or empirical models.28,29 Ab initio methods such as density functional theory (DFT) are more accurate but require vastly more computational resources. One empirical model used to approximate the molecular polarizability is the bond polarizability model.28 It assumes the polarizability of a molecule is the sum of independent contributions from its chemical bonds. The bond polarizability constants are sometimes called electro-optical parameters in literature.30−32 Values for bond polarizability constants can be obtained either by fitting to experimental Raman peak intensities of molecules containing the bonds of interest16 or by fitting molecular polarizabilities computed using DFT.30 The bond polarizability model has been used successfully to represent the molecular polarizability of hydrocarbons,16,21,22,33 small molecules,30 and silica polymorphs.34,35 In our work, we will use the bond polarizability model to compute the Ramanactivity of molecular vibrations in condensed phases. The main objective of this paper is to develop a comprehensive approach to predict Raman spectra of bulk phases from molecular dynamics simulation. In developing this method, we will simulate the Raman spectra of various condensed phases of PE and n-alkanes and compare our

■

METHODS FOR PREDICTING RAMAN SPECTRA

Normal Mode Analysis. In our work, we use two different approaches to compute Raman spectra. The first is normal-mode analysis (NMA), which starts with the set of vibrational normal modes, obtained by diagonalizing the Hessian matrix. The Hessian is a square matrix containing the second-order partial derivates of the potential energy with respect to atomic positions. We examine the Raman activity of each mode numerically by computing the time-dependent molecular polarizability along the vibrational trajectory. The thermal vibrational amplitude of each mode is determined through equipartition. At discrete time points along the trajectory, we use the bond polarizability model (discussed in the next section) to compute the polarizability tensor from the molecular configuration. It will oscillate at the same frequency as the mode, with some amplitude that determines its Raman intensity. The total Raman intensity is proportional to the signal power ϕ(ω), given by the sum of the squared Fourier amplitudes [Pαβ(ω)] from all polarizability tensor elements (eq 1).

ϕ(ω) =

∑ Pαβ2 (ω) αβ

(1)

A normal mode can be Raman inactive, if the polarizability contributions of different copies of the same moiety oscillate out of phase and cancel. (See the Appendix for a short discussion.) This can happen if the molecule has some particular symmetry, such as inversion symmetry. For such molecules, vibrational modes will be either even or odd under inversion. The polarizability is a rank-two symmetric tensor (hence even under inversion), while a dipole is a vector (hence odd under inversion). Thus, for an even mode, contributions to the oscillating polarizability of moieties related by inversion contribute in phase, while contributions to an oscillating dipole contribute out of phase and cancel. (The reverse is true for an odd mode.) Hence even modes are Raman-active and IR silent, and odd modes are the reverse. For highly symmetric molecules, group theoretical methods are useful in identifying Raman active and inactive modes. But more generally, Raman selection rules are imposed automatically by the bond polarizability model, which amounts to a constitutive relation that by construction has all the symmetry properties of the molecule itself. The normal mode approach for computing Raman spectra works well for highly ordered phases, for which harmonic normal modes are relatively well-defined. It does not work well for disordered phases such as polymer melts. For disordered phases, large system sizes are required to obtain a suitable average over local modes, and accurate local energy minima must be found so that the Hessian has only positive frequencies. More fundamentally, typical configurations of disordered phases may not be well-described as small oscillations about a minimum-energy state. Another shortcoming of the normal mode approach is that the line widths of Raman peaks are not predicted. Each vibrational mode only corresponds to a single frequency because no damping effects are considered. Often in practice, the line width of Raman spectral features is an important clue for phase identification. To compute Raman spectra with finite line widths, we must use the time-dependent polarizability. B


Article

Macromolecules Time-Dependent Polarizability. Line widths in experimental Raman spectra are a consequence of vibrational damping, which in turn arises from fluctuations in molecular conformations. In a real polymer system, vibrational motion of a molecule is perturbed by random buffeting from its neighbors, caused by thermal fluctuations. The more disordered the system, the more molecular vibrations are damped. Damping shifts the mode frequency downward and broadens the peak in the Raman spectrum. Thermal fluctuations also cause molecules to deviate from their lowest-energy configurations, so that the potential that describes the mean restoring force is no longer the energy but the free energy. This leads to shifts in the frequency of vibrational modes, which are more severe for more disordered phases and higher temperatures. Molecular dynamics simulation can capture the time-varying motion of molecules, if appropriate force fields are used. We can thus generate a time-fluctuating signal of the system polarizability directly from a molecular dynamics simulation trajectory, using the bond polarizability model. In this approach, fluctuations and damping caused by thermal noise are naturally included in the simulation trajectory and, hence, in the predicted Raman spectrum. The Raman spectrum is proportional to the time-averaged power spectrum of time-dependent system polarizability Pαβ(t). To compute it, we divide a long trajectory into multiple short time series, and for each time series compute the Fourier transform Pαβ(ω). At each frequency, the Raman intensity averaged over all orientations of the system is proportional to the sum of the squared Fourier amplitude of all tensor components (eq 1). We average the power spectrum over the multiple short trajectories. Computing Raman spectra from the time-dependent polarizability has several important advantages. First, simulated peaks now have line widths, which contain information about damping and mode lifetimes. Second, this approach can be used for disordered systems such as the rotator and melt phases, which are not well-described by small fluctuations around a minimum energy state, with a set of well-defined harmonic modes. Of course, because this approach does not separate the timedependent signal into contributions from normal modes, we cannot use it to make correspondences between spectral features and vibrational motions. Instead, we apply the normal mode approach to more ordered systems to identify spectral features with known vibrational motions. It is well-known experimentally that spectral features tend to persist in some form near the same frequency for progressively more disordered phases because the phase structure has relatively modest effects on local molecular vibrations. Hence we can identify peaks in the Raman spectrum obtained for disordered phases using the time-dependent polarizability, by correspondence with peaks in the Raman spectrum of more ordered systems obtained using normal modes.

■

The bond polarizability model is intuitively appealing because polarization in an external field ultimately results from the perturbation of electron wave functions. The electrons most easily perturbed are those most loosely held, which are the valence electrons that form chemical bonds. Also, bonds are naturally anisotropic; valence electrons may polarize differently along a bond than transverse to the bond. This anisotropy naturally gives rise to configuration-dependent polarizability. Interactions between induced dipoles are neglected in the bond polarizability model. This is certainly convenient but also can be defended in applications to condensed phases, with the same arguments that appear in the derivation of the ClausiusMossotti relation. In short, the local fields from other induced dipoles acting at the site of a given polarizable entity tend to cancel, either for ordered phases with inversion symmetry or for disordered liquids. We represent the polarizability tensor P as a sum over bonds, as P=

∑ bonds i

[biL L(ni) + biT T(ni)]

(2)

Here ni is the unit tangent vector for bond i, Lαβ(n) = nαnβ and Tαβ(n) = δαβ − nαnβ are the longitudinal and transverse projection operators, respectively, and bLi and bTi are the bond longitudinal and transverse polarizability coefficients for bond i, respectively. For saturated hydrocarbons such as n-alkanes and polyethylene, we have only CC and CH bonds. Then the coefficients {bLi } are either bL,CC or bL,CH depending on the type of bond i, and likewise {bTi } are either bT,CC or bT,CH. Thus, we have four polarizability coefficients in the model, that must be fitted somehow to experiments or calculations of molecular polarizability. In the simplest version of bond polarizability models, the polarizability coefficients bLi and bTi are taken to be constants. In more detailed models, they can depend on bond length, in recognition that stretching a bond can change its polarizability. In our work, we neglect length dependence of the CH bond polarizabilities but retain a linear dependence on bond length for the CC bond polarizabilities. We neglect length dependence of CH polarizabilities because the CH bond stretching frequency (around 3100 cm−1) lies well outside the region where the Raman spectra are affected by the state of the system (crystal, rotator, or liquid). In contrast, we write the bond polarizabilities for the CC bonds as b(r) = b + b′Δr, in which Δr = r − r0 is the shift of the bond length r from its minimum energy value r0, and the coefficient b′ makes the polarizability depend linearly on bond length. The CC stretching modes have frequencies around 1100 cm−1 and lie within the region where structural features of condensed phases affect the Raman spectra. It turns out that retaining the length dependence of CC bond polarizabilities is necessary to predict reasonable Raman amplitudes of certain prominent features in the experimental spectra. We fit the four parameters in our bond polarizability model for saturated hydrocarbons to polarizabilities computed for a set of n-alkane conformers, for which density functional theory (DFT) calculations are feasible. We generated n-alkane configurations in two different ways. First, we created a small set of six conformers each for C10H22 and C12H26 n-alkanes, using the conformer library tool in Spartan (Wave function Inc.). Later, we created a larger set of 21 configurations of C10H22, by simulating a single molecule in vacuum using Gromacs37 (version 4.6.1). We used Gaussian0938 to compute

BOND POLARIZABILITY MODEL

To simulate the Raman spectra for condensed polymer phases, we need an efficient method to compute the system polarizability from the molecular configurations. Typically, molecular polarizability can be computed with a high degree of accuracy using ab initio methods. However, these calculations are very expensive and thus limited to small systems with few molecules and short simulation trajectories. For polarizability calculations to be practical for the scale of our simulations, we use an empirical approach called the bond polarizability model, which assumes that individual bond polarizabilities are additive. Each contribution depends upon the bond type, orientation, and length. Bond polarizability models have been widely employed to represent the polarizability of single molecules, polymers, silica glasses, and nanotubes.16,28,30,36 For polymeric materials, only a few studies have focused on computing the Raman spectra of crystals and uses rather small molecular systems.22 C


Article

Macromolecules

for the different conformers. The best fit values corresponding to (•) data points in the figure are (in units of Bohr3) bT,CC = 7.98, bL,CC = 14.6, bT,CH = 2.49, and bL,CH = 0.040. A closer investigation into the goodness of fit reveals that trade-offs are allowed between the magnitudes of CC and CH bond contributions, without compromising the model accuracy. In other words, a nearly as good a fit can be obtained by shifting either the isotropic or anisotropic CC bond coefficient by some amount, while compensating by increasing the corresponding CH bond coefficient proportionally. A more detailed investigation of the fitting parameters can be found in the Appendix. It turns out that the “best fit” bond coefficient values do not reproduce well the expected features in the Raman spectra. The best fit value of the bond anisotropy bL,CC − bT,CC is too small to describe well the Raman intensity for CC asymmetric stretch and CH bend vibrations, both of which depend largely on polarizability fluctuations caused by CC bond anisotropy. In a similar fashion, the Raman intensity for CH twist and CH bend vibrations mainly involve motion of CH bonds and so depend primarily on the CH bond anisotropy. A too-small value of the bond anisotropy leads to very small amplitude polarizability oscillations, which in turn gives too low a Raman intensity. To remedy this shortcoming, we modified the bond polarizability coefficients slightly, to increase both CC and CH bond anisotropy. The modified bond coefficients are bT,CC = 1.87, bL,CC = 16.55, bT,CH = 2.42, and bL,CH = 4.49, all in units of Bohr3. Figure 2 (panels a and b) compares the predicted polarizability components using the modified polarizability coefficients (in × markers) to DFT results. We find that compared to the original bond coefficients (in • markers), polarizability tensor predictions using the modified coefficients resulted in only a slight increase (about 3%) in the sum of square errors. The adjusted values fall slightly outside the formal confidence ellipsoid for the least-squares fit, but that ellipsoid is probably too optimiztic, as the errors in the fit are probably not Gaussian noise. We refer readers to the Appendix for details of the fitting process, confidence ellipsoids, adjustment of coefficients, and sensitivity of the Raman intensity to the coefficient values. To give a good account of Raman amplitudes for the CC stretching modes, it turns out the CC polarizability coefficients must depend on bond length. (We present a comparison between the final predicted Raman spectra with and without stretching contributions in the Appendix.) To fit the bond length dependence of the CC polarizability coefficients, we generate a new training set of conformers in which CC bond lengths vary. Specifically, we increased the length of each CC bond by 0.01 Å in turn, for the Spartan-generated set of six C10H22 and six C12H26 conformers, resulting in a total of 120 new conformers. We computed the shift in polarizability tensor ΔP̃ for this new training set using DFT, subtracting the polarizability of the modified structure from that of the original conformer. We then used least-squares fitting to minimize the difference between the predicted and calculated polarizability shift:

the polarizability tensor of each molecular configuration, using the B3LYP functional and the 6-311+G(3df,2p) basis set. For fitting purposes, we regard the result from Gaussian09 as the “true” polarizability tensor P̃ for each molecule. We determine the four polarizability coefficients in our model by minimizing the square errors χ2 in each element of the polarizability tensor, summed over molecular configurations c. That is, we minimize χ2 =

∑ ∑ [Pαβ(c) − Pαβ̃ (c)]2 configs c αβ

(3) T,CC

Because P is a linear function of the four coefficients {b , bL,CC, bT,CH, bT,CH}, the minimization is a linear least-squares fit, which can be performed analytically and leads to a unique minimum. Figure 2 displays scatter plots of the predicted polarizability tensor versus the “true” values, element by element. Data points

Figure 2. Scatter plots of predicted polarizability tensor elements versus DFT calculations, for training set of conformers of C10H22 and C12H26 n-alkanes: (a) diagonal elements and (b) off-diagonal elements. (•) best-fit bond polarizability constants; (×) modified bond polarizability constants. All polarizability components are in units of Bohr3.

shown by (•) markers predict the polarizability tensor using the linear best fit values as the bond polarizability coefficients. Because the molecular polarizability is only mildly anisotropic, the polarizability tensor tends to be proportional to the identity tensor. Thus, the diagonal elements tend to be in the same range of values, while the off-diagonal elements are smaller and scattered about zero. Therefore, we present scatter plots of the diagonal and off-diagonal elements separately, in Figure 2 (panels a and b). From the figure, it is apparent that the bond polarizability model gives a reasonable account of both the diagonal and off-diagonal elements of the polarizability tensor

χs2 =

∑ ∑ [ΔPαβ(c) − ΔPαβ̃ (c)]2 configs c αβ

(4)

Here ΔP̃ is the change in the “real” polarizability tensor calculated using DFT for each conformer after increasing the CC bond length. The shift ΔP in the predicted polarizability is given by D


Article

Macromolecules ΔP =

∑ CC bonds i

[b′iL L(ni) + b′Ti T(ni)]Δri

and GGA-PBE functionals, with inferior results). We start with geometry optimization using the LBGFS algorithm and apply an “ultra fine” energy cutoff (750.0 eV) and an “ultra fine” SCF tolerance (5.0e−7 eV/atom). The K-points mesh used for the single chain and crystal systems are 1 × 1 × 6 and 2 × 3 × 6, respectively. After optimization, the resulting lattice constants are a = 8.88 Å, b = 5.96 Å, and c = 2.57 Å. To determine the Raman-active vibrational modes, a phonon calculation at the Γpoint is performed (option CALCULATE_Raman set to true in the .param file). The resulting output file contains Raman frequencies, intensities, and vibrational trajectories. Figure 4 (panels a and b) compares the Raman modes of a single periodic chain and a periodic orthorhombic unit cell to

(5)

in which the change in the length of bond i is Δri (in units of Å), and b′Li and b′Ti are the coefficients of length dependence of the longitudinal and transverse polarizability coefficients bLi and bTi . We obtain best-fit values b′Li = 20.31 and b′Ti = 2.625, in units of Bohr3/Å. Figure 3 presents a scatter plot of the change in all polarizability tensor elements (both diagonal and off-diagonal),

Figure 3. Scatter plot of predicted versus calculated change in all polarizability tensor elements (diagonal and off-diagonal), from increasing length of each CC bond in turn for all training set conformers. Polarizabiltiy components are in units of Bohr3.

predicted from the bond polarizability model versus computed from DFT. Evidently, the bond polarizability model gives a reasonable account of polarizability changes arising from fluctuations in the CC bond length.

■

RAMAN MODES FROM DFT We first investigate the Raman-active modes of a single polyethylene chain and a periodic polyethylene crystal using density functional theory. Ab initio calculations offer the greatest accuracy in representing real molecules, although with high computational cost. Hence, the systems we used for these DFT calculations are small, consisting of only a few atoms. Our goal is to first relate peaks in the experimental spectra to a few key normal modes and then examine each mode by visualizing its motions in a single molecule or unit cell. For our calculations, we use the CASTEP module in Material Studio,39 which uses the plane-wave pseudopotential method suitable for periodic systems. For a single chain, we first import into CASTEP the structure of a single all-trans alkane in Protein DataBank format. We take advantage of symmetry by connecting the molecule to itself across the periodic boundary, and representing the molecule with the minimum unit cell of a single ethylene monomer (C2H4), to maximize computational efficiency. We likewise construct a periodic “herringbone” alkane crystal, by importing a structure file of an orthorhombic unit cell containing two molecules. Each molecule is made up of a single ethylene monomer (C2H4) in the trans configuration. The molecules are oriented 90° to each other in the plane normal to the chain axis. The initial crystal lattice constants are a = 7.51 Å, b = 4.89 Å, and c = 2.47 Å. All CASTEP calculations are performed using the GGAPW91 exchange correlation functional, which gives good agreement with the experiment (we also tried LDA-CAPZ

Figure 4. Raman modes (Γ points) determined using the plane-wave pseudopotential method for (a) a single periodic PE chain and (b) a periodic PE orthorhombic crystal; the red bracket highlights peak splitting. The labeled peaks corresponds to 1, CC asymmetric stretch; 2, CC symmetric stretch; 3, CH rock; 4, CH twist; and 5, CH bend. Dashed lines are experimental peak positions for the orthorhombic crystal. Intensity values are in arbitrary units.

experimental Raman frequencies of an PE orthorhombic crystal. We find that Raman mode frequencies remain mostly unchanged between the single chain versus the crystal, with only minor peak shifts. The main qualitative difference is the splitting in the crystal spectrum of the CH bend peak above 1400 cm−1 (discussed below). This suggests that modes in single chain and crystal systems are closely related. We exploit this when we optimize existing classical force fields to give better agreement with single-chain vibrational frequencies computed using DFT (discussed in the next section), arguing that improved agreement for single-chain E


Article

Macromolecules

CH bending mode is lower when the chains move in phase than when they move out of phase. The Raman mode frequencies for the periodic crystal computed using DFT match closely to experimental results. In particular, the appearance of peak splitting for the CH bending mode is consistent with the experimental signature used to identify the orthorhombic phase. However, the peak separation is smaller in the DFT results (15 cm−1) than in the experiment (25 cm−1).

modes will likely also lead to better agreement for modes in crystals and other condensed phases. From the single chain spectra in Figure 4a, we pick the top five modes with the greatest Raman activity and use CASTEP to generate vibrational trajectories for each mode. We label these modes in order of increasing wave numbers. We find that the observed molecular motions match expected vibrations and frequencies of real PE molecules. The vibrational motions of each labeled peak in the figure are identified as 1, CC asymmetric stretch; 2, CC symmetric stretch; 3, CH rock; 4, CH twist; and 5, CH bend. Figure 5 shows a schematic of each type of motion. The red and blue arrows indicate motions of carbon and hydrogen atoms, respectively.

■

RAMAN MODES FROM CLASSICAL FORCE FIELDS Before we generate the Raman spectra from molecular dynamics simulations, we must verify that our classical MD force fields correctly describe vibrational motions, with frequencies comparable to vibrations in real PE molecules. For this purpose, we examine Raman-active modes of single chain and crystalline PE through normal-mode analysis using a classical force field. We compute the entire set of eigenmodes of the system and then the Raman intensity for each mode using the bond polarizability model. We compare the modes with the highest Raman activity to our DFT results and to values for real polyethylene. Single-Chain Normal Modes. We represent a single decane (C10) molecule in vacuum, bonded to itself across the periodic boundary, using the OPLS all-atom force field. We place the molecule in a box for which the transverse dimensions (x and y) extend far beyond the chain, so the molecule does not interact with its periodic image. Covalently bonding the first and last carbon atom through the periodic boundary effectively simulates an “infinite” PE chain. Calculations are performed using Gromacs37 version 4.6.1. Normal mode analysis begins with a stringent energy minimization, using the LBFGS algorithm in double precision with energy tolerance 10 kJ mol−1 nm−1. The lowest energy conformation of the molecule is the trans configuration. We compute the mass-weighted Hessian matrix using md_run. Eigenvalues and eigenvectors are obtained using g_nmeig and g_anaeig. We use the eigenvectors to generate a vibrational trajectory for each mode and compute the Raman activity with the bond polarizability model. We select the five eigenmodes with the highest Raman activity and compare them to modes computed with DFT for the single chain, as well as the experimental Raman frequencies for a PE crystal. Figure 6a displays the results, in which classical frequencies appear as blue solid lines, DFT values as red dashed lines, and experimental results as gray dashed lines with corresponding labels above the figure. Here we have omitted from the DFT results the CH wagging mode (∼1340 cm−1). This mode is not used to characterize PE; it has very weak intensity in the experiment and also shows very weak Raman activity in our normal-mode analysis using the classical force field. We visualize the motion of each classical mode in Figure 6a and match them to the motions observed in DFT. The OPLS force field produces the correct vibrational motions; however, the frequencies are not very accurate compared to DFT, so much so that the modes come in the wrong order. As remarked previously, we expect the frequencies for single chain modes to be quite close to the experimental frequencies for a crystal. Evidently, the unmodified OPLS force field is not very accurate in representing the modes of interest. Classical force fields are developed as approximate representations of a quantummechanical system, in which some set of structural and

Figure 5. Vibrational motions of PE chain in each Raman active mode. 1, CC asymmetric stretch; 2, CC symmetric stretch; 3, CH rock; 4, CH twist; and 5, CH bend. Red arrows indicate carbon atom motions, blue arrows indicate hydrogen atom motions.

Except for the CH rocking mode,40 these five modes are used experimentally as spectral signatures to characterize PE phases.5,6,8,9 The CH rocking mode has a very low experimental intensity and is difficult to observe. In addition to the labeled peaks in Figure 4a, there is a sixth peak in the figure which corresponds to a CH wagging motion. Although it is Ramanactive in principle (not ruled out by symmetry), it has even lower experimental intensity and so is likewise not used to characterize PE. One difference between the single-chain and the crystal spectra is that the crystal spectrum has two peaks for each type of Raman vibration because the crystal system contains two chains instead of one. For the most part, these doublets are degenerate in frequency, or nearly so. Visualizing the molecular motions in the crystal, we realize that for each doublet, the two chains in the unit cell can oscillate either in phase or out of phase with respect to each other. If the chains interact as they move, this can serve to split the doublet, as for the CH bending mode. Because of the orthorhombic packing of the chains in the crystal, the restoring force and thus the frequency of the F


Article

Macromolecules

Table 1. Compare Original and Modified Spring Constants kb of the Classical OPLS All-Atom Force Field bonded force (kb)

OPLS-AA

modified OPLS-AA

CC bond (kJ mol−1 nm−2) HCH angle (kJ mol−1 rad−2) CCH angle (kJ mol−1 rad−2) CCC angle (kJ mol−1 rad−2)

224262.4 276.144 313.800 488.273

291541.0 275.438 472.830 500.565

Next, we use this modified force field to investigate the spectra of bulk systems. Normal Modes for Crystals. The characteristic Raman signature of the PE orthorhombic crystal is the peak splitting occurring around 1416 cm−1 for the CH bending mode. The splitting is unique to the orthorhombic packing structure; it disappears when chains become more disordered, transitioning to the rotator phase or to the melt. Thus, when simulating the spectra of a PE crystal it is essential that the peak splitting is captured. Previously, we computed the Raman vibrations of an orthorhombic crystal using DFT and observed that the crystal unit cell generates two degenerate modes for each type of vibrational motion. To verify that this degeneracy is also observed with our modified OPLS classical force field, we performed normal-mode analysis on a periodic PE crystal. We constructed a system composed of 8 × 6 periodically connected C40 molecules packed in the orthorhombic crystal structure. Figure 7b shows a schematic of the system. Periodic boundary conditions are applied to simulate an infinite PE crystal. To obtain normal modes, a stringent energy minimization is performed, followed by diagonalization of the mass-weighted Hessian matrix. We apply the bond polarizability model to the resulting eigenmodes to compute the Raman intensities. Performing normal-mode analysis on a large system results in a large number of eigenmodes, among which the Γ point is often difficult to identify. Normal mode analysis identifies 17280 eigenmodes for our 8 × 6 C40 crystal “supercell”. These eigenmodes form vibrational bands. Figure 8a displays these vibrational bands within the relevant frequency range, in light gray. At the Γ point, all periodic images of the minimal unit cell (containing two C2H4 monomers) oscillate in phase. In other modes in the same vibrational band, the local motion is the same as at the Γ point but with different temporal phase at different locations (i.e., modulated by plane waves). For direct comparison to our previous DFT results for modes at the Γ point only, we need a way to extract the Γ point from each vibrational band. Diagonalizing the Hessian for a single unit cell composed of two C2H4 monomers would give the Γ point eigenmodes of the orthorhombic crystal. Indeed, this is the same unit cell we used in our DFT calculations. We want to use Gromacs to compute the eigenmodes because of its convenient routines for minimizing the energy and computing the Hessian of a very complicated Hamiltonian. Unfortunately, we cannot directly represent such a small system in Gromacs because of restrictions on the minimum system size. The linear dimensions of the simulation volume must all be greater than twice the cutoff distance for nonbonded interactions. To obtain the Γ point frequencies from the Hessian for a large system, we develop a “Hessian projection” method that can compute the Hessian matrix for the unit cell from the much larger supercell Hessian matrix. In essence, the method partitions the supercell into identical copies of the unit cell

Figure 6. Raman-active modes of a single C10 periodic chain in vacuum from normal-mode analysis using (a) unmodified OPLS allatom force field and (b) modified OPLS force field. Classical frequencies are blue solid lines, DFT values are red dashed lines; experimental results for PE crystal are dashed gray lines.

dynamical properties are chosen for force matching. Ordinarily, matching the properties of infrared or Raman vibrational modes would not be a priority. Force Field Optimization. For our purposes, we must optimize the OPLS force field so MD simulations can reproduce accurate Raman frequencies. On the basis of the motions of the relevant modes, we assert that the most important force field parameters to adjust are the intrachain spring constants for the CC bond, HCH angle, CCH angle, and CCC angle. We tune these four parameters to minimize the sum of square errors between OPLS and DFT calculated frequencies for a single chain. We have verified that these four parameters have the most pronounced effect in altering peak positions. When adjusted individually, each parameter shifts the frequencies of multiple peaks. Hence, we minimize the error tuning all four parameters simultaneously, using the conjugate gradient method. The minimization was done “by hand”, in the sense that each step required a new set of simulations, which were set up and launched manually. We stopped iterations when the classical spectrum no longer made visible improvements. Table 1 shows the original and the modified force field parameter values, which are generally rather close. The largest change is an increase in about 50% for the C−C−H angle parameter. These changes do not appear to affect other properties of the system relevant to our work. Figure 6b shows Raman frequencies obtained from the modified OPLS force field. The resulting peak positions are now in reasonable agreement with both DFT and experiment. G


Article

Macromolecules

Figure 8. Raman activity of eigenmodes from (a) C40 periodic crystal and (b) C40 periodic chain in vacuum. Normal modes obtained from original Hessian shown by gray, solid lines; Γ points are shown by black, solid lines. Experimental frequencies for PE crystal indicated by dashed, gray lines. Red bracket indicates peak splitting.

Figure 8a show that eigenmodes from the original supercell spectrum displayed the greatest Raman activity near the Γ frequencies. Using the same method, we also calculated the Γ points of a C40 single chain system in vacuum, as shown in Figure 8b. Because this system is small, there is no difficulty in identifying the Γ points directly from the supercell spectrum. However, to further demonstrate the validity of the Hessian projection method, we projected the supercell Hessian onto a single C2H4 monomer unit cell. The resulting eigenmode frequencies shows an exact match to the Γ frequencies of the supercell. With the spectra greatly simplified, we can now compare Raman frequencies obtained using the modified OPLS force field, to DFT and experimental results. As expected, there are two degenerate modes associated with each type of Ramanactive vibration. Using the eigenvectors to construct vibrational trajectories for each mode, we confirm that the degeneracy corresponds to in-phase and out-of-phase motion of the two chains in the unit cell. For a single chain system, the degeneracy is no longer present. The characteristic peak splitting for the CH bending mode is also observed. However, the splitting (10 cm−1) is smaller than DFT results (15 cm−1) and smaller still than the experiment (25 cm−1). In searching for an explanation of the discrepancy in the CH bend peak splitting, we found that the splitting is larger when the system is more dense. The peak splitting results from the

Figure 7. Size of orthorhombic unit cell (red) compared to system simulated with molecular dynamics. (a) Single periodic chain simulated using MD; a unit cell contains 2 carbon atoms per chain. (b) 8 × 6 C40 orthorhombic crystal; a single unit cell contains 2 chains.

and likewise partitions the matrix elements of the supercell Hessian into smaller submatrices that correspond to copies of the unit cell. At the Γ point, all these unit cell copies move in the same way. Hence we obtain the unit cell Hessian by summing all the submatrices and normalizing by the number of copies in the supercell. The eigenmodes of the projected Hessian are the modes of the unit cell at the Γ point. We present details of the Hessian projection method in the Appendix. We obtain the Γ frequencies of our simulated 8 × 6 C40 crystal by projecting the Hessian matrix down to a single C2 orthorhombic unit cell. Figure 7 (panels a and b) shows the unit cell structure in red, viewed from different directions. We expect Γ points to have the greatest Raman intensity compared to other modes in the vibrational band because all monomers in the unit cell are vibrating in phase. Conversely, out of phase motions decrease Raman intensity because it creates asymmetric oscillations in the molecule, thus reducing the overall net change in the molecular polarizability. Indeed, our results in H


Article

Macromolecules

fs. During this run, we switch to the Parrinello−Rahman barostat (τp = 100) to maintain the correct thermodynamic ensemble. We divide the total signal into 10 segments of 25 ps, compute the power spectrum of each segment, and average the ten spectra to obtain the Raman spectra of Figure 10.41 A 25 ps

additional restoring force on the CH bending mode when adjacent chains in the unit cell move out of phase. Tighter packing increases the repulsive interactions between adjacent chains and hence increases the vibrational frequency for out of phase motion compared to motion in phase. At finite temperature, when chains are fluctuating thermally about their minimum energy configurations, the effective repulsion between adjacent chains will increase, and with it the CH bend peak splitting.

■

RAMAN SPECTRA FROM MD TRAJECTORIES We have shown that with modest modifications to the OPLS classical force field, we obtain Raman-active normal-mode frequencies reasonably close to experimental values. With the normal mode approach, we can analyze and visualize independent vibrational modes. However, this approach is limited to ordered phases and cannot predict line widths. Now, we compute Raman spectra of both crystalline and disordered phases from the time-dependent polarizability. The bond polarizability model makes this a practical calculation for large systems, disordered or not, since the contribution of each bond is additive, and readily computed from the bond orientation and length. Figure 9 shows an example of the time-dependent polarizability tensor components, calculated from a short (25 ps)

Figure 10. Time-averaged Raman spectra of orthorhombic crystal simulated at various temperatures under NσT conditions. Dashed lines show experimental Raman peak positions for PE orthorhombic phase.

time series with 2 fs timesteps has a maximum frequency of 16600 cm−1 and a frequency increment of 1.33 cm−1, which can represent the spectra of interest with ample line shape resolution. We average multiple spectra to reduce the noise; we do not apply any Gaussian smoothing of the resulting spectra. Figure 10 shows distinct peaks near the experimental Raman frequencies of the polyethylene crystal. The peaks in the figure have line widths in reasonable semiquantitative agreement with experimental spectra. For example, the full width at halfmaximum (fwhm) at 200 and 300 K for the CH twisting peak in Figure 10 are about 8 cm−1. The experimental fwhm for the same peak in Figure 1 at 303 and 310 K for the crystal and rotator phases are about 6 cm−1. The line width reflects damping and finite mode lifetimes as the result of thermal fluctuations. Note that a line width of about 10 cm−1 corresponds to a mode lifetime of around 1 ps, which is a typical collision time in a condensed phase of organic molecules at standard temperature and pressure. Collisions between molecules under these conditions are rather efficient in disrupting local vibrational modes and redistributing energy between them. The experimental intensity for the CH rocking peak at 1171 cm−1 is extremely weak and is not used as a spectral feature to characterize PE phases. In our results, the CH rocking peak is likely obscured by the intense CC symmetric stretch peak, since previous normal-mode analysis (Figure 8) showed that the CC symmetric stretch and CC rocking frequencies are quite close in our simulations, more so than in the experiment. The predicted spectra reproduce the characteristic peak splitting of the CH bending mode (around 1450 cm−1). Remarkably, the spectrum at 200 K has a splitting of over 20 cm−1, closer to the experimental value of 25 cm−1 than results from DFT (15 cm−1) and normal-mode analysis (10 cm−1). This appears to be an effect of finite temperature; fluctuating chains interact more strongly with their neighbors, which increases the difference between in-phase and out-of-phase mode frequencies, discussed further below.

Figure 9. Time-dependent polarizability components calculated from a 25 ps trajectory of PE crystal system. (black) Pxx, (red) Pyy, (blue) Pzz.

simulation trajectory for a PE crystal system. (Only diagonal components are shown, but off-diagonal components behave similarly.) The signals clearly show strong oscillations, which suggests strong contributions from relatively few Raman-active vibrational modes. To accurately represent such time-dependent data, it is crucial to sample the trajectory frequently. In our work, we save system “snapshots” every 2 fs (i.e., every two MD timesteps). The highest frequency mode of interest here is the CH bending mode with a period of about 23 fs, which is wellrepresented by a time series with 2 fs between snapshots. We first investigate the time-averaged Raman spectrum of the orthorhombic crystal. We use the same periodic array of C40 molecules and modified OPLS force field as in our normalmode analysis. For the simulations, we start with energy minimization using steepest descents to a tolerance of 10 kJ/(mol·nm). Then we equilibrate for 500 ps under constant stress and temperature NσT conditions. A velocity-rescaling thermostat (τT = 0.1) and Berendsen barostat (τp = 1) maintain system temperature and pressure at 200 K and 1 bar. The time step in all MD runs is 1 fs. To obtain the input trajectory for computing the time-dependent polarizability, we simulate after equilibration for another 250 ps, saving configurations every 2 I


Article

Macromolecules

Figure 11. Computed Raman spectra of various PE phases simulated under NσT conditions. All spectra averaged over 100 independent trajectories of 25 ps. Dashed lines are experimental peak positions for PE orthorhombic phase.

or −90°. The lattice constants are also increased sightly to a = 0.768 nm and b = 0.498 nm, consistent with previous simulations and experiment.3 For the R2 rotator phase which is more rotationally disordered than R1, we expanded the orthorhombic system onto a regular hexagonal lattice then rotated each chain randomly about its axis by an arbitrary angle. Each unit cell still contains two chains, with lattice constants a = 0.821 nm and b = 0.477 nm. From previous simulation work, we know that the R1 rotator phase is not stable in alkanes using the OPLS force field.3 However, drastic structural rearrangements are strongly inhibited by using periodically self-bonded molecules. The CC bond across the periodic boundary in the z direction restricts movement of the chain ends. Although this rotator phase is artificially constructed, the resulting Raman spectrum does resemble the spectrum of an alkane rotator phase. Polyethylene can form a metastable monoclinic phase, under conditions of elongation at atmospheric pressure.43 The monoclinic unit cell consists of a single chain, with lattice constants a = 0.451 nm, b = 0.482 nm, and γ = 122.3°.3 We construct a monoclinic array of 6 × 8 periodically connected C40 chains, for consistency with our crystal and rotator phase systems. To construct a melt phase, nonperiodic molecules must be used, since chain ends must move freely for the system to melt. Starting with a 8 × 6 crystalline array of C40 nonperiodic chains, we equilibrate at 500 K and 1 bar under NσT conditions until the crystal melts. Figure 11 compares the simulated Raman spectra of all PE phases. All spectra are obtained from the averaging of 100 independent spectra computed from 25 ps trajectories performed under NσT conditions at atmospheric pressure, at temperatures appropriate to each phase (200 K for orthorhombic, monoclinic, and R1; 300 K for R2; and 500 K for melt). Our simulated Raman spectra reproduces the qualitative features of experimental spectra for crystal, rotator, and melt phases. The peak splitting for the CH bending mode is indeed a signature of the orthorhombic phase only. In the monoclinic phase all chains are structurally identical, so there is no analog

We investigate the effect of temperature on predicted Raman spectra by simulating the same crystal at 100, 200, and 300 K. Figure 10 shows that with decreasing temperature, line widths become narrower and even modest splittings between nearly degenerate modes become visible. The spectrum at 100 K shows well-defined doublets for every type of Raman vibration, while peak splittings are less distinct at 300 K. An additional reason for this may be that at lower temperatures the system is more dense, and crystal splittings become stronger with increasing density, as we saw in our normal-mode analysis. Figure 10 shows that at very low temperature (100 K), multiple low intensity peaks begin to emerge from the background. These peaks do not seem to be “noise” (i.e., cannot be removed by averaging over more independent spectra). They appear to correspond to vibrational modes with weak Raman intensity and become visible at low temperatures because modes are so weakly damped that lineshapes become narrow. At higher temperatures, thermal noise broadens these weak lines into obscurity. We remark that the predicted Raman activity is plotted in log scale. This allows better visualization and side-by-side comparison between peaks with disparate intensities. However, a drawback is that the peak line widths appear artificially broad in log scale compared to the experimental spectra which is plotted in linear scale. To verify, we measured the full widths at half-maximum of our predicted Raman peaks and found them to be roughly consistent with experimental observations. Spectra of Disordered Phases. An advantage of using time-dependent polarizability to compute Raman spectra is that the method applies to disordered systems as well as crystals. Here, we compute the Raman spectra of various bulk PE phase structures, including the metastable monoclinic phase, R1 and R2 rotator phases, and the melt phase. Figure 11 presents the spectra, along with a schematic of each phase structure. Rotator phases are stable in polyethylene only at elevated pressures. At atmospheric pressure, rotator phases are stable only in alkanes within a narrow range of molecular weights.42 To build the R1 rotator phase, we start with an orthorhombic crystal composed of 8 × 6 periodically connected C40 chains and then randomly rotate each chain about its axis by 0°, + 90°, J


Article

Macromolecules

calculations and experimental Raman peak positions. To use MD simulation for computing spectroscopic properties, it becomes important to check vibrational frequencies against DFT calculations and experiment as part of MD force field parametrization. Our method is useful for its ability to directly link an assumed microstructure with its Raman signature. Using nalkanes and polyethylene as a test case, we have successfully captured the main features of experimental Raman spectra for crystalline and disordered phases, including the position and lineshapes of peaks in the frequency range of 1000−1500 cm−1. This range includes the CC stretch, CH twist, and CH bend modes, used experimentally to identify alkane condensed phases. In our evaluation of computed Raman spectra, we focus on predicting frequencies, representing lineshapes, and determining modes as we assess the potential for Raman spectroscopy to identify phases. We avoid relying on quantitative spectral intensities. Predicted peak intensities are sensitive to details of the bond polarizability model, as well as to its simplifying assumptions. Fundamentally, molecular polarizability is a manybody problem; polarization of a bond can be affected by its immediate surroundings, which is not accounted for in the model. However, the assumption that bonds contribute independently to the polarizability is what makes the model simple enough to use in repeatedly computing the polarizability for systems of thousands of atoms. Another reason for caution against analysis of the predicted spectral intensities is that the spectrum is generated from simulation trajectories using classical force fields. As a result, higher frequency vibrations which are quantum-mechanically “frozen out” at lower temperatures contribute too much in simulation trajectories and ultimately result in too much Raman intensity in the predicted spectrum. We deal with this problem in part by avoiding the CH bond stretching region of the Raman spectrum, which anyway is not useful in identifying phases. However, other vibrational motions which may only be partially classical in reality are treated as fully classical in our simulation. Fortunately, over the relatively narrow frequency range of interest, the ratio of the true quantum-mechanical vibrational amplitudes to the classical amplitudes is nearly constant on a log scale, so that the shape of the computed spectra is little affected. In future work, we can use this method to compute the Raman spectra of hypothetical structures, as an aid to their identification. For example, we can compute the Raman signatures of our best guess for the structure of proposed rotator phases in isotactic polypropylene, which presently are not well-characterized.

of in-phase versus out-of-phase modes for two chains in the unit cell. This is likewise true of the R1 and R2 phases, which randomize the herringbone packing of the orthorhombic phase and thereby render all chains at least statistically identical in structure. As the system becomes increasingly more disordered, transforming to the rotator phases and ultimately the melt, the peaks become broader, again consistent with experimental spectra. This is likely the result of two effects: first, at higher temperatures, stronger fluctuations give increased damping and shorter mode lifetimes. Second, peaks become much broader in disordered phases presumably because heterogeneous local environments seen by identical moieties on different chains lead to random shifts in local vibration frequencies. The monoclinic phase is distinguishable from the R1 and R2 phases, mainly by its sharper peaks. In contrast, our results suggest that there are no qualitative features in the Raman spectrum that can be used to differentiate between the R1 and R2 phases. Only slight differences in peak width can be observed, with somewhat narrower peaks in the R1 phase than in the more disordered R2 phase. Yet the two phases are definitely structurally different, which illustrates a limitation of Raman spectroscopy alone in distinguishing phases. The single CH bending peak in both rotator phases (and monoclinic) exhibit a slight upward shift in wavenumber from the crystal CH bending peak. Although no corresponding shift is observed in the experimental rotator phase spectra, the experimental peak intensity near this frequency is complicated by contributing signals from Fermi resonance.5 No frequency shift is observed in peaks other than CH bending, suggesting that the majority of chains remain in the trans conformation. The melt phase spectra shows that the previously distinct peaks for CC asymmetric and CC symmetric stretch modes combine into one broad feature. This is consistent with experimental observations on melting, of a broad amorphous CC stretching peak that emerges at a wavenumber between the two crystalline CC stretching modes. We remark that because our melt phase is simulated at such a high temperature, thermal noise leads to a high background intensity compared to spectra of other phases, which causes the relative peak intensity to appear lower. However, this artifact is unavoidable in our simulations because lowering the temperature sometimes causes our chains to crystallize.

■

CONCLUSIONS We have developed a comprehensive approach to predict Raman spectra from molecular dynamics simulation, which we have applied to polyethylene and alkanes. We test our approach by comparing to Raman signatures used by experimentalists to distinguish PE bulk phases. We compute Raman spectra in two ways: (1) using normalmode analysis, which works well for ordered phases with welldefined vibrational modes, but does not predict lineshapes and (2) by directly computing the power spectrum of the timedependent polarizability of the entire simulation. Both methods rely on a bond polarizability model, which approximates the total polarizability as a sum of contributions from each bond. Parameters of the bond polarizability model are fitted to density functional theory (DFT) results for a set of alkane conformers. Our molecular simulations are performed with a modified OPLS force field. We have tuned several bond-angle and bondstretching parameters to improve agreement of normal-mode frequencies from classical simulations with results from DFT

■

APPENDIX A: SIMPLIFIED TREATMENT OF RAMAN SPECTROSCOPY Raman scattering is most often presented in a quantummechanical description but can be described for the most part classically, as follows. The induced dipole moment (μ) of a molecule is proportional to the electric field, μ = αE (6) where E is the incident electric field and α is the molecular polarizability. The molecular polarizability (α) is a tensor, which describes the deformability of the electron cloud about the molecule by an external field. We can expand the molecular polarizability K


Article

Macromolecules

Figure 12. Confidence ellipse (95% confidence level) for CC and CH bond polarizabilities: (a) isotropic and (b) anisotropic. • best fit bond polarizability constants, × modified bond polarizability contants, and □ bond polarizability constants from Smirnov and Bougeard.30

the mode is Raman active. Likewise, when identical moieties vibrating with a given normal mode lead to dipoles that oscillate in phase, the mode is infrared active: the oscillating dipole radiates, scattering light elastically. Because polarizability is a tensor and dipoles are vectors, the two objects transform with opposite signs under mirror symmetry operations. Hence if identical moieties on a vibrating molecule contribute polarizabilities oscillating in phase, the same moieties would have dipoles oscillating out of phase, and vice versa. Thus, modes of a given parity under a reflection symmetry are either Raman active and IR silent, or vice versa.

with respect to the amplitude of normal modes. For a normal mode k with frequency νk, the time-dependent molecular polarizability αk is expressed as αk = α0 + αk′ × Q k 0 × cos(2πνkt + ϕk )

where αk′ =

( ), Q ∂α ∂Q k

k

(7)

is the normal coordinate along the

0

vibrational mode, Qk0 is the mode amplitude, ϕk is the phase shift, and α0 is the unperturbed polarizability. In Raman spectroscopy, an incident electric field with amplitude E0 and frequency ν0 is applied to probe the vibrational motions of the molecules. E = E0 × cos(2πν0t )

■

APPENDIX B: BOND POLARIZABILITY MODEL FITTING The bond polarizability model predicts the polarizability tensor from molecular geometry. By fitting this model to a training set of known molecular polarizability calculated from DFT, we obtained four bond coefficients that describe longitudinal and transverse contributions of CC and CH bonds. We are interested not only in these best fit values but also in how well they are known. Using the machinery of the linear least-squares fit, we determine confidence intervals for each parameter. It turns out these confidence intervals are strongly correlated. The “confidence hyperellipsoid”, a 4D surface enclosing the region of acceptable parameter values, is very far from spherical, so that changes in one parameter can be compensated by changes in another. To gain insight into the confidence region, it is helpful to express the bond polarizabilities in terms of isotropic and traceless anisotropic contributions. (We emphasize that the results found by fitting isotropic and anisotropic coefficients are completely equivalent to those found by fitting transverse and longitudinal coefficients.) For each bond, we can write

(8)

The induced dipole follows from substituting eq 7 and 8 into eq 6: 1 αk′ × E0Q k 0 cos[2π (ν0 + νk)t + ϕk ] 2 1 + αk′ × E0Q k 0 cos[2π (ν0 − νk)t − ϕk ] 2

μ = α0 × E0 cos(2πν0t ) +

(9)

Eq 9 shows that the induced dipole is the sum of three components with frequencies ν0, ν0 + νk, ν0 − νk. The first term represents elastic scattering of light (Rayleigh scattering), and remaining two terms represent inelastic scattering, in which the frequency is either the sum (anti-Stokes) or difference (Stokes) of the incident field and vibrational mode frequencies. The Raman spectrum comprises Stokes and anti-Stokes scattering, in which the frequency of scattered light is either increased or decreased by “beating” together with the oscillating polarization. Molecular conformations can have symmetry, for example with respect to mirror planes and inversion centers. Symmetric molecules have normal modes which are either even or odd under the symmetry operation, which leads to selection rules for Raman activity, by applying eq 9: ⎛ ∂α ⎞ ⎟⎟ ≠ 0 αk′ = ⎜⎜ ⎝ ∂Q k ⎠0

bL L(n) + bT T(n) = bI 1 + bU Q(n)

(11)

Here 1 is the identity tensor and Q(n) is the traceless uniaxial “quadrupolar” tensor [i.e., Qαβ(n) = nαnβ − (1/3) δαβ]. Equating terms proportional to δαβ and proportional to nαnβ, we connect bI and bU to bL and bT, with the relations bT = bI + bU and bL = bI − 2bU. Physically, bI is the isotropic polarizability of the bond and bU the traceless anisotropic polarizability. Fitting bI and bU turns out to be separate problems because the isotropic and

(10)

When identical moieties on a molecule vibrating with a given normal mode contribute polarizabilities that oscillate in phase, L


Article

Macromolecules

Figure 13. Time-dependent Raman spectra for orthorhobmic crystal at 200 K computed using (a) modified bond coefficients (black) vs best fit bond coefficients (blue) and (b) modified bond coefficients with CC stretching contriubtions (black) vs no CC stretching (red).

are slightly outside. The same figures also present bond coefficients reported by Smirnov and Bougeard.30 These values were obtained by fitting a similar linear model to a training set of various small molecules containing CC and CH bonds. The total sum of square errors computed using the new bond coefficients is only about 3% larger than the best fit value, and the resulting scatter plot of predicted versus DFT polarizability values are scarcely distinguishable by eye. The confidence ellipses are probably too optimiztic, since they are based on the assumption that the errors result from Gaussian noise rather than systematic errors from shortcomings of the model. We compare the resulting Raman spectra predicted using both the best fit bond coefficient values and the modified values. Figure 13a shows the spectra for the orthorhombic crystal simulated at 200 K. We observe that use of the best fit coefficients falsely amplifies some lower intensity peaks around 1060 cm−1, giving the appearance of a lower intensity CC asymmetric stretch peak relative to the background at around 1040 cm−1. In addition, the intensity of the CH bending peak at around 1440 cm−1 also appears much lower when best fit bond coefficients are used, suggesting that certain vibrational motions are more sensitive to changes in bond anisotropy. Figure 13b illustrates the sensitivity of Raman intensities with respect to CC bond stretching. We find that the CC asymmetric stretch peak around 1060 cm−1 is barely visible when bond stretching contributions are not incorporated into the bond polarizability model. The intensities for all other Raman peaks remain unaffected, suggesting that these vibrational motions depend largely on fluctuations in bond orientation.

anisotropic polarizabilities transform differently under rotations. Thus, the confidence ellipse for bI,CC and bI,CH is independent of the ellipse for bU,CC and bU,CH. Figure 12a displays the confidence ellipse for the isotropic polarizabilities, centered on the best fit values bI,CC = 10.18 and bI,CH = 1.673. The ellipse is extremely long and narrow, with a slope of −1/2. This indicates that nearly as good a fit can be obtained, by shifting the best fit bI,CC upward by some amount Δ and compensating by reducing bI,CH by Δ/2. This is because our training conformers are long n-alkanes, made mostly of CH2 groups; the isotropic polarizability can come either from the CC bonds or from the CH bonds, which are about twice as numerous. Indeed, it was essential that we include alkanes of different length in our training set, or we could not parse the isotropic polarizability into CC and CH contributions. The difference between C10H22 and C12H26 from this perspective is that the latter has a slightly smaller fraction of CH bonds (26/37 = 0.7 vs 22/31 = 0.71). We could probably tighten this ellipse by including a greater range of chain lengths in the training set. However, this would eventually push against another implicit assumption, that all CH bonds contribute the same way, and all CC bonds likewise. The difference between C10H22 and C12H26 in the model depends on the different fractions of methyl CH bonds (6/31 = 0.19 vs 6/37 = 0.16) in the two molecules. The approach assumes these methyl CH bonds polarize the same as CH bonds in CH2 moieties, which must break down for sufficiently short alkanes. Figure 12b displays the confidence ellipse for the anisotropic polarizabilities, centered on the best fit values bU,CC = −2.197 and bU,CH = 0.8168. The ellipse is likewise long and narrow, with a slope of about 0.54. Again, nearly as good a fit can be obtained by shifting the best fit value bU,CC upward and compensating by increasing bI,CH. Because the CC and CH bonds meet at tetrahedrally coordinated carbons with nearly fixed bond angles, anisotropic polarization of the overall alkane regardless of conformer can arise either from anisotropic CC bonds or from anisotropic CH bonds. Again, we must include chains of different length in the training set to parse the contribution from CC and CH bonds. The × data point in Figure 12 (panels a and b) displays a new set of bond coefficients, adjusted to increase the bond anisotropy and better reproduce experimental spectral features. The new bond coefficients fall on the major axis of the confidence ellipses. The new isotropic bond parameters reside within the confidence region, while the anisotropic parameters

■

APPENDIX C: HESSIAN PROJECTION In simulations of crystals, we are often obliged to simulate a “supercell” consisting of many copies of the unit cell, to keep the linear dimension of the simulation volume longer than twice the cutoff and avoid self-interactions with periodic images. We can obtain the Γ point modes of such a supercell, by projecting its Hessian onto a smaller matrix corresponding to a single unit cell. Suppose the original supercell consists of N copies of the unit cell, which contains n atoms. The Hessian matrix then consists of an N × N array of blocks, each of which is 3n × 3n. The j,j diagonal blocks are the self-interactions of the jth copy of the unit cell; the off-diagonal i,j blocks are the interactions of the ith copy of the unit cell, with the nearest periodic image of the jth copy. M


Article

Macromolecules In a normal mode at the Γ point, atoms in every copy of the unit cell move in the same way. The system displacement vector V, of dimension N × 3n, then consists of N copies of the same 3n-dimensional vector v. The system energy in harmonic approximation is then a quadratic form consisting of the double sum over the blocks, contracted with the vector v. (To obtain the energy per unit cell, divide by N.) We call the double sum over blocks the “projected Hessian”, since it is the Hessian for a single unit cell interacting with all its multiple periodic images in the supercell. Figure 14 illustrates the Hessian projection method for a one-dimensional system, containing three copies of a unit cell

(6) Naylor, C. C.; Meier, R. J.; Kip, B. J.; Williams, K. P. J.; Mason, S. M.; Conroy, N.; Gerrard, D. L. Raman Spectroscopy Employed for the Determination of the Intermediate Phase in Polyethylene. Macromolecules 1995, 28, 2969−2978. (7) Chai, C. K.; Dixon, N. M.; Gerrard, D. L.; Reed, W. Rheo-Raman studies of polyethylene melts. Polymer 1995, 36, 661−663. (8) Migler, K. B.; Kotula, A. P.; Hight Walker, A. R. Trans-Rich Structures in Early Stage Crystallization of Polyethylene. Macromolecules 2015, 48, 4555−4561. (9) STROBL, G. R.; Hagedorn, W. Raman spectroscopic method for determining the crystallinity of polyethylene. J. Polym. Sci., Polym. Phys. Ed. 1978, 16, 1181−1193. (10) Androsch, R.; Di Lorenzo, M. L.; Schick, C.; Wunderlich, B. Mesophases in polyethylene, polypropylene, and poly(1-butene). Polymer 2010, 51, 4639−4662. (11) Yuan, S.; Li, Z.; Hong, Y.-l.; Ke, Y.; Kang, J.; Kamimura, A.; Otsubo, A.; Miyoshi, T. Folding of Polymer Chains in the Early Stage of Crystallization. ACS Macro Lett. 2015, 4, 1382−1385. (12) Wang, Z.-G.; Hsiao, B. S.; Srinivas, S.; Brown, G. M.; Tsou, A. H.; Cheng, S. Z. D.; Stein, R. S. Phase transformation in quenched mesomorphic isotactic polypropylene. Polymer 2001, 42, 7561−7566. (13) Zia, Q.; Mileva, D.; Androsch, R. Rigid Amorphous Fraction in Isotactic Polypropylene. Macromolecules 2008, 41, 8095−8102. (14) Olley, R. H.; Bassett, D. C. On the development of polypropylene spherulites. Polymer 1989, 30, 399. (15) Gatos, K. G.; Minogianni, C.; Galiotis, C. Quantifying crystalline fraction within polymer spherulites. Macromolecules 2007, 40, 786. (16) Snyder, R. G.; Kim, Y. Conformation and low-frequency isotropic Raman spectra of the liquid n-alkanes C4-C9. J. Phys. Chem. 1991, 95, 602−610. (17) Snyder, R. G. Chain conformation from the direct calculation of the Raman spectra of the liquid n -alkanes C 12 − C 20. J. Chem. Soc., Faraday Trans. 1992, 88, 1823−1833. (18) Meier, R. J. Studying the length of trans conformational sequences in polyethylene using Raman spectroscopy: a computational study. Polymer 2002, 43, 517−522. (19) Hallmark, V. M.; Bohan, S. P.; Strauss, H. L.; Snyder, R. G. Analysis of the low-frequency isotropic Raman spectrum of molten isotactic polypropylene. Macromolecules 1991, 24, 4025−4032. (20) Galabov, B. S.; Dudev, T. In Vibrational Intensities; Durig, J. R., Ed.; Vibrational Spectra and Structure; Elsevier Science B.V.: Amsterdam, The Netherlands, 1996; Vol. 22. (21) Henssge, E.; Dumont, D.; Fischer, D.; Bougeard, D. Analysis of infrared and Raman spectra calculated by molecular dynamics. J. Mol. Struct. 1999, 482−483, 491−496. (22) Dumont, D.; Henssge, E.; Fischer, D.; Bougeard, D. Simulation of vibrational spectra of polymers by molecular dynamics calculations. Macromol. Theory Simul. 1998, 7, 373−379. (23) Berens, P. H.; White, S. R.; Wilson, K. R. Molecular dynamics and spectra. II. Diatomic Raman. J. Chem. Phys. 1981, 75, 515−529. (24) Noid, D. W.; Koszykowski, M. L.; Marcus, R. A. A spectral analysis method of obtaining molecular spectra from classical trajectories. J. Chem. Phys. 1977, 67, 404−408. (25) Horn, P.; Kertesz, M. Conformational Preferences of βCarotene in the Confined Spaces inside Carbon Nanotubes. J. Phys. Chem. C 2010, 114, 12139−12144. (26) Putrino, A.; Parrinello, M. Anharmonic Raman Spectra in HighPressure Ice from Ab Initio Simulations. Phys. Rev. Lett. 2002, 88, 176401. (27) Horníček, J.; Kaprálová, P.; Bouř, P. Simulations of vibrational spectra from classical trajectories: Calibration with ab initio force fields. J. Chem. Phys. 2007, 127, 084502. (28) Long, D. A. Intensities in Raman Spectra. I. A Bond Polarizability Theory. Proc. R. Soc. London, Ser. A 1953, 217, 203−221. (29) Thole, B. T. Molecular polarizabilities calculated with a modified dipole interaction. Chem. Phys. 1981, 59, 341−350. (30) Smirnov, K. S.; Bougeard, D. Quantum-chemical derivation of electro-optical parameters for alkanes. J. Raman Spectrosc. 2006, 37, 100−107.

Figure 14. 1D Hessian matrix of supercell containing 6 atoms.

of two atoms. Each block is a 2 × 2 matrix of second derivatives of the total potential V, with respect to atom positions {xi} in the corresponding unit cells. (Shaded in gray is one “reference” block, corresponding to self-interactions of one unit cell.)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Scott T. Milner: 0000-0002-9774-3307 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Kalman Migler for helpful conversations and permission to reproduce Figure 1, and support from NSF-DMR 1507980 and NRT 1449785.

■

REFERENCES

(1) Strobl, G. A thermodynamic multiphase scheme treating polymer crystallization and melting. Eur. Phys. J. E: Soft Matter Biol. Phys. 2005, 18, 295−309. (2) Sirota, E. B.; Herhold, A. B. Transient phase-induced nucleation. Science 1999, 283, 529−532. (3) Wentzel, N.; Milner, S. T. Crystal and rotator phases of nalkanes: A molecular dynamics study. J. Chem. Phys. 2010, 132, 044901. (4) Sirota, E. B. Remarks concerning the relation between rotator phases of bulk n-alkanes and those of Langmuir monolayers of alkylchain surfactants on water. Langmuir 1997, 13, 3849−3859. (5) Kurelec, L.; Rastogi, S.; Meier, R. J.; Lemstra, P. J. Chain Mobility in Polymer Systems: On the Borderline between Solid and Melt. 3. Phase Transformations in Nascent Ultrahigh Molecular Weight Polyethylene Reactor Powder at Elevated Pressure As Revealed by in Situ Raman Spectroscopy. Macromolecules 2000, 33, 5593−5601. N


Article

Macromolecules (31) Abbate, S.; Gussoni, M.; Masetti, G.; Zerbi, G. Infrared and Raman intensities of polyethylene and perdeuteropolyethylene by electro-optical parameters. Single chain. J. Chem. Phys. 1977, 67, 1519−1531. (32) Gussoni, M.; Abbate, S.; Zerbi, G. Raman intensities: Transferability of electrooptical parameters. J. Raman Spectrosc. 1977, 6, 289−298. (33) Cates, D. A.; Strauss, H. L.; Snyder, R. G. Vibrational Modes of Liquid n-Alkanes: Simulated Isotropic Raman Spectra and Band Progressions for C5H12-C20H42 and C16D34. J. Phys. Chem. 1994, 98, 4482−4488. (34) Liang, Y.; Miranda, C. R.; Scandolo, S. Infrared and Raman spectra of silica polymorphs from an ab initio parametrized polarizable force field. J. Chem. Phys. 2006, 125, 194524. (35) Smirnov, K. S.; Bougeard, D. Raman and infrared spectra of siliceous faujasite. A molecular dynamics study. J. Raman Spectrosc. 1993, 24, 255−257. (36) Bougeard, D.; Smirnov, K. S. Calculation of off-resonance Raman scattering intensities with parametric models. J. Raman Spectrosc. 2009, 40, 1704−1719. (37) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (38) Frisch, M. J. et al. Gaussian 09; Gaussian Inc.: Wallingford, CT, 2009. (39) Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J. First principles methods using CASTEP. Z. Kristallogr. - Cryst. Mater. 2005, 220, 567−570. (40) Gall, M. J.; Hendra, P. J.; Peacock, O. J.; Cudby, M. E. A.; Willis, H. A. The laser-Raman spectrum of polyethylene: The assignment of the spectrum to fundamental modes of vibration. Spectrochimica Acta Part A: Molecular Spectroscopy 1972, 28, 1485−1496. (41) Stoica, P.; Moses, R. Spectral Analysis of Signals; Prentice-Hall: Upper Saddle River, NJ, 2005. (42) Kraack, H.; Deutsch, M.; Sirota, E. B. n-Alkane Homogeneous Nucleation: Crossover to Polymer Behavior. Macromolecules 2000, 33, 6174−6184. (43) Peacock, A. J. Handbook of Polyethylene: Structures, Properties, and Applications; Marcel Dekker Inc.: New York, 2000.

O


Predicting Raman Spectra of Condensed Polymer Phases from MD

Recommend Documents