Equation of State Based Slip Spring Model for Entangled Polymer

Article pubs.acs.org/Macromolecules

Equation of State Based Slip Spring Model for Entangled Polymer Dynamics Georgios G. Vogiatzis, Grigorios Megariotis, and Doros N. Theodorou* School of Chemical Engineering, National Technical University of Athens, 9 Heroon Polytechniou Street, Zografou Campus, GR-15780 Athens, Greece ABSTRACT: A mesoscopic, mixed particle- and field-based Brownian dynamics methodology for the simulation of entangled polymer melts has been developed. Polymeric beads consist of several Kuhn segments, and their motion is dictated by the Helmholtz energy of the sample, which is a sum of the entropic elasticity of chain strands between beads, slip springs, and nonbonded interactions. Following earlier works in the field [Phys. Rev. Lett. 2012, 109, 148302], the entanglement effect is introduced by the slip springs, which are springs connecting either nonsuccessive beads on the same chain or beads on different polymer chains. The terminal positions of slip springs are altered during the simulation through a kinetic Monte Carlo hopping scheme, with rate-controlled creation/destruction processes for the slip springs at chain ends. The rate constants are consistent with the free energy function employed and satisfy microscopic reversibility at equilibrium. The free energy of nonbonded interactions is derived from an appropriate equation of state, and it is computed as a functional of the local density by passing an orthogonal grid through the simulation box; accounting for it is necessary for reproducing the correct compressibility of the polymeric material. Parameters invoked by the mesoscopic model are derived from experimental volumetric and viscosity data or from atomistic molecular dynamics simulations, establishing a “bottom-up” predictive framework for conducting slip spring simulations of polymeric systems of specific chemistry. Initial configurations for the mesoscopic simulations are obtained by further coarse-graining of well-equilibrated structures represented at a greater level of detail. The mesoscopic simulation methodology is implemented for the case of cis-1,4-polyisoprene, whose structure, dynamics, thermodynamics, and linear rheology in the melt state are quantitatively predicted and validated without a posteriori fitting the results to experimental measurements.

■

INTRODUCTION One of the fundamental concepts in the molecular description of structure−property relations of polymer melts is chain entanglement. When macromolecules interpentrate, the term entanglements intends to describe the topological interactions resulting from the uncrossability of chains. The fact that two polymer chains cannot go through each other in the course of their motion changes their dynamical behavior dramatically, without altering their equilibrium properties. Anogiannakis et al.1 have examined microscopically at what level topological constraints can be described as a collective entanglement effect, as in tube model theories, or as certain pairwise uncrossability interactions, as in slip link models. They employed a novel methodology, which analyzes entanglement constraints into a complete set of pairwise interactions (links), characterized by a spectrum of confinement strengths. As a measure of the entanglement strength, these authors used the fraction of time for which the links were active. The confinement was found to be mainly imposed by the strongest links. The weak, trapped, uncrossability interactions cannot contribute to the lowfrequency modulus of an elastomer or the plateau modulus of a melt. In tube model theories,2−5 it is postulated that the entanglements generate a confining mean-field potential, © 2017 American Chemical Society

which restricts the lateral monomer motion to a tube-like region surrounding each chain. In polymer melts the confinement is not permanent but leads to a one-dimensional diffusion of the chain along its tube, called reptation.6 An alternative, discrete, localized version of the tube constraint is utilized in models employing slip links.7−13 The tube is replaced by a set of slip links along the chain, which restrict lateral motion but permit chain sliding through them. The real chain is represented by its primitive path, which is a series of strands of average molar mass Me connecting the links. Hua and Schieber11,12 considered that the molecular details on the monomer or Kuhn-length level are smeared out in the slip link model, while the segmental network of generic polymers is directly modeled, by introducing links between chains, which constrain the motion of segments of each chain into a tubular region. The motion of segments is updated stochastically, and the positions of slip links can either be fixed in space or mobile. When either of the constrained segments slithers out of a slip link constraint, they are considered to be disentangled, and the slip link is destroyed. Conversely, the end Received: August 5, 2016 Revised: March 6, 2017 Published: March 23, 2017 3004

DOI: 10.1021/acs.macromol.6b01705 Macromolecules 2017, 50, 3004−3029

Article

Macromolecules

by incorporating the tension force from chain segments and an osmotic force caused by density fluctuations. The entanglement points, modeled as slip links, can also fluctuate spatially (or three-dimensionally). The creation and annihilation of entanglements are controlled by the number of monomers at chain ends. The longest relaxation time and the self-diffusion coefficient scaling, as predicted from the model, were found in good agreement with experimental results. Later on, the primitive chain network model was extended to study the relationship between entanglement length and plateau modulus.20−23 It was also extended to study star and branched polymers,24 nonlinear rheology,25−27 phase separation in polymer blends,28 block copolymers,29 and the dynamics of confined polymers.30 Recently, Masubuchi has compiled an excellent review of simulations of entangled polymers.31 Coarse graining from an atomistic molecular dynamics (MD) simulation to a mesoscopic one results in soft repulsive potentials between the coarse-grained blobs, causing unphysical bond crossings influencing the dynamics of the system. The Twentanglement uncrossability scheme proposed by Padding and Briels32,33 was introduced in order to prevent these bond crossings. Their method is based on considering the bonds as elastic bands between the bonded blobs. As soon as two of these elastic bands make contact, an “entanglement” is created which prevents the elastic bands from crossing. In their method, entanglements were defined as the objects which prevent the crossing of chains. Only a few of them were expected to be entanglements in the usual sense of long-lasting topological constraints, slowing down the chain movement. Simulations of a polyethylene melt32 showed that the rheology can be described correctly, deducing the blob−blob interactions and the friction parameters from MD simulations. Chappa et al.34 proposed a model in which the topological effect of noncrossability of long flexible macromolecules is effectively taken into account by slip springs, which are local, pairwise, translationally and rotatationally invariant interactions between polymer beads that do not affect the equilibrium properties of the melt. The slip springs introduce an effective attractive potential between the beads; in order to eliminate this effect, the authors introduced an analytically derived repulsive compensating potential. The conformations of polymers and slip springs are updated by a hybrid Monte Carlo (MC) scheme. At every step, either the positions of the beads are evolved via a short dissipative particle dynamics (DPD) run or the configuration of the slip springs is modified by MC moves, involving discrete jumps of the slip springs along the chain contour and creation or deletion at the chain ends. The number of slip springs can vary during the simulation, obeying a prescribed chemical potential. That model can correctly describe many aspects of the dynamical and rheological behavior of entangled polymer liquids in a computationally efficient manner, since everything is cast into only pairwiseadditive interactions between beads. The mean-square displacement of the beads evolving according to this model was found to be in favorable agreement with the tube model predictions. Moreover, the model exhibited realistic shear thinning, deformation of conformations, and a decrease of the number of entanglements at high shear rates. Müller and Daoulas were the first to introduce the hopping of slip springs by means of discrete MC moves.35 These authors thoroughly investigated the ability of MC algorithms to describe the single-chain dynamics in a dense homogeneous melt and a lamellar phase of a symmetric diblock copolymer.

of one segment can hop toward another segment and create another new entanglement or slip link. The governing equations in the slip link model can be split into two parts: the chain motion inside its tube is governed by Langevin equations and the tube motion is governed by deterministic convection and stochastic constraint release processes. Based on the tube model,4 it is assumed that the motion of the primitive path makes the primary contribution to the rheological properties of entangled polymer melts. Therefore, from the microscopic information given by the slip link model, these authors could precisely access the longest polymer chain relaxation time. Moreover, by employing an elegant formulation, the macroscopic properties of polymer melts, i.e., stress and dielectric relaxation, can be extracted from the ordering, spatial location, and aging of the entanglements or slip links in the simulations. Later, Schieber and co-workers studied the fluctuation effect on the chain entanglement and viscosity using a mean-field model.14 Shanbhag et al.15 developed a dual slip link model with chain-end fluctuations for entangled star polymers, which explained the observed deviations from the “dynamic dilution” equation in dielectric and stress relaxation data. Doi and Takimoto16 adopted the dual slip link model to study the nonlinear rheology of linear and star polymers with arbitrary molecular weight distribution. Likhtman17 has shown that the standard tube model cannot be applied to neutron spin-echo measurements because the statistics of a one-dimensional chain in a three-dimensional random-walk tube become wrong on the length scale of the tube diameter. He then introduced a new single-chain dynamic slip link model to describe the experimental results for neutron spin echo, linear rheology, and diffusion of monodisperse polymer melts. All the parameters in this model were obtained from one type of experiment and were applied to predict other experimental results. The model was formulated in terms of stochastic differential equations, suitable for Brownian dynamics (BD) simulations. The results were characterized by some systematic discrepancies, suggesting possible inadequacy of the Gaussian chain model for some of the polymers considered and possible inadequacy of the time−temperature superposition (TTS) principle. Del Biondo et al.18 extended the slip link model of Likhtman17 in order to study an inhomogeneous system. These authors studied the dependence of the relaxation modulus on the slip link density and stiffness and explored the nonlinear rheological properties of the model for a set of its parameters. The crucial part of their work was the introduction of excluded volume interactions in a mean field manner in order to describe inhomogeneous systems and application of this description to a simple nanocomposite model. They concentrated on an idealized situation where the fillers were well dispersed, with a simple hard-core interaction between the fillers and the polymer matrix. However, addressing real situations where the fillers are poorly dispersed and partially aggregated was clearly possible within their framework. Masubuchi et al.19 performed several multichain simulations for entangled polymer melts by utilizing slip links to model the entanglements. These authors proposed a primitive chain network model, in which the polymer chain is coarse-grained into a set of segments (strands) going through entanglement points. Different segments are coupled together through the force balance at the entanglement node. The Langevin equation is applied to update the positions of these entanglement nodes 3005


Article

Macromolecules

Laradji et al.46 originally introduced the concept of conducting particle-based simulations with interactions described via collective variables and also two ways of calculating these collective variables: grid-based and continuum weightingfunction-based. Along the same lines, Soga et al.47 investigated the structure of an end-grafted polymer layer immersed in poor solvent through MC simulations based on Edwards’s Hamiltonian, incorporating a third-order density functional. Pagonabarraga and Frenkel48 carried out dissipative particle dynamics (DPD) simulations of nonideal fluids. In their work the conservative forces needed for the DPD scheme were derived from a free energy density obtained from an equation of state (EOS) (e.g., the van der Waals EOS). Later on, Daoulas and Müller49 have also employed the density functional representation of nonbonded interactions in order to study the single chain structure and intermolecular correlations in polymer melts and fluctuation effects on the order−disorder transition of symmetric diblock copolymers. In the present work we introduce a grid-based density-field-based scheme for dealing with nonbonded interactions in coarsegrained simulations which takes into account a specific equation of state. The purpose of this work is twofold. First, our main goal is to develop a consistent and rigorous methodology capable of simulating polymeric melts over realistic (i.e., 10−2 s) time scales. In our previous works,50,51 we have developed a methodology in order to generate and equilibrate (nanocomposite) polymer melts at large length scales (on the order of 100 nm). This coarse-grained representation is based on the idea that the polymer chains can be described as random flights at length scales larger than that defined by the Kuhn length of the polymer. Now, we develop a methodology to track the dynamics of the system at a coarse-grained level by invoking a free energy functional motivated by an EOS capable of describing real polymeric materials, where both conformational and nonbonded contributions are taken into account. The EOS-based nonbonded free energy allows us to simulate a compressible model which, in principle, can deal with equation of state effects, phase transitions (e.g., cavitation), and interfaces. It is computed as a functional of the local density by passing an orthogonal grid through the simulation box. The voxels of the grid employed are large enough that their local thermodynamic properties can be described well by the chosen macroscopic EOS. The validity of this assumption is tested. At our level of description, i.e., beads consisting of 50 chemical monomers, forces among beads are not pairwise-additive, and many-body effects are likely to be important.52−54 Pairwise additive effective potentials between beads55,56 may not work well under conditions of high stress, where phase separation phenomena such as cavitation take place, or near interfaces, where long-range attractions are important in shaping structure. Thus, resorting to drastically coarse-grained models based on collective variables is a promising route.49 Introducing a nonbonded compensating potential like Chappa et al.34 eliminates the problem of chain conformations, caused by the artificial attraction of the slip springs. Here, we refrain from explicitly including a term of this kind, in order to check whether the use of our nonbonded energy functional can provide an alternative way for avoiding unphysical agglomerations. Once the free energy is known, BD simulations driven by the free energy functional can be used in order to obtain thermodynamic averages and correlation functions. When

Three different types of single-chain dynamics were considered: local, unconstrained dynamics, slithering-snake dynamics, and slip link dynamics. Ramı ́rez-Hernań dez et al.36,37 replaced the discrete hopping of a slip spring along the chain contour by a one-dimensional continuous Langevin equation. In their method, slip springs consist of two rings that slip along different chain contours and are connected by a harmonic spring. The rings move in straight lines between beads belonging to different chains, scanning the whole contour of the chains in a continuous way. Comparison with experimental results was made possible by rescaling the frequency and the modulus in order to match the intersection point of storage and loss moduli curves of polystyrene melts. More recently, Ramı ́rez-Hernań dez et al.38 presented a more general formalism in order to qualitatively capture the linear rheology of pure homopolymers and their blends as well as the nonlinear rheology of highly entangled polymers and the dynamics of diblock copolymers. The number of slip springs in their approach remained constant throughout the simulation, albeit their connectivity changes. These authors have later presented a theoretically informed entangled polymer simulation approach wherein the total number of slip springs is not preserved, but instead, it is controlled through a chemical potential that determines the average molecular weight between entanglements.39 The idea of slip springs was in parallel used by Uneyama and Masubuchi,40 who proposed a multichain slip spring model inspired by the single chain slip spring model of Likhtman.17 Differently from the primitive chain network model of the same authors, they defined the total free energy for the new model and employed a time evolution equation and stochastic processes for describing its dynamical evolution. All dynamic ingredients satisfy the detailed balance condition and are thus capable of reproducing the thermal equilibrium which is characterized by the free energy. The number of slip springs varies. Later, Langeloth et al.41,42 presented a simplified version of the slip spring model of Uneyama and Masubuchi,40 where the number of slip springs remains constant throughout the simulation. Working on a different problem than melt rheology, Terzis et al.43−45 have invoked the microscopic description of entanglements and the associated processes envisioned in slip link models, in order to generate entanglement network specimens of interfacial polymeric systems and study their deformation to fracture. The specimens were created by sampling the configurational distribution functions derived from a selfconsistent field (SCF) lattice model. The specimens generated were not in detailed mechanical equilibrium. To this end, these authors developed a method for relaxing the network with respect to its density distribution and thereby imposing the condition of mechanical equilibrium, without changing the network topology.45 The free energy function of the network was minimized with respect to the coordinates of all entanglement points and chain ends. Contributions to the free energy included (a) the elastic energy due to stretching of the chain strands and (b) the free energy due to the repulsive and attractive (cohesive) interactions between segments. The latter was calculated by superimposing a simple cubic grid on the network and taking into account contributions between cells and within each cell. A density functional was used for nonbonded interactions depending on a segment density field derived from coarse-grained particle positions tracked in a simulation. 3006


Article

Macromolecules macromolecules interpenetrate, the term entanglement intends to describe the interactions resulting from the uncrossability of chains. At this level of description, we introduce the entanglements as slip springs connecting beads belonging to different chains. In the course of a simulation, the topology of the entanglement network changes through the introduction of elementary kinetic Monte Carlo (kMC) events governed by rate expressions which are based on the reptation picture of polymer dynamics and the free energy defined. Our second objective concerns parametrizing this methodology in a bottom-up approach, avoiding the introduction of nonmeaningful parameters and ad-hoc fitting of the results. All ingredients of our model are based on either more detailed (atomistic) simulations or experimental findings. For every parameter we introduce, we thoroughly explain its physical meaning and provide rigorous guidelines for estimating its exact value or its range of variation. Our work is different from the extensive current literature on the subject in a number of ways. All observables stem out of an explicit Helmholtz energy functional. The stress tensor of the model, a necessary ingredient for rheology calculations, is rigorously derived from the free energy functional including all contributions (bonded and nonbonded). Its special characteristics are presented and discussed. In contrast to previous relevant studies, a rigorous kinetic MC scheme (with rate constants defined in terms of the free energy) is used for tracking slip spring hopping dynamics in the system. A hybrid integration scheme allows for rate-controlled discrete slip spring processes to be tracked in parallel with the integration of equations of motion of the polymeric beads. A new densityfield-based scheme for dealing with nonbonded interactions in coarse-grained simulations, reproducing a specific equation of state, is presented, and its numerical application is thoroughly reviewed. All necessary links to more detailed levels of simulation and experimental observables are established, without introducing a posteriori tuning of the parameters. Finally, extensive quantitative validation against existing experimental findings is performed.

Figure 1. Network representation of a polymer melt or rubber. End beads, internal beads, and cross-link beads (if present) are shown as red, blue, and black spheres, respectively. Entropic springs along the chain contour are shown in blue, while slip springs are shown in red.

stochastically by free chain ends in the polymer network. Along with the BD simulation of the bead motion in the periodic simulation box, a parallel simulation of the evolution of the slip springs present is undertaken. The rates used for the kMC procedure are described in detail in the following paragraphs. The initial contour length between consecutive entanglement points (slip spring anchoring beads) on the same chain is commensurate with the entanglement molecular weight, Me, of the simulated polymer. With the above mesoscopic model, the relative importance of reptation, constraint release (CR), and contour-length fluctuation (CLF) mechanisms depends on the specific melt of interest. As expected, in a monodisperse sample of long molecules, reptation plays a dominant role. In contrast, in a bidisperse sample composed of short and long macromolecules, CR and CLF may dominate the relaxation process.57−60 Model Free Energy. We postulate that the entangled melt specimen of given spatial extent, defined by the edge vectors of our periodic simulation box (Lx, Ly, Lz), under temperature T, is governed by a Helmholtz energy function, A, which has a direct dependence on the set of local densities {ρ(r)}, the temperature, T, and the separation vectors between pairs of connected polymer beads, {rij ≡ rj − ri}:

■

MODEL AND METHOD Polymer Description. Our melt consists of linear chains, represented by specific points or beads (i.e., internal nodal points, end points, or permanent cross-link points) along their contour, connected by entropy springs. Each coarse-grained bead represents several Kuhn segments of the polymer under consideration. Our construction results in a set of nodes for each chain, where each node i has a specific contour position along the chain, positional vector in three-dimensional space, ri, and pairing with other nodes, as shown in Figure 1. The piece of chain between two nodes (blue spring in Figure 1) is referred to as a strand. This paper will be concerned with melts of linear chains, so permanent cross-link points will not be considered further. We introduce the effect of entanglement by dispersing slip springs,34 which are designed to bring about reptational motion of chains along their contours, as envisioned in theories and simulations of dynamics in entangled polymer melts and as observed by topological analysis of molecular configurations evolving in the course of MD simulations.1 In more detail, a slip spring connects two internal nodal points or chain ends on two different polymer chains and is stochastically destroyed when one of the nodes it connects is a chain end. To compensate for slip spring destruction, new slip springs are created

A({rij}, {ρ(r)}, T ) = Ab({rij}, T ) + A nb({ρ(r)}, T )

(1)

The first term on the right-hand side of eq 1, Ab, is the contribution of bonded interactions, whereas the second one, Anb, is the contribution of the nonbonded interactions to the Helmholtz energy. It should be noted that we employ the full free energy (including ideal gas contribution) and not the excess free energy in our formulation. In the discussion of the nonbonded contribution to the Helmholtz energy, Anb, we will elaborate further on our choice of employing the full Helmholtz energy, rather than the excess one (with reference to an ideal gas of beads). We start by considering the bonded contribution to the Helmholtz energy, which can be written as a sum over all bonded pairs, (i, j), where i is connected to j with either intramolecular springs or slip springs: 3007


Article

Macromolecules Ab({rij}, T ) = Ab({rij}, T ) =

translational entropy associated with allowing the beads which reside in one voxel to visit other voxels. However, we have validated our results against those obtained by employing the excess Helhmoltz energy. They are practically identical because, in our implementation, the cells of the discretization grid are so large that they behave almost as macroscopic domains. The translational entropy associated with exchange of a bead between voxels is negligible in comparison with the entropy associated with roaming of the bead within a voxel. In this work we invoke the Sanchez−Lacombe equation of state,62 which gives for the Helmholtz energy:

∑ A pair (rij , T ) (i , j)

(2)

The sum runs over all pairs, where each pair is thought to interact via a distance-dependent Helmholtz energy, Apair(rij,T). The elastic force depends on the coordinates of the nodes, but it does not depend on the local network density. The elastic force between connected beads arises due to the retractive force acting to resist the stretching of a strand. This force originates in the decrease in entropy of a stretched polymer strand. In this approximation, the force on bead i, due to its connection to bead j is Fijb

= −∇rji A pair (rij , T )

ASL (ρ , T ) = − nrkBT *ρ ̃ + ⎤ ⎡⎛ 1 ⎞ 1 nrkBT ⎢⎜ − 1⎟ × ln(1 − ρ ̃) + ln(ρ ̃)⎥ − nkBT ln(w) r ⎠ ⎦ ⎣⎝ ρ ̃

(3)

The Gaussian approximation to Apair (rij) can be used for most extensions.61 The conformational entropy of strands is taken into account via a simple harmonic expression: intra A pair (rij ,

rij 2 3 T ) = kBT 2 NijbK 2

(8)

where ρ̃ = ρ/ρ* is the reduced density with ρ being the mass density of the melt and ρ* the close-packed mass density. The temperature of the melt is denoted by T, and T* = ϵ/kB corresponds to an equivalent temperature defined in terms of ϵ which is the nonbonded, mer−mer interaction energy (kB being Boltzmann’s constant). The n chains of the melt, whose molecular weight per chain is M, are thought to be composed of r Sanchez Lacombe segments each (“r-mers”). Finally, w is a parameter quantifying the number of configurations available to the system. The last term in eq 8 does not depend on the density of the system, being an ideal gas contribution (please see discussion in Appendix B). All Sanchez−Lacombe parameters are presented in Table B1 of Appendix B. The Helmholtz energy density, avol(ρ,T), can be calculated as

(4)

where Nij is the number of Kuhn segments assigned to the strand, bK the Kuhn length of the polymer, and kB the Boltzmann constant. The summation is carried over all pairs (i, j) which lie along the contour of the chains. The Helmholtz energy of slip springs, which are included to account for the entanglement effect, is described by the equation sl A pair (rij , T ) =

rij 2 3 kBT 2 2 lss

(5)

where lss is an adjustable parameter (i.e., root-mean-square equilibrium slip spring length) which should be larger than the Kuhn length, bK, and smaller than or equal to the tube diameter of the polymer under consideration, app. In order to account for nonbonded (excluded volume and van der Waals attractive) interactions in the network representation, we introduce a Helmholtz energy functional: A nb({ρ(r)}, T ) =

∫box d3r avol[ρ(r), T ]

a vol(ρ , T ) =

= ρamass(ρ , T )

rRT *ρ ̃ rRT + M M ⎡ 1 ⎤ RT 1 × ⎢( − 1) ln(1 − ρ ̃) + ln(ρ ̃)⎥ − ln(w) ⎣ ρ̃ ⎦ r M (10)

amass(ρ , T ) = −

(6)

and by using r = (Mp*)/(ρ*RT*), the above expression becomes amass(ρ , T ) = −

Ncells

acc ∑ V cell, ka vol (ρcell, k , T ) k=1

(9)

where the amass(ρ,T) denotes the Helmholtz energy per unit mass of the system and NA is Avogadro’s number. On the basis of eq 8, we can calculate the above:

where avol is a free energy density (free energy per unit volume) and ρ(r) is the local mass density at position r. Expressions for avol(ρ,T) may be extracted from equations of state. Local density, ρ(r), will be resolved only at the level of entire cells, defined by passing an orthogonal grid through the simulation domain. Thus, eq 6 will be approximated by a discrete sum over all cells of the orthogonal grid: A nb({ρ(r)}, T ) =

A (ρ , T ) A (ρ , T ) A (ρ , T ) = = ρNA nM V nM /(ρNA )

p* 2

ρ*

ρ+

⎞ p*T̃ ⎡⎛ ρ* ⎢⎜ − 1⎟ ρ * ⎣⎝ ρ ⎠

⎛ ρ⎞ ρ*RT * ⎛ ρ ⎞⎤ × ln⎜1 − ln⎜ ⎟+ ⎟⎥ ρ* ⎠ Mp* ⎝ ⎝ ρ*w ⎠⎦

(7)

where Vacc cell,k is the accessible volume of cell k. The cell density, ρcell,k, must be defined based on the beads in and around cell k. The spatial discretization scheme employed for nonbonded interactions is thoroughly presented in Appendix A. In deriving nonbonded forces (see eq 13), it would have been more correct to use the excess Helmholtz energy density (relative to an ideal gas of beads), rather than the total Helmholtz energy density and to rely on the particle-based simulation for effectively including the contribution of translational entropy. Our formulation uses the total, rather the excess, Helmholtz energy density, thereby missing the

(11)

The Helmholtz energy density (free energy per unit volume), avol(ρ,T), is a vol(ρ , T ) = ρamass(ρ , T ) = −p*ρ 2̃ + ⎡⎛ 1 ⎤ ⎞ ρ*RT * ρ*RT * ln(ρ ̃) − ln(w)⎥ p*T̃ρ ̃⎢⎜ − 1⎟ × ln(1 − ρ ̃) + Mp* Mp* ⎠ ⎣⎝ ρ ̃ ⎦

(12)

where everything is cast in terms of the reduced variables T*, p*, ρ*, and the molecular weight of a chain, M. Upon 3008


Article

Macromolecules

parameter for the Gibbs energy minimization. In that case, the volume is expressed as V = V9(1 + ϵxx + ϵyy + ϵzz). The Gibbs energy function, eq 15, is, of course, valid only for small deformations away from the reference state. In our case, we let the system relax under atmospheric hydrostatic pressure, starting from an equilibrated configuration at the average density for the temperature under consideration. Slip springs have not been introduced yet at this point. Slip springs can either be placed randomly in the initial configuration, or they can be allowed to be created, following the fluctuating kMC scheme described in Appendix E. In the former way, the number of slip springs is chosen so as to be consistent with the molar mass between entanglements, Me, of the polymer under consideration. Every possible pair of beads in the melt can be coupled by a slip spring, given that their separation distance is less than lss (eq 5). Otherwise, the initially slip spring free configuration is used as an input, and slip springs are generated on the go by utilizing the fluctuating kMC scheme (Appendix E) with a hopping pre-exponential factor, νhop, at the upper limit of the allowed range (Appendix D). In this way, slip springs are introduced in an unbiased and rigorous way while the melt is in parallel equilibrated. The procedure is interrupted at the point where the number of entanglements (slip springs) reaches the desired value, and the frequency factor νhop is fixed at the value that ensures conservation of the number of slip springs. Once the slip springs have been introduced, the Gibbs energy of the system (eq 15) can be minimized again in their presence, at fixed topology of the entanglement network, in order to ensure the proper simulation box dimensions. However, slip springs contribute little (on the order of 1%) to the stress tensor of the system. Brownian Dynamics. When simulating a system of coarsegrained particles, some degrees of freedom are treated explicitly, whereas others are represented only by their stochastic influence on the former ones. In our model the effect of the surrounding melt on the motion of the coarsegrained beads is mimicked by introducing a stochastic force plus a frictional force into the equations of motion of the beads. When the stochastic force contains no correlations in space or time, one obtains the simplest form of stochastic dynamics, called Brownian dynamics (BD).65−67 The theory of Brownian motion was developed to describe the dynamic behavior of particles whose mass and size are much larger than those of the host medium particles. In this case the position Langevin equation becomes

integration over the domain of a system with given mass, the last term within the brackets of eq 12, which is strictly proportional to density, yields a density-independent contribution. Thus, it does not affect the equations of motion or thermodynamic properties. All needed parameters (i.e., T*, p*, and ρ*) can be obtained from experimental studies.63 The force due to nonbonded interactions on a bead i is given by Finb = −∇ri A nb({ρcell }, T ) = − =−

∑ k ∈ cells

acc V cell, k

⎤ ∂ ⎡⎢ acc ⎥ ∑ V cell, ka vol (ρcell, k , T ) ∂ri ⎢⎣ k ∈ cells ⎦⎥

∂a vol(ρ , T ) ∂ρ

∂ρcell, k ρ = ρcell, k

∂ri (13)

with the derivative of avol(ρ,T) with respect to density being ⎛ ∂a vol(ρ , T ) p*T ⎡ 2p* ρ ⎞⎤ ⎢1 + ln⎜1 − = − 2ρ − ⎟⎥ ∂ρ ρ*T * ⎣ ρ * ⎠⎦ ⎝ ρ* +

⎛ ρ ⎞⎤ RT ⎡ ⎢1 + ln⎜ ⎟⎥ M⎣ ⎝ ρ*w ⎠⎦

(14)

and the derivative ∂ρcell,k/∂ri given by eq A6 in Appendix A. The terms of the Helmholtz energy density containing the parameter w (or the thermal wavelength as those derived in Appendix B), which are linear in the density, do not contribute to the forces. Upon integration over the entire domain of the primary simulation box, these terms give a constant contribution to the Helmholtz energy, and their gradient with respect to any bead position is zero. Generation of Initial Configurations. Initial configurations for the linear melt are obtained by field theory inspired Monte Carlo (FT-i MC) equilibration of a coarse-grained melt, wherein chains are represented as freely jointed sequences of Kuhn segments subject to a coarse-grained Helfand Hamiltonian which prevents the density from departing from its mean value anywhere in the system.50,51 The coarse-graining from the freely jointed chain model to the bead−spring model involves placement of beads at regular intervals along the contour of the chains obtained after the MC equilibration. As already discussed in the previous section, at the new (coarser, mesoscopic) level of description, the polymer is envisioned as a network of strands connecting internal nodal points and end points. Starting from a well equilibrated configuration 9 , obtained from a FT-i MC simulation, we determine the box size and shape for which the Gibbs energy function:64 1 G(T , τ ) = A({rij}, {ρ(r)}, T ) − V9 Tr(τ ) 3

vi(t ) =

(16)

The systematic force Fi(t) is the explicit mutual force between the N particles, and -i(t ) represents the effect of the medium on the particles. Each particle is characterized by its mass mi and the friction coefficient ζi, measured in kg/s. The systematic force, Fi, is to be derived from the free energy following eqs 3 and 13:

− V9 ∑ ταβ ϵαβ = G({ri}, τ ) αβ

1 1 Fi({ri(t )}) + -i(t ) ζi ζi

(15)

becomes minimal under the given, externally imposed τ. By prescribing the stress tensor, τ, we minimize the free energy with respect to the positions of the nodal points, {ri}, which we assume to follow the macroscopic strain in an affine way. The presence of any symmetry element in the stress tensor reduces the number of minimization parameters, and in the special case where only hydrostatic pressure is applied on the system, the sum of the last two terms in eq 15 is equivalent to −pV, letting the volume V of the deformed configuration be the only

Fi({ri(t )}) = −∇ri(t )A({rij(t )}, {ρ(r, t )}, T )

(17)

with the expression of the Helmholtz energy as given in eq 1. The stochastic force -i is assumed to be stationary, Markovian, and Gaussian with zero mean and to have no correlation with prior velocities nor with the systematic force. For large values of (ζi/mi)Δt in the diffusive regime, when the friction is so strong 3009


Article

Macromolecules that the velocities relax within Δt, the BD algorithm of van Gunsteren and Berendsen,68 eq 18 can be used:

assigned to every possible transition, p6 → 5 = k 6 → 5 ΔtkMC. The transition is undertaken only if p6 → 5 exceeds a pseudorandom number sampled from a uniform distribution in the interval [0,1). The use of small kMC time steps, ΔtkMC, ensures that ∑5 p6 → 5 ≪ 1 at every step kMC is attempted. In the following we will present a microscopically reversible and physically meaningful formulation for calculating all rates connected with slip spring jumps, slip spring destruction, and formation. Our formulation is further elaborated in Appendices D, E, and F, where detailed considerations are presented. In order to develop a formalism for elementary events of slip spring hopping, creation, or destruction, we need expressions for the rate of slippage along the chain backbone. We concentrate on the slip spring connecting beads a0 and b0. ss) The Helmholtz energy of the system is denoted as A(N a0∧b0, where in our notation we have incorporated both the number of slip springs present and the connectivity at which the particular slip spring finds itself. An individual jump of one end of a slip spring along the chain backbone, e.g., from bead b0 to bead b+, takes place with rate

ri , α(tn + Δt ) = ri , α(tn) + 1⎡ 1 2⎤ ⎢Fi , α(tn)Δt + Fi ,̇ α(tn)(Δt ) ⎦⎥ + 9 i , α(Δt ) ⎣ ζi 2 (18)

where the random variable 9 i, α(Δt ) is sampled from a Gausssian distribution with zero mean and width: ⟨9 i , α 2(Δt )⟩ =

2kBTref Δt ζi

(19)

The time derivatives of the forces, Ḟ i,α, are estimated by finite differences. However, more sophisticated schemes, like smart Monte Carlo, have been developed.35 Slip Spring Kinetic Monte Carlo. Chain slippage through entanglements is modeled as hopping of slip springs along the contours of the chains they connect, without actual displacements of their beads. Let the ends of a slip spring connect beads a0 and b0 along chains a and b (see Figure 2). In order to track

⎛ A‡ ,(Nss) − A (Nss) ⎞ a0 ∧ b0 ⎟ k hop = ν0 exp⎜⎜ − ⎟ kBT ⎝ ⎠

(20)

conforming to a transition state theory (TST) picture of the slippage along the backbone as an infrequent event, which involves a transition from state 6 ≡ a0 ∧ b0 to state ss) 5 ≡ a0 ∧ b+ over a barrier of free energy, A‡,(Nss) − A(N a0∧b0 (cf. Figure 3). Since the connectivity of the Nss − 1 slip springs Figure 2. Illustration of the “hopping” scheme of slip springs along the chain contour.

hopping events in the system, we use a discretized version of the kinetic Monte Carlo (kMC) algorithm, which runs in parallel with the BD integration of beads’ equations of motion. The dynamics of slip spring jumps is envisioned as consisting of infrequent transitions from one state to another, with long periods of relative inactivity between these transitions. We envision that each state corresponds to a single free energy basin, and the long time between transitions arises because the system must surmount an energy barrier, A‡, to get from one basin to another. Then, for each possible escape pathway from an original to a destination basin (e.g., 6 → 5 ), there is a rate constant k 6 → 5 that characterizes the probability, per unit time, that an escape from state 6 to a new one, 5 , will occur. These rate constants are independent of what state preceded state 6 . We sample slip spring transitions at regular predefined time intervals ΔtkMC which are multiples of the time step used for the BD integration, ΔtkMC = nkMCΔt. Every nkMC steps of the BD integration along the trajectory of the system, we freeze the coordinates of all beads and calculate the transition probabilities, khopΔtkMC, for all possible transitions of the slip springs. The kMC time interval, ΔtkMC, is chosen such that the number of transitions performed at every kMC step is much smaller than the number of slip springs present in the system and preferably on the order of unity (i.e., infrequent events governed by Poisson process statistics). The kMC step consists in visiting all slip springs of the system and computing the probabilities for them to perform a jump, considering all possible states accessible from their current one. A probability is

Figure 3. Schematic representation of free energy levels defined in the slip spring hopping formulation. The transition considered is the slippage of the slip spring initially connecting the beads a0 and b0, along chain b, i.e., b0 → b+. The height of the free energy barrier in our simulations was in the range of (3−8)kBT, as indicated in the figure. Please see text for details.

does not change during the transition, we can define A 6‡ → 5 as the difference between the total Helmholtz energy of the configuration with a slip spring at the transition state, A‡,(Nss), and the total Helmholtz energy of a configuration that is missing the particular slip spring but is otherwise identical, A(Nss−1), i.e., A‡ ,(Nss) = A(Nss − 1) + A 6‡ → 5 . Accordingly, the free ss) energy of the starting configuration can be split as A(N a0∧b0 = (Nss−1) + Aa0∧b0. All configurations depicted in Figure 3 are A identical but for the presence and end point positions of the 3010


Article

Macromolecules Nss-th slip spring. We use the logical conjunction operator “∧” to denote the connectivity of the system. As depicted in Figure 2, in the case considered none of the nearest-neighbor beads of the ends of the slip spring are chain ends. Double jumps (e.g., a0 → a+ and b0 → b+) are disallowed, which is a reasonable approximation for small time steps, ΔtkMC. Based on the proper (following Appendix D) selection of the pre-exponential factor νhop, the number of steps between two successive kinetic Monte Carlo steps, nkMC, should be chosen such that the overall probability, phop = khopΔtkMC = khop nkMCΔt, is less than 1 for every slip spring in the system. On the basis of the selected νhop and Δt, we set nkMC = 100. Any number smaller than that can be used and any number larger than that renders hopping probabilities per kMC step of the order of unity and thus has to be avoided. Our results are insensitive to nkMC within the aforementioned range. Assuming that all slip springs attempting to jump face the same free energy barrier, A 6‡ → 5 , on top of A(Nss−1), eq 20 can be more conveniently rewritten as

Figure 4. Reduced rate of slip spring hopping, khop, divided by νhop, as a function of slip spring length (in units of lss, assuming Gaussian springs). The theoretical Gaussian distribution of the slip spring length is also presented in the second ordinate of the figure.

⎛ A‡ ⎞ ⎛A ⎞ 6 → 5 − A a0 ∧ b0 ⎟ = νhop exp⎜ a0 ∧ b0 ⎟ k hop = ν0 exp⎜⎜ − ⎟ kBT ⎝ kBT ⎠ ⎝ ⎠ (21)

where the rate of hopping, khop, depends directly on the energy stored in the slip spring at its initial state (not of the entire configuration), Aa0∧b0, while the dependence on the height of the free energy at the barrier (i.e., A 6‡ → 5 ) has been absorbed into the pre-exponential factor νhop = ν0 exp( −βA 6‡ → 5 ), with β = 1/(kBT). Equation 21 leaves us with only one adjustable parameter to be determined, the pre-exponential frequency factor νhop. In Appendix D, a simple theoretical guideline in order to arrive at an estimate of νhop is presented, given the fact that the average hopping rate, ⟨khop⟩, can be accessed. The representation of slip springs in terms of entropic Gaussian springs allows for a direct estimation of the expected hopping rate, ⟨khop⟩. The distribution of the length of the slip springs, Pe(r), and the corresponding hopping rate, khop(r), are presented in Figure 4. Given the probability distribution, the expected average rate of hopping per slip spring can be estimated as ⟨khop⟩ = ∫ 0∞Pe(r)khop(r)dr/∫ 0∞Pe(r)dr. By assuming a maximum extension of slip springs up to 3lss, a distance excluding less than 10−4 of the cumulative probability distribution, the expected hopping rate is roughly equal to ⟨khop⟩/νhop ≃ 38. Let us assume that the ends of the slip spring connect beads a0 and b0 as those are depicted in Figure 5. One of the beads connected by the slip spring (e.g., b0) is a chain end. One end of the slip spring can hop along the chain, while the second can either hop moving away from the chain end or be destroyed. The rate of destruction of a slip spring is equal to the rate of hopping along the chain, khop. This assumption keeps the necessary adjustable parameters of our model to an absolute minimum. Through the intrinsic dynamics of the system, it is possible for a slip spring end to leave its chain. This process mimics the disentanglement at the chain ends and the process of CR. The process of slip spring destruction is introduced in the model in order to represent the chain disentanglement, as that is envisioned by the polymer tube theories. In the case that both ends of a slip spring are attached to chain ends, based on the scheme we have introduced, there are equal probabilities for

Figure 5. Illustration of the “hopping” scheme of a slip springs whose one end is a chain end.

any of them sliding across any of the chain ends and getting destroyed. As far as the event of creation of a slip spring is concerned, two schemes have been developed; the former preserves the number of slip springs, while the latter allows a fluctuating number of slip springs obeying microscopic reversibility. For the sake of clarity, we will present the constant scheme here, while the fluctuating scheme is presented in Appendix E. Since both schemes yield identical trajectories for the systems studied, we will refrain from further distinction between them in the following. In order to preserve the number of slip springs constant throughout the simulation, we assume that a new slip spring is created whenever an existing slip spring is destroyed through slippage past a chain end. In creating a new slip spring, a free end of a chain can be paired to another bead if they are separated by a distance less than αattempt. Thus, in order to ensure microscopic reversibility, a slip spring can be destroyed only if its extent is less than αattempt at the time the destuction is decided to take place. If a slip spring is considered for destruction, all possible connections of a0 (i.e., the bead that is not a candidate for slippage along its chain contour, see Figure 5) with other beads b lying inside a sphere of radius αattempt centered at a0 are identified. If no such beads can be found, the destruction/ creation move is rejected, due to violation of the microscopic reversibility, since the system would not be able to access its present state through the kMC moves. The Rosenbluth weight 3011


Article

Macromolecules

where C∞ is the characteristic ratio of PI, N is the length of the chain measured in monomers, l is the root-mean-squared average carbon−carbon bond length over an isoprene monomer, and b is the Kuhn length of polyisoprene. The factor of 4 is required because each isoprene contains four backbone bonds. Using Mark’s70 values for C∞ and l of 4.7 and 0.1485 nm, respectively, we obtain b = 0.958 nm and 2.21 isoprene units per Kuhn length.71 The effective length of an isoprene unit has been shown to have a slight temperature dependence;72 however, we ignore this effect in our model. Each bead in our representation consists of nKuhns/bead = 10 PI Kuhn segments. Thus, its mass is mbead = 1506 g/mol, and its characteristic mean-square end-to-end distance (if considered as a random walk) should be (nKuhns/beadb2) = 9.1776 × 10−18 m2. As far as the parametrization of the slip springs is concerned, Fetters et al.69 have estimated the entanglement molecular weight of PI, Me, to be 5430 g/mol, thus allowing us to define the average required number of slip springs present in the system as

Figure 6. Illustration of the slip spring creation procedure.

of the old configuration is accumulated, by iterating over all neighboring beads b, W6 = ∑b exp( −βA a0 ∧ b). A chain end a0′ in the system is randomly selected with probability 1/(2n), with n being the total number of chains. All possible connections of a0′ with beads b′ lying inside a sphere of radius αattempt centered at it are identified (see Figure 6) and the Rosenbluth weight W = ∑′ exp( −βA ) is accumulated. From all candidate 5

b

a0′∧ b ′

anchoring points, b′, for the slip spring emanating from chain end a′0, one is chosen to serve as the end of the new slip spring, b′0 with probability Pa0′ ∧b0′ = exp( −βA a0′∧ b0′ )/W5 . The new connectivity of the system is finally accepted with probability ⎡ ⎛ A a0 ∧ b0 − A a0′∧ b0′ ⎞ W5 ⎤ ⎢ ⎥ P6accept ⎟ → 5 = min 1, exp⎜ − ⎢⎣ kBT ⎝ ⎠ W6 ⎥⎦

⟨Nss⟩ =

(22)

(23)

where khop,a0∧b0 as defined in eq 21 and ΔtkMC is the elapsed time from the previous kMC step. By design the overall destruction/creation scheme satisfies microscopic reversibility. The corresponding proof can be found in Appendix F. Simulation Details. In this section we will discuss the parametrization of our methodology in a bottom-up approach. During the development of the model we have tried to keep its adjustable parameters to an absolute minimum, without introducing new parameters which cannot be mapped to physically relevant observables. We will elaborate on how values of the parameters are obtained with respect to structure, thermodynamics (e.g., equation of state), entanglement density, and friction/hopping rates, either from experimental evidence or from more detailed simulation levels. The systems considered were cis-1,4-polyisoprene (natural rubber) melts in the range of molar masses, M = 21 kg/mol to M = 120 kg/mol. All simulations were carried out in the canonical statistical ensemble, at the temperature of T = 400 K, where rheological measurements are available.69 The simulation box was cubic with varying edge length from 30 to 100 nm. The mean-squared end-to-end distance of a PI chain can be either expressed in terms of the number of isoprene units or Kuhn lengths: ⟨R e,0 2⟩ = C∞4Nl 2 = NKuhnb2

(25)

with n being the number of chains present and Z = M/Me. Based on the chain discretization we have introduced above, the average distance between the slip springs is roughly four coarsegrained beads. It is not at all clear or well-defined what a slip spring is. In this work we try to relate the parameters of the slip springs to quantities obtainable from atomistic simulations or experiments, so as to reduce the number of adjustable parameters. However, one slip spring may not correspond to exactly one entanglement. Other interpretations can be invoked, i.e., thinking of the slip spring as a soft constraint of the motion of the chain perpendicular to its contour and not as a “kink” in the primitive path. The stiffness of the slip springs can be tuned by the parameter lss, introduced in eq 5. We set this equal to a reasonable estimate of the entanglement tube diameter of PI, lss = αpp ≃ 6.2 nm.73 At this point, we should stress an important design rule for the degree of coarse-graining to be chosen. In order for the chains to preserve their unperturbed fractal nature in the melt, entropic springs along them should be stiffer than slip springs, i.e., lss > nKuhns/bead b. If this is not the case, nKuhns/bead should be tuned (i.e., the degree of coarse-graining should be adjusted) in order to ensure that entropic attraction between beads belonging to the same chain is at least an order of magnitude stronger than attraction due to slip springs (connecting mostly beads belonging to different chains). We have used the parameters of the Sanchez−Lacombe equation of state, given by Rudolf et al.,63 who have conducted pVT measurements on several polymers under isothermal conditions, above their glass transition temperature. These authors suggest p* = 383.0 MPa, T* = 631.2 K, and ρ* = 0.961 g/cm3 for molten polyisoprene of M = 2594 g/mol. Based on these values, the density of polyisoprene at T = 300 K is estimated as ρ(300 K) = 0.908 g/cm3. Klopffer et al.74 have characterized the rheological behavior of a series of polybutadienes and polyisoprenes over a wide range of temperatures. The viscoelastic coefficients resulting from the time−temperature superposition principle were determined. A Rouse theory modified for undiluted polymers

If the creation of a new slip spring a0′ ∧b0′ is rejected, based on the criterion of eq 22, the slip spring a0∧b0 remains as is, i.e., anchored to beads a0 and b0. Otherwise, the actual probability for a specific slip spring destruction/creation event to take place at any point in the simulation where changes in slip spring ends are considered is Pa0 ∧ b0 → a0′∧ b0′ = k hop, a0 ∧ b0ΔtkMCP6accept →5

n(Z − 1) 2

(24) 3012


Article

Macromolecules was used to calculate the monomeric friction coefficient, ζ0. The monomeric friction coefficient, ζ0, characterizes the resistance encountered by a monomer unit moving through its surroundings. It was concluded that, within experimental error, a single set of Williams−Landel−Ferry (WLF) parameters75 at Tg was adequate to characterize the relaxation dynamics irrespective of the vinyl content of the polybutadienes and polyisoprenes. These authors proposed that the variation of the monomeric friction coefficient with temperature can be given by log ζ0(T ) = log ζ∞ +

C1gC2g T − Tg + C2g

functions presented in Figure 7 include a number of salient points which should be discussed. Based on our picture of slip

(26)

with the parameters = 13.5 ± 0.2, = 45 ± 3 K, log ζ∞ = −10.4 dyn s cm−1, and Tg = 211.15 K. At a temperature of 298 K, ζ0(298 K) = 1.61 × 10−6 dyn s cm−1, while at a temperature of 400 K, ζ0(400 K) = 1.1508 × 10−11 kg/s. If we think of the friction coefficient as being proportional to the mass of the entity it refers to, we can estimate the friction coefficient of our coarse-grained beads as mbead ζbead = ζ0 mmonomer (27) Cg1

Cg2

Figure 7. Probability density functions of the residence time between two successive slip spring jumps, for two different molar masses of the chains (50 and 100 kg/mol). Two different values of the hopping frequency factor, νhop = 10−1 s−1 and νhop = 101 s−1, are considered in the main figure and in the inset to the figure, respectively.

where mmonomer = 68.12 g/mol refers to the molar mass of a PI monomer. The friction coefficient of a coarse-grained bead at the temperature of 413 K is ζbead = 2.54 × 10−10

kg s

spring hops as infrequent transitions, the hopping process should follow Poisson statistics, which is evident from the inset to the figure, where the distributions of residence times are clearly straight lines in semilogarithmic axes. There are minor deviations at long time scales, which are expected, since the population of slip springs being stationary until that time is extremely small. We have studied two molecular weights, namely 50 and 100 kg/mol. It is evident that the distribution of residence times is independent of the molecular weight of the chains. The use of the same pre-exponential frequency factor, νhop, yields indistinguishable results, as far as the individual motion of the slip springs is concerned. This is crucial for the methodology we have developed, since we will employ the same νhop across different molecular weights, and we expect the differences in chain dynamics to arise from the difference of their sizes (as they should) and not from tuning the slip spring dynamics. In order to evaluate the effective hopping rate of the slip springs, khop, we invoke a hazard-plot analysis79−81 of the hopping events which have taken place along the BD trajectory. The hazard rate, h(t), is defined such that h(t)dt is the probability that a slip spring which has survived a time t in a certain state since its last transition will undergo a transition (i.e., jump) at the time between t and dt. The cumulative hazard is defined as H(t) = ∫ t0h(t′) dt′. By assuming a Poisson process consisting of elementary transitions with first-order kinetics from one state to the other, the Poisson rate can be extracted as the slope of a (linear) plot of the cumulative hazard versus the residence time. Thus, the effective rate constant of the slip spring hopping can be obtained from the hazard plot of the residence time of the slip springs at a specific topology (i.e., the elapsed time between consecutive slip spring jumps), which is presented in Figure 8. At first, we can confirm that the hopping processes follow Poisson statistics, since the main parts of the plots are straight lines. However, there is a significant error

(28)

Moreover, Doxastakis et al.76,77 have estimated the selfdiffusion coefficient of unentangled PI chains consisting of 115 carbon atoms at T = 413 K to be DC115 = 4.4 × 10−11

m2 s

(29)

The length of these chains corresponds to 23 monomers (or 11 PI Kuhn segments). If we use the Rouse model78 to predict the diffusivity of these chains, based on the parameters we have chosen above, that would be DRouse,C115 =

kBT Nmonomersζ0

= 2.15 × 10−11

m2 s

(30)

where Nmonomers is the chain length measured in monomers and ζ0 the monomeric friction coefficient. Our estimation of the self-diffusivity of Rouse chains, based on the model parameters we have introduced, coheres with what Doxastakis et al. have measured both experimentally and by all-atom MD simulations. Finally, as far as the implementation of the model is concerned, we integrate the equations of motion, eq 18, by employing a time step Δt = 10−11 s and up to at least 10 trun with trun = 10−2 s. The kMC scheme is carried out every nkMC = 100 steps, thus resulting in ΔtkMC = 10−9 s. The frequency factor needed for the hopping rates of the kMC scheme is allowed to vary in the range 10−1 ≤ νhop ≤ 10 s−1 and is kept constant for the whole range of molecular weights studied. In the following, its optimal value will be thoroughly discussed.

■

RESULTS AND DISCUSSION Hopping Dynamics. At first we study the distribution of residence times between successive slip spring jumps, in the course of a BD/kMC simulation. The probability density 3013


Article

Macromolecules

Figure 8. Hazard plot of the time between slip spring transitions for systems consisting of two different chain lengths. Linear fits are also included which provide an estimate of the slope, i.e., the effective hopping rate of the slip springs. In the inset to the figure, the dependence of the effective hopping rate, keff, on the hopping preexponential frequency factor, νhop, is presented for the 100 kg/mol system.

Figure 9. Distribution of the distance between successive beads (strands) along the contour of the chains for three different molar masses. In the inset to the figure, the distribution of the length of the slip springs is presented.

springs along the chain are considered Gaussian with an unperturbed length which equals the total length of the Kuhn segments represented by a coarse-grained bead. The conformational features, at least at the strand level, continue to be respected during the BD simulation. Moreover, the distribution of strand lengths during the simulation coincides with a Gaussian distribution centered at the unperturbed strand length, as is theoretically expected. Our simulation scheme produces trajectories of the system consistent with the imposed Helmholtz energy function, eq 2. In the inset to Figure 9 we examine the distribution of slip spring lengths. Slip springs represent entanglements of a chain with its surrounding chains. The polymer tube model considers a tube formed around the primitive path of the chain, which fluctuates in time. Recent simulations have shown that the probability of finding segments of the neighboring chains inside the tube of the chain under consideration is Gaussian.82 This is also the case in our simulations. The use of Gaussian entropic springs for describing the free energy of the slip springs results in a Gaussian distribution of slip spring lengths, conforming to the picture obtained by more detailed simulations and theoretical arguments. The distributions obtained from the simulations are slightly shifted to higher distances. However, they are independent of the molecular weight of the chains under consideration. Finally, we examine the chain dimensions, as those can be quantified by the end-to-end length, Re. In Figure 10 we examine the distribution of end-to-end distances for two molecular weights, lying at the extremes of our range of study. Shorter chains, of M = 21 kg/mol, behave extremely close to their unperturbed, Gaussian statistics. This is not the case for larger chains which become conformationally stiffer compared to the theoretical prediction. End-to-end departures occur at distances commensurate or larger than the average distance between slip spring anchoring beads. At this length scale, slip springs bring about an extension of the chains. However, even this extension does not affect the scaling law ∼ N1/2 of the endto-end length of the chains, as can be observed in the inset to Figure 10, where the end-to-end distance is calculated as a

related to fitting the hazard plots due to the fact that at short time scales a lot of short-lived events are monitored. Given eq 20, knowledge of νhop which is a controlled parameter and the effective rate of slip spring jumps from the hazard plots (i.e., khop) allows us an estimation of the average free energy of the slip springs, Aa0−b0, entering eq 20. The effective hopping rate, khop, is barely sensitive to the frequency factor, νhop, as can be observed in the inset to Figure 8. The dependence of the former on the latter seems logarithmica rule of thumb that can be used if fine-tuning of the dynamics of the model is needed. This observation seems to be in apparent contradiction to the discussion connected with Figure 4. However, it is fully justified, since when deriving a theoretical simplified estimate of the average hopping rate on the pre-exponential frequency factor, a number of assumptions were made. The entanglement creation and destruction processes were not taken into account. If, for example, a slip spring to be destroyed does not find the proper candidates for its new anchoring points, it remains as is, perturbing the statistics. Moreover, slip spring jumps are not fully independent events. The jump of a slip spring at a specific time step may help (or hinder) the remaining ones to relax in the course of the simulation or vice versa. Finally, the slope of the hazard plots that we have used in order to estimate the hopping rate is sensitive to the recrossing events taking place at short time scales (e.g., a specific slip spring jumping back and forth). Summarizing, the effective hopping rate presented in Figure 8 incorporates cooperative effects between different slip springs and reduction of slip spring jumps due to slip spring creation rules employed. However, we expect that the linear dependence may be recovered for higher νhop frequency factors. Structural Features. We start by examining the structure of the melt chains at the level of individual strands, where we refer to a strand as the distance between successive beads along the contour of a chain. The distribution of the length of the strands (entropic springs) that connect the beads along the contour of the chains is depicted in Figure 9. The entropic 3014


Article

Macromolecules κT ≡ −

∂ ln ρ ̃ ⎞ 1 ∂V ⎞ ⎟ = ⎟ ∂p ⎠ V ∂p ⎠ T T −4

(32) −1

which is equal to 1.27 × 10 bar . The small difference (lower compressibility in our simulation) is probably attributable to the additional contributions to the Helmholtz energy, i.e., entropy springs and slip springs. Our model can faithfully reproduce the compressibility of polymeric melts with reasonable accuracy. Realistic compressibility is missing from the majority of coarse-grained models and methods. Polymer Dynamics. A rigorous way of studying the mobility of a polymeric melt is to calculate the mean-squared displacement (MSD) of the structural units of the chains. In order to avoid chain end effects,83,84 only the innermost beads along the chain contribute to the calculations g1(t ) = Figure 10. Distribution of the end-to-end distances of the chains for two different molecular weights. In the inset to the figure, the scaling of the end-to-end distance is presented, with respect to the number of beads, Nbeads, of subchains along the original chains.

2n inner + 1

i = N /2 + n inner

∑

⟨(ri(t0 + t ) − ri(t0))2 ⟩

i = N /2 − n inner

(33)

with the value of the parameter ninner quantifying the number of innermost atoms, on each side of the middle segment of each chain that are monitored. In our case, ninner is set in such a way that we track more than half of the chain, excluding one-fifth of the chain close to one end and one-fifth close to the other end. Moreover, the motion of the chains can be described in terms of the mean-squared displacement of their centers of mass, g3(t):

function of the number of beads of subchains. Polymer chains, in both cases, behave as random walks, with a larger “Kuhn” step, though. Our minimalistic approach has significantly reduced, but not eliminated, perturbations in chain conformations caused by the unphysical attraction introduced by the slip springs. This attraction is not causing obvious problems under the considered conditions but is always present in our simulations. The only way of restoring chain conformations is by introducing a compensating potential as that of Chappa et al.,34 which fully compensates the effect of slip springs on equilibrium properties. Thermodynamics. Thermodynamics is naturally introduced in our model via the nonbonded energy which is dictated by a suitable equation of state (e.g., the Sanchez− Lacombe equation of state). Several thermodynamic properties can be calculated. As a proof of concept, the compressibility of our PI melts can be calculated. We can calculate the time evolution of local densities at the cells of our computational grid, as presented in Appendix A. By treating each cell of the grid as a system capable of exchanging mass with its surroundings, compressibility can be estimated by the following fluctuation formula: ⟨(δnKuhns/cell)2 ⟩ κT = V kBT ⟨nKuhns/cell⟩2

1

g3(t ) = ⟨(rcm(t0 + t ) − rcm(t0))2 ⟩

(34)

Figure 11 contains results for the time dependence of the mean-square displacements g1(t) and g3(t) for two PI melts of different molecular weights. To a first approximation, the dynamics of short chains in a melt can be described by the Rouse model.78 As far as the lower molecular weight chains are concerned, it can be seen in Figure

(31)

where the averages ⟨...⟩ are taken over all cells of the grid and along the trajectory of the simulation. The averages appearing in the numerator and denominator of eq 31 are the variance and the mean number of Kuhn segments assigned to each cell, respectively. It should be noted that local densities of grid cells are obtained via the smearing scheme introduced in Appendix A. The isothermal compressibility estimated from our simulations varies from 1.09 × 10−4 to 1.04 × 10−4 bar−1 for samples from 21 to 120 kg/mol, respectively. On the other hand, starting from the macroscopic equation of state, the compressibility is

Figure 11. Time evolution of the mean-square displacements of the inner beads (g1(t)) and center of mass of the chains (g3(t)) for M = 21 and 100 kg/mol. Slip springs are present. Alongside, the scaling regimes of the tube model are also included. 3015


Article

Macromolecules 11 that the mean-square center-of-mass displacement, g3(t), remains almost linear at all times; this means the intermolecular forces between polymers are too weak to affect diffusive behavior and play a minor role compared to the bonded interactions. However, small departures are fully justified due to the few slip springs present even it the short-chain system. The bead mean-squared displacements, g1(t), exhibit a subdiffusive behavior that arises from chain connectivity and are characterized by a power law of the form g1(t) ∼ t1/2. After an initial relaxation time where a change in g1(t) occurs, a regular diffusive regime is entered, where g1(t) ∼ t. This sequence of scaling trends is predicted by the Rouse theory. The limiting behavior of the chains’ center-of-mass displacement can yield an estimate of the diffusivity of the chains: limt→∞ g3(t) = 6Dcmt, which has been found in excellent agreement with the diffusivity predicted by all-atom MD simulations of PI chains of the same molecular weight by Doxastakis et al.76 The introduction of slip springs in those short-chain systems does not seem to affect the scaling laws of the unentangled melt significantly. Moreover, Figure 11 includes results for the mean-square displacements, g1(t) and g3(t), as a function of time, obtained from the simulation of a 100 kg/mol PI melt with slip springs present. At short time scales, the segmental MSD curves, g1(t), for both melts coincide. However, at longer time scales, the characteristic signature of an entangled polymer system can be observed in the longer-chain system. It can be seen that at short times the bead mean-square displacement, g1(t), shows a scaling regime with a power law t1/2; at intermediate times, a regime with a power law t1/4 appears; eventually we observe a crossover to regular diffusion at long times. The mean-square displacement of the chain center-of-mass, g3(t), also exhibits subdiffusive behavior at intermediate times, with a scaling behavior t1/2, as predicted by the tube model; at long times, regular diffusion is achieved. From the long time behavior of g3(t), we can estimate the longest relaxation time, τd. It should be noted that all transitions between scaling regimes are very smooth in our model. Despite the fact that crossovers between the scaling regions are not accurately discerned in our results, the scaling of polymer melts, as expected from tube theory, is observed. From the long time behavior of the mean-square displacement of the center of mass of the chains, g3(t), we compute the diffusion coefficient and the longest relaxation time (disengagement time), τd, as a function of molecular weight. The results are presented in Figure 12. As expected, D is a monotonically decreasing function of M (or N) with a power scaling law of −2 (tube model) or −2.3 (experimental observations85) for wellentangled melts. Our slip spring simulations can recover the aforementioned scaling behavior of the diffusion coefficient, and their predictions lie in the correct order of magnitude, cohering with the results from detailed atomistic simulations of Doxastakis et al.77 In the inset to Figure 12, the molecular weight dependence of the disengagement time, τd, is presented. For large M (or N), the longest relaxation time should grow as N3 (as predicted from the tube model) or with slightly higher exponents (around 3.586). Our results fit well to a scaling law of 3.3, which is reasonable. Experimental results at a significantly lower temperature86 are also included for comparison. The exponent of the experimental measurements is slightly larger (around 3.5) for PI measurements at a temperature of 300 K. Rheology. Linear rheological properties can be characterized through the shear relaxation modulus, G(t) = τxy(t)/γ,

Figure 12. Diffusion coefficient, D, plotted against the molecular weight of the chains, M. The dashed line provides a guide to the eye showing the M−2.3 scaling law.85 Results from earlier atomistic simulations are also included.77 In the inset, the disengagement time, τd, is plotted as a function of the molecular weight, M. The dashed line is a power law fit with an exponent 3.3, and experimental estimates of τd at a lower temperature are also included.86

with γ being a small shear deformation and xy two orthogonal axes. In computer simulations, the most convenient way of evaluating G(t) is to use the fluctuation−dissipation theorem:87 G (t ) =

V ⟨ταβ(t0 + t )ταβ(t0)⟩ kBT

(35)

where αβ stand for any two orthogonal directions. One can show that the stress relaxation after any small deformation will be proportional to G(t), and thus G(t) fully characterizes the linear rheology of a polymeric system with given external parameters. However, the stress autocorrelation function (acf) is notoriously difficult to calculate due to huge fluctuations at long times. In order to improve the accuracy of our calculations, we average ⟨ταβ(t0 + t)ταβ(t0)⟩ over all possible ways of selecting a pair of perpendicular axes α and β.87 We compute the time correlation functions in eq 35 by using the multipletau correlator algorithm of Ramı ́rez et al.88 Let us first consider Figure 13, which displays the time evolution of the shear relaxation modulus, as defined in eq 35, for PI melts of different chain lengths. The initial value of the modulus is ρkBT, with ρ being the mean number density of beads in the polymer melt (slightly higher for longer chains). The relaxation of the shortest chains (21 kg/mol) considered is typical of unentangled polymer melts and is close to the shear relaxation modulus predicted by the Rouse model. In particular, the intermediate time scale behavior of the stress acf is consistent with the Rouse model scaling, G(t) ∼ t−1/2, while the long time decay is exponential, with a time constant characterizing the longest stress relaxation time in the system. It is clear that the longest relaxation time from the stress autocorrelation function follows the same power law scaling with N as the end-to-end vector relaxation time. Upon increasing the chain length, a plateau starts to appear; this indicates that the viscoelastic character of polymer melts is captured by our model. However, we should note the log−log 3016


Article

Macromolecules

balance between all relaxation mechanisms involved in the rheology of linear melts.91,92 Thus, compared to the recent experimental results, we may reasonably underestimate the exponent (slope in the log−log plot of η versus N), given that the slope in the log−log plot of D versus N is reasonable. Furthermore, our simulation results are in good agreement with previous slip spring simulations of other research groups,18,40 which exhibit the same or even lower scaling exponent. In experimental practice, the cosine and sine Fourier transforms of G(t) can be measured:

∫0

G″(ω) = ω

∫0

m

G (t ) =

(37)

∞

cos(ωt )G(t ) dt

⎛ t⎞ ⎟ ⎝ τi ⎠

∑ Gi exp⎜− i=1

(38)

(39)

where Gi and τi are the amplitude and the relaxation time of the mode number i (often called the Maxwell modes). Our results for the storage and loss moduli of the PI melts under study are presented in Figure 14.

scale of Figure 13, which implies that the rubbery plateau still involves significant relaxation, even for M = 120 kg/mol. The shear viscosity, η, can be computed from the stress relaxation function through a Green−Kubo relation:

∫0

sin(ωt )G(t ) dt

by applying small oscillatory deformation and recording the inphase and out-of-phase responses G′ (ω) and G″ (ω), respectively, which are called storage and loss moduli. A straightforward way to obtain the Fourier transform of G(t) is to fit it with a sum of exponentials:

Figure 13. Temporal evolution of the shear relaxation modulus, G(t), for a series of PI melts at T = 400 K. In the inset to the figure, the zero shear viscosity of the melts is presented, alongside experimental results, at a slightly lower temperature (T = 373 K).89

η=

∞

G′(ω) = ω

∞

Gαβ (t ) dt

(36)

In the inset to Figure 13 we present the viscosity as a function of the molecular weight. The viscosity, η, is a monotonically increasing function of N, which exhibits a scaling behavior N3.4, as experimental data suggest.89 It should be noted that the tube model predicts a N3 dependence of viscosity. Our simulation results seem to cohere with the N3.4 scaling. Moreover, there is very good quantitative agreement with available experimental data for the same system. Experimental results have been obtained at a slightly lower temperature, that being the reason for higher viscosity values for the same molecular weight. It is very promising that our methodology, parametrized in a bottom-up approach, captures the correct values of viscosity without any adjustment or rescaling of the shear relaxation function. The precise scaling of viscosity of polymer melts is still elusive and remains an open question for theoretical, simulation, and experimental studies.90 The results presented by Auhl et al.86 imply a M3.6 scaling law for polyisoprene melts, in contrast to the established M3.4 scaling of the early experimental studies. A slip spring simulation scheme can capture modes of entangled motion beyond pure reptation (M3 scaling). In linear response contour length fluctuation (CLF), the Brownian fluctuation of the length of the entanglement path through the melt modifies early time relaxation. Similarly, the process of constraint release (CR), by which the reptation of surrounding chains endows the tube constraints on a probe chain with finite lifetimes, contributes to the conformational relaxation of chains at longer times. Both processes of CLF and CR exist in slip spring models, and their relative importance contributes to the quantitative prediction of linear rheology. In our case, achieving a M3.4 scaling law can be considered as a

Figure 14. Loss and storage moduli for monodisperse polyisoprene samples of different molecular weights. Measured plateau modulus from Fetters et al.69 is also included.

The forms of the G′(ω) and G″ curves for PI are typical of terminal and plateau responses for monodisperse linear polymers in the entangled region.89,93 For short chains of M = 21 kg/mol, relaxation is rather fast; the longest relaxation time being on the order of 10−4 s. For longer chains, several important features arise. At low frequencies, G′ ∼ ω2 and G″ ∼ ω, providing values of the zero-shear viscosity η and steadystate recoverable compliance, J0e :94 3017


Article

Macromolecules G″(ω) ω

(40)

G′(ω) 1 lim η2 ω→ 0 ω 2

(41)

η = lim

ω→ 0

Je0 =

When moving to higher frequencies, G″(ω) passes through a shallow maximum, Gmax ″ at ωmax. The terminal loss peak becomes progressively more well-defined (especially for 100 and 120 kg/mol), i.e., further separated from the transition response, with increasing chain length. The plateau modulus, G0,N sim can be determined by an integration over a fully resolved terminal loss peak:89 G N0,sim =

2 π

∞

∫−∞ G″(ω) d ln ω = 0.39 MPa

(42)

which is in excellent agreement with the measured plateau modulus, G0N = 0.42 MPa. The majority of the plateau modulus values reported in the literature have been obtained by this method.69,95,96 Even for the longest chains, the loss modulus (open symbols in Figure 14) shows that some relaxation is nevertheless present, revealing a power law close to ω−1/4. It is often observed that rheological response measured at different temperatures can be reduced to a single master curve at the reference temperature T0 if one shifts the time (or frequency) appropriately. According to the time−temperature superposition (TTS) principle, the complex relaxation modulus of rheologically simple polymers, defined as G*(ω) = G′(ω) + iG″(ω), measured at different temperatures, obeys G*(ω , T ) = )(T )G*(((T )ω , T0)

Figure 15. Loss and storage moduli for monodisperse polyisoprene samples of 100 kg/mol. Experimental results have been taken from Auhl et al.86 and have been properly shifted to the simulation temperature. Measured plateau modulus from Fetters et al.69 is also included.

(43)

where T0 is the reference temperature, and ((T ) and )(T ) are horizontal and vertical shift factors for this particular reference temperature.94 The range of validity of TTS includes only those processes whose rates are proportional to the rate of one particular basic process (e.g., monomer diffusion), which in turn depends on temperature. In polymer melts, the shift factors are usually given by empirical equations: log ((T ) = )(T ) =

C1(T − T0 + C2/M ) T − T∞ + C2/M

ρ (T )T ρ(T0)T0

Figure 16. Schematic representation of a grid cell and a nodal point with its surrounding cube.

(44)

in quantitative agreement with the experimental observations. The minor discrepancies observed can be attributed to several reasons. In the case of our simulations, polyisoprene melts are strictly monodisperse, which does not hold for the experimental studies, where a finite polydispersity exists. The plateau value obtained from the simulation curve around 1/τe is the same as the one measured by rheological experiments as well as the shape of the whole storage modulus curve in the frequency range of 102−106 s−1. The use of the kMC scheme for the slip spring hopping introduces a relaxation time scale on the order of 10−6 s, which accounts for the different slope of the simulated storage modulus curve with respect to the experimental in the range ω > 106 s−1. As far as the loss modulus is concerned, discrepancies start to appear at high frequencies (close to the shortest relaxation time, τ0 ≃ ζbK2/ (3π2kBT)), which correspond to glassy behavior that cannot be captured by a coarse-grained model which omits truly microscopic, atomistic processes. Experimental loss modulus curves are monotonically increasing with a power law ω1/2 at high frequencies, which is not the case for our soft, coarsegrained model.

(45)

with the former being the Williams−Landel−Ferry (WLF)75 equation (with two material constants C1 and C2) and the latter incorporating the dependence of polymer density on temperature. In our case, we use the experimentally derived curves of Auhl et al.,86 horizontally and vertically shifted to our temperature of interest (400 K). The values of the relevant parameters are T0 = 25 °C, T∞ = −114.03 °C, C1 = 4.986, and C2 = 14.65 with M measured in kg/mol.86 The relation for the temperature dependence of the PI density was obtained from the Sanchez−Lacombe equation of state used for describing the nonbonded interactions of the our coarse-grained model. In Figure 15, the predicted moduli of the polyisoprene sample of 100 kg/mol are compared against the experimental measurements of Auhl et al.86 shifted to the temperature of the simulations following the TTS principle described in the previous paragraph. Despite the fact that experimental measurements have been conducted for a slightly lower molecular weight, M = 94.9 kg/mol, the predicted moduli are 3018


Article

Macromolecules

functional of the local density by passing an orthogonal grid through the simulation box. Slip springs are envisioned as connecting nodes on different polymer chains or on the same chain. Equations for a stochastic description of the dynamics are derived from the coarse-grained Helmholtz energy function. All beads execute Brownian motion in the high friction limit. The ends of the slip springs execute thermally activated hops between adjacent beads along chain backbones, these hops being tracked by kMC simulation. In addition, creation/ destruction processes are included for the slip springs. A slip spring is destroyed when one of its ends slips past the free end of a polymer chain. A new slip spring is formed when a chain end captures a bead of another chain lying within a certain radius from it, according to a prescribed rate constant. Parameters needed in the model are derived from experimental volumetric and melt viscosity data or from atomistic molecular dynamics simulations. Initial configurations for the network are obtained from MC simulations of linear melts. Tests of the simulation code on molten linear (non-cross-linked) cis-1,4polyisprene of high molar mass at equilibrium have given satisfactory results for the mean-square displacement of beads and for the shear relaxation modulus. By following the rigorous procedures described before, a consistent methodology has been set in place. To implement this methodology, we have developed a general purpose inhouse software program to simulate polymer melts of arbitrary geometry. Predicting the rheology of binary blends of chains of different geometries poses formidable challenges for tube and slip link models.97 Nonlinear rheology can also be studied with the formulation developed in this work, and strain rate experiments can be conducted.98 The methodology has already been implemented to study tensile deformations of cross-linked cis-1,4-polyisoprene, with very encouraging results.99 Morover, it can be readily extended to systems containing solid surfaces, particles, or cavities.100 Therefore, when all relevant relaxation processes are correctly treated, the viscoelastic response of the melt should be fully obtained from the simulation findings.

Figure 17. Total and intermolecular pair distribution functions of polymeric beads accumulated in 0.5 nm thick bins. The edge length of the nonbonded discretization cells, lx = ly = lz, is 5 nm. In the inset to the figure, the compressibility is presented as a function of the observation length scale. Compressibility is estimated from the static structure factor, Stot(q), of a 80 kg/mol polyisoprene melt. The experimental value, κT = 1.27 × 10−4 bar−1, predicted from the Sanchez−Lacombe EOS, is also included.63

Table 1. Sanchez−Lacombe Notation symbol n M r V ρ ρmol = n/V v* rv* V* = n(r v*) ρ* = M/(rv*NA) T*= ε*/kB p*= ε*/v* T̃ = T/T* p̃ = p/p* ρ̃ = ρ/ρ* = ρmolM/ (ρ*NA) ṽ = 1/ρ̃ = V/V*

■

explanation

units

number of molecules molecular weight Sanchez−Lacombe segments per molecule volume of the system mass density of the system molecular density close packed volume of a mer close packed volume of an r-mer molecule close packed volume of the n r-mers close packed mass density characteristic temperature characteristic pressure reduced temperature reduced pressure reduced density

− kg −

reduced volume

−

m3 kg/m3 m−3 m3 m3

■

APPENDIX A. DISCRETE NONBONDED INTERACTIONS SCHEME FOR SLIP SPRING SIMULATIONS The simplest option for relating the positions and masses of the beads to ρcell is to envision each bead i as a cube containing mass mi, of edge length hi, centered at position ri = (xi, yi, zi), as shown in Figure 16. Our scheme was designed to permit fast analytical calculation of nonbonded forces. It ensures continuity of the nonbonded forces as functions of bead position, but not of the derivatives of these forces (which are not needed in our simulation setup). It has proven satisfactory for use in the Brownian dynamics simulations conducted in this work and can be readily extended to systems containing solid surfaces, particles, or cavities. Other schemes for determining the instantaneous density field from the particle positions,51 including gridless ones,101 seem not to be straightforwardly usable in the aforementioned cases. The cell dimensions along the x, y, and z directions will be denoted as lx, ly, and lz, respectively. We will focus on a cell extending between xcell − lx and xcell along the x-direction, between ycell − ly and ycell along the y-direction, and between zcell − lz and zcell along the z-direction. In the regular grid considered, if (0,0,0) is taken as one of the grid points, xcell, ycell, and zcell will be integer multiples of lx, ly, and lz, respectively. In the following, we assume that hi < min(lx, ly, lz).

m3 kg/m3 K kg/(ms2) − −

SUMMARY AND CONCLUSIONS The first steps toward a consistent coarse-grained model capable of reproducing the rheological properties of polymer melts have been presented. The methodology and the corresponding computer code have been developed for the case of a pure polymer melt. Chains are modeled as sequences of beads, each bead encompassing approximately 10 Kuhn segments. The Helmholtz energy of the system is written as a sum of three contributions: entropy springs, representing the entropic elasticity of chain strands between beads; slip springs, representing entanglements; and nonbonded interactions. The Helmholtz energy of nonbonded interactions is estimated by invoking an arbitrary equation of state and is computed as a 3019


Article

Macromolecules The mass contributed by the node to the cell is mi ,cell

V = mi cubei ∩ cell Vcubei

ρcell = (A1)

⎧⎡ ⎛ ⎛ ⎞⎤ ⎫ ⎞ h h Vcube i ∩ cell = max⎨⎢min⎜xi + i , xcell⎟ − max⎜xi − i , xcell − lx⎟⎥, 0⎬ ⎝ ⎠ ⎝ ⎠⎦ ⎭ ⎣ 2 2 ⎩

(A3)

i

hi < xcell + hi 2 h − ly < yi + i < ycell + hi 2 h − lz < zi + i < zcell + hi 2

xcell − lx < xi +

⎧⎡ ⎛ ⎞ ⎛ ⎞⎤ ⎫ h h × max⎨⎢min⎜yi + i , ycell ⎟ − max⎜yi − i , ycell − ly⎟⎥, 0⎬ ⎠ ⎝ ⎠⎦ ⎭ 2 2 ⎩⎣ ⎝

ycell

⎧⎡ ⎛ ⎞⎤ ⎫ ⎛ ⎞ h h × max⎨⎢min⎜zi + i , zcell⎟ − max⎜zi − i , zcell − lz⎟⎥, 0⎬ ⎠⎦ ⎭ ⎝ ⎠ 2 2 ⎩⎣ ⎝

zcell

(A2)

(A4)

According to the above approach, the force on node i due to nonbonded interactions is

As defined by eq A2, Vcube i∩cell is a linear function of the bead coordinates. Clearly, if cube i lies entirely within the cell, Vcube i∩cell = Vcube i and, consequently, mi,cell = mi. If, however, the borders of the cube i intersect the borders of the considered cell, then node i will contribute a mass mi,cell < mi to the cell. The total mass contributed by bead i to all cells in which it participates will always be mi. The density ρcell in the considered cell is estimated as

∑

Fi = −∇ri A nb = −

acc V cell

cells ∩ cube i

∂a vol(ρ , T ) ∂ρ

∇ri ρcell ρ = ρcell

(A5)

From eqs A1, A2 and A3 one can obtain the components of ∇ri ρcell. Along the x direction,

⎧ ⎪0 ⎪ ⎪ m Vcubei ∩ cell 1 i ⎪ acc hi ⎪ V cell Vcubei xi + 2 − xcell + lx ⎪ ⎪ = ⎨0 ⎪ ⎪ m V 1 ⎪− i cubei ∩ cell acc ⎪ V cell Vcubei xcell − xi + hi 2 ⎪ ⎪ ⎪0 ⎩

and similarly for ∂ρcell/∂yi and ∂ρcell/∂zi. The derivatives are bounded but not continuous. To make them continuous, an extension of the “smearing scheme” for beads which uses a continuous density distribution for the contribution of each bead would be required. Two nodes whose cubes lie entirely within a cell experience the same nonbonded force (zero). This is not true, however, of nodes whose cubes intersect cell borders. The edge length of the cube assigned to node i, hi, can be set approximately equal to the root-mean-square end-to-end distance of the strands assigned to a node: ⎞1/2 ⎛m hi = ⎜ i bK 2⎟ ⎠ ⎝ mK

∑ mi ,cell

Clearly, only beads i whose cubes have a nonzero overlap with the considered cell will contribute to the summation of eq A3. The position vectors ri of these nodal points will necessarily lie within the considered cell or its immediate neighbors. The precise conditions for cube i to have common points with the considered cell are

with Vcube i∩cell being the volume of the intersection of cube i, associated with bead i, and the considered cell, while Vcube i = hi3 is the volume of cube i. Under the condition hi < min(lx, ly, lz), Vcube i∩cell is obtainable as

∂ ρ ∂xi cell

1 acc V cell

if xi ≤ xcell − lx − if xcell − lx −

hi 2

hi h < xi < xcell − lx + i 2 2

hi h ≤ xi ≤ xcell − i 2 2 h h − i < xi < xcell + i 2 2

if xcell − lx + if xcell

if xi ≥ xcell +

hi 2

(A6)

between polymeric beads. On very short length scales, the total distribution function diverges as 1/r because of the intramolecular contribution dictated by the Gaussian chain model. This spike formed is an artificial effect, caused by the fact that chain self-intersection is not prevented in the Gaussian chains used in our calculations. If an atomistic description were used, this spike would be replaced by a series of intramolecular peaks reflecting the bonded geometry and conformational preferences of polyisoprene molecules. In the regime of small distances, segments of other molecules are expelled from the volume of a reference chain and this gives rise to a slight “correlation hole” effect in the intermolecular g(r). For large length scales, the intermolecular distribution function approaches unity, as the intramolecular distribution function approaches zero. The use of the nonbonded scheme employed in this work leads to a uniform density profile throughout the volume of the simulation box (the nonbonded discretization length employed for Figure 17 was 5 nm).

(A7)

where mK and bK are the mass and the length of a Kuhn segment of the polymer under consideration. Figure 17 shows the variation of the total and the intermolecular pair distribution functions with the distance 3020


Article

Macromolecules ⎡ ⎛ 1⎞ ⎤ ρ 2̃ + p ̃ + T̃ ⎢ln(1 − ρ ̃) + ⎜1 − ⎟ρ ̃⎥ = 0 ⎝ ⎣ r⎠ ⎦

In a dense polymer melt, density fluctuations are strongly suppressed on large length scales, resulting in the screening of excluded volume interactions among the polymer beads. The correlations of local segment density in homopolymer melts are characterized by the structure factor Stot(q), i.e., the Fourier transform with respect to distance r − r′ of the density correlation function ⟨ρ(r) ρ(r′)⟩, where brackets denote averaging in the canonical ensemble. Since our simulations deal with explicit particle coordinates, Stot(q) can be calculated according to 1 Stot(q) = nN

∑e

It should be noted that, in the isothermal−isobaric ensemble, p̃ and T̃ are the independent variables, while ρ̃ is the dependent one. Therefore, eq B3 defines the value of ρ̃ at given (T̃ , p̃) that minimizes the free energy. Equations B1 and B3 contain the complete thermodynamic description of the model fluid. Since the Sanchez−Lacombe theory is not a corresponding states theory, the state parameters are predicted to vary with chain length. It has been shown that thermal expansion coefficients, compressibilities, and free volumes are predicted by the Sanchez−Lacombe EoS to decrease with increasing degree of polymerization, in agreement with experiment. The decrease of the free volume with increasing molecular weight is supported by the observation that the glass temperature increases with increasing molecular weight.102 If we wish to obtain some physical insight into the parameter w, we can compare the Sanchez−Lacombe EoS with the ideal gas EoS. The former should reduce to the latter, in the limit of zero molecular density, ρmol = n/V:

2

nN −iq·rj

j=1

(A8)

where rj stands for the position vector of bead j. The limiting value of the total structure factor at large length scales, q → 0, is directly connected to the compressibility of the system: lim Stot(q) =

|q|→ 0

1 kBTρ0 κT

(A9)

lim ASL = lim + Aid.g.

where ρ0 is the mean segment density and κT the isothermal compressibility. By calculating Stot(q) for a range of q-vectors, the compressibility of the system can be estimated as a function of the observation length scale. In the inset to Figure 17 the compressibility of a polyisoprene system is plotted as a function of the scattering vector q used for the calculation of the total structure factor. Coarse-grained simulations with soft interactions are performed with potential expressions and parameters which yield more compressible systems than their atomistically simulated counterparts. This applies to our simulations, too. It is evident that the compressibility of the system deviates from the macroscopic one (∼10−4 bar−1), as q increases (i.e., as one moves to smaller length scales). However, at large length scales compressibility values lie in the range expected for experimental systems.

ρmol → 0+

(B4)

ρmol → 0

where we have used ASL and Aid.g. to denote the Helmholtz energy obtained from the Sanchez−Lacombe and the ideal gas EoS, respectively. The Helmholtz energy can be obtained from eq B1: ASL (ρ , T ) = G − pV = G − p ̃

nrε* ρ̃

(B5)

⎡ 1 ⎛ ρ ̃ ⎞⎤ =−nrkBT *ρ ̃ + nrkBT ⎢(v ̃ − 1) ln(1 − ρ ̃) + ln⎜ ⎟⎥ ⎣ r ⎝ w ⎠⎦

Hence, the limit appearing in the left-hand side of eq B4 reduces to

■

⎡ ⎤ ⎛ρ M⎞ lim + ASL = nkBT ⎢ lim + ln⎜ mol ⎟ − r ⎥ ⎢⎣ ρmol → 0 ⎝ NAρ*w ⎠ ⎥⎦ ρmol → 0

APPENDIX B. EQUATION OF STATE CONSIDERATIONS One of the most widely used equations of state (EOS) for polymer fluids is the one derived by Sanchez and Lacombe.62 These authors have employed a lattice formulation, wherein the polymer chains are occupying discrete lattice sites, while there also exist vacant lattice sites (holes). The Gibbs energy, based on the Sanchez−Lacombe equation of state, can be expressed in dimensionless variables as

(B6)

On the other hand, the Helmholtz energy of an ideal gas can be written as Aid.g.(ρ , T ) =

Gid.g.(ρ , T ) − pid.g. V

=

nμid.g. (ρ , T ) − nkBT

⎡ ⎛ ρN ∏ Λ 3 ⎞ ⎤ A i i ⎟ = nkBT ⎢ln⎜⎜ − 1⎥ ⎟ ⎢⎣ ⎝ M A0 ⎠ ⎥⎦

⎡ G 1 ⎛ ρ ̃ ⎞⎤ G̃ ≡ = − ρ ̃ + pṽ ̃ + T̃ ⎢(v ̃ − 1) ln(1 − ρ ̃) + ln⎜ ⎟⎥ ⎣ nrε* r ⎝ w ⎠⎦ (B1)

(B7)

with Λi being the thermal wavelength of atom i of a molecule of the fluid and A0 the configurational integral of a single molecule with respect to all but three translational degrees of freedom. The molar mass of the molecule is denoted by M. At this point we can introduce the density of molecules, ρmol, in eq B7 and consider its limit at ρmol → 0+:

where T̃ , p̃, ṽ, and ρ̃ are the reduced temperature, pressure, volume, and density. The parameter w is connected to the number of different configurations available to a system of n rmers. It will be examined in detail later. The Sanchez−Lacombe parameters are presented in Table 1. The corresponding equation of state (EoS) can be extracted, by minimizing G: ∂G̃ ⎞ ⎟ =0 ∂v ̃ ⎠T̃ , p̃

(B3)

⎡ ⎛ ρ ∏ Λ 3 ⎞⎤ i ⎥ ⎟⎟ lim + Aid.g. = nkBT lim +⎢ln⎜⎜ mol i ρmol → 0 ρmol → 0 ⎢ A e ⎠⎥⎦ 0 ⎣ ⎝

(B2)

(B8)

Finally, by setting the right-hand side term of eq B6 equal to the right-hand side term of eq B8, we obtain

which yields the equation of state for the system: 3021


Article

Macromolecules ⎡ ⎛ ρ ∏ Λ 3 ⎞⎤ ⎡ ⎛ ρ M ⎞⎤ i ⎥ ⎟⎟ lim +⎢ln⎜ mol r ⎟⎥ = lim +⎢ln⎜⎜ mol i ⎥ ⎢ ρmol → 0 ⎢ ρ → 0 N we e A ρ * ⎠ ⎝ ⎦ ⎣ mol ⎠⎥⎦ A 0 ⎣ ⎝

ταβ = (B9)

b ({rij}, T ) = ταβ

⎛ ∂(A /m) ⎞T ⎟ τ = ρ · ⎜ ⎝ ∂ ⎠

(C7)

nb ({ρcell }, ταβ

3

b ταβ ({rij},

∂A ∂βγ

=

⎧

⎫ ⎡ ∂A pair (rij , T ) ∂rij ⎤⎪ ⎥⎬ ∂rij ∂βγ ⎥⎦⎪ ⎩ (i , j) ⎢⎣ ⎭

3

⎪

γ=1

∑ (i , j)

1 V9

∑

1 V9

∑

∂A pair (rij , T ) ∂rij

γ=1

∂A pair (rij , T ) ∂rij

(i , j)

3

∑

3

∑ γ=1

∂A pair (rij , T ) rij , β

(i , j)

∂rij

rij

∂rij ∂βγ

αγ

∂rij ∂rij , β ∂rij , β ∂βγ 3

∑ γ=1

∂rij , β ∂βγ

αγ

αγ (C9)

where we have made use of the fact that the partial derivative of the Euclidean norm of a vector a = (a1, a2, ..., an) with respect to one of its components, aj, is ⎛ ∂∥a∥ ∂ ⎜ = ∂aj ∂aj ⎜⎝

(C4)

n

⎞

i=1

⎠

∑ ai 2 ⎟⎟ =

1 2∥a∥

n

⎡ ∂ ⎤ aj (ai 2)⎥ = ∥a∥ ⎥⎦ ⎣⎢ ∂aj

∑⎢ i=1

(C10)

It should be noted that eq C9 is valid for any kind of deformation (both affine and nonaffine). We envision that, at a certain time, a homogeneous deformation is applied on the polymer that displaces bond ends affinely, in the sense that their positions are changed in the same way as material points in a macroscopic continuum description. Straight parallel lines in the reference configuration map to straight parallel lines in the deformed configuration Let ri and rj be the positions of the start and the end of the bonded beads before the deformation and ri′ and rj′ the positions of the same beads after the deformation; then

∂(A({rij}, {ρcell }, T )/m) ∂βγ

⎧ ⎡ ∂A pair (rij , T ) ⎤⎫ ⎪ ⎪ ⎥⎬ ⎨ ∑ αγ⎪∑ ⎢⎢ ⎥ ∂βγ ⎦⎪ ⎩ (i , j) ⎣ ⎭ γ=1

∑ αγ⎪⎨∑ ⎢

1 V9

=

where A is the total Helmholtz energy (eq 1) and m is the total mass in the system. It should be noted that density, ρ9 , refers to the reference, 9 , configuration of the system.  denotes the deformation gradient tensor defined through a mapping of an infinitesimal vector dx of the initial configuration onto the infinitesimal vector dx′ after the deformation. Equation C4 can be written in component form as

γ=1

1 V9

=

(C3)

⎛ ∂(A({rij}, {ρ }, T )/m) ⎞T cell ⎟⎟ τ = ρ9 ·⎜⎜ ∂ ⎝ ⎠

(C8)

3

1 V9

T) =

=

where V9 is the volume of the simulation box in the reference state 9 . In general the prefactor should be 1/V with V being the current volume of the system. The density, ρ, generally changes with  and V = V9 det(). In our model, the free energy per unit mass, A({rij}, {ρ(r)},T)/m, is a function of the separation vectors between all connected beads, {rij}, local densities, {ρ(r)} ≡ {ρcell}, and temperature, T. The thermodynamic stress tensor, τ, is given by eq C1:

∑ αγ

γ=1

(C2)

If we assume that our simulation box does not exchange mass with its environment, the stress tensor can further be simplified:

3

∑ αγ

is the nonbonded contribution to the stress tensor. We start by calculating the bonded contribution to the stress tensor, which depends only on the distances between connected pair of atoms, rij, and temperature, T:

where A is the Helmholtz energy, m is the total mass in the system, and  denotes the deformation gradient tensor. The same equation can be written in component form as

ταβ = ρ9

3

1 T) = V9

⎧ ⎡ ⎤⎫ ⎪ ∂ ⎪ acc ⎢ ∑ V cell, ⎥⎬ a ( , T ) ×⎨ ρ k vol cell, k ⎪ ⎦⎥⎪ ⎩ ∂βγ ⎢⎣ k ∈ cells ⎭

(C1)

∂(A /m) ταβ = ρ ∑ αγ ∂βγ γ=1

⎧ ⎫ ⎤⎪ ⎪ ∂ ⎡ ⎢ ⎥ ⎨ ⎬  A ( r , T ) ∑ αγ ⎪ ⎢∑ pair ij ⎥⎪ ⎦⎭ γ=1 ⎩ ∂βγ ⎣ (i , j) 3

1 V9

is the bonded contribution to the stress tensor, and

APPENDIX C. DERIVATION OF THE MODEL STRESS TENSOR We consider the free energy per unit mass, A(T,V)/m of the system. The thermodynamic stress tensor, τ, is103,104

γ=1

b nb ταβ ({rij}, T ) + ταβ ({ρcell }, T )

(B10)

■

∑ αγ

∂A nb({ρcell }, T ) ⎤ ⎥ ⎥⎦ ∂βγ

where

that connects the parameter w with the thermal wavelengths of the molecules and their intramolecular configurational integral, A0 . The molecular weight, M, enters eq B10 because ρ* refers to the close-packed mass density of the equation of state.

3

∂βγ

⎣

γ=1

+

(C6)

NAρ* ∏i Λi3 (r − 1) 1 e = w M A0

1 V9

m

⎡ ∂Ab({rij}, T )

3

∑ αγ⎢⎢

=

which results in

ταβ =

ρ9

(C5)

Invoking the functional dependence of our Helmholtz energy (eq 1): 3022


Article

Macromolecules

r′i = ri

(C11)

r′j = rj

(C12)

=

By subtracting eq C12 from eq C11, we get: +

r′ij = rij

1 V9

∑ αγ ∑

+

∑ βγrij ,γ

=

(C14)

γ=1

+

and the derivative appearing in eq C9 can now be calculated: ∂rij , β ∂βγ

1 V9 =

= rij , γ

(C15) −

Equation C9 takes the form b ταβ ({rij}, T ) =

=

1 V9

∑

∂A pair (rij , T ) rij , β ∂rij

(i , j)

1 V9

∑ (i , j)

rij

=

3

∑ rij ,γ αγ

γ=1

1 V9

⎡ ∂a vol(ρcell, k , T ) ⎤ acc ⎢Vcell, ⎥ k ⎢ ⎥⎦ ∂βγ k ∈ cells ⎣

3 γ=1

1 V9 1 V9

3

∑

a vol(ρcell, k , T ) ∑ αγ

k ∈ cells 3 k ∈ cells

∂βγ

γ=1


k ∈ cells 3

k ∈ cells


k ∈ cells

ρcell, k

∂βγ

∂ρcell, k 3

∑

∂βγ

∂a vol(ρcell, k , T ) ∂ρcell, k

γ=1

∑

acc ∂Vcell, k

γ=1

acc Vcell, k ∑ αγ

∑ 1 V9

∂a vol(ρcell, k , T ) 3

∑

∂βγ

γ=1 acc Vcell, k ∑ αγ

∑

acc ∂Vcell, k

acc ∂Vcell, k

γ=1

∂a vol(ρcell, k , T ) ∂ρcell, k

k ∈ cells

∂βγ

3

∑ αγ γ=1

acc ∂Vcell, k

∂βγ

⎡ ∂a vol(ρcell, k , T ) ⎤ 3 ∂V acc ⎢a vol(ρ , T ) − ρ ⎥ ∑ αγ cell, k cell, k cell, k ⎢ ⎥⎦ γ = 1 ∂ρcell, k ∂βγ k ∈ cells ⎣

∑

(C17)

γ=1

where avol(ρcell,k,T) is the nonbonded Helmholtz energy density in cell k. The term in brackets is the negative of the pressure, −p(ρcell,k,T), as that is predicted by the equation of state (eq B3) under given density ρcell,k and temperature, T. At this point we have to calculate the derivatives expressing the variation of the volume of a cell with respect to an element of the acc deformation gradient tensor, ∂V cell, k / ∂βγ . The determinant of the deformation gradient tensor is the ratio of volumes or densities of the deformed and initial configurations:

∂A pair (rij , T ) rij , βrij , α ∂rij

1 V9

⎡ ∂V acc ⎤ ⎢ cell, k a vol(ρcell, k , T )⎥ ⎢ ∂βγ ⎥⎦ k ∈ cells ⎣

∑ αγ ∑

3

rij , β′ =

γ=1

1 V9

1 = V9

the components of the undeformed through the relation:

⎤⎫ ⎪ ∂ ⎡⎢ acc ⎥⎬ V a T ρ ( , ) ∑ cell, vol k cell, k ⎪ ⎪ ⎢ ⎥ ∂  ⎣ ⎦ ⎩ βγ k ∈ cells ⎭

∑ αγ ⎨

3

(C13)

Thus, one component of the deformed vector is connected to

⎧ ⎪

3

1 V9

nb ταβ ({ρcell }, T ) =

rij (C16)

where the virial theorem is recovered.105 We now move to the estimation of the nonbonded contribution to the stress tensor, τnb αβ. We calculate the stress by taking the derivative of the equation-of-state free energy density, avol(ρcell,k,T), with respect to the deformation gradient

det() =

tensor,  , considering the set of space discretization voxels

acc ′ ρcell, k V cell, ρ V′ k = acc = = ′ k ρ′ ρcell, V9 V cell, k

(C18)

where the use of primes denotes the deformed configuration. The derivative of the determinant of  with respect to the tensor  itself is calculated by the following equation:106

(introduced in Appendix A) as a permanent scaffold into which beads would be reassigned upon deformation. In that approach, a tensile deformation in a periodic system would amount to

∂(det()) = det()(−1)T ∂

extending the system such that it occupies an additional layer of voxels, for example. This would indeed be infinitesimal for a

(C19)

Finally, the nonbonded contribution to the stress tensor, substituting the above terms in eq C17 is

very large model system. There are practical problems with this, however, so we consider our voxels as a means of spatial

nb ταβ ({ρcell }, T ) = −

discretization in integrating the free energy density over 3dimensional space, allowing them to deform affinely following

1 =− V9

the macroscopically applied deformation. This may not be too

p(ρcell, k , T ) ∑ αγ

k ∈ cells

γ=1

acc ∂V cell, k

∂βγ

acc ∂V cell, k ∂(det()) ∑ p(ρcell,k , T ) ∑ αγ ∂det() ∂βγ k ∈ cells γ=1

= −δαβ

nonbonded contribution to the stress tensor can be cast as

3

∑

3

bad if the voxel size exceeds the correlation length of density fluctuations in the polymer. Following our assumptions, the

1 V9

1 V9

∑

acc p(ρcell, k , T )V cell, k

k ∈ cells

(C20) 3023


Macromolecules

■

where all volumes, V9 and Vacc cell,k refer to the undeformed (reference) state of the system. Equation C20 could be fully anticipated. The contribution of the equation of state to the stress tensor of the system is limited to the diagonal components and its magnitude is the negative of the weighted average pressure over the cells of the grid, with volumes Vacc cell,k being the weights multiplying the individual contributions.

APPENDIX E. MICROSCOPIC REVERSIBILITY IN SLIP SPRING FORMATION AND DESTRUCTION (FLUCTUATING NUMBER OF SLIP SPRINGS SCHEME) Consider a specific configuration bearing a slip spring connecting an end b0 of a chain contained in it with a bead a0 of the same or another chain (which may be either an interior or an end bead). Let us call this configuration a0∧b0, for simplicity. Consider also a configuration which is identical to a0∧b0, the only difference being that the slip spring connecting a0 and b0 is missing. Let us call that configuration a0∧b× (see Figure 5). If a0 is not a chain end, the only way in which one can go from a0∧b0 to a0∧b× is slippage of the slip spring past the chain end b0. The probability per unit time for observing the transition is

■

APPENDIX D. INITIAL ESTIMATION OF THE HOPPING FREQUENCY FACTOR The mean-square displacement of a slip spring in time Δt is khopΔt (nKuhns/beadb2), where the mean-square displacement of the chain along the primitive path due to the slip springs is mean‐square displacement of slip springs nss/chain

⟨Δrcm 2⟩ =

k hopΔt

=

k hopΔt

=

nKuhns/beadb2

Nss

2N

8 a0 ∧ b0 → a0 ∧ b× = k hop, a0 ∧ b0Pa0 ∧ b0

nss/chain 2Nss n

N

= 2Dcmalongcontour Δt

beads

(D1)

8 a0 ∧ b× → a0 ∧ b0 = k form, a0 ∧ b× → a0 ∧ b0Pa0 ∧ b×

with Nbeads and Nss being the total number of beads and slip springs in the system, respectively, and n the number of chains. The factor 2 is included because a slip spring is attached to two chains. N is the chain length measured in number of Kuhn segments. The chain diffusivity by means of slip spring jumps is Dcmalong contour =

Nss

N

beads

8 a0 ∧ b0 → a0 ∧ b× = 8 a0 ∧ b× → a0 ∧ b0

kBT Nζ

k form, a0 ∧ b× → a0 ∧ b0 = k hop, a0 ∧ b0

(D2)

4kBT Nss 1 ζ Nbeads (nKuhns/beadb2)

(D3)

(D4)

The last expression provides an upper limit to the hopping rate, in the case the diffusion of the polymer chains was based only on their motion along their primitive paths. Following the discussion of the main text, an estimate of the average hopping rate, in the case of Gaussian slip springs, is ⟨k hop⟩ ≃ 38νhop

(D5)

which, in view of eq D4, allows us for an upper limit for νhop: νhop

Equation of State Based Slip Spring Model for Entangled Polymer

Recommend Documents