General Methodology for Estimating the Stiffness of Polymer Chains

Article pubs.acs.org/Macromolecules

General Methodology for Estimating the Stiffness of Polymer Chains from Their Chemical Constitution: A Single Unperturbed Chain Monte Carlo Algorithm Panagiotis-Nikolaos Tzounis, Stefanos D. Anogiannakis, and Doros N. Theodorou* School of Chemical Engineering, National Technical University of Athens, GR 15780 Athens, Greece S Supporting Information *

ABSTRACT: The spatial dimensions and the stiffness (characteristic ratio, C∞) of polymer chains are intimately related to key macroscopic properties such as the plateau modulus and the melt viscosity. Furthermore, these molecular features are very important in the selection and design of copolymer species used in directed self-assembly lithographic processes. We have developed a general methodology for predicting the chain dimensions of any polymer chain in the unperturbed state starting from its detailed atomistic structure. The methodology is based on performing Metropolis Monte Carlo (MC) simulation, leading to equilibration of the conformational distribution of a single unperturbed polymer chain, subject only to local interactions along its backbone. To define what constitutes local interactions, the maximal topological distance of repeat units between which nonbonded forces are active is varied systematically, until a maximum in stiffness is achieved. Our methodology was validated by comparing the predicted characteristic ratios for a series of polymers against the corresponding values estimated from MC simulations of the same polymers in the melt state based on the same force field. Furthermore, we have predicted the characteristic ratios for three polymers used in directed self-assembly lithographic processes and shown that they are in good agreement with reported experimental values.

■

mean-square end-to-end distance ⟨R2⟩ scales linearly with the number of backbone bonds n. Flory’s characteristic ratio1 is defined as

INTRODUCTION Polymer chain dimensions, usually expressed by the meansquare end-to-end distance ⟨R2⟩ or the mean-square radius of gyration ⟨Rg2⟩, and conformational stiffness, usually quantified via Flory’s characteristic ratio, C∞, are related to a plethora of macroscopic properties and play a significant role in the design of new polymeric materials as well as in the development of coarse-grained models for such materials. Chains under theta (Θ) conditions behave as if they are free of nonlocal interactions between topologically distant segments along their backbones and subject only to local interactions.1 We will call such chains unperturbed. Conformation in the amorphous bulk (polymer melts, polymer glasses) can be considered as unperturbed to a very good approximation (Flory’s random coil hypothesis),2 although recent theoretical and simulation work on soft repulsive models has uncovered weak long-range intrachain correlations which lead to a power law, rather than exponential, decay of the orientational correlation function of bonds with distance along the backbone.3 For a sufficiently long unperturbed chain, the © XXXX American Chemical Society

C∞ ≡ lim

n →∞

⟨R2⟩ nl 2

(1)

2

where l is an average squared backbone bond length. Furthermore, ⟨R2⟩ and ⟨Rg2⟩ are related via ⟨R g 2⟩ = ⟨R2⟩/6

(2)

A long unperturbed chain can be mapped onto a random flight of equal steps (equivalent freely jointed chain) that has the same mean-squared end-to-end distance and contour length as the original chain. For an unperturbed chain whose backbone consists of n identical bonds of length l, with successive bonds Received: March 28, 2017 Revised: May 16, 2017

A

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

Article

Macromolecules forming an equilibrium angle π−θ with each other, the contour length is L = nl cos(θ/2). The equivalent freely jointed chain has NK steps (Kuhn segments), each of length b (Kuhn segment length) such that4−7 NK =

cos2(θ /2) L2 = n, C∞ ⟨R2⟩

b=

generator matrix techniques. The conformational properties of a plethora of polymeric systems such as polypropylene (PP),29 p o l y ( m e t h y l m e t h a c r y la t e) ( P M M A ) , 3 0 − 3 2 p o l y (dimethylsiloxane) (PDMS), 33,34 and poly(lactic acid) (PLA)35,36 were predicted using RIS theory under unperturbed conditions. In more recent years, Monte Carlo (MC) sampling has been used in lieu of the analytical computation of averages of RIS theory. Sampling can be performed either in a discretized conformation space37 or in continuous space.38 Computationally, the stiffness of polymer chains is usually extracted from atomistic simulations of either polymer melts or single unperturbed chains. Efficient equilibration of long-chain polymer melt simulations is readily achieved using connectivity altering Monte Carlo (MC) simulations of multichain models with three-dimensional boundary conditions at melt density; by analyzing the resulting equilibrated melt configurations, ⟨R2⟩, ⟨Rg2⟩, and the characteristic ratio C∞ are estimated. All interactions, bonded and nonbonded, local and nonlocal, intra- and intermolecular, are included in these melt simulations. The chain dimensions for many polymer melts, including melts of chains containing polar groups, have been estimated via connectivity altering Monte Carlo (MC) simulations.39−47 In some cases the atomistic model was coarse-grained into one involving fewer degrees of freedom and a smoother effective energy hypersurface, which can be equilibrated with smaller computational effort.48 Metropolis Monte Carlo simulations in combination with the RIS model, the so-called rotational isomeric state Metropolis Monte Carlo (RMMC) simulations, have been used to sample the unperturbed configurations of single polymer chains.49,50 Typical moves employed in such simulations are rotations of randomly chosen bonds along the backbone, which are accepted or rejected by a Metropolis selection criterion.51 Taking into account only local interactions, the configuration of a single polymer chain is finally equilibrated and the chain dimensions and the stiffness of the chain are estimated. RMMC simulations were widely used to estimate the unperturbed dimensions of many polymers, some of them possessing a complex chemical structure.52−59 Single unperturbed chain simulations are much less demanding computationally than melt simulations and therefore easier to employ. A major problem with single unperturbed chain simulations, however, is the very definition of “local” interactions. In other words, how can one define a range within which all interactions in the polymer chain can be considered as local, such that single chains subject to these interactions adopt unperturbed conformations representative of a melt? In addition, in RMMC single unperturbed chain simulations all bond lengths and bond angles of a chain are kept fixed. Thus, in RMMC simulations fully flexible force fields cannot be implemented. The work of Destrée et al.60 was a useful contribution in the direction of defining local interactions and Θ conditions. In the Metropolis Monte Carlo simulations of single polyethylene (PE) chains which were undertaken in that work, and in particular in the attractive part of the Lennard-Jones potential used to describe nonbonded interactions, a multiplicative factor β was introduced, which represents the effect of a solvent on the solvent-mediated nonlocal interactions of the single PE chain in dilute solution. The dependence of Flory’s characteristic ratio, Cn, on the number of backbone bonds, n, was examined for different values of β. The optimal value of the

C∞ ⟨R2⟩ = l L cos(θ /2) (3)

The Kuhn length is a key parameter in many coarse-grained models and theories, a notable example being field theory of inhomogeneous polymers.8 For example, in directed self assembly (DSA) lithography of block copolymers, the Kuhn length or, equivalently, the mean-square end-to-end distance for given chain length, is a very important input parameter to “single chain in mean field” calculations of the ordered structures expected under given conditions and of the thermodynamics and the kinetics of defect formation and annihilation.9−12 Chain dimensions are intimately related to key viscoelastic properties such as the plateau modulus and the entanglement molar mass Me. Through empirical equations unperturbed chain dimensions have been directly related to several viscoelastic properties. These relations were validated through the pioneering work of Fetters et al.13−15 An important parameter for both rheological properties and miscibility of athermal polymer blends16 is the packing length p, defined as the ratio of the volume occupied by the segments of a chain in a melt over the mean-square end-to-end vector of the chain: p=

υ M = K2 2 ρNA⟨R ⟩ b

(4)

where NA is Avogadro’s number and M, ρ, and υK are the chain molar mass, the mass density of the melt, and the volume corresponding to a Kuhn segment, respectively. The packing length has been successfully correlated to the entanglement molar mass through the Fetters equation Me = p3 ρNAnt 2

(5)

where NA is Avogadro’s number and nt is the number of entanglement strands in a cubed tube diameter.17 The parameter nt was found to have the universal value of 20.6 ± 8% for a series of polymers. Through the same model the entanglement tube diameter d can be related to the packing length via d = ntp.18 Experimentally, the stiffness and the unperturbed dimensions of polymer chains are mainly estimated by two methods. The first is intrinsic viscosity measurements in Θ-solutions, based on the relation between ⟨R2⟩ and the intrinsic viscosity [η].19−21 The second is small-angle neutron scattering (SANS), applicable directly in the melt or glassy state by isotopic labeling.22−24 A breakthrough in the theoretical prediction of chain conformational stiffness from chemical constitution was the development of the rotational isomeric state (RIS) theory.25−28 In this theory, conformation space is discretized by definition of a number of rotational states for each backbone bond. Statistical weights for pairs of successive bonds to assume specific pairs of states are determined through force field-based conformational analysis of small molar mass analogues. The partition function for the unperturbed chain is written as a product of statistical weight matrices and averages over all conformations, such as ⟨R2⟩ and ⟨Rg2⟩, are obtained through B


Article

Macromolecules parameter β, which corresponds to Θ conditions, was chosen as the one which leads to an asymptotic limit of Cn for high n values, i.e., the one for which the unperturbed scaling ⟨R2⟩ ∼ n1 is obeyed. At this value, βΘ, PE chains can be considered as unperturbed, subject only to “local” interactions. The multiplicative factor β was also introduced in the work of Sariban et al.,61 where Monte Carlo simulations of single polymethylene chains were performed, applying a force field with rigid bonds and bond angles, in order to successfully simulate Θ conditions. The dependence of Flory’s characteristic ratio, Cn, on the number of backbone bonds, n, was also examined for different values of β and a similar criterion, where Cn becomes constant as a function of n, was used to identify Θ conditions. In this paper we address the problem of estimating the conformational stiffness of unperturbed chains of given chemical constitution at prescribed temperature through single chain simulations based on an atomistic force field. We first describe a general methodology we have developed for estimating the chain dimensions and stiffness based on Metropolis Monte Carlo simulations of single unperturbed chains. Furthermore, we describe a variable range calculation scheme we have implemented for nonbonded interactions (both Lennard-Jones and Coulombic) in order to define which interactions can be considered as “local”. Next, we validate our methodology by comparing the stiffnesses of six polymer chains sampled with our single chain algorithm against the corresponding estimates from atomistic or coarse-grained multichain simulations in the melt state based on the same force field. Finally, we discuss the ability to apply different kinds of force fields within our methodology, provide some details on the scheme we developed to identify nonlocal interactions, and assess the performance of our algorithm in comparison with Monte Carlo simulations of polymer melts.

Figure 1. Flow diagram of the Metropolis Monte Carlo algorithm we developed to sample single unperturbed chains.

In Figure 2, we show the three basic MC moves we perform during the equilibration of a single unperturbed chain by our

■

METHODOLOGY The general methodology we have developed to estimate the dimensions and the stiffness of any polymer chain is based on performing Metropolis Monte Carlo (MC) simulations leading to equilibration of the conformational distribution of a single unperturbed polymer chain. In this section we will first describe the basic steps of the Monte Carlo algorithm we have developed. Then, we will describe how we managed to define what constitutes “local” interactions, by systematically examining the dependence of the chain stiffness on the range over which nonbonded interactions are active. In Figure 1, we show the flow diagram of the algorithm we have developed. Note that first we need to generate an initial configuration for the single chain. In our work we used the amorphous and the polymer builders of the MAPS suite to generate all the initial configurations of the polymers we studied.62 Having an initial configuration, we enter the main MC loop; in each step of this loop, one of the five available MC moves described below is randomly selected and attempted, leading to a new trial configuration of the chain. The energy difference ΔU, between the trial and the old configurations, is calculated and used in a standard Metropolis Monte Carlo acceptance criterion51 to decide if the trial configuration will be accepted or rejected. Subsequently, we move to the next iteration of the algorithm. Eventually, an ensemble of equilibrated unperturbed chains is created from which we can estimate the conformational properties of the chain such as the Flory’s characteristic ratio, the Kuhn length, etc.

Figure 2. Three basic Monte Carlo moves performed for equilibration of a single unperturbed chain: (a) initial configuration of the chain, (b) single atom displacement move, (c) flip atom move, and (d) rotate strand (pivot) move.

algorithm. In this work, as also described in the Results section, the force fields we apply are united-atom; this means that hydrogen atoms are fused into the carbon atoms with which are connected, forming pseudoatoms known or “united atoms”. Furthermore, the force fields we apply are fully flexible, which means that there are no constraints imposed on bond stretching, bond bending, or torsion around bonds. In this our models differ from RMMC simulations, where bonds and bond angles are kept fixed and a continuous set of torsional states is invoked. However, we have designed our MC procedure for use with both fully flexible chain models, free of any geometric constraints, and models with infinitely stiff bond stretching force constants, in which backbone bond lengths are fixed. The first move (see Figure 2a), applicable to fully flexible models, is a single atom displacement move, where an atom of the chain is randomly selected and displaced to a new position which is also selected randomly inside a cube with edge length Δr centered at the original atom position. The edge length of the cube is the first free parameter of our algorithm. In C


Article

Macromolecules general, higher Δr values lead to longer atom displacements, which are more difficult to be accepted, and thus to lower acceptance rates for the move. In this work we set Δr = 0.5 Å. Changes in chirality which may be a result of this move count as rejections. The second move (see Figure 2b) is called f lip atom move.63 It is applicable to both fully flexible chain models and models with fixed bond lengths. In this move a single backbone atom to which no branches are attached is randomly selected and rotated around the axis which connects the previous and the next backbone atoms by a randomly selected angle Δω, with −10° ≤ Δω ≤ 10°. During this move all bond lengths are kept fixed, and thus the bond stretching energy of the chain before and after the move remains the same. The third and the most effective move we perform is the rotate strand or pivot move.64−66 In this move a backbone bond is randomly selected, and the part of the chain which either follows or precedes the chosen bond (this is also selected randomly) is rotated around the bond by a randomly chosen angle Δφ, with −180° ≤ Δφ ≤ 180°. One can easily see that this move can alter drastically the conformation of the chain. Furthermore, during this move no bond lengths or angles are changed. Thus, during the pivot move both the bond stretching and the bond angle bending energies remain the same. The three basic moves described above were performed in all cases, for both linear and branched polymer chains. For branched polymer chains we have introduced two additional moves, which are schematically shown in Figure 3. The first one is the f lip branch move. In this move a branch is randomly selected and is flipped around the axis that is defined by the two backbone atoms flanking the backbone atom to which the branch is connected. The flip angle Δω is picked from a uniform distribution within the interval −10° ≤ Δω ≤ 10°. Following the move, the local geometry is checked to ensure that the chirality of the backbone atom bearing the branch has not changed. Changes in chirality count as rejections. The second branch move is the rotate branch move, where a branch is randomly selected and a branch segment is rotated around a randomly selected rotatable bond along the branch. One could imagine this move as the pivot move which acts along branches. These two moves can effectively relax the conformations of the branches of both short and long branched polymer chains. As implemented, all moves described above satisfy the condition of microscopic reversibility.67 Perhaps the most important contribution of our methodology is the scheme we have implemented to determine what constitutes local interactions. Figure 4 shows schematically the conceptual difference between local and nonlocal interactions. By definition, an unperturbed chain is subject only to local interactions, such as the ones denoted in Figure 4a. If all nonlocal, nonbonded interactions were included in the Hamiltonian of a long chain, the chain would shrink, departing from ⟨R2⟩ ∼ n1 scaling, vacuum being a very bad solvent for it. Below a certain temperature, the chain would collapse and become a globule or even crystallize.68 Including only local interactions in the Hamiltonian yields ⟨R2⟩ ∼ n scaling. As amply demonstrated in past simulation work, however,46,69 the characteristic ratio C∞ (eq 1) is a function of the range assumed for the local interactions. How should the range of local interactions be defined in the single chain atomistic simulation for the sampled conformations to be representative of actual unperturbed chains encountered in melts?

Figure 3. Additional two Monte Carlo moves performed during the equilibration of a single unperturbed chain with branches: (a) initial configuration of the chain, (b) flip branch move, and (c) rotate branch move.

Figure 4. Schematic representation of local and nonlocal interactions in a single polymer chain.

In order to quantify the range of nonbonded forces included in our atomistic single-chain model, we introduce the parameter Δnpair. This parameter corresponds to the maximum topological distance along the backbone between electrically neutral groups (usually repeat units) for which nonbonded interactions, both van der Waals and Coulomb, are active with respect to the force field we have chosen to represent the polymer. Atoms of a D


Article

Macromolecules given group i interact with all atoms belonging to groups with index numbers j ∈ [i − Δnpair, i + Δnpair] along the chain. To illustrate this, we show in Figure 5a a schematic representation

are in very good agreement with the corresponding estimates extracted from bulk melt simulations based on the same force field. Further discussion of the criterion is provided in the Results section. The choice of a maximum topological distance along the chains (Δnpair) beyond which nonlocal interactions are switched off corresponds to an abrupt cutoff of the longrange interactions. One could introduce a more gradual screening of the long-range interactions. For example, in the work of Destrée et al.60 a modification of the attractive term of the Lennard-Jones potential, which is used to describe the van der Waals interactions, is introduced by a multiplicative factor, β, as we mentioned in the Introduction. The value of β was systematically varied until an optimal value was found which corresponds to Θ conditions. A similar strategy was followed by Sariban et al.61 where the same multiplicative factor, β, was used to continuously screen long-range interactions. However, in these works, a simple, linear, nonpolar polymer (PE) was studied. The extension of this methodology to polymer chains with more complicated chemical structure, and especially when partial charges are included, would not be trivial. On the other hand, in terms of a topological criterion, one could introduce a switching function: ⎧1 n ≤ n0 ⎪ f (n) = ⎨ ⎛ n ⎞ ⎟ n > n0 ⎪ exp⎜ − ⎩ ⎝ nmax ⎠

(6)

Figure 5. (a) Schematic representation of a single PDMS chain segment. With red dashed circles we mark the groups (neutral repeat units) with which the ith group is allowed to interact, given that Δnpair = 2. (b) Dependence of the characteristic ratio, C∞, on Δnpair for PDMS at 298 K. The insets display snapshots of the PDMS chains corresponding to each of the Δnpair values considered. All snapshots are captured with VMD software (v.1.9.3).70

where n is the topological distance of two monomers along the polymer chain, n0 is the maximal topological distance where Lennard-Jones and Coulomb interactions are fully taken into account, and nmax is a screening distance. Thus, the pairwise potential could be modified as

of a PDMS chain where we mark with red dashed circles the neutral groups with which group i is allowed to interact for Δnpair = 2. The higher the value of Δnpair, the longer the range of nonlocal interactions taken into account. Baschnagel et al.71 introduced a topological parameter similar to Δnpair, called Δmax, in order to study the effect of the range of nonbonded interactions on the unperturbed dimensions of single PE chains. By generating various configurations of single PE chains using a simple Monte Carlo sampling technique, and applying a fully flexible force field, Baschnagel et al. estimated the characteristic ratio of PE and its temperature dependence, but unfortunately without proposing a criterion for choosing the optimal Δmax value. In order to determine the optimal value of Δnpair for a given polymer chain, we first calculate the stiffness of the chain as a function of Δnpair. As we observe in Figure 5b, different values of Δnpair lead to different estimates for the characteristic ratio C∞ of the PDMS chain. As we will also show in the Results section, however, there is a value of Δnpair where the stiffness becomes maximal. For higher Δnpair values more attractive nonlocal interactions become active and the chain begins to shrink until it collapses on itself, as shown in the inset PDMS chain snapshots of Figure 5b. We assert that the optimal value of Δnpair is that where the stiffness becomes maximal. This is a rather empirical criterion; nevertheless, as we will show in the Results section, estimates of the characteristic ratio for a series of polymers based on single chain sampling using this criterion

where r is the Euclidean distance between two interacting atoms and ULJ and UCoul are the Lennard-Jones and Coulomb energy, respectively. With this modification, a continuous screening of long-range interactions is achieved. However, using this approach, two new free parameters are introduced, making a systematic study of their optimal values more complicated. In addition, the need for a gradual screening of the long-range interactions increases the computer time, making a single chain simulation more computationally demanding. The approach introduced here, wherein a Heaviside-like step function is assumed to describe the screening of long-range interactions, is very simple to implement, less computationally demanding and, as we will show and discuss in the following sections, able to reproduce successfully the unperturbed dimensions of linear or branched polymer chains, whether they contain polar groups or not.

Upair(r ) = f (n)[ULJ(r ) + UCoul(r )]

(7)

■

RESULTS In this section we present our estimates of the stiffness (characteristic ratio), as obtained by applying our single chain Monte Carlo methodology to six different polymers. Note that in all figures where the characteristic ratio (Cn) is plotted versus the number of backbone bonds (n), error bars are very small, comparable to the width of the lines, and therefore are omitted. i. Polyethylene (PE). The first system we studied was a linear polyethylene (PE). The main reason we used this system is to validate our methodology by comparing our estimate of E


Article

Macromolecules

isotactic polypropylene (i-PP). Using the MAPS platform suite,62 the initial configuration of a single i-PP chain consisting of 500 propylene repeat units was generated. The force field we applied was the one proposed by Logotheti and Theodorou.79,80 Skeletal carbon atoms and pendant hydrogens are represented explicitly, while the methyl branch groups are represented as united atoms. Harmonic bond stretching and bond angle bending potentials and a 3-fold symmetric torsional potential are used. Local interactions between atoms which are one or two bonds apart are forbidden, while for atoms which are three or four bonds apart, local nonbonded interactions are described by a Lennard-Jones 12−6 potential using a different set of σij and εij from the one used for interactions between atoms which are more than four bonds apart or belong to different chains. We have performed long single chain MC runs of about 1 billion steps, at T = 450 K, with Δnpair from 2 to 6. Here, Δnpair refers to interactions between groups which consist of a single propylene unit. As we observe in Figure 7, for all Δnpair, Cn increases with n until it reaches a plateau value. Here, the plateau value is

the characteristic ratio against available estimates from Monte Carlo simulations of linear PE melts. The force field we applied was TraPPE-UA, developed by Martin and Siepmann,72,73 in which there are interactions between united-atoms, each united atom consisting of a skeletal carbon fused with the hydrogens attached to it. For the torsional potential we applied the expression invoked in the anisotropic united-atom (AUA) model proposed by Toxvaerd,74 as the original TraPPE-UA torsional potential has been shown to overestimate the stiffness.75 Furthermore, instead of fixing the bond length, we applied a harmonic bond stretching potential of the form 1/2Kb(l − l0)2 with a large stiffness parameter Kb = 1000 kcal/ (mol Å2) and l0 = 1.54 Å. Therefore, the force field we used is fully flexible. Using the MAPS platform suite,62 the initial configuration of a single linear PE chain was generated. The chain consisted of 500 united atoms (498 CH2 and 2 CH3 united atoms located at the ends). We performed long single chain MC runs of about 1 billion steps, at T = 450 K, with varying Δnpair from 1 to 5. Here, Δnpair refers to interactions between groups, each group consisting of an ethylene unit −CH2−CH2−. For example, Δnpair = 1 means that for a given ethylene group, i, nonlocal interactions are active only between ethylene groups i, i − 1, and i + 1. In addition, nonbonded interactions are forbidden by the force field between united atoms that are connected by up to three bonds.76 In Figure 6, we plot the characteristic ratio of subchains of the simulated chain, Cn, as a function of the number of

Figure 7. Characteristic ratio, Cn, of isotactic PP as a function of the number of backbone bonds n, for different values of Δnpair at 450 K.

maximal for Δnpair = 3, unlike PE, where the maximum is observed for Δnpair = 2. A possible explanation of this effect could be that the presence of the methyl branches in PP leads to the selection of a wider range in which interactions can be considered as “local”. In other words, the range of Δnpair = 2 which corresponds to PE is not enough to fully describe all “local” interactions in PP. Following the same empirical rule we used to estimate the stiffness of PE, we obtain the characteristic ratio of isotactic PP, C∞ = 5.90, which is in very good agreement with the estimate C∞ = 5.49 ± 0.33 from molecular dynamics melt simulations based on the same force field79 and a little underestimated but still in reasonable agreement in comparison with the experimental value C∞ = 6.2.13,81 Moving further, we would like to check the ability of our algorithm to capture the stiffness dependence on different tacticities. To this end, we have also generated the initial configurations of a single syndiotactic s-PP chain and a single atactic (50% meso and 50% racemo dyads, Bernoullian tacticity sequence) a-PP, consisting of 500 propylene repeat units each. Applying the same force field, we performed long single chain MC simulations of about 1 billion steps. For Δnpair = 3, which is also the optimum Δnpair value for both a-PP and s-PP, we compare, in Figure 8, the estimated stiffness of the three PP chains we studied at a temperature of 650 K. We find that s-PP is the stiffest chain with C∞ = 12.3, while a-PP (C∞ = 7.9) is

Figure 6. Characteristic ratio, Cn, of linear PE as a function of the number of backbone bonds, n, for different values of Δnpair at 450 K. The error bars in this plot and in all other similar figures are comparable to the thickness of the lines.

backbone bonds in the subchain, n, for different Δnpair values. We observe that for all Δnpair Cn increases with n, until it reaches a plateau value. We also note that for Δnpair = 2 the plateau value reaches a maximum, while for Δnpair > 2 the plateau value decreases with Δnpair. We choose to estimate the characteristic ratio, C∞, from the maximum plateau value, C∞ = 8.25, before the chain begins to shrink. Our estimate is in excellent agreement with the estimated value of 8.26 ± 0.15, obtained from Monte Carlo simulations of linear PE melts,77,78 applying the same force field as we did. Therefore, the empirical rule we use to estimate C∞ is seen to work in this case. The estimated C∞ values from simulations are not so close to the experimental measurements, which yield C∞ = 7.3.13 ii. Polypropylene (PP). In order to further validate our Monte Carlo methodology, we chose to estimate the stiffness of a branched polymer bearing simple short branches, namely F


Article

Macromolecules

Figure 8. Comparison of the characteristic ratios, Cn, of the syndiotactic, atactic, and isotactic forms of PP as a function of the number of backbone bonds, n, for Δnpair = 3 at 650 K.

Figure 9. Characteristic ratio, Cn, of atactic PS as a function of the number of backbone bonds, n, for different values of Δnpair at 450 K.

large side groups one should better use a wider range over which nonbonded interactions are active.50,52 Finally, from Δnpair = 9, where the plateau reaches its maximum, we can estimate the characteristic ratio of a-PS, C∞ = 8.6, which is in very good agreement with the estimate C∞ = 8.7 ± 0.1 from coarse-grained MC simulations83 and in excellent agreement with the experimental value of C∞ = 8.584 at the same temperature. iv. Poly(ethylene oxide) Dimethyl Ether (PEODME). Having validated our methodology by comparing our estimates for the characteristic ratio of PE, i-PP, and a-PS against various estimates from melt MC, MD, and coarse-grained MC simulations, respectively, we turn to a linear polymer which contains polar groups that require inclusion of Coulombic interactions in the force field for its proper description. The polymer we chose was poly(ethylene oxide) dimethyl ether (PEODME), i.e., poly(ethylene oxide) (PEO) terminated with methyl groups. The force field we applied in our simulations is the same as in the work of Wick and Theodorou.45 Electrostatic Coulomb interactions are included but are prohibited between atoms separated by less than three bonds. For atoms separated by exactly three bonds Coulomb interactions are taken into account, but reduced in magnitude by a scaling factor of 0.5. On the other side, Lennard-Jones interactions are forbidden between atoms separated by fewer than four bonds. A single PEODME chain consisting of 250 monomers was generated using the MAPS platform suite.61 We performed long single chain MC simulations (up to 1 billion steps) at T = 450 K and estimated the characteristic ratio as a function of the number of backbone bonds for Δnpair = 1 to Δnpair = 4. Here, Δnpair refers to interactions between neutral ethylene oxide repeat units, −CH2−CH2−O−. We note here that the choice of using electrically neutral groups for the definition of local interactions is very important. Trial runs with different definitions of the interacting groups that violated electroneutrality led to unrealistic configurations for various values of Δnpair, e.g., caused the PEODME to collapse upon itself or become very stiff with C∞ > 20. In Figure 10, we observe that the plateau value of the characteristic ratio plots is maximized for Δnpair = 2, as in the case of PE. This behavior supports further the argument that the optimal Δnpair depends strongly on the presence of side groups, as PEODME is linear and has no side groups. From our single chain simulation of PEODME with the optimal value Δnpair = 2 we estimate C∞ = 5.25, which is very close to the estimate C∞ = 5.4−5.6 obtained from connectivity-altering MC

stiffer than i-PP (C∞ = 5.53 at this temperature). Qualitatively, the results are reasonable as s-PP is far stiffer than both i-PP and a-PP. Quantitatively, the stiffness of both s-PP and a-PP is overestimated in comparison with the experimental results which are 9.1 and 6.0,13,81 respectively, possibly due to the force field and in particular due to the potential which is used to describe the dihedral interactions. Unfortunately, there are no available estimations for the stiffness of s-PP and a-PP from atomistic melt simulations using the same force field. However, we may attempt a comparison with melt simulations using a similar force field69 in which s-PP is once again the stiffest, with C∞ = 8.5 estimated from Monte Carlo melt simulations and C∞ = 9.6 estimated from continuous unperturbed chains, while the characteristic ratio of a-PP is C∞ = 6.2 estimated from Monte Carlo melt simulations and C∞ = 5.2 estimated from continuous unperturbed chains. Once again, qualitatively, our results are in agreement with past simulations, but quantitative differences are most probably attributable to the different force field we applied. iii. Atactic Polystyrene (a-PS). The next polymer for which we chose to estimate stiffness was atactic polystyrene. A single atactic polystyrene chain consisting of 750 styrene units (50% meso and 50% racemo dyads, Bernoullian tacticity sequence) was generated using the MAPS platform suite62 and used as initial configuration for a single chain MC run of 7 billion steps at 450 K. The force field developed by Harmandaris et al.82,83 was applied in our Monte Carlo simulations. It is a united atom force field in which improper torsion contributions are included to keep the phenyl ring planar and also to maintain the tetrahedral structure between the phenyl ring and its connection with the backbone. In Figure 9, we observe that the characteristic ratio reaches a plateau value for all Δnpair at higher n (n > 1000) in comparison with PE and i-PP. Furthermore, the maximum of the plateau is reached at Δnpair = 9, which is higher than the corresponding values in PE and i-PP. We believe that this effect is mainly due to the presence of the large phenyl groups which leads to the need for a wider range, and thus higher Δnpair, to capture all local interactions. In general, as we observe from our results on PE, i-PP, and a-PS the presence of large side groups in a polymer results in an increase of the Δnpair value at which the characteristic ratio plateau is maximized. A similar observation seems to have been made in RMMC simulations. Although no systematic study of the effect of the range where the interactions are considered as “local” has been presented in connection with RMMC, it is proposed that for polymers with G


Article

Macromolecules

Unfortunately, there are no estimates for the C∞ of PDMS from either atomistic or coarse-grained simulations. However, we can compare our estimates of Cn for n = 40, where C40 = 4.9, with the corresponding estimations from MD simulations performed by Curro et al.86 using the same force field, in which C40 = 5.28. We note that the two predictions are reasonably close. We can also compare our predictions against the experimental estimate C∞ = 5.81.91 The deviation from experiment could be due to the force field we applied in our simulations and its possible difficulty to realistically predict the conformational properties of PDMS. vi. Atactic Poly(methyl methacrylate) (a-PMMA). The last polymer we have studied was atactic PMMA, which is widely used in directed self-assembly lithography and is one of the best known polyacrylates. The force field we applied to simulate a-PMMA was TraPPE-UA. Thanks to the transferability of the force field, we have managed to collect all the parameters we needed (bond stretching, bond bending, torsional, etc.) from a series of papers, where simulations of carboxyl esters, ethers, and carboxyl acids were performed.92−96 A single atactic (50% meso and 50% racemo dyads, Bernoullian tacticity sequence) PMMA chain consisting of 250 monomers was generated using the MAPS platform suite.62 Long single chain MC simulations (up to 1 billion steps) were performed at T = 413 K. The characteristic ratio versus the number of backbone bonds for Δnpair = 2 to Δnpair = 6 was estimated. In PMMA, we chose Δnpair to refer to interactions between neutral MMA repeat units. In Figure 12, we observe that the plateau values of characteristic ratio plots become maximal at Δnpair = 3. For

Figure 10. Characteristic ratio, Cn, of PEODME as a function of the number of backbone bonds, n, for different values of Δnpair at 450 K.

simulations of PEODME melt performed by Wick and Theodorou45 using the same force field and to the experimental value C∞ = 5.5.85 Thus, our methodology seems applicable to polymers containing polar groups for which Coulombic interactions must be included in the force field. v. Polydimethylsiloxane (PDMS). Another polymer with a wide range of applications which we chose to estimate conformational stiffness by applying our Monte Carlo algorithm is PDMS. It is a silicon-based organic polymer containing Si−O dipoles along its backbone. Each Si is connected to two methyl groups, making the PDMS structure a very interesting candidate to test our methodology. We chose to apply a united-atom force field developed by Curro et al.86−89 in our simulations. This is a hybrid force field which mixes class I (e.g., Lorentz−Berthelot) with class II90 mixing rules and 9−6 with 12−6 Lennard-Jones equations to describe nonbonded interactions. A single PDMS chain consisting of 250 monomers was generated using the MAPS platform suite.62 Long single chain MC simulations (up to 1 billion steps) were performed at T = 298 K, and the characteristic ratio versus the number of backbone bonds for Δnpair = 1 to Δnpair = 4 was estimated. In PDMS, we chose Δnpair to refer to neutral dimethylsiloxane repeat units −Si(CH3)2O− following the same strategy as in PEODME. In Figure 11, we observe that the plateau value of the characteristic ratio plots is maximized for Δnpair = 2, as in PE and PEODME, from which we estimate C ∞ = 5.28.

Figure 12. Characteristic ratio, Cn, of PMMA as a function of the numbers of backbone bonds, n, for different values Δnpair at 413 K.

these values of Δnpair we estimate C∞ = 7.35. Unfortunately, there are no estimates for C∞ from either atomistic or coarsegrained simulations. However, we can compare our estimate directly against the corresponding experimental measurement at 413 K, C∞ = 8.22,91 and we note reasonable agreement.

■

TEMPERATURE DEPENDENCE OF STIFFNESS In this section, we study the temperature dependence of the characteristic ratio C∞ for the six polymers we have simulated. To this end, we have performed single chain MC simulations at three different temperatures for each polymer at the optimal Δnpair values we found previously and calculated the temperature coefficient

Figure 11. Characteristic ratio, Cn, of PDMS as a function of the number of backbone bonds, n, for different values of Δnpair at 298 K. The green line for Δnpair = 4 corresponds to a PDMS which has collapsed on itself and is a globule. A snapshot of this chain is available in the Methodology section, in the inset of Figure 5b.

κ= H

d ln(C∞) d ln⟨R2⟩ = dT dT

(8) DOI: 10.1021/acs.macromol.7b00645 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

force fields we apply. In any case, our method could prove to be a valuable tool which can help in the direction of force field development and parametrization.

In Figure 13, we plot the logarithm of the characteristic ratio ln(C∞) versus the temperature for the six polymers we studied.

■

SUMMARY AND DISCUSSION In this section, we will first summarize the results of this work by presenting tables and figures which compare the estimated properties for all systems studied. Then, we will discuss the ability of our methodology to support any force field, united atom or not. Furthermore, we will focus on the scheme we have introduced to systematically examine the effect of nonlocal interactions on the stiffness of polymer chains. Finally, we will compare the simulation times needed to estimate the stiffness of a long PE chain by applying our single chain MC algorithm and by performing a MC simulation of a PE melt. In Table 2 we present the estimated stiffness for the six polymer chains we studied, and we first compare our estimates against the results from the corresponding melt simulations (either atomistic or coarse-grained), where available, using the same force field. The agreement between the results is very good and strongly validates our single chain methodology. Next, we make an attempt for a comparison against experimental estimates of stiffness from SANS measurements for all the polymer systems we studied. With the exception of PE and possibly PDMS, the general agreement is good. The overestimated and underestimated results for PE and PDMS, respectively, can be explained by the force fields and mainly by the torsional potentials which were applied. In Figure 14, we show the dependence of the characteristic ratio C∞ of the six polymer systems we studied on the range

Figure 13. Logarithm of the characteristic ratio as a function of the temperature for all six polymers we studied. Lines represent linear fits on the data.

We observe that there is a linear behavior between ln(C∞) and temperature over the ranges studied. We estimate κ values from the slopes of the linear fits performed on this plot. In Table 1, we present the estimated temperature coefficients for the six polymer chains we studied, and we compare our Table 1. Temperature Parameters κ Estimated from Our Single Chain Metropolis Monte Carlo Methodology and Comparison with Reported Experimental Results polymer

κ × 103 (K−1) (our methodology)

κ × 103 (K−1) (exptl, SANS)

PE i-PP a-PS PEODME PDMS a-PMMA

−0.95 −0.56 0.10 −0.70 −1.09 −0.96

−1.06 ± 0.0797 0.081 0.0 ± 0.198 −0.324 N/A 0.1 ± 0.198

estimates against the results from the corresponding SANS measurements. We note that for PE and a-PS there is a very good agreement between our simulations and SANS measurements. Also, there is reasonable agreement for PEODME. On the contrary, there is a deviation between our estimates and the experimental ones for i-PP and a-PMMA. Unfortunately, we have not found any experimental results for PDMS. In the case of a-PMMA, differences from experiment may be due to the tacticity in the experimental polymer not being exactly Bernoullian as postulated here. In general, the discrepancies we observe between our simulations and experiments concerning the magnitude and temperature coefficient of stiffness for specific polymers are probably mainly due to the

Figure 14. Characteristic ratio, C∞, as a function of the different values of Δnpair for the six polymer we studied.

over which nonbonded interactions are considered as “local”, expressed by the Δnpair parameter. It is evident that for all polymers there is a value of Δnpair where the stiffness is maximized. According to our criterion, this value corresponds

Table 2. Characteristic Ratio (Stiffness) Values Estimated from Our Single Chain Metropolis Monte Carlo Methodology and Comparison with Reported Experimental Results and Results from Atomistic Simulations of the Corresponding Melt Systems polymer

T (K)

⟨l⟩ (Å)

Δnpair

C∞ (our methodology)

C∞ (melt simulations)

C∞ (exptl)

PE i-PP a-PS PEODME PDMS a-PMMA

450 450 450 450 298 413

1.54 1.51 1.54 1.46 1.64 1.54

2 3 9 2 2 3

8.25 5.90 8.60 5.25 5.28 7.35

8.26 ± 0.1577,78 5.49 ± 0.3379 8.70 ± 0.183 5.4045 N/A N/A

7.313 6.213,81 8.584 5.524,85 5.891 8.291

I


Article

Macromolecules

effectively screened. In RMMC simulations a parameter similar to Δnpair, called Maxbonds, was introduced. Unfortunately, in many works52−59 the value of Maxbonds was chosen arbitrarily. Furthermore, no systematic analysis of the effect of this value on the conformational properties of the sampled polymer chains was undertaken. Interestingly, although the optimal Maxbonds value was empirically chosen in RMMC simulations, it was proposed that for chains with bulky side groups attached to their contours Maxbonds should be greater than for chains with less complex structure. We also noted this effect when comparing the optimum Δnpair values found in PE, i-PP, and PS systems. As expected, the simulation times for the estimation of the characteristic ratio using our single chain sampling methodology are much lower than the ones needed to perform full MC simulations of the corresponding polymer melts. For example, we need approximately 2 min CPU time on an Intel(R) Core(TM) i7-4770 CPU at 3.40 GHz to estimate the characteristic ratio of a PE chain consisting of 500 CH2 backbone units, as is shown in the video provided in the Supporting Information, while we need approximately 1 week on the same machine to fully equilibrate a corresponding PE melt consisting of 16 chains with 500 CH2 backbone units per chain and then estimate the characteristic ratio.

to the correct range where all interactions of the chain which are taken into account are indeed “local”, and the chain can be considered as unperturbed. Higher values of Δnpair lead to shrinkage of the chains until they finally collapse into a globule, as was observed in PDMS. Another interesting result from this plot is that the optimal Δnpair for PS is substantially higher than for the other polymers. As we already discussed in the Results section, this effect could be explained by the presence of the bulky phenyl groups which leads to the need for a wider range, and thus higher Δnpair, of local interactions. As we described in the previous section, the force fields we applied to validate our methodology by predicting the characteristic ratios of six polymer chains were mainly united atom force fields (with the exception of the hybrid i-PP force field), cast in terms of pseudoatoms defined by fusing hydrogen atoms with the carbon atoms to which they are connected. Our choice of united-atom representations was mainly motivated by our intention to compare against full melt simulations carried out with the same force field. However, our methodology can also support any explicit atom (all-atom) force field (OPLSAA,99−101 DREIDING,102 etc.) or other more sophisticated force fields, such as class II force fields (COMPASS,103 TEAMFF,104 etc.) without any restrictions, performing the same MC moves as described in the Methodology section. The application of a more detailed force field inevitably leads to an increase in the simulation times needed to fully equilibrate the conformational distributions of the unperturbed chains. Simulation times for single chains, however, are much lower than those needed to perform well-equilibrated simulations of the corresponding polymer melts (see below); thus, increasing the level of detail in the representation of single chains is certainly viable. The range over which nonlocal interactions are active in our simulations, expressed by the Δnpair parameter, refers to interactions between groups of atoms, usually repeat units, which should be electrically neutral. A different choice may lead to unrealistic chain conformations, e.g., chains which collapse onto themselves or very stiff chains. These unrealistic conformations result from attractive or repulsive Coulomb forces, respectively, between uncompensated charges (monopoles). Such forces become dominant if charged groups are defined for the calculation of nonbonded interactions. Therefore, one should be very careful and partition the unperturbed chains in neutral groups in order to calculate correctly the local interactions. A similar strategy, where interactions between neutral groups are involved in single chain MC simulations, was also followed in Blomqvist et al.’s105,106 single chain RIS Metropolis MC simulations, designed for studying the configurational properties of polylactic acids, poly(glycolic acid)s, and side group polyesters. As mentioned in the Methodology section, the criterion we propose for choosing the optimal range of nonlocal interactions in estimating the stiffness of unperturbed polymer chains is empirical. It could be imagined as a limiting case where the contour of the polymer chain “unfolds” to a maximal extent subject to local interactions along its contour, right before attractive interactions between topologically more distant segments begin to shrink it until it collapses for high Δnpair values. It turns out that this criterion of maximal unfolding yields conformations comparable to those encountered in a melt, where nonlocal intramolecular interactions are compensated by interactions with other, surrounding chains and excluded volume effects on the chain conformation are

■

CONCLUSIONS In this paper we have described a general methodology for predicting the stiffness and the unperturbed chain dimensions of any polymer chain by conducting single chain Metropolis Monte Carlo simulations. Five Monte Carlo moves were introduced to effectively relax both linear and branched polymer chains. The most effective moves were found to be the rotate strand (pivot) and the rotate branch moves, which can drastically change the conformations of long strands and branches along the polymer chains. Furthermore, we have proposed and implemented a variable range calculation scheme to analyze nonbonded interactions (both van der Waals and electrostatic). Our scheme is based upon defining electrically neutral groups (e.g., monomer units) along the chain and progressively increasing the range of nonbonded interactions Δnpair between pairs of such groups along the chain. Through this scheme we have managed to precisely define “local” nonbonded interactions as corresponding to that Δnpair value which maximizes the predicted chain stiffness and use them to sample single unperturbed chain conformations comparable to those encountered in melts and glasses. Finally, we have validated our methodology by comparing the predicted characteristic ratios of a series of polymers, both linear and branched, polar and nonpolar, against the ones estimated from atomistic simulations carried out in the melt state with the same force field, from coarse-grained melt simulations, or from SANS experiments. The agreement was generally very good. Apart from rapidly providing estimates of the characteristic ratio and the Kuhn length based on atomic-level chemical structure, our single chain simulations are useful in evaluating different force fields in terms of their ability to reproduce polymer chain dimensions in the amorphous bulk. In current work, our methodology is used in the direction of quickly screening among already synthesized or hypothetical polymers of various chemical constitutions and architectures in terms of the size of their coils in three-dimensional space, without the need to perform full-blown condensed phase atomistic or coarsegrained simulations. J


Article

Macromolecules

■

(12) Li, W.; Müller, M. Thermodynamics and Kinetics of Defect Motion and Annihilation in the Self-Assembly of Lamellar Diblock Copolymers. Macromolecules 2016, 49, 6126−6138. (13) Fetters, L. J.; Lohse, D. J.; Garcia-Franco, C. A.; Brant, P.; Richter, D. Prediction of Melt State Poly(α-olefin) Rheological Properties: The Unsuspected Role of the Average Molecular Weight per Backbone Bond. Macromolecules 2002, 35, 10096−10101. (14) Fetters, L. J.; Lohse, D. J.; Graessley, W. W. Chain Dimensions and Entanglement Spacings In Dense Macromolecular Systems. J. Polym. Sci., Part B: Polym. Phys. 1999, 37, 1023−1033. (15) Fetters, L. J.; Lohse, D. J.; Richter, D.; Witten, T. A.; Zirkel, A. Connection between Polymer Molecular Weight, Density, Chain Dimension, and Melt Viscoelastic Properties. Macromolecules 1994, 27, 4639−4647. (16) Fredrickson, G. H.; Liu, A. J.; Bates, F. S. Entropic Corrections to the Flory-Huggins Theory of Polymer Blends: Architectural and Conformational Effects. Macromolecules 1994, 27, 2503−2511. (17) Lin, Y.-H. Number of Entanglement Strands per Cubed Tube Diameter, a Fundamental Aspect of Topological Universality in Polymer Viscoelasticity. Macromolecules 1987, 20, 3080−3083. (18) Theodorou, D. N. Equilibration and coarse-graining methods for polymers. In Computer Simulations in Condensed Matter: From Materials to Chemical Biology; Ferrario, M., Ciccotti, G., Binder, K., Eds.; Springer: Berlin, 2006; pp 419−448. (19) Xu, Z.; Hadjichristidis, N.; Carella, J. M.; Fetters, L. J. Characteristic Ratios of Atactic Poly(vinylethylene) and Poly(ethylethylene). Macromolecules 1983, 16, 925−929. (20) Kawaguchi, S.; Imai, G.; Suzuki, J.; Miyahara, A.; Kitano, T.; Ito, K. Aqueous solution properties of oligo - and poly(ethylene oxide) by static light scattering and intrinsic viscosity. Polymer 1997, 38, 2885− 2891. (21) Zioga, A.; Ekizoglou, N.; Siakali-Kioulafa, E.; Hadjichristidis, N. Characteristic ratio of poly(tetrahydrofurfuryl acrylate) and poly(2ethylbutyl acrylate). J. Polym. Sci., Part B: Polym. Phys. 1997, 35, 1589− 1592. (22) Horton, J. C.; Squires, G. L.; Boothroyd, A. T.; Fetters, L. J.; Rennie, A. R.; Glinka, C. J.; Robinson, R. A. Small-Angle Neutron Scattering from Star-Branched Polymers in Molten State. Macromolecules 1989, 22, 681−686. (23) Fetters, L. J.; Graessley, W. W.; Krishnamoorti, R.; Lohse, D. J. Melt Chain Dimensions of Poly(ethylene - 1-butene) Copolymers via Small Angle Neutron Scattering. Macromolecules 1997, 30, 4973− 4977. (24) Smith, G. D.; Yoon, D. Y.; Jaffe, R. L.; Colby, R. L.; Krishnamoorti, R.; Fetters, L. J. Conformations and Structures of Poly(oxyethylene) Melts from Molecular Dynamics Simulations and Small-Angle Neutron Scattering Experiments. Macromolecules 1996, 29, 3462−3469. (25) Flory, P. J. Statistical Mechanics of Chain Molecules; WileyInterscience: New York, 1969. (26) Flory, P. J. Foundations of Rotational Isomeric State Theory and General Methods for Generating Configurational Averages. Macromolecules 1974, 7, 381−392. (27) Mattice, W.; Suter, U. W. Conformational Theory of Large Molecules: The Rotational Isomeric State Model in Macromolecular Systems; Wiley-Interscience: New York, 1994. (28) Helfer, C. A.; Mattice, W. L. The Rotational Isomeric State Model. In Physical Properties of Polymers Handbook; Mark, J. E., Ed.; Springer: New York, 2007; pp 43−57. (29) Suter, U. W.; Flory, P. J. Conformational Energy and Configurational Statistics of Polypropylene. Macromolecules 1975, 8, 765−776. (30) Sundararajan, P. R.; Flory, P. J. Configurational Characteristics of Poly(methyl methacrylate). J. Am. Chem. Soc. 1974, 96, 5025−5031. (31) Yoon, D. Y.; Suter, U. W.; Sundararajan, P. R.; Flory, P. J. Conformational Characteristics of Poly(methyl acrylate). Macromolecules 1975, 8, 784−789. (32) Sundararajan, P. R. Conformational Analysis of Poly(methyl methacrylate). Macromolecules 1986, 19, 415−421.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.7b00645. United-atom force field applied in single-chain MC PMMA simulations; (b) characteristic ratio, Cn(n) plots for different values of Δnpair for i-PP, a-PP, and s-PP at 650 K (PDF) Video displaying the convergence of the C n (n) calculation with Δnpair = 2 as obtained from our single chain MC algorithm applied to C500 PE at 450 K (AVI)

■

AUTHOR INFORMATION

Corresponding Author

*(D.N.T.) E-mail: [email protected]. ORCID

Doros N. Theodorou: 0000-0002-4763-9739 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/ 2007−2013) under grant agreement no. 619793 CoLiSA.MMP. Scienomics SARL is gratefully acknowledged for making their software available to us through the Scienomics Group of Scientific Excellence (SGSE) program. We are grateful to Eleftherios Koulierakis for his help with some calculations. We also thank Dr. Nikolaos A. Romanos for his help with the force field we applied to i-PP simulations.

■

REFERENCES

(1) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press: New York, 1986. (2) Flory, P. J. Principles of Polymer Chemistry; Cornell University Press: Ithaca, NY, 1953. (3) Wittmer, J. P.; Beckrich, P.; Meyer, H.; Cavallo, A.; Johner, A.; Baschnagel, J. Intramolecular long-range correlations in polymer melts: The segmental size distribution and its moments. Phys. Rev. E 2007, 76, 011803. (4) Flory, P. J. Statistical Mechanics of Chain Molecules; Wiley: New York, 1969. (5) Xu, Z.; Hadjichristidis, N.; Fetters, L. J.; Mays, J. W. Structure/ chain-flexibility relationships of polymers. In Physical Properties of Polymers; Springer: Berlin, 1995; pp 1−50. (6) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press: New York, 2003. (7) Wu, S. Chain Structure and Entanglement. J. Polym. Sci., Part B: Polym. Phys. 1989, 27, 723−741. (8) Fredrickson, G. F. The Equilibrium Theory of Inhomogeneous Polymers; Oxford University Press: New York, 2013. (9) Li, W.; Nealey, P. F.; de Pablo, J. J.; Müller, M. Defect Removal in Course of Directed Self-Assembly is Facilitated in the Vicinity of the Order-Disorder Transition. Phys. Rev. Lett. 2014, 113, 168301. (10) Hur, S.-M.; Khaira, G. S.; Ramírez-Hernández, A.; Müller, M.; Nealey, P. F.; de Pablo, J. J. Simulation of Defect Reduction in Block Copolymer Thin Film by Solvent Annealing. ACS Macro Lett. 2015, 4, 11−15. (11) Li, W.; Müller, M. Directed self-assembly of block copolymers by chemical or topographical guiding patterns: optimizing molecular architecture, thin-film properties, and kinetics. Prog. Polym. Sci. 2016, 54−55, 47−75. K


Article

Macromolecules (33) Flory, P. J.; Crescenzi, V.; Mark, J. E. Configuration of the Poly(dimethylsiloxane) Chain. III. Correlation of Theory and Experiment. J. Am. Chem. Soc. 1964, 86, 146−152. (34) Zhang, L.-X.; Xia, A.-G.; Wang, X.-H; Xu, J.-M. Configurationdependent properties of poly(dimethylsiloxane) chain. Chin. J. Polym. Sci. 1999, 17, 347−354. (35) Tonelli, A. E.; Flory, P. J. The Configurational Statistics of Random Poly(lactic acid) Chains. I. Experimental Results. Macromolecules 1969, 2, 225−227. (36) Brant, D. A.; Tonelli, A. E.; Flory, P. J. The Configurational Statistics of Random Poly(lactic acid) Chains. II. Theory. Macromolecules 1969, 2, 228−235. (37) Theodorou, D. N.; Suter, U. W. Shape of Unperturbed Linear Polymers: Polypropylene. Macromolecules 1985, 18, 1206−1214. (38) Dodd, L. R.; Theodorou, D. N. Atomistic Monte Carlo Simulation and Continuum Mean Field Theory of the Structure and Equation of State Properties of Alkane and Polymer Melts. In Atomistic Modeling of Physical Properties; Monnerie, L. M., Suter, U. W., Eds.; Advances in Polymer Science No. 116; Springer-Verlag: Berlin, 1994; pp 249−281. (39) Moorthi, K.; Kamio, K.; Ramos, J.; Theodorou, D. N. Monte Carlo Simulation of Short Chain Branched Polyolefins: Structure and Properties. Macromolecules 2012, 45, 8453−8466. (40) Alexiadis, O.; Mavrantzas, V. G.; Khare, R.; Beckers, J.; Baljon, A. R. C. End-Bridging Monte Carlo Simulation of Bulk and Grafted Amorphous Polyethylene Above and Below the Glass Transition. Macromolecules 2008, 41, 987−996. (41) Peristeras, L. D.; Rissanou, A. N.; Economou, I. G.; Theodorou, D. N. Novel Monte Carlo Molecular Simulation Scheme Using identity-Altering Elementary Moves for the Calculation of Structure and Thermodynamic Properties of Polyolefin Blends. Macromolecules 2007, 40, 2904−2914. (42) Pandey, Y. N.; Doxastakis, M. Detailed atomistic Monte Carlo simulations of a polymer melt on a solid surface and around a nanoparticle. J. Chem. Phys. 2012, 136, 094901. (43) Baig, C.; Alexiadis, O.; Mavrantzas, V. G. Advanced Monte Carlo Algorithm for the Atomistic Simulation of Short- and LongChain Branched Polymers: Implementation for Model H-Shaped, A3AA3Multiarm (Pom-Pom), and Short-Chain Branched Polyethylene Melts. Macromolecules 2010, 43, 986−1002. (44) Mulder, T.; Harmandaris, V. A.; Lyulin, A. V.; van der Vegt, N. F. A.; Michels, M. A. J. Molecular simulation via connectivity-altering Monte Carlo and scale-jumping methods: Application to amorphous polystyrene. Macromol. Theory Simul. 2008, 17, 393−402. (45) Wick, C. D.; Theodorou, D. N. Connectivity-Altering Monte Carlo Simulations of the End Group Effects on Volumetric Properties for Poly(ethylene oxide). Macromolecules 2004, 37, 7026−7033. (46) Gestoso, P.; Nicol, E.; Doxastakis, M.; Theodorou, D. N. Atomistic Monte Carlo Simulation of Polybutadiene Isomers: cis-1,4Polybutadiene and 1,2-Polybutadiene. Macromolecules 2003, 36, 6925−6938. (47) Kamio, K.; Moorthi, K.; Theodorou, D. N. Coarse Grained End Bridging Monte Carlo Simulations of Poly(ethylene terephthalate) Melt. Macromolecules 2007, 40, 710−722. (48) Vogiatzis, G. G.; Theodorou, D. N. Local Segmental Dynamics and Stresses in Polystyrene−C60 Mixtures. Macromolecules 2014, 47, 387−404. (49) Subramanian, V.; Asirvatham, P. S.; Balakrishnan, R.; Ramasami, T. Molecular mechanics studies on polypropylene and polymethylmethacrylate polymers. Chem. Phys. Lett. 2001, 342, 603−609. (50) Honeycutt, J. D. A general simulation method for computing conformational properties of single polymer chains. Comput. Theor. Polym. Sci. 1998, 8, 1−8. (51) Metropolis, N.; Rosenbluth, A.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087−1092. (52) Sulatha, M. S.; Purushotham, S.; Natarajan, U. RMMC simulations of the chain properties of polyester homopolymers from

1,4-cyclohexanedimethanol and 2,2,4,4-tetramethyl- 1,3- cyclobutanediol. Polymer 2002, 43, 6295−6305. (53) Li, W.; Cui, S. W.; Wang, Q.; Yada, R. Y. Study of conformational properties of cereal β-glucans by computer modeling. Food Hydrocolloids 2012, 26, 377−382. (54) Pavel, D.; Yarovsky, I.; Shanks, R. Prediction of liquid crystalline properties of poly(1,4-phenylene sebacate-oxybenzoate) by Monte Carlo simulations. Polymer 2005, 46, 2003−2010. (55) Karyappa, R. B.; Natarajan, U. Molecular Simulations of the Conformational Properties of Atactic Poly(2-ethylbutyl methacrylate). J. Appl. Polym. Sci. 2012, 125, 1586−1591. (56) Karyappa, R. B.; Natarajan, U. Monte Carlo Simulations of Chain Dimensions and Conformational Properties of Various Poly(nalkyl methacrylates) in Solution. J. Macromol. Sci., Part B: Phys. 2008, 47, 1075−1086. (57) Patil, R. D.; Mark, J. E. Evaluations of forcefields for aromatic polysiloxanes, and some applications to poly(diphenylsiloxane). Comput. Theor. Polym. Sci. 2000, 10, 189−195. (58) Cail, J. I.; Stepto, R. F. T.; Taylor, D. J. R.; Jones, R. A.; Ward, I. M. Further computer simulation studies of the orientational behaviour of poly(ethylene terephthalate) chains. Phys. Chem. Chem. Phys. 2000, 2, 4361−4367. (59) Wu, Y.; Li, W.; Cui, W.; Eskin, N. A. M.; Goff, H. D. A molecular modeling approach to understand conformation-functionality relationships of galactomannans with different mannose/galactose ratios. Food Hydrocolloids 2012, 26, 359−364. (60) Destrée, M.; Lyulin, A.; Ryckaert, J.-P. Monte Carlo Prediction of the Structure Factor of Polyethylene. Macromolecules 1996, 29, 1721−1727. (61) Sariban, A.; Brickmann, J.; van Ruiten, J.; Meier, R. J. Simulating Θ-Conditions in Polymers by Modifying the Long-Range Interactions: A Monte Carlo Approach to Polymethylene. Macromolecules 1992, 25, 5950−5956. (62) Scienomics, MAPS platform: version 3.4.1, France, 2014. (63) Mavrantzas, V. G.; Theodorou, D. N. Atomistic Simulation of Polymer Melt Elasticity: Calculation of the Free Energy of an Oriented Polymer Melt. Macromolecules 1998, 31, 6310−6332. (64) Lal, M. ‘Monte Carlo’ computer simulation of chain molecules. I. Mol. Phys. 1969, 17, 57−64. (65) Madras, N.; Sokal, A. D. The pivot algorithm: A highly efficient Monte Carlo method for the self-avoiding walk. J. Stat. Phys. 1988, 50, 109−186. (66) Fröhlich, M. G.; Sewell, T. D. Pivot Algorithm and Push-off Method for Efficient System Generation of All-Atom Polymer Melts: Application to Hydroxyl-Terminated Polybutadiene. Macromol. Theory Simul. 2013, 22, 344−353. (67) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (68) Binder, K.; Baschnagel, J.; Müller, M.; Paul, W.; Rampf, F. Simulation of phase transitions of single polymer chains: Recent advances. Macromol. Symp. 2006, 237, 128−138. (69) Theodorou, D. N. Variable Connectivity Monte Carlo Algorithms for the Atomistic Simulation of Long-Chain Polymer Systems. In Bridging Time Scales: Molecular Simulations for the Next Decade; Lecture Notes in Physics 605; Nielaba, P., Mareschal, M., Ciccotti, G., Eds.; Springer-Verlag: Berlin, 2002; pp 67−127. (70) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (71) Baschnagel, J.; Qin, K.; Paul, W.; Binder, K. Monte Carlo Simulation of Models for Single Polyethylene Coils. Macromolecules 1992, 25, 3117−3124. (72) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (73) Martin, M. G.; Siepmann, J. I. Novel Configurational-Bias Monte Carlo Method for Branched Molecules. Transferable Potentials for Phase Equilibria. 2. United-Atom Description of Branched Alkanes. J. Phys. Chem. B 1999, 103, 4508−4517. L


Article

Macromolecules (74) Toxvaerd, S. Equation of state of alkanes II. J. Chem. Phys. 1997, 107, 5197−5204. (75) Karayiannis, N. C.; Mavrantzas, V. G.; Theodorou, D. N. A Novel Monte Carlo Scheme for the Rapid Equilibration of Atomistic Model Polymer Systems of Precisely Defined Molecular Architecture. Phys. Rev. Lett. 2002, 88, 105503. (76) Martin, M. G.; Siepmann, J. I. Novel Configurational-Bias Monte Carlo Method for Branched Molecules. Transferable Potentials for Phase Equilibria. 2. United-Atom Description of Branched Alkanes. J. Phys. Chem. B 1999, 103, 4508−4517. (77) Karayiannis, N. C.; Giannousaki, A. E.; Mavrantzas, V. G.; Theodorou, D. N. Atomistic Monte Carlo simulation of strictly monodisperse long polyethylene melts through a generalized chain bridging algorithm. J. Chem. Phys. 2002, 117, 5465−5479. (78) Foteinopoulou, K.; Karayiannis, N. C.; Laso, M.; Kröger, M. Structure, Dimensions, and Entanglement Statistics of Long Linear Polyethylene Chains. J. Phys. Chem. B 2009, 113, 442−455. (79) Logotheti, G. E.; Theodorou, D. N. Segmental and Chain Dynamics of Isotactic Polypropylene Melts. Macromolecules 2007, 40, 2235−2245. (80) Romanos, N. A.; Theodorou, D. N. Crystallization and Melting Simulations of oligomeric α1 Isotactic Polypropylene. Macromolecules 2010, 43, 5455−5469. (81) Zirkel, A.; Urban, V.; Richter, D.; Fetters, L. J.; Huang, J. S.; Kampmann, R.; Hadjichristidis, N. Small-Angle Neutron Scattering Evaluation of the Temperature Dependence of Atactic Polypropylene and Poly(1-butene) Chain Dimensions in the Melt. Macromolecules 1992, 25, 6148−6155. (82) Harmandaris, V. A.; Adhikari, N. P.; van der Vegt, N. F. A.; Kremer, K. Hierarchical Modeling of Polystyrene: From Atomistic to Coarse-Grained Simulations. Macromolecules 2006, 39, 6708−6719. (83) Mulder, T.; Harmandaris, V. A.; Lyulin, A.; van der Vegt, N. F. A.; Kremer, K.; Michels, M. A. J. Structural Properties of Atactic Polystyrene of Different Thermal History Obtained from a Multiscale Simulation. Macromolecules 2009, 42, 384−391. (84) Mark, J.; Ngai, K.; Graessley, W.; Mandelkern, L.; Samulski, E. Physical Properties of Polymers; Cambridge University Press: Cambridge, 2003. (85) Annis, B. K.; Kim, M.-H.; Wignall, G. D.; Borodin, O.; Smith, G. D. A Study of the Influence of LiI on the Chain Conformations of Poly(ethylene oxide) in the Melt by Small-Angle Neutron Scattering and Molecular Dynamics Simulations. Macromolecules 2000, 33, 7544−7548. (86) Frischknecht, A. L.; Curro, J. G. Improved United Atom Force Field for Poly(dimethylsiloxane). Macromolecules 2003, 36, 2122− 2129. (87) Habenschuss, A.; Tsige, M.; Curro, J. G.; Grest, G. S.; Nath, S. K. Structure of Poly(dialkylsiloxane) Melts: Comparison of WideAngle X-ray Scattering, Molecular Dynamics Simulations, and Integral Equation Theory. Macromolecules 2007, 40, 7036−7043. (88) Sides, S. W.; Curro, J.; Grest, G. S.; Stevens, M. J.; Soddemann, T.; Habenschuss, A.; Londono, J. D. Structure of Poly(dimethylsiloxane) melts: Theory, Simulation, and Experiment. Macromolecules 2002, 35, 6455−6465. (89) Nath, S. K.; Frischknecht, A. L.; Curro, J. G.; McCoy, J. D. Density Functional Theory and Molecular Dynamics Simulation of Poly(dimethylsiloxane) Melts near Silica Surfaces. Macromolecules 2005, 38, 8562−8573. (90) Maple, J. R.; Hwang, M.-J.; Stockfisch, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. Derivation of Class II Force Fields. I. Methodology and Quantum Force Field for the Alkyl Functional Group and Alkane Molecules. J. Comput. Chem. 1994, 15, 162−182. (91) Fetters, L. J.; Lohse, D. J.; Colby, R. H. Chain dimensions and entanglement spacings. In Physical Properties of Polymers Handbook; Mark, J. E., Ed.; Springer: New York, 2007; pp 447−454. (92) Kamath, G.; Robinson, J.; Potoff, J. J. Application of TraPPE-UA forcefield for determination of vapor-liquid equilibria of carboxylate esters. Fluid Phase Equilib. 2006, 240, 46−55.

(93) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 6. United-Atom Description for Ethers, Glycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108, 17596−17605. (94) Maerzke, K. A.; Schultz, N. E.; Ross, R. B.; Siepmann, J. I. TraPPE-UA Force Field for Acrylates and Monte Carlo Simulations for Their Mixtures with Alkanes and Alcohols. J. Phys. Chem. B 2009, 113, 6415−6425. (95) Clifford, S.; Bolton, K.; Ramjugernath, D. Monte Carlo Simulation of Carboxylic Acid Phase Equilibria. J. Phys. Chem. B 2006, 110, 21938−21943. (96) Kamath, G.; Cao, F.; Potoff, J. J. An Improved Force Field for the Prediction of the Vapor-Liquid Equilibria for Carboxylic Acids. J. Phys. Chem. B 2004, 108, 14130−14136. (97) Boothroyd, A. T.; Rennie, A. R.; Boothroyd, C. B. Direct Measurement of the Temperature Dependence of the Unperturbed Dimensions of a Polymer. Europhys. Lett. 1991, 15, 715−719. (98) Boothroyd, A. T.; Rennie, A. R.; Wignall, G. D. Temperature coefficients for the chain dimensions of polystyrene and polymethylmethacrylate. J. Chem. Phys. 1993, 99, 9135−9144. (99) Jorgensen, W. L.; Tirado-Rives, J. The OPLS Potential Functions for Proteins. Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (100) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (101) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. L. OPLS All-Atom Force Field for Carbohydrates. J. Comput. Chem. 1997, 18, 1955−1970. (102) Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III DREIDING: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897−8909. (103) Sun, H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applicationss–Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338−7364. (104) Sun, H. Prediction of fluid densities using automatically derived VDW parameters. Fluid Phase Equilib. 2004, 217, 59−76. (105) Blomqvist, J. RIS Metropolis Monte Carlo studies of poly(Llactic), poly(L,D-lactic) and polyglycolic acids. Polymer 2001, 42, 3515−3521. (106) Blomqvist, J.; Pietilä, L.-O.; Mannfors, B. RIS Metropolis Monte Carlo studies of some aliphatic main chain and side group polyesters. Polymer 2001, 42, 109−116.

M


General Methodology for Estimating the Stiffness of Polymer Chains

Recommend Documents