Coarse grained simulation approaches provide powerful tools for the prediction of the equilibrium properties of polymeric systems. Recent efforts have sought to develop coarse-graining strategies capable of predicting the non-equilibrium behavior of entangled polymeric materials. Slip-link and slip-spring models, in particular, have been shown to be capable of reproducing several key aspects of the linear response and rheology of polymer melts. In this work, we extend a previously proposed multi-chain slip-spring model in a way that correctly incorporates the effects of the fluctuating environment in which polymer segments are immersed. The model is used to obtain the equation of state associated with the slip-springs, and the results are compared to those of related numerical approaches and an approximate analytical expression. The model is also used to examine a polymer melt confined into a thin film, where an inhomogeneous distribution of polymer segments is observed, and the corresponding inhomogeneities associated with density fluctuations are reflected on the spatial slip-spring distribution.

A central challenge in development of coarse-grained modeling approaches for soft matter is that of capturing the wide spectrum of spatial and temporal scales associated with the underlying collective behavior.1–4 In general, coarse-graining reduces the number of degrees of freedom and finer details are only kept in an average sense. Effective interactions between surviving interacting sites become softer, and the interaction potential between two interacting sites, separated by a distance r, remains finite in the limit r0. Therefore, by ignoring the finer details regarding the local non-crossability between molecules, one can achieve faster equilibration of materials over mesoscopic length scales. There are, however, systems for which microscopic details related to the non-crossability between molecules have consequences that manifest themselves at long time scales. Examples include high molecular weight polymer melts and concentrated solutions. For such materials, topological constraints play a central role on their mechanical and rheological properties.3,5 Several phenomenological top-down approaches have therefore been proposed to reincorporate the effect of topological interactions on coarse grained polymer models. In this work, we focus on a particular class of models that offers the potential to describe the dynamics of entangled polymer melts subject to structural or compositional inhomogeneities. These models are generally referred to as slip-spring models.6–11 

In recent work, we have presented a simulation approach that is able to describe the linear and non-linear rheology of entangled polymer melts.6 This approach relies on a multi-chain representation of the macromolecular fluid, and effective fluctuating interactions, mediated by the so-called slip-springs (SS),7,10 between neighboring pairs of macromolecules are introduced to represent the topological effects that arise from the non-crossability constraint. In general, there is not a correspondence between a physical entanglement and a slip-spring, and the mapping between these two concepts should be elucidated by more detailed comparisons between atomistic and coarse-grained simulations. Depending on the implementation, the total number of slip-springs can be kept constant, or it can be allowed to fluctuate through a chemical potential that determines the average molecular weight between slip-springs. In both cases, our results have been found to be in quantitative agreement with experimental data and suggest that the proposed formalism may also be used to describe the dynamics of inhomogeneous systems, such as composites, block copolymers, or polymers near surfaces. In previous work, we also made a connection between our proposed multi-chain model and the existing single-chain mean-field formalisms6,12 and suggest that it may be possible to develop systematic coarse-graining formalisms that start with atomistic, multi-chain representations of a material and end up with predictive constitutive models for thermodynamics and rheology.

As encouraging as our past results and those of others have been, in efforts to compare the results of our multi-chain model with those of a single-chain theory, and to do so at a modest computational expense, several assumptions were made regarding the process of creation and destruction of slip-springs between a given pair of molecules. In particular, the local environment around a specific polymer segment was considered to be constant (i.e., non-fluctuating). That assumption is inherent to single-chain descriptions, but it need not be made in particle-based simulation approaches. Importantly, for materials that exhibit inhomogeneities (e.g., copolymers, composites, or polymers near a surface), it is necessary to include the effects of a fluctuating local environment in order to capture the equation of state of the material and its coupling to dynamics. With that goal in mind, in this work, we revisit our previous model for entangled polymers and incorporate the effects of the fluctuating environment in which a polymer segment is immersed. We use this new version to obtain the equation of state associated with the slip-springs and compare it to results from an analytical expression proposed by Chappa and co-workers.7 We then apply the model to explore the spatial behavior of slip-springs in inhomogeneous systems, where confinement and phase segregation may play an important role on the dynamics.

General details pertaining to the “Theoretically-Informed Entanglement Polymer Simulations” or TIEPOS approach have been presented in the literature. Here we provide a short summary of the model and methods, and highlight important differences between past implementations and the formalism employed here. Without loss of generality, we consider an AB copolymer melt consisting of Nc macromolecules in volume V at temperature T. The polymers are represented by the discretized Gaussian chain model, and each macromolecule is composed of Nα beads, α[1,Nc]. The system’s configuration is specified by the coordinates of all polymer beads {ri(s)} at different times, where ri(s) is the position of the sth bead in the ith chain. Intra-molecular interactions are given by HbkBT=32i=1Ncs=1Ni1[𝐫i(s+1)𝐫i(s)]2/b2, where b2 is the mean-squared bond length of an ideal chain and kB is Boltzmann’s constant. Generalizations to more complex situations, including blends of polymers with different bond lengths or chain stiffness, are straightforward. The inter-molecular interactions are expressed as a functional of the local densities ϕA,B(𝐫),13–17 

nbkBT=ρVd𝐫[χϕAϕB+κ2(ϕA+ϕB)2],
(1)

where ρ is the bead number density. The repulsion between unlike monomers, quantified by the Flory-Huggins parameter χ, is represented by the first term in the functional. A finite melt compressibility is prescribed by the second term, which corresponds to Helfand’s quadratic approximation,18 where the parameter κ is related to the inverse isothermal compressibility. By regularizing the scalar density fields, ϕα(𝐫), employing a density cloud w(r) associated with each interacting segment of the polymer chains, it is possible to transform the field-based description in Eq. (1) into a particle-based representation with a pairwise interaction potential, HnbkBT=i<jU(𝐫ij), where the sum runs over all pairs of polymer segments in the system. The interaction range, rint, should be large enough such that individual beads interact with a sufficiently large number of particles, in the simulations reported in this work, we have chosen rint0.12Re. These intra- and inter-molecular interactions are soft and fail to capture entangled dynamics. One must therefore resort to slip-springs (or other forms of effectively describing non-crossability constraints) to represent such topological interactions.

As mentioned above, entanglement effects are represented in our work by fluctuating interactions, mediated by slip-springs (SS). A SS consist of two anchoring units that move along different chain contours, i.e., self-entanglements are excluded in the present approach, and connect the respective polymer segments via a harmonic spring: usskBT=32[𝐑k𝐑nk]2/bss2, where Rk denotes the position of the kth anchoring unit and nk denotes its partner. The strength of the spring is quantified by bss, which is chosen here as bss = b. In the implementation adopted here, the SS can only hop between neighboring polymer beads, and at most one SS occupies a given segment. Here we note that previous implementations have also considered “dynamic” springs that evolve according to their own equation of motion;10,11,19 however, the short-time dynamics on the length scale below a coarse-grained segment has no particular physical significance.20 In the grand-canonical implementation considered here, the total number of SS fluctuates around an average controlled by a chemical potential μ. The average number of slip-springs per chain is Z. To initialize the SS configuration, for each SS pair, a polymer bead is randomly chosen and the first anchoring unit is attached to that segment. Then, a nearby segment is chosen to place the partner anchoring unit, in such a way as to generate a Boltzmann distribution of SS lengths.

The time evolution of the system is performed by relying on a hybrid scheme, where polymer configurations are updated via dynamic simulations over a short time interval, τss, and SS configurations are updated via Monte Carlo moves at the end of such a time interval. The slip-springs can be updated in a variety of ways; one could use a stochastic variable, τ^ss, with an arbitrary probability distribution Ψ(τ^ss), and use this function to generate a set of τ^ss’s and update slip-springs after these different waiting times. In the interest of simplicity, here we have opted for a scheme with Ψ(τ^ss)=δ(τ^ssτss). This choice provides good agreement with experiments. In this work, also for simplicity, we adopt a Brownian dynamics approach to evolve polymer segments. Thus, the time evolution of configuration {rj} is governed by

γd𝐫j(t)dt=𝐟jb(t)+𝐟jnb(t)+𝐟jss(t)+𝝃j(t),
(2)

where 𝐟jb and 𝐟jnb are the intra- and intermolecular forces acting on the polymer segment j, and 𝐟jss is the force associated with the slip-springs. Here 𝝃j represents the stochastic influence of the coarse-grained microscopic degrees of freedom and satisfies the fluctuation-dissipation relation: ξiα(t1)ξjβ(t2)=2γkBTδijδαβδ(t1t2), where α,β=x,y,z, and γ is the friction coefficient.

During the SS configuration update, one performs Nss+2Nc attempts to update SS where Nss(t) are the number of SS present at that time t, as well as to create new SS at chain ends (see the supplementary material for details). For each attempt, an existing SS is selected with probability PSS, or a chain-end is chosen randomly with probability Pch = 1 − PSS. In our implementation, we have chosen PSS=Z2+Z. Then, one of the following MC moves is attempted:

  1. A hopping move, which involves the transfer of one end of a slip-spring from one polymer segment to a neighboring one, along the same polymer. We propose to move in one or the other direction along the chain with equal probability. This move is accepted with probability p=min(1,eΔuss/kBT), where Δuss is the difference on the slip-spring energy between the trial and old configuration.

  2. A destruction move involves the deletion of an existing slip-spring pair if one of the rings reaches the chain end and moves out of the contour with probability 1/2. This move is accepted with probability p = min(1,PchPSS2NssNcβeW1). Here, W = euss(r1)/kBT+k=2Meuss(rk)/kBT, where r1 is the bond length associated with the slip-spring to be destroyed.

  3. A creation move creates a new slip-spring pair by adding a new slip-spring end to a randomly selected unoccupied chain end. The location of the ring’s partner is selected from all M unoccupied polymer beads within a radius Rc with probability euss(rj)/kBT/W, where W=k=1Meuss(rk)/kBT. This move is accepted with probability p = min(1,PSSPchNc2(Nss+1)βe1W).

where βe=eμ/kBT. Thus, the value of the chemical potential dictates the average number of slip-springs per chain. Note that, as in previous implementations, the main relaxation mechanisms of entangled polymers are incorporated naturally in this approach: reptation, contour length fluctuations, and constraint release.

The important point to highlight here is the following: in this multi-chain model, at each instant of time there is a fluctuating discrete number of possible neighbors to anchor the other end of a slip-spring. Also the microstates associated with those neighboring beads can change from one time to another. All these microscopic details are encoded in W, which is, therefore, a fluctuating quantity. Note that in the absence of fluctuations in the local environment, as is the case of single-chain models, then W is constant, let us say Wo, and could be adsorbed into βe, thus renormalizing the SS chemical potential, i.e., μsc=μ+kBTlnWo. That assumption was made in a previous report.6 As discussed below, while its consequences are minor for homogeneous systems, e.g., homopolymer melts, it could lead to incorrect predictions for inhomogeneous materials, including block copolymers, thin films, or solutions. The implementation of the model outlined above correctly incorporates such microscopic fluctuations, and thus, W is a fluctuating quantity, varying both in space and time.

In this section, we explore the effects of incorporating a fluctuating W into a multi-chain description of homopolymers, as indicated in Section II. Note that W depends on particle density ρ; the higher the number density, the smaller become the relative fluctuations of W. If W is assumed to be a constant, it plays no role on the acceptance probabilities and, therefore, the average number of slip-spring per chain Z is no longer sensitive to changes in particle density. In the acceptance rules outlined above, W appears explicitly, and number density therefore affects the resulting Z.

We fix the value of the chemical potential in both schemes and vary the polymer bead density (the results reported in this section use N = 32, and the length is measured in units of Re). It is important to remark that the chemical potential in both schemes does not have the same numerical value and does not depend on ρ either. Figure 1 shows Z as a function of density ρ. As expected, for constant W, the number of slip-springs per chain does not depend on the bead density for large enough values (green dashed curve). In contrast, for a fluctuating W, the new set of acceptance rules gives rise to a monotonic dependence on ρ. In fact, for relatively small densities, a linear dependence is found, which agrees with the prediction of Ref. 7 for a model similar to that considered here: This behavior can be explained by invoking an analytical prediction for the SS equation of state, originally derived in Ref. 7, which reads34:

Z=ρβeNc(NcN3)d𝐫g(𝐫)euss(𝐫)/kBT,
(3)

where g(𝐫) is the slip-spring pair correlation function, which only includes pairs of segments that can be connected by a slip-spring and βe=1/ζ=exp(μ/kBT) denotes the inverse fugacity of SS. Note, in particular, that if the melt’s microscopic structure did not depend on density, i.e., if g(r) were not affected by changes of ρ, then for a fixed chemical potential Z would linearly grow with ρ. However, if by changing ρ there is a change of the slip-spring pair correlation function, then a non-linear dependence should be expected.

FIG. 1.

Average number of entanglements per chain, Z, as a function of ρ. Green circles are obtained using the previous acceptance rules for the creation-destruction processes, while black circles are simulation results obtained by using the correct rules described on the text. The chemical potential differs in numerical value in both cases, but the rest of physical parameters are the same as those used in Ref. 1. The red curve is a linear fit to the small density data.

FIG. 1.

Average number of entanglements per chain, Z, as a function of ρ. Green circles are obtained using the previous acceptance rules for the creation-destruction processes, while black circles are simulation results obtained by using the correct rules described on the text. The chemical potential differs in numerical value in both cases, but the rest of physical parameters are the same as those used in Ref. 1. The red curve is a linear fit to the small density data.

Close modal

One key difference between the model in Ref. 7 and ours is that in the model of Chappa et al.7 a single segment can be occupied by any number of SS, only SS-connections to the identical bead or nearest neighbors along the same chain molecule were excluded. In our case, only one SS can be attached to a given polymer segment and the two ends of a SS have to reside on different polymers. We stress, however, here that at this point it is still unclear whether one or multiple SS per bead provide a better description of entangled polymers.

Interestingly, the curve Z vs ρ departs from a linear relationship for large densities. We computed the slip-spring pair correlation function for different particle densities. In Figure 2, we show g(r) for two extreme cases: ρ = 400/Re3 and ρ = 8000/Re3. As can be seen in the figure, the local microscopic structure has a strong dependence on the particle density, thereby affecting the number of SS per chain. Since an increase in density gives rise to a larger g(r) at short distances, the number of potential SS-partners of a segment increases faster than ρ, and Eq. (3) predicts an enhanced increase of Z with density, contrary to what is observed in Fig. 1. Thus, the deviations from the linear behavior in Figure 1 indicate a saturation in the number of SS due to the additional constraints of single-segment occupancy of the intermolecular SS. Also, we note again that the fact that we are excluding more than a single SS per segment induces effective contact interactions (and associated correlations) between SS that were not considered in the analysis leading to Eq. (3).

FIG. 2.

Slip-spring pair correlation function g(r) for two different number densities, as indicated on the labels.

FIG. 2.

Slip-spring pair correlation function g(r) for two different number densities, as indicated on the labels.

Close modal

If the number density is fixed at a given value (ρ = 1798/Re3), we can explore how Z is affected by the chemical potential. As reported in Figure 3, a monotonic dependence of the chemical potential is obtained. Again the growth of Z with SS fugacity ζexp(μ/kBT) agrees with that expected from Eq. (3). However, in this case, a linear relationship was reported in Ref. 7. The origin of this discrepancy between both approaches relies on one key difference: first, as mentioned above, in the approach considered here each polymer bead can be occupied by only one slip-spring and the two ends of a SS are connected to different polymers, whereas in Ref. 7 there are no such restrictions on the SS connectivity. At a fixed density, the incorporation of slip-springs has no effect on the microscopic structure of the melt, even for large values of the SS chemical potential (see Figure 4), i.e., for this density the compensating potential that eliminates the effect of SS on equilibrium properties7 is not necessary. Therefore, as the chemical potential is increased, the number of available unoccupied beads decreases, giving place to a slower growth of Z with μ. In the model introduced in Ref. 7, any number of SS on a single bead is permitted, and increasing ζexp(μ/kBT) leads to a linear increase of the number of SS that can be accommodated between any neighboring pairs of segments except nearest neighbors along a chain.

FIG. 3.

Average number of entanglements per chain, Z, as a function of exp(βμ). Results obtained by using the algorithm described in this work. Density is fixed at ρ=1798/Re3.

FIG. 3.

Average number of entanglements per chain, Z, as a function of exp(βμ). Results obtained by using the algorithm described in this work. Density is fixed at ρ=1798/Re3.

Close modal
FIG. 4.

Slip-spring pair correlation function g(r) for two different chemical potential values at ρ=1798/Re3, Z 2 and 14, respectively.

FIG. 4.

Slip-spring pair correlation function g(r) for two different chemical potential values at ρ=1798/Re3, Z 2 and 14, respectively.

Close modal

We conclude this section by comparing the analytical equation of state obtained for the multi-chain system (Eq. (3)) with that predicted by the single-chain theory of Schieber et al.21–24 In the latter theory, for discrete chains, the equation of state reads

Z=Nβesc+1,
(4)

where βesc=eμsc/kBT, the superscript “sc” refers to the single-chain prediction, and, βesc is chosen to be equal to the average number of beads between entanglements, Ne. Note that the theory of Schieber et al. excludes double SS occupancy by construction, whereas Eq. (3) captures the multi-chain correlations. Our model contains both ingredients, and the following considerations provide a reasonable approximation for dense systems that allows a comparison. For weakly entangled melts, the single-chain equation of state can be approximated by ZNβesc. On the other hand, ignoring that in Ref. 7 SS must not connect the same or nearest-neighbor beads along a chain, we can approximate Eq. (3) by

Z=NβeW(ρ),
(5)

where

W(ρ)=d𝐫ρg(𝐫)euss(𝐫)/kBT,
(6)

is the average of the Rosenbluth factor defined in the creation-destruction rules described above. In the case of homogeneous systems, W(ρ) is independent of space and time. Thus, for weakly entangled homopolymers in bulk, we can identify βesc = βeW(ρ)1, i.e., the chemical potentials are related by μsc = μ+kBTlnW(ρ). To corroborate this approximate relationship, we first plot the data in Figure 3 in terms of βe = exp(μ/kBT). This is shown in the main panel of Figure 5. As can be seen, a similar dependence of Z on βe is obtained (see Figure 2 in Ref. 1), but the latter has quite different numerical values. Then, we use the time average of W sampled during the simulations, W and plot Z vs μ+kBTlnW. The resulting curve is presented in the inset of Figure 5 along with the equation of state predicted by the single-chain theory, Eq. (4). As can be seen, the simulations follow the analytical predictions, even for systems with small βesc. Therefore, for homogeneous systems, the results presented in Ref. 1 provide a reasonable approximation to the results obtained using the fluctuating-environment creation-destruction rules adopted here. Note, however, that for inhomogeneous systems the spatiotemporal character of W(r,t;μ) will affect the creation-destruction rules locally, possibly leading to inhomogeneous slip-spring distributions, as expected, for example, in phase segregated block copolymers. A specific example is considered in Section III B.

FIG. 5.

Main panel: Average number of entanglements per chain, Z, as a function of βe. Inset: Average number of entanglements per chain, Z, as a function of an effective chemical potential. Line is the theoretical prediction and symbols are numerical simulation results. N = 32 and ρ=1798/Re3.

FIG. 5.

Main panel: Average number of entanglements per chain, Z, as a function of βe. Inset: Average number of entanglements per chain, Z, as a function of an effective chemical potential. Line is the theoretical prediction and symbols are numerical simulation results. N = 32 and ρ=1798/Re3.

Close modal

In what follows, we consider two inhomogeneous a homopolymer melt confined into a thin film, and a phase segregated block copolymer thin film. In both cases, instead of evolving the polymer configurations by dynamical equations (Eq. (2)), we have chosen to perform local Monte Carlo moves to relax polymer chains.28 The reason for this has to do with the fact that hard boundary conditions are easier to handle with an MC approach, and we are only interested in the polymer and slip-spring density distributions, as opposed to dynamic properties. As above, polymer segments are relaxed by a certain number of MC steps; the SS configuration is updated using the rules described in Section II. The coarse grained model we use requires that we use a large chain discretization, N, in order to describe the narrow interface of a nearly incompressible polymer melt in the vicinity of a solid.25 That is computationally demanding, and we have therefore opted here for an effective bonded interaction which helps to reproduce the configurational statistical properties of Gaussian chains under confinement with a low value of N.25 Note that an alternative technique, particularly useful when the statistical segment length is smaller than the polymer-solid interface width, would be to use an external potential obtained via Self-Consistent Field Theory (SCFT)26 or use a higher-order discretization of the Gaussian paths.27 Explicitly, we add the extra term to the bonded contribution25 

VwallkBT=ln(1e6z1z2/b2),
(7)

where z1 and z2 are the perpendicular coordinates, with respect to the hard surface, of two connected segments along the backbone. This equation is obtained by modeling the hard walls as infinite step functions and solving the diffusion equation with Dirichlet boundary conditions on the impenetrable surface (see Ref. 25 for details). This term compensates for the lost of accessible configurations close to the hard interface. The films have a thickness Lz and the hard walls are located at z = 0 and z = Lz. In the directions parallel to the surfaces, we impose periodic boundary conditions. Using these effective interactions, it is possible to recover the density profiles predicted by the field theory with N as low as N = 32. In Figure 6, we present the density profile for polymer segments (black curve) of a melt confined into a thin film of thickness Lz = 1.25Re. As can be seen in the figure, there is a depletion layer close to the surfaces, where most of the segments are chain-ends (red dashed curve), which are entropically favored to be localized near the walls. Thus, this breaking of translational symmetry induces changes to the polymer configurations and it becomes stronger as the thickness of the films gets smaller. As a consequence of the inhomogeneous distribution of polymer segments and chain configurations, the slip-springs also distribute inhomogeneously inside the film (green dashed curve); as can be appreciated in the figure, most of the entanglements are located in the middle of the film. We should note here that recent work29–31 addressing the effect of strong confinement on polymer configurations and entanglements in microscopic models has reported an entanglement spatial distribution that displays the same qualitative behavior as the one reported in here. Importantly, not only does the distribution of slip-springs become inhomogeneous, but their number decreases as the films become thinner (see Figure 7). This effect is consistent with experimental observations showing that entangled melts have a lower viscosity as they undergo nanoscale confinement,32,33 which is an indirect indication that chains become disentangled under such circumstances. As mentioned above, in our current approach, we have not considered self-entanglements, which are likely to arise in these thin film systems.

FIG. 6.

Density profiles of polymer thin films confined by hard, impenetrable solid walls. N = 32 and ρ=1798/Re3.

FIG. 6.

Density profiles of polymer thin films confined by hard, impenetrable solid walls. N = 32 and ρ=1798/Re3.

Close modal
FIG. 7.

Average number of entanglements per chain, Z for different film thickness, Lz. Symbols are numerical simulation results and the red curve is a linear fit to the data. N = 32 and ρ=1798/Re3.

FIG. 7.

Average number of entanglements per chain, Z for different film thickness, Lz. Symbols are numerical simulation results and the red curve is a linear fit to the data. N = 32 and ρ=1798/Re3.

Close modal

An analogous situation arises in lamellae-forming block copolymers. For these materials, we chose a discretization parameter N = 64, and the Flory-Huggins parameter is chosen as χN=37, with a lamellar periodicity of Lo = 1.79Re. The melt is confined into a thin film (thickness Lz = Lo, and we have also made use of Eq. (7)) and the bottom substrate is patterned with stripes of width w = Lo/2, and periodicity equal to Lo. The stripes are preferential for block A, and the space between the stripes is preferential for block B. Thus, the translational symmetry is broken in two dimensions: perpendicular to the substrate and perpendicular to the A-B interfaces. As expected, there is an inhomogeneous distribution of slip-springs in both directions. More importantly, for the parameters we have chosen, there is a minimum in the slip-spring spatial distribution located at the A-B interface, and a shallow minimum in the middle of the domains (Figure 8). Entanglements are in this case localized in the middle of the blocks, as opposed to the chain ends. The latter results are quite promising, as they suggest some correlations between molecular conformations and slip-springs, as expected to occur for real entanglements. We should note here that it could be possible that the local chain orientation would play a stronger role defining the local entanglement length; thus SS update rules could require a term associated with the local orientation. To address this issue, more comparative studies between microscopic models and the current approach are needed.

FIG. 8.

Density profiles of block copolymers and slip-springs perpendicular to the A-B interface. χN=37 and Lo = 1.79Re.

FIG. 8.

Density profiles of block copolymers and slip-springs perpendicular to the A-B interface. χN=37 and Lo = 1.79Re.

Close modal

A previously reported multi-chain slip-spring formalism has been extended to incorporate the effects of a fluctuating local number of entanglements. The simulation algorithm was used to obtain the equation of state associated with the slip-springs, and its results were found to be in good agreement with predictions from an analytical expression.7,21 The formalism was also used to examine the effects of confinement and phase separation on the entanglement distribution. In both cases, these inhomogeneities are clearly reflected on the spatial distribution of slip-springs and are expected to play a major role in the dynamic response of the material. Results for the linear behavior and non-linear rheology of these materials will be reported in a future publication.

See supplementary material for a detailed description of the algorithm to update SS configurations.

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. L.S. and M.M. thank the German Science Foundation for financial support (No. DFG Mu 1674/16-1). Marat Andreev gratefully acknowledges the support from NIST through a CHiMaD postdoctoral award. We gratefully acknowledge the computing resources provided on Blues, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.

1.
P. G.
de Gennes
,
J. Chem. Phys.
55
,
572
(
1971
).
2.
P. G.
de Gennes
,
Scaling Concepts in Polymer Physics
(
Cornell University Press
,
1986
).
3.
R. G.
Larson
,
The Structure and Rheology of Complex Fluids
(
Oxford
,
1999
).
4.
R. B.
Bird
,
R. C.
Armstrong
, and
O.
Hassager
,
Dynamics of Polymeric Liquids
, Vol. 1 (
John Wiley & Sons, Inc.
,
1987
).
5.
M.
Doi
and
S. F.
Edwards
,
The Theory of Polymer Dynamics
(
Oxford
,
1986
).
6.
A.
Ramírez-Hernández
,
B. L.
Peters
,
M.
Andreev
,
J. D.
Schieber
, and
J. J.
de Pablo
,
J. Chem. Phys.
143
,
243147
(
2015
).
7.
V.
Chappa
,
D.
Morse
,
A.
Zippelius
, and
M.
Müller
,
Phys. Rev. Lett.
109
,
148302
(
2012
).
8.
T.
Uneyama
and
Y.
Masubuchi
,
J. Chem. Phys.
137
,
154902
(
2012
).
9.
M.
Langeloth
,
Y.
Masubuchi
,
M. C.
Böhm
, and
F.
Müller-Plathe
,
J. Chem. Phys.
138
,
104907
(
2013
).
10.
A.
Ramírez-Hernández
,
M.
Müller
, and
J. J.
de Pablo
,
Soft Matter
9
,
2030
(
2013
).
11.
A.
Ramírez-Hernández
,
F. A.
Detcheverry
,
B. L.
Peters
 et al,
Macromolecules
46
,
6287
(
2013
).
12.
J. D.
Schieber
and
M.
Andreev
,
Annu. Rev. Chem. Biomol. Eng.
5
,
367
(
2014
).
13.
F. A.
Detcheverry
,
D. Q.
Pike
,
P. F.
Nealey
,
M.
Müller
, and
J. J.
de Pablo
,
Phys. Rev. Lett.
102
,
197801
(
2009
).
14.
D. Q.
Pike
,
F. A.
Detcheverry
,
M.
Müller
, and
J. J.
de Pablo
,
J. Chem. Phys.
131
,
084903
(
2009
).
15.
M.
Müller
and
K. Ch.
Daoulas
,
J. Chem. Phys.
128
,
024903
(
2008
).
16.
M.
Müller
,
J. Stat. Phys.
145
,
967
(
2011
).
17.
M.
Müller
and
J. J.
de Pablo
,
Annu. Rev. Mater. Res.
43
,
1
(
2013
).
18.
E.
Helfand
,
J. Chem. Phys.
62
,
999
(
1975
).
19.
A. E.
Likhtman
,
Macromolecules
38
,
6128
(
2005
).
20.
M.
Müller
and
K. Ch.
Daoulas
,
J. Chem. Phys.
129
,
164906
(
2008
).
21.
J. D.
Schieber
,
J. Chem. Phys.
118
,
5162
(
2003
).
22.
D. M.
Nair
and
J. D.
Schieber
,
Macromolecules
39
,
3386
(
2006
).
23.
R. N.
Khaliullin
and
J. D.
Schieber
,
Macromolecules
42
,
7504
(
2009
).
24.
M.
Andreev
and
J. D.
Schieber
,
Macromolecules
48
,
1606
(
2015
).
25.
A.
Ramírez-Hernández
,
F. A.
Detcheverry
, and
J. J.
de Pablo
,
J. Chem. Phys.
133
,
0
64905
(
2010
).
26.
M.
Müller
,
B.
Steinmüller
,
K. Ch.
Daoulas
,
A.
Ramírez-Hernández
, and
J. J.
de Pablo
,
Phys. Chem. Chem. Phys.
13
,
10491
(
2011
).
27.
M. H.
Müser
and
M.
Müller
,
J. Chem. Phys.
142
,
174105
(
2015
).
28.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulations
, 2nd ed. (
Academic Press
,
2002
).
29.
D. M.
Sussman
,
W.-S.
Tung
,
K. I.
Winey
,
K. S.
Schweizer
, and
R. A.
Riggleman
,
Macromolecules
47
,
6462
(
2014
).
30.
W.-S.
Tung
,
R. J.
Composto
,
R. A.
Riggleman
, and
K. I.
Winey
,
Macromolecules
48
,
2324
(
2015
).
31.
32.
H. D.
Rowland
,
W. P.
King
,
J. B.
Pethica
, and
G. L. W.
Cross
,
Science
322
,
720
(
2008
).
33.
J. M.
Rathfon
 et al,
Macromolecules
44
,
5436
(
2011
).
34.

Assuming that there is not a constraint on how many SS a single bead can hold, the physical interpretation of Eq. (3) is as follows: the probability of a given pair of beads, i and j, to have a SS is eμ/kBTuss(r)/kBT, where r is their separation. Now, the probability of having such separation is given by g(r)/V, thus the average number of SS between this pair is: Zij=eμ/kBTVd3rg(r)euss(r)/kBT. The total number of SS is obtained by summing over all possible pairs, and from that Eq. (3) follows (see Ref. 7 for details).

Supplementary Material