Communication: One size fits all: Equilibrating chemically different polymer liquids through universal long-wavelength description

Predicting properties of polymeric materials with computer simulations often require their description with microscopic detail. This is, however, challenging due to slow dynamics of entangled “spaghetti-like” polymers¹ and the need to address systems with dimensions significantly larger than the coil size (to avoid finite-size effects). To circumvent these difficulties, hierarchical strategies gradually equilibrating the material on different observation scales are attractive (see Ref. 2 and references therein). First, the crudest level is addressed with models, representing a large amount of microscopic degrees of freedom by a single effective particle. Details on shorter wavelengths are then gradually reinserted until all chemical details are recovered. Computations remain tractable since backmapping requires only local sampling.

Implicitly hierarchical approaches assume that the long-wavelength structure is correctly described by the crude model and is insensitive to chemical details. Scale-decoupling in polymers can be rigorously justified and linked to universality of the long-wavelength behavior, often described by generic laws adsorbing chemistry-specific details into few parameters.³ Since different polymers can be mapped to the same point of parameter space, classes of systems with identical long-wavelength properties can be defined. This allows simplifying significantly the hierarchical modeling of polymeric materials by creating “material-genomic” libraries of morphologies with the correct long-wavelength behavior of different classes. For an entire class, such morphology must be generated only once, preferably using the most simple microscopic representation. This system can be even a chemistry non-specific model (e.g., bead-spring) which nevertheless maps on the same point of parameter space. Atomistic representations of any class-member can be recovered reinserting chemical details into the common long-wavelength description with standard techniques.² 

Here, the implementation of this concept is demonstrated for amorphous homopolymer melts. Despite their simplicity, they present significant interest for basic polymer physics (e.g., as a framework for studying rheology) and industry. Polymer melts are systems where long strongly interdigitating molecules fill space in a random walk fashion. The number of molecules crossing the volume of a test chain in a melt of type⁴ γ is proportional to $N ̄ ( γ ) 1 / 2 = ρ ( γ ) R ( γ ) 3 / N ( γ )$⁠, where ρ_(γ), R_(γ), and N_(γ) are the monomer density, root mean-square end-to-end distance, and polymerization degree, respectively. The invariant degree of polymerization, $N ̄ ( γ )$⁠, presents a natural choice for arranging homopolymer melts into groups with common static long-wavelength properties. Indeed, mesoscale conformations^5,6 and liquid structure^3,7,8 are known to be universal functions of this quantity (see below). Generally for dynamical properties, non-universalities can appear even in the long-wavelength limit;⁹ however, there is evidence that the invariant degree of polymerization of subchains between two consecutive entanglements is similar for all polymers.^10–12 $N ̄ ( γ )$ plays a key role in more complex systems, e.g., it controls in a universal way thermodynamic properties in symmetric block-copolymers.^13,14

Let us consider different melts (indexed by α), having the same $N ̄ ( l )$ as the lth library melt (see Fig. 1). For each melt, the regime of universal long-wavelength behavior can be defined introducing a coarse-graining length scale, ΔL_(γ) (γ = α, l), in the spirit of renormalization group theories.^3,15 In principle, the observation scale separating chemistry-specific and universal behavior cannot be exactly defined, thus, ΔL_(γ) presents at this stage an arbitrary “large” scale. Introducing ΔL_(γ) renormalizes the microscopic coordinate space as $r ̃ =r/Δ L ( γ )$⁠. The shortest subchains that are resolved during coarse-graining contain a number of monomers, N_b(γ), such that their root mean-square end-to-end distance, R_b(γ), equals ΔL_(γ). ΔL_(γ) can be simply related to N_b(γ) within the Flory hypothesis (FH).^3,16 Linking monomers into long chains in combination with incompressibility reduces intermolecular monomer-monomer contacts on the scale of chain-size creating a “correlation hole.” FH proposes that intramolecular interactions tending to swell a chain are canceled by interactions with the other molecules, forming the “soft walls” of the correlation hole. Hence, the polymer conformations on the mesoscale are ideal random walks, so that $R b ( γ ) 2 = N b ( γ ) b e ( γ ) 2$⁠. Thus, $N b ( γ ) =Δ L ( γ ) 2 / b e ( γ ) 2$⁠. The chemistry-specific coefficient, $b e ( γ ) 2$⁠, is the squared effective bond length.¹ 

All melts are represented in the renormalized space by ensembles of chains of spheres (blobs) with the same average diameter $R ̃ b ( γ ) =1$⁠. Combining known theoretical results, we argue that under certain conditions such an ensemble constitutes a common long wavelength description for melts at given $N ̄ ( l )$⁠. In practice, stored (library) configurations will have a fixed number of chains, n_(l), and volume, V_(l). Hence, to identify these conditions, it is natural to consider the canonical ensemble. The α-type and the library melt will have the same conformations and structure in the renormalized space when they (a) are represented by the same amount of blobs, $N CG = N ̄ ( l ) / ρ ( γ ) 2 b e ( γ ) 4 Δ L ( γ ) 2$⁠,¹⁷ (b) contain an equal number of chains, n, and (c) have the same rescaled volume, $V ̃ = V ( γ ) /Δ L ( γ ) 3$ (in fact, since $N ̄ ( γ )$ is the same, the last condition follows from the other two).

When ideal random walks are coarse-grained into a blob representation, the statistics (a) of the distance, $l ̃$⁠, between the center-of-mass (COM) of two sequential blobs and (b) of the angle, θ, between two vectors joining the COM of a blob with the COMs of the preceding and the succeeding blob are universal in renormalized space.¹⁸ Since the distributions of internal coordinates $l ̃$ and θ are universal, the conformations of coarse-grained chains in the library and α-type melts with an identical amount of blobs will be the same. The average mean-square distance of two blobs with ranking numbers s₁ and s₂ along the chain contour is a representative quantifier of these conformations, given by

$〈 l ̃ 2 〉= 2 3$⁠, p = 〈cos(π − θ)〉 ≃ 0.22, and s = |s₂ − s₁|.

FH was corrected,^5,6 demonstrating that the volume spanned by a chain presents a hierarchy of nested correlation holes of all possible subchains. This leads to weak self-avoidance and deviations from the ideal random walk behavior. The power-law dependence is replaced by^5,6

c₁ and c₂ are constants which, in general, are chemistry-specific (though for the end-to-end distance of long chains, the leading order correction in Eq. (2) reduces¹⁴ to a universal form).

Since the melts have the same number of chains, it is straightforward to show that they have the same invariant degree of polymerization also on the level of subchains, where $N ̄ b ( γ ) 1 / 2 =n N CG / V ̃$⁠. Due to the same $N ̄ b ( γ )$⁠, the liquid structure on the scale of blobs will be the same. This conclusion can be traced back to de Gennes predicting³ that the depth of the correlation hole for entire chains scales as $∼ N ̄ ( γ ) − 1 / 2$⁠. Modeling studies highlighted the universality of the shape of the correlation hole on intermediate scales.⁷ Blob-packing has been considered within integral equation theory predicting for their intermolecular pair distribution function,⁸ 

where A_o = 6^3/2/2π² and X₀ is a universal function.

In summary, theoretical arguments suggest that configurations of homopolymer melts with the same invariant degree of polymerization can be obtained by backmapping a common blob-based representation. The latter is defined in spaces renormalized by the coarse-graining scale and originates from a library melt. In all practical applications, chain length is defined in terms of N_(α); thus, the matching library configuration is chosen such that $N ̄ ( l ) 1 / 2 = ρ ( α ) b e ( α ) 3 N ( α ) 1 / 2$⁠. The coarse-graining scales of the modeled and library melts are arbitrary, provided that (a) ΔL_(γ) are sufficiently large for the subchains to be approximated by ideal random walks and (b) $Δ L ( α ) = ( ρ ( l ) b e ( l ) 2 / ρ ( α ) b e ( α ) 2 ) Δ L ( l )$⁠. The last condition establishes the requirement that the number of the blobs be the same (equal to N_CG). Identifying $N ̄ ( l )$ and matching ΔL_(γ) requires chemistry-specific information, i.e., $b e ( γ ) 2$⁠. Thus, conventional microscopic simulations of small samples of melts with moderate chain-lengths are performed. $R b ( γ ) 2$ is calculated and fitted with the RHS of Eq. (2), to extract $b e ( γ ) 2$⁠, c₁, and c₂. These data provide also an estimate of ΔL_(γ) at which the ideal random walk approximation is accurate. To perform the reference simulations, several techniques are available, including configuration-assembly^19,20 and rebridging Monte Carlo algorithms.²¹ 

Prior to fine-graining, the microscopic length scale is recovered back-transforming (cf. Fig. 1) the coordinate space of the blob-based description as $r=Δ L ( α ) r ̃$⁠. In principle ΔL_(α) (equivalently N_b(α)) can be large, hampering direct reinsertion of microscopic details into the common blob-based representation. This difficulty can be circumvented²² through a hierarchical scheme where a low-resolution blob-based description undergoes a sequence of fine-graining steps. Each step increases resolution substituting every blob with a pair of smaller ones. On these intermediate scales, conformations and liquid structure are affected substantially by microscopic features, so that the blob-based description becomes chemistry-specific. Nevertheless, fine-graining is feasible²² considering the blobs as soft spheres²³ with simple interactions parameterized from the reference microscopic simulations. Alternatively, somewhat more complex chemistry-specific blob-based models can be developed from integral equation theories,^8,24,25 reducing in the future the need for such calibration data. Once the blobs become sufficiently small, microscopic details can be introduced. The general strategy² is to reinsert the underlying subchains so that they comply, e.g., with the size and location of the COM of the blobs, and relax the reintroduced degrees of freedom.

As a proof-of-principle, we consider here homopolymer melts described with a generic microscopic model.²⁶ Linear chains are represented with monomers linked by finitely extensible nonlinear elastic (FENE) springs augmented by an angular potential¹⁹ U(θ) = κ_θ(γ)(1 − cosθ_o), where θ_o is the angle between two sequential springs. Nonbonded interactions are captured through a Weeks-Chandler-Andersen (WCA) potential. For FENE and WCA interactions, the standard parameterization²⁶ is employed. We mimic chemical diversity by varying chain stiffness and/or monomer number density in the range 0 ≤ κ_θ(γ) ≤ 1.5 and 0.60 ≤ ρ_(γ) ≤ 0.85, respectively. All lengths and energies are expressed in units of σ (the characteristic WCA length) and thermal energy, k_BT.

For the considered κ_θ(γ) and ρ_(γ), the required microscopic reference samples contained 400–700 chains with N_(γ) = 500, at most. These moderate sized systems were efficiently equilibrated through a configuration-assembly procedure^19,20 implemented within the ESPResSo+ + package.²⁷ From the reference configurations, $b e ( γ ) 2$ is extracted as previously described (see Figure S1 and Tables S1 and S2 in the supplementary material²⁸) and plotted as a function of κ_θ(γ) in the inset of Fig. 2. The main panel of Fig. 2 presents $R b ( γ ) 2 / N b ( γ ) b e ( γ ) 2$ as a function of N_b(γ) for κ_θ(γ) = 0 and 1.5. Fig. 2 demonstrates that for all considered κ_θ(γ), the deviations from the ideal random walk limit drop below 2% roughly at the same threshold $N b ( γ ) th =100$⁠. For practical purposes, for such deviations, we accept the ideal random walk statistics as a valid approximation. Thus, $Δ L ( γ ) th = b e ( γ ) N b ( γ ) th$ is taken as the smallest possible coarse-graining scale.

All library melts represent the same “chemical substance” described by κ_θ(l) = 1.5 and ρ_(l) = 0.85. Creating the library with stiffer homopolymers is advantageous since less monomers per chain are required to achieve the desired $N ̄ ( l )$⁠, as follows from $N ( γ ) = N ̄ ( γ ) / ρ ( γ ) 2 b e ( γ ) 6$⁠. The library samples cover a broad range of invariant degrees of polymerization $7.5×1 0 3 ≤ N ̄ ( l ) ≤15×1 0 3$ which are representative of values in experiments. The samples were generated employing the hierarchical scheme²² mentioned before. Since the universal long-wavelength description is not available at this stage, the melt configuration for fine-graining is obtained²² from direct Monte Carlo equilibration of a chemistry-specific blob-based model. The hierarchical strategy allows the preparation of large library melts, e.g., for $N ̄ ( l ) =7.5×1 0 3$ and 15 × 10³, the samples contain n_(l) = 1000 chains with N_(l) = 500 and 1000 monomers, respectively.

In the library melts, the universal description is reproduced with high accuracy already on the threshold scale, $Δ L ( l ) th ≃16.5$⁠. As an illustration, Fig. 3(a) presents for two melts with $N ̄ ( l ) =7.5×1 0 3$ and 15 × 10³ the internal distance plot $R ̃ CG 2 ( s ) /s$⁠. Within error bars, the plots (solid and open circles) follow each other, reproducing closely the universal behavior of Eq. (1) (dashed line). The magnitude of deviations is quantified²⁸ in Figure S2. The blob-blob intermolecular pair distribution function is presented in Fig. 3(b) (solid and open circles) and is well described by the universal form in Eq. (3).

Configurations of all other systems can be obtained by fine-graining the blob-based representation characterizing a library melt with the same $N ̄ ( l )$⁠, on the scale ΔL_(l) = 16.5. As an example, we discuss two chemically different melts with κ_θ(α) = 0.75, ρ_(α) = 0.60 (melt I) and κ_θ(α) = 0, ρ_(α) = 0.85 (melt II). From the definition of invariant degree of polymerization and the data on $b e ( α ) 2$ from Fig. 2, it follows that both melts for N_(melt I) = 1500 and N_(melt II) = 2000 map on $N ̄ ( l 1 ) =7.5×1 0 3$⁠. For N_(melt I) = 3000 and N_(melt II) = 4000, both correspond to $N ̄ ( l 2 ) =15×1 0 3$⁠. Because of scale matching, the universal representation corresponds to ΔL_(melt I) = 26.9 and ΔL_(melt II) = 26.4, representing an underlying number of monomers N_b(melt I) = 300 and N_b(melt II) = 400, respectively. The universal description is hierarchically fine-grained²² (cf. Fig. 1) until reaching blob-based descriptions corresponding to N_b(melt I) = 75 and N_b(melt II) = 50. Microscopic details are introduced into these small blobs, employing constrained reinsertion of underlying subchains and gradual elimination of overlaps between monomers.^19,20

To demonstrate equilibration of melts obtained from this procedure, Fig. 3(c) presents their internal distance plots, $R ( α ) 2 ( s ) /s$⁠, calculated on microscopic basis. s denotes the difference of the ranking numbers of monomers along chain contour and $R ( α ) 2 ( s )$ is their mean-square end-to-end distance. In all cases, $R ( α ) 2 ( s ) /s$ in the backmapped melts follow closely their counterparts in independent reference simulations (they deviate at most by 2%, see Ref. 28 Figure S3). These observations strongly verify our approach, since internal distance plots are extremely sensitive¹⁹ to the quality of equilibration. Fig. 3(d) compares the total monomer-monomer pair distribution functions, g_(α)tot(r), in the backmapped and reference melts, demonstrating that the local structure is correctly reproduced. It is instructive to convert the equilibrated melts I and II to blob-based representations in the renormalized space, corresponding to the initial coarse-graining scales ΔL_(melt I) = 26.9 and ΔL_(melt II) = 26.4. The calculated $R ̃ CG 2 ( s ) /s$ and $g CG ( r ̃ )$ are plotted in Figs. 3(a) and 3(b) and reproduce the universal behavior of the library melts.

These results confirm that a single melt can serve as a blueprint for configurations of different homopolymers described with microscopic resolution. Describing all systems through bead-spring models does not compromise the applicability of the approach to real materials. It has been already demonstrated that such models offer microscopic (albeit coarse-grained) description of real polymers when appropriately parameterized.^29–31 Moreover, the atomistic description can be recovered by backmapping the bead-spring representation through standard schemes.^29,31 So far, large samples with experimentally relevant invariant degrees of polymerization, as those in the current work, have been generated only on drastically coarse-grained level.³² Here, these systems are described microscopically with minimum computational costs. For example, equilibrating a sample of melt II with N_(melt II) = 4000 at $N ̄ ( l ) =15×1 0 3$ (4 × 10⁶ monomers in total) required only six days on 32 processors (3.0 GHz). The concept was elaborated and verified for melts of linear chains but should be straightforward to extend to certain cases of non-linear molecules, e.g., branched or concatenated ring-polymers (named³ “Olympic gels”). In contrast, the implementation to nonconcatenated ring-polymers is challenging since coarse-graining into blobs relaxes microscopic non-crossability of subchains. It is important to explore whether similar strategies are applicable to materials with complex long-wavelength structure, depending on more control parameters,^13,14 e.g., microphase-separated block-copolymers.

It is a pleasure to thank Livia Moreira for fruitful discussions and Carlos Marques for carefully reading our manuscript. The computing time granted by the John von Neumann Institute for Computing (NIC) on the supercomputer JUROPA at Jülich Supercomputing Centre (JSC) is gratefully acknowledged.