A method for calculating the nonequilibrium entropy of a flowing polymer melt via atomistic simulation

Nonequilibrium thermodynamics of polymeric liquids is a demanding practical and philosophical challenge for scientists and engineers. Unlike small molecular isotropic liquids, the macromolecular configurational microstructure of long chain-like molecules allows for immense changes in the liquid’s entropy under nonequilibrium conditions, such as an imposed flow field, in which the configurational microstate of the fluid can change from tightly coiled to highly stretched; even crystalline states can be induced under flow.^1–5 Powerful theoretical frameworks of nonequilibrium thermodynamics of polymeric systems have been developed;^6–12 however, these rely on physically appropriate expressions for the free energy and entropy of the nonequilibrium liquid. Unfortunately, there is no way to determine the nonequilibrium entropy changes experimentally. On the other hand, one might expect that entropy, a purely statistical property of the polymeric liquid based on counting the occupancy of configurational microstates at various levels of coarse-graining of the constituent macromolecules, could be readily calculated via atomistic simulation of nonequilibrium flow systems, similar to the internal energy. This expectation is not justified in practice, however, because of two major problems: (1) a theory is necessary (a) to relate atomic positions to configurational microstates at the chosen level of coarse-graining and (b) to relate these microstates to entropy, and (2) immense computational resources are required to perform the intensive integrations involved in determining the three-dimensional probability distribution functions that relate the configurational mesostates to the entropy.

Problem 1(a) boils down to the determination of the proper variable(s) necessary to quantify the relevant configurational information for the application at hand within the observational window of the experiment. Polymeric liquids can be described at various levels of coarse-graining based on the time and length scale of the observational window; at each level, the relevant variable represents the state space of configurational mesostates, and each point in that state space represents a distinct mesostate of the system. Depending on the level of coarse-graining, each mesotate can correspond to many configurational microstates of the constituent polymer chains whose internal configurations, although degenerate, are all compatible with the assigned mesostate. Hence, different choices for the relevant variables may result in different values for the entropy; this is because more coarse-grained variables have fewer available configurations (mesostates) than less coarse-grained variables, and therefore, each mesotate (of a more coarse-grained variable) has a higher degree of degeneracy of associated configurational microstates.

Problem 1(b) was addressed by Wall 80 years ago in his classical (3-page!) treatise on statistical thermodynamics of rubbers.¹³ Starting with Boltzmann’s entropy expression, $S=kBlnΩ$⁠, where $Ω$ is the number of distinct configurational mesostates and $kB$ is Boltzmann’s constant, Wall derived an exact expression for the change in entropy under nonequilibrium conditions as

where $ni$ is the number of macromolecules occupying configurational mesostate i and $neq,i$ is the same but at equilibrium conditions, echoing Gibbs’ entropy expression. This equation would solve the problem of calculating nonequilibrium entropy, at least philosophically if not practically, were it not for problem 1(a). Problem 1(a) arises because in order to determine $ni$ and $neq,i$ one must first have an idea of what is being counted, and as it turns out, there is not a unique way to quantify the configurational mesostates of a macromolecular liquid in terms of microstructural variables. Nevertheless, assuming that the ensemble is large enough, Eq. (1) will provide a proper entropy because counting the number of molecules occupying each mesostate would cover essentially all of the microstates available to that particular mesostate. Different variables (associated with different degrees of coarse-graining) might be relevant for different purposes, and hence, it is necessary to understand the relevant observation window for any particular application. For example, to calculate the stress state of the liquid under flow conditions, at large length and time scales, the configurational state of the entire chain is relevant to the entropy calculation, whereas to calculate the entropy of flow-induced crystallization at high strain rates, a more refined level of coarse-graining would be necessary to count the configurational mesostates. In the former case, the end-to-end vector might be chosen as the coarse-grained variable to classify the configurational microstates of the polymer chains into countable mesostates, whereas in the latter case, a more highly detailed counting would be required with highly nondegenerate mesostates comprised of relatively few microstates.

At present, consider a situation in which the configurational mesostate of the fluid is described by the end-to-end vector (⁠$R$⁠) of each polymer molecule. Each mesostate at this level of coarse-graining is highly degenerate, being composed of an ensemble of many possible internal configurations (microstates) that so happen to have the same chain endpoints. The various populated configurational mesostates accessible to this quantity are labeled $Ri$⁠, and hence, $ni=Nψ(Ri)d3Ri$⁠, where $N$ is the total number of macromolecules and $ψ$ is the associated probability distribution function, with a similar expression for $neq,i$⁠. Equation (1) can then be expressed in terms of the probability distributions associated with the various configurational states represented by the $Ri$ as

Now that the first issue has been resolved [assuming that one has made a proper choice for problem 1(a)], the calculation of entropy reduces to the determination of $ψ$⁠. However, as alluded to above, there are no reproducible experimental methods to accomplish this. Realistic atomistic simulations of polymeric solutions and melts offer an alternative, virtual platform that provides detailed information regarding the time evolution of phase space consisting of positions and momenta of each constituent atomic particle. Once the variables of problem 1(a) have been defined, the atomic positions can be used to determine the instantaneous values of the relevant variables, such as $R$ for each individual macromolecule, and then to compute their 3d probability distribution functions (PDFs), such as $ψ$⁠. This is when problem 2 arises: these simulations of entangled polymeric liquids with relevant time and lengths scales spanning many orders of magnitude require immense computational resources to perform. This is true even if one is not interested in calculating the entropy, but the problem becomes even more acute when one wishes to evaluate the relevant 3d PDFs required for application of Eq. (2). These calculations require an accurate methodology for extrapolating limited amounts of data from a large number of subcells within the simulation box, as described below. One of the goals of this Communication is to introduce such a methodology and apply it to steady-state shear and elongational flow of a polyethylene (PE) liquid.

The simulated fluid is a monodisperse, linear, entangled C₁₀₀₀H₂₀₀₂ PE melt undergoing steady-state planar Couette (shear) flow (PCF) and planar extensional flow (PEF) over a wide range of deformation rates. Specifically, the PCF simulations were performed in the range $0<Wi≤1170$⁠, where the Weissenberg number, $Wi$⁠, is defined as $Wi≡γ̇τd$ with $γ̇$ being the shear rate and $τd=5270$ ns being the disengagement time of the polymer under quiescent conditions.¹⁴ The PEF simulations were performed in the range $0<De≤15$⁠, where $De≡ε̇τR$ is the Deborah number, $ε̇$ is the extension rate, and $τR=137$ ns is the Rouse relaxation time of the polymer at equilibrium.¹⁴ Nonequilibrium molecular dynamics (NEMD) simulations were performed in the NVT statistical ensemble at 450 K and the experimental density of 0.766 g/cm³. The modified Siepmann–Karaborni–Smit (SKS) potential model^15,16 was used to quantify the energetic interactions between the atomistic constituents of the polyethylene macromolecules, which were viewed as united atoms composed of a carbon atom and its associated bound hydrogen atoms.¹⁵ Full information about the simulation methodology and details including force field parameters, simulation box sizes, and the number of particles are provided in the supplementary material.

Computation of the entropy change from NEMD simulation via Eq. (2) was performed by evaluating the integrand, $I=ψlnψeq/ψ$⁠, which requires the calculation of 3d equilibrium, $ψeq$⁠, and nonequilibrium, $ψ$⁠, PDFs of mesostates based on chain end-to-end vectors. Density functions were calculated by binning the ensemble of end-to-end vectors consisting of all chains at steady-state and then normalizing the resultant histograms—see the supplementary material for details regarding the choice of a suitable number of bins. Note that the integration interval is theoretically $(−∞,+∞)$ for all components of $R$⁠; however, it is practically performed in the range $[−Rmax,Rmax]$⁠, where $Rmax≡‖R‖max$ is the theoretical maximum chain end-to-end distance, i.e., a chain with all dihedral angles occupying trans conformations. Even within this narrow range, the probability densities for many values of $R$ are negligibly small. In particular, $ψeq$ and $ψ$ at low and intermediate deformation rates effectively go to zero as the components of $R$ approach $Rmax$⁠. This situation is aggravated further when the PDFs are calculated from molecular dynamics simulations because of the limited accuracy due to the restricted ensemble available. Such zero probabilities could potentially be problematic when evaluating the integrand of Eq. (2), as they appear as arguments of a logarithm function. These numerical problems can be partially avoided using algebraic manipulations of the integrand $I$⁠,

The problem of $ln(ψ)|ψ=0$ could be handled by using L’ Hôpital’s rule to show that $limψ→0ψ⁡ln(ψ)=0$⁠. Hence, the second term of Eq. (3) vanishes for $ψ→0$⁠; however, the same does not apply to the first term. Specifically, $limψeq→0ψ⁡ln(ψeq)$ is undefined. Therefore, zero values of the equilibrium PDF, $ψeq$⁠, should be avoided altogether.

As discussed above, zero values of $ψeq$ are inevitable if one solely relies on values calculated from simulation. However, the equilibrium distribution configuration of flexible polymeric liquids can be reliably assumed to be Gaussian. Indeed, the C₁₀₀₀H₂₀₀₂ liquid simulated in the present study was shown to exhibit Gaussian PDFs under quiescent and weak flow conditions.^14,17,18 In Eq. (3), $ψeq$ is the PDF for the chain end-to-end vector (not its magnitude), and hence, it is a multivariate distribution, with $x$⁠, $y$⁠, and $z$ components of the end-to-end vector as independent variables. $ψeq$ can be estimated very well with a Gaussian distribution function expressed as

where $μT=(⟨Rx⟩,⟨Ry⟩,⟨Rz⟩)$ is the mean vector composed of ensemble averages of components of the end-to-end vector $Rα$ and $Σαβ=Cov(Rα,Rβ)$ is the covariance matrix. Note that this estimate is only necessary when $ψeq$ happens to be zero in particular bins on the extreme edges of the variable range, as discussed above. $ψeq$ obtained from Eq. (4) can become indefinitely close to zero but never exactly zero; therefore, this equation mitigates the numerical error of $ψ⁡ln(ψeq)$ when $ψeq=0$ within certain bins of the simulation. See the supplementary material for additional remarks regarding the Gaussian PDF, Eq. (4), under equilibrium conditions. After the evaluation of $I$⁠, the triple integral of Eq. (2) can be computed via numerical integration. Trapezoidal rule was employed using the same number of sample points as the number of bins of the PDFs used in the calculation of $I$ in each direction.

The 3d distributions required to calculate entropy via Eq. (2) are difficult to display graphically. To illustrate conceptually the large variations in these distributions with respect to type of flow and flow strength, 1d plots of the PDFs of the magnitude of the end-to-end vector of the simulated PE melt under both steady-state PCF and PEF are displayed in Fig. 1 at multiple values of dimensionless strain rate. It is apparent that the Gaussian distribution of the end-to-end vector at equilibrium remains approximately Gaussian only for very low strain rates in both types of flow. In PEF, the chains initially extend slightly with increasing strain rate (⁠$De<0.3$⁠) until a critical $De$ (⁠$≈0.3$⁠) is attained where a coil–stretch transition occurs,¹⁹ which induces a configurational microphase separation composed of sheet-like regions of highly stretched chains enveloping ellipsoidal domains of tightly coiled macromolecules (⁠$0.3≤De≤1.5$⁠) that produces bimodal PDFs.²⁰ At high strain rates (⁠$De≥3$⁠), the chains are almost exclusively highly stretched, providing tall, narrow distributions approaching the fully extended chain length of 1290 Å, with flow-induced crystallization eventually occurring at $De≥15$ as induced by random nucleation events.^4,21 In PCF, the approximately Gaussian distributions at low values of strain rate broaden significantly as $Wi$ increases, ultimately becoming essentially bimodal at high strain rates (⁠$Wi>1000$⁠). This broadening of the PDFs is caused by the semiperiodic rotational cycles of the macromolecules induced by the vorticity of the imposed shear field after the entanglement network has begun to degrade (⁠$Wi>10$⁠).^14,18 Under both types of flow, plots of the PDFs of the magnitude of the end-to-end vector reveal wide variations in the shape of the underlying nonequilibrium configurational distribution functions with respect to the type of flow and the flow strength that ultimately determine the entropy of the liquid. Any a priori notion of applying Gaussian chain statistics is eroded at all but the lowest strain rates lying within the linear viscoelastic flow regime.

The entropy change, $ΔS$ (relative to the quiescent state), as calculated from the PDFs discussed above using Eq. (2), is displayed in Fig. 2 for both PCF and PEF (blue data points). The entropy decreases monotonically with increasing flow strength in both types of flow, as expected, since the quiescent state represents the most probable ensemble distribution and hence the highest entropy state. The decrease is mild at low strain rates within the linear viscoelastic regime, where the individual macromolecules are still relatively coiled and their PDFs approximately Gaussian. However, at intermediate strain rates under both types of flow, the entropy drops substantially, particularly under PEF where the distributions become highly peaked at high chain extensions, which implies a large reduction in the number of available configurations that the chain molecules can inhabit. The broad PDFs of PCF, on the other hand, result in less significant entropy changes because the wide range of available chain configurations results in higher entropy.

The entropy as calculated above is exact for the choice of variable used to describe the system, within applicable statistical limitations; however, from a practical perspective, it is very computationally intensive to evaluate. First, the simulations must be highly detailed (in the present case, 180 000 atoms) and cover both very small (fs) and relatively long (ms) time scales, ultimately achieving steady-state flow conditions and maintaining these until replicable ensemble averages can be obtained. The computational requirements necessary to perform the required simulations are therefore enormous by present day’s standards. Consequently, it would be very beneficial to calculate the entropy without having to introduce additional details into the simulation, either by increasing the size of the simulation cell or by increasing the run time. Unfortunately, calculation of the 3d PDFs necessary to evaluate the entropy via Eq. (2), as described above, press the extreme limits of contemporary computational resources.

It is common to express distribution functions in terms of their moments, such as the second-order dimensionless conformation (second moment) tensor, $C=⟨3RR⟩/Req2$⁠, where $Req2$ is the mean-square magnitude of the end-to-end vector in the unperturbed state and the angular brackets denote the ensemble average. This tensor then provides an ensemble-average representation of the distribution function from which it was derived. Although there is necessarily a loss of information during averaging, generally some defining characteristics of the underlying PDFs are represented. For example, the trace and determinant of the conformation tensor describe the mean-square extension of the macromolecules and directional shape of their spatial distribution. Despite the requisite loss of information, the second moment has the advantage of requiring far less statistical detail to compute than the 3d PDFs discussed above.

The nonzero components of the symmetric conformation tensor under both PCF and PEF are displayed in Fig. 3 as functions of strain rate. In both flows, $trC$ is the sum of the diagonal components of $C$⁠, each representing the mean-square extension of the macromolecules in the respective direction. In PCF, an additional independent nonzero component, $C12$⁠, represents the angular orientation of the macromolecules induced by shearing forces in the flow (⁠$x1$⁠)-gradient (⁠$x2$⁠) plane (shown in the inset). The average molecular extension increases gradually at low strain rates in both PCF and PEF and then significantly at intermediate rates before plateauing, with the PEF extension substantially higher than under PCF at high rates. In PCF, the $C12$ component increases linearly at low strain rate (see the inset: the plot is on a semi-log scale) up to $Wi≈10$ where the rotational cycles begin to develop and then attains a maximum and decreases thereafter as the distribution widens substantially and the molecules tend to orient along the flow direction.

Booij derived an entropy expression based on the dynamics of a dilute polymer solution, as quantified by the Rouse model.²² This expression was based on a derivation of van Wiechen and Booij, which demonstrated that the Rouse model could be re-expressed exactly in terms of second moments of the nonequilibrium Gaussian distribution functions developed under flows in the linear viscoelastic regime.²³ These second moments, one for each relaxational mode of the polymer macromolecule, followed the independent evolution equations corresponding to individual upper-convected Maxwell modes, as discussed in the supplementary material. Limited to the longest relaxation mode corresponding to the end-to-end vector of the molecule, the predictions of $trC$ and $C12$ (for PCF only) are also shown in Fig. 3. The average extension of the molecules rapidly diverges as the system moves out of the linear viscoelastic flow regime, in obvious disagreement with the simulated values. Furthermore, in PCF, the $C12$ component increases linearly with $Wi$ indefinitely, whereas the simulated value attains a maximum near $Wi≈10$⁠. Booij’s expression for the entropy, based on second moments of the underlying Gaussian PDFs (herein limited to the longest relaxation mode), is given by

A fundamentally similar expression was also derived by Sarti and Marrucci for the same dilute solution model.²⁴ The predictions for the entropy provided by this expression are also displayed in Fig. 2. As one would have expected, these entropies diverge rapidly outside of the linear viscoelastic regime.

It might appear pointless to discuss Booij’s entropy expression, given that the dilute solution for which it was originally developed is far removed from the entangled polyethylene melt simulated herein. However, the form of this expression turns out to be very useful if the second moments of the Gaussian distributions appearing in Eq. (5) are replaced with the values of the second moments of the C₁₀₀₀H₂₀₀₂ liquid taken directly from the NEMD simulations. The entropies computed using this method, $ΔSsm$⁠, are displayed in Fig. 2 as calculated from an expression similar to Booij’s entropy,

where $Csim$ is the dimensionless conformation tensor computed directly from the simulations corresponding to the second moment tensor of the 3d distribution function of $R$⁠, as displayed in Fig. 3. As evident, this entropy is generally statistically indistinguishable from the exact, fully 3d calculation according to Eq. (2) under both PCF and PEF at all strain rates. There appears to be a small systematic error at high strain rates, but this drift is within the uncertainty associated with the 3d calculation according to Eq. (2).

The implications of the preceding results are quite remarkable. Although developed for the Gaussian distributions within the linear viscoelastic flow regime, from a physical perspective Eq. (5) appears to give a very good approximation of the exact system entropy, as computed using Eq. (2). This is very surprising since the distributions developed in both types of flows are highly non-Gaussian (see Fig. 1) and their second moments are effectively averaged over a wide range of dynamical phenomena, including flow-induced phase separation (PEF: $0.3≤De≤1.5$⁠), flow-induced crystallization (PEF: $De≥15$⁠), and flow-induced rotation (PCF: $Wi≥20$⁠).^4,14,18–20 It is also remarkable that higher moments of the distributions are not necessary to quantify the entropy, even though the distributions are highly non-Gaussian and at high flow rates can be far from the linear viscoelastic flow regime from which Eq. (5) was originally derived.^22,23 The practical implication is just as remarkable since the second moments required to attain this approximation are drastically simpler to calculate than the fully 3d PDFs required for Eq. (2), and they require much less computational wherewithal to compute. Indeed, these second moments are readily obtained, with substantial statistical certainty, directly from the NEMD simulations. Consequently, use of the approximate expression in terms of second moments greatly increases the practical utility of computing the entropy in highly nonlinear, high strain rate flows of polymeric liquids.

The discussion above was limited to a description of the polymer molecules based on the end-to-end vector. Results for other possible segmental discretizations of the overall macromolecule are presented and discussed in the supplementary material.

See the supplementary material for additional methodology details and a discussion of the alternative chain discretizations.

This research is part of the Frontera Computing Project at the Texas Advanced Computing Center (TACC) at the University of Texas at Austin. Frontera is made possible by National Science Foundation Award No. OAC-1818 253. Financial support was provided by the NSF under Award No. CBET-1602 890.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.