Do Transport Properties of Entangled Linear Polymers Scale with

Oct 16, 2013 - Thus, the influence of entanglements on entropy scalings can be quantified. We make use of thermodynamic integration to determine the e...
0 downloads 10 Views 2MB Size
Article pubs.acs.org/Macromolecules

Do Transport Properties of Entangled Linear Polymers Scale with Excess Entropy? Evangelos Voyiatzis,* Florian Müller-Plathe, and Michael C. Böhm Eduard-Zintl-Institut fur Anorganische und Physikalische Chemie and Center of Smart Interfaces, Technische Universität Darmstadt, Petersenstrasse 20, D-64287Darmstadt, Germany S Supporting Information *

ABSTRACT: We analyze the range of validity of the empirical excess entropy scalings proposed independently by Rosenfeld and Dzugutov for the viscosity and the diffusion coefficient of monodisperse Lennard-Jones chains. They are long enough to be considered as entangled. Thus, the influence of entanglements on entropy scalings can be quantified. We make use of thermodynamic integration to determine the exact excess entropy. Different ways of approximating the excess entropy either based on conformational information or by employing the self-associating fluid theory (SAFT) are explored. The correlations between the various excess entropy estimates are investigated in detail. The conformational route appears to be the most suitable way for the studied Lennard-Jones systems. Its prediction is in qualitative agreement with those obtained by thermodynamic integration, which is computationally more demanding. The SAFT results are not always consistent with the simulation data. The main qualitative feature of both the Rosenfeld and Dzugutov scaling is the coexistence of two linear regimes in the correlation of different entropy expressions. The sharpness of the transition regime depends on the way the excess entropy is approximated. The Dzugutov scaling for the viscosity seems to have certain advantages in comparison to the Rosenfeld one. The influence of the chain length on the parameters appearing in the excess entropy scaling relationships is examined. Despite the pronounced chain length dependence of the scaling parameters, a master curve relating the diffusion coefficient to the excess entropy is derived when performing a suitable renormalization. The relation between the scaling parameters of the viscosity and the diffusion coefficient arewith one exceptionin line with the predictions of the Stokes−Einstein law.

1. INTRODUCTION A long-standing aim in the field of soft matter physics is the formulation of bridging relations between the transport properties of a system and its structural or thermodynamic parameters. The connection between the viscosity, diffusion coefficient, and thermal conductivity, on the one hand, and the excess entropy, on the other, has received particular attention.1−8 There have been several experimental and computational studies that tried to address this issue and to provide empirical evidence for the validity of such a correlation.9−13 They intended to express the transport properties as single-valued functions of the excess entropy. The main motivation of these approaches is easy to explain. If such relations are valid, the transport coefficients, which are difficult to estimate directly, can be determined indirectly from the excess entropy by scaling relations. The excess entropy, on the other hand, can be computed with sufficient accuracy by liquid state theories. Rosenfeld1,2 was among the first to propose a scaling law applicable to fluids and to study these correlations systematically. The original Rosenfeld scaling1,2 was an empirical approach. Because of its importance, there have been various attempts to rationalize its theoretical background.10−13 They aimed to rationalize the scaling behavior of the diffusion coefficient of pure and binary mixtures of atomic and molecular liquids. Especially for binary mixtures, the scaling behavior of the mutual diffusion coefficient received much interest.9 In the © XXXX American Chemical Society

same spirit, an exact scaling relation for polydisperse systems of hard spheres was derived quite recently.10 The Rosenfeld scaling was tested in several “molecular” fluids such as ionic liquids3,4 and hydrocarbons5,14 as well as in real fluids with low molecular weight.15 Modifications in the formalism to treat compounds with anisotropic interactions were derived in ref 11. Apart from the Rosenfeld approach, two other scaling laws have been quite successful in connecting structural properties of a system to its dynamics, namely the formalisms of Dzugutov6 and Bretonnet.7 The main characteristics of the Dzugutov scaling is the possibility to relate dynamical properties to the two-body part of the excess entropy as estimated from the pair correlation function rather than the whole excess entropy of the system. Bretonnet’s scaling law is biparametric, relating the diffusion coefficient to the overall excess entropy of the sample and to the packing density. The latter is defined as the product of the number density of particles and the third power of the position of the first peak in the pair correlation function. On the basis of arguments from mode-coupling theory,12,13 Das and co-workers have shown that scaled diffusivityexcess entropy diagrams are characterized by two intersecting linear regimes.16 The transition between them is more or less sharp, although there is no discontinuity in the curves. Received: August 1, 2013 Revised: September 18, 2013

A

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

studies so far, the two-body part of the excess entropy, which is estimated via the pair correlation function, was assumed to approximate the total excess entropy quantitatively. In the present work a systematic analysis of the influence of radial inhomogeneities in the pair correlation function on the twobody entropy is included. The importance of the three-body entropy relative to the two-body part is also discussed. The present study is the first attempt to estimate the loss of information that has to be expected when truncating the excess entropy expansion after the two-body term. The three-body entropy is computed using Kirkwood’s superposition approximation.27 A direct calculation of the three-body correlation function is thereby avoided. Thus, its numerical accuracy is not plagued by insufficient sampling. Correlations between the various approaches to approximate the excess entropy are established and quantified. The two most frequently employed scaling laws for the diffusion coefficient, the “macroscopic” procedure proposed by Rosenfeld1,2,28 and the “microscopic” one by Dzugutov,6 are examined. We systematically analyze their range of applicability and determine the scaling constants as a function of the chain length. A similar analysis is carried out for the viscosity. The constants from the viscosity scaling are compared with the corresponding constants of the diffusion coefficient scaling by employing the Stokes−Einstein law.

A less rigorous phenomenological relationship between dynamical properties and the configurational arrangement of the particles, as probed by the product Tvγ, where T is the temperature, v the specific volume, and γ a material constant, was established in ref 17. The scaled diffusion coefficient and viscosity were defined by the same macroscopic length-scale and characteristic thermal velocity as proposed by Rosenfeld. The dynamical properties of more than one hundred real substances were tested.17 A similar density scaling was also proposed by Lopez et al.18 The concept of excess entropy scaling of dynamical properties was further extended to local dynamic quantities, such as characteristic times for molecular reorientations and intramolecular motions.8 A scaling law based on Rosenfeld arguments was proposed for the rotational diffusion coefficient of ionic liquids.3 Most of the attempts to explain characteristic times for intramolecular motions by a scaling law focused on either collective or dipole relaxation times. They were extracted by analyzing the decay of coherent scattering functions or the decay of the second Legendre polynomial of the angles formed by polarization vectors.5 Some attempts were made to estimate whether an excess entropy scaling is also valid for confined samples.19−22 These studies focused on simple systems such as confined small molecules in porous media,22 hard spheres,19 and LennardJones particles.20 In the context of scaling concepts, it is of great interest to test whether quantities other than the excess entropy, which also reflect the conformation and the arrangement of the particles, can be employed in the scaling of dynamical properties. The effective packing density, which was defined as the product of the number density of particles and the third power of the equivalent hard-sphere diameter, is one of such possible quantities. The determination of the effective packing density requires the mapping of the actual model onto an equivalent hard-sphere one. Excess entropy scaling for homogeneous Lennard-Jones chains was tested in ref 23. The chains consisted of less than 10 monomers and were represented as freely rotating linear strings of interacting sites. They were short enough to prevent any impact of entanglements on the system dynamics. Although two different approaches were used to estimate the excess entropy in ref 23, no attempt was made to determine the excess entropy directly. The excess entropy contribution due to spatial inhomogeneities in the pair correlation functions was ignored. The aims of the present study are as follows: (i) to quantify the accuracy of the various approaches to estimate the excess entropy of an entangled system and (ii) to analyze whether the excess entropy scaling of the diffusion coefficient and the viscosity is still valid in a system which cannot be classified as a molecular fluid. To this end, the model of long monodisperse Lennard-Jones chains is selected as it captures the essential physics for the present study. Extensive molecular dynamics (MD) simulations are carried out to derive the viscosity and diffusion coefficient as well as the excess entropy of the systems under consideration. Results for the viscosity of Lennard-Jones chains with more than 50 monomers have not been presented in the literature up to now. In this context it should be mentioned that entanglements have a great influence on the dynamics as well as on the conformation of the chains.24 The exact excess entropy of the systems is calculated by thermodynamic integration. An estimate of the excess entropy is obtained by using an equation of state (EOS) which is based on the self-associating fluid theory (SAFT).25,26 In most of the

2. MODEL AND METHODS The monodisperse systems under investigation are formed by 100 freely jointed linear chains. Their length, n, varies; samples with n = 20, 30, 40, 50, 60, or 70 monomers per chain have been simulated. The systems with a chain length greater than 20 are in the entangled regime.23,29 The density, ρ, is defined as ρ = N/V where N is the total number of monomers in the system and V the volume. We have chosen ρ = 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0. In this manuscript all quantities are reported in reduced units.30 Each system is studied at four reduced temperatures, T = 2, 3, 4, or 5. A list of the most important symbols employed in the present work is provided in Table 1. The nonbonded interactions between monomers are modeled by a 12−6 Lennard-Jones potential with a spherical cutoff of 3.30 Longrange corrections to the pressure and the total energy30 are taken into account while the nearest-neighbor nonbonded interactions are excluded. The potential energy for the bond stretching is described via a harmonic potential. The equilibrium bond length is equal to 1. The spring constant has a value of 3000 which is in line with values proposed previously in the literature.23,25 The initial configurations at the chosen densities were created using a simplified “bond growth” procedure. The only restriction imposed was that the distance between any pair of monomers should exceed a value of 1. The configuration was subjected to an energy minimization using the Polak−Ribiére version of the conjugate gradient method.31 Further equilibration was carried out by MD simulations with 13 × 106 time steps of length 0.001 in the canonical (N, V, T) ensemble. The equations of motion were integrated using the velocity−Verlet algorithm and a Nose− Hoover thermostat.30 The production runs were performed in the microcanonical (N, V, E) ensemble. The time for the determination of the diffusion coefficient covered 40 × 106 time steps. The minimum time for the viscosity determination was 80 × 106 time steps, depending primarily on the chain length. All computations were performed using the LAMMPS software.32 The ratio between the standard deviation of the total energy and its average value never exceeded 5 × 10−5. The self-diffusion coefficient of the centers of mass (COM) of the chains, DCOM, was calculated using the Einstein formalism

DCOM = lim

t →∞

B

1 ⟨[rCOM, i(t ) − rCOM, i(0)]2 ⟩ 6t

(1)

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

Table 1. List of Symbols in Alphabetical Ordera A

Helmholtz free energy

n

Aex,SAFT

Helmholtz free energy by SAFT parameter in soft-core potential diffusion coefficient scaled Dzugutov center-ofmass diffusion coefficient scaled Rosenfeld center-ofmass diffusion coefficient pair correlation function spherically averaged pair correlation function angle-dependent pair correlation function

P

number of monomers per chain pressure tensor

r

position vector

S Sex

entropy excess entropy

Sex,SAFT

excess entropy by SAFT ideal gas entropy entropy due to nbody correlations two-body translational entropy three-body translational entropy temperature time nonbonded potential volume

a D D̃ COM D̅ COM g2 g2,ave g2,conf

Sideal Sn S2,trans

g3

three-body correlation function

H kB m N

Hamiltonian of the system T Boltzmann constant t monomer mass unb number of monomers in the V system Greek Symbols collision frequency η̅

ΓE λ η η̃

S3,trans

coupling parameter zero shear viscosity scaled Dzugutov zero shear viscosity

The method proposed by Beutler et al.34 is adopted in the present work to couple the Hamiltonian of the system, H, to the transition parameter λ. The dispersion interactions between two nonbonded monomers, unb, are derived by the formula ⎛ ⎞ 1 1 ⎟ − unb(rij ; λ) = 4λ 4⎜⎜ 2 6 2 a(1 − λ)2 + rij 6 ⎟⎠ ⎝ [a(1 − λ) + rij ] (4)

For λ = 1, the Hamiltonian of the original system is recovered while for λ = 0 the dispersion interactions have been switched off completely and the system has been “transformed” to an ensemble of noninteracting freely jointed chains. The terms containing the prefactor a in the denominators ensure that, in the ideal chain limit, the nonbonded interactions become singularity-free. The noninteracting system is described by the Rouse model whose Helmholtz free energy and entropy can be computed analytically.24 The numerical computation of the integral appearing in eq 3 is carried out employing the GaussKronrod integration rule in a mesh of 15 points.35 The averages appearing in eq 3 have been estimated from MD simulations over 106 time steps. They converged to a relative error of at most 3 × 10−2. In the denominators of eq 4, we selected a = 0.3. 3.2. Configurational Route. The entropy of a molecular system is closely linked to its internal structure as well as to conformational properties which can be quantified by manybody correlation functions. A formal expansion of the entropy, S, in terms of these correlation functions has been derived on the grounds of statistical mechanics:36

scaled Rosenfeld zero shear viscosity density chain density

ρ ρn

a

Greek symbols appear at the end of the table. All indices employed have been explained in the text.

S = Sideal − 1 2 ρ 6 1 + ρ2 6



where ⟨...⟩ denotes an ensemble average and rcom,i(t) the position vector of the COM of the ith chain at time t. The zero-shear viscosity, η, of the systems is also computed by the Einstein formula

η=

V d lim ⟨∑ (Gab(t ) − Gab(0))2 ⟩ 20kBT t →∞ dt ab



∫0

H (λ )

1⎛

⎜⎜⟨H(λ)⟩λ ⎝

∂H(λ) ∂λ

⎞ ⎟⎟ dλ λ⎠

∂H(λ) ∂λ

∫ ∫ ∫ g3(r1, r2, r3) ln(δg3(r1, r2, r3)) dr1 dr2 dr3 ∫ ∫ ∫ (g3(r1, r2, r3) − 3g2(r1, r2)g2(r2, r3) + 3g2(r1, r2) (5)

Sideal is the entropy of an ideal-gas ensemble of monomers at the same density and temperature as the system under consideration, g2 the pair, i.e. two-body, correlation function, g3 the three-body correlation function and ri the position vector of the ith monomer. The function δg3(r1,r2,r3) is defined as δg3(r1,r2,r3) = (g3(r1,r2,r3))/(g2(r1,r2)g2(r1,r3)g2(r2,r3)). Thus, the entropy, S, and the excess entropy, Sex, of the system can be written as

3. ENTROPY ESTIMATION The direct calculation of the exact excess entropy in polymeric systems is not possible as it cannot be estimated by a suitable averaging formalism. This difficulty has triggered the development of several methods to obtain appropriate estimates of the excess entropy.33 This section describes three different routes to address this issue: by employing thermodynamic integration, an entropy series expansion in terms of many-body correlation functions or via the SAFT based EOS formalism. 3.1. Thermodynamic Integration. The excess entropy in a molecular system, Sex, can be determined in the canonical (N, V, T) ensemble by thermodynamic integration.34 1 kBT 2

∫ ∫ (g2(r1, r2) ln(g2(r1, r2)) − g2(r1, r2) + 1) dr1 dr2

− 1) dr1 dr2 dr3 + higher order terms

(2)

with Gab(t) = ∫ 0t (Pab(t′)−(1)/(3)δab∑yPyy(t′))dt′ and kB the Boltzmann constant. The double sum in eq 2 runs over the Cartesian directions. The pressure tensor P was computed using the atomic virial expression;30 it was recorded every 100 time steps.

Sex =

1 ρ 2

N

S = Sideal +

∑ Si i=2

N

and

Sex = S − Sideal =

∑ Si i=2

(6)

where Si is the entropy due to i-body correlations. We expect that Sex will be of negative sign as the entropy of the ideal reference system with complete randomness exceeds the entropy terms of a finite series expansion as well as other entropy estimates encountered in the different entropy definitions. It is quite often assumed that the dominant term in a series expansion is the two-body contribution, S2, and that the remaining Si terms can be neglected.3−5,14,23 Moreover, it is supposed that the system is homogeneous, i.e., that there are no angular fluctuations in the pair correlation function. In this case all contributions to S2with the exception of the translational onecan be neglected. Thus, eqs 5 and 6 are reduced to

λ

(3) C

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules Sex ≅ S2,trans = −2πρ

Article

∫0

− g2,ave(r ) + 1)r 2 dr



(g2,ave(r ) ln(g2,ave(r )) (7)

where S2,trans is the translational two-body entropy; g2,ave(r) is a spherically averaged pair correlation function of the monomers. It includes both intra- and intermolecular pairs, but excludes the contribution due to the bonded first neighbors of each monomer.23 As a result of the homogeneity assumption, the integration over all pairs of atoms is replaced by integration over angularly uniform pair correlation functions. Figure 1a

Figure 2. Chain length and density dependence of the translational pair entropy, S2,trans, (parts a and b) and the excess entropy, Sex,SAFT, based on the SAFT theory (parts c and d) for temperatures between 2 and 5. The density ρ in parts a−c is equal to 1.0 while the chain length n in parts b−d is 50.

the observed changes in the pair correlation function. The empirical formula S2,trans = −3.1877(1 − 0.2046 ln(T))(1 − 0.3739ρ−1)(1 + 0.0905n−1) provides an accurate correlation of S2,trans with T, ρ, and n. The mean square error is equal to 0.0240, which shows that a factorization is possible. In Figure 3, we have correlated S2,trans and the excess entropy from thermodynamic integration, Sex, for six Lennard-Jones

Figure 1. Structure of g2,ave(r) as a function the distance r measured relative to a chosen reference monomer. First neighbors on the same chain are excluded. In part a, g2,ave(r) is shown for four different densities ρ at a temperature T of 5. In part b, g2,ave(r) is presented for four different temperatures at a density ρ of 0.8. The chains have 50 monomers. The bond contribution was not considered in the evaluation of g2,ave(r).

presents the variation of the g2,ave(r) pair correlation function for 50-mer chains of different densities ρ at a temperature of 5.0. The position of the first peak of g2,ave(r) has a weak dependence on ρ and shifts outward with decreasing densities. The remaining weak minima and maxima of the pair correlation function are located at roughly the same distances for all densities studied. Higher peaks and lower valleys are formed with increasing density. Figure 1b presents the g2,ave(r) pair correlation function as a function of temperature for a density of 0.8. A reduction of the temperature results in an increased cohesiveness of the system; the first peak becomes higher. The pair correlation function is more sensitive to density than to temperature variations. The temperature dependence of the translational pair entropy, S2,trans, calculated from g2,ave(r) is presented in Figure 2, parts a and b, as a function of the chain length and density, respectively. The density ρ in Figure 2a is equal to 1.0 while the chain length n in Figure 2b is equal to 50. S2,trans reacts more sensitive to density and temperature variations than to length variations. As mentioned above, the estimated excess entropy is always negative. The bond-stretching and nonbonded interactions in the systems under consideration reduce the number of accessible configurations compared to the ideal reference state. The translational pair entropy increases slightly with increasing chain length. The variation of S2,trans is in line with

Figure 3. Correlation between the excess entropy as predicted by thermodynamic integration, Sex, and the translational pair correlation entropy S2,trans for chain lengths between 20 and 70 monomers. The data have been derived within temperature and density intervals considered in Figure 1. In the inset we have plotted the Sex and S2,trans correlation for the 70mer alone to show the standard deviation for each entry. The error bars correspond to the standard deviation. The dashed line refers to the best linear fit between S2,trans and Sex.

chains of different lengths. The entropies for each chain have been derived for the temperatures and densities employed above. We have prepared an inset with the entropy correlation for the chain with 70 monomers to emphasize the almost linear behavior of the curve and to show the error bars of our D

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

atures considered in the plot are the same as employed in the foregoing entropy correlations. We can identify two different linear regimes for the systems under consideration. Irrespective of the chain length, the transition region is found for S2,trans values between −1.6 and −1.7. Above the “transition domain” all curves have the same slope. Only the constant term depends slightly on the chain length. The correlations can be approximated by a single master curve. Below the transition domain, the sum S2,trans + S2,conf depends somehow more on the chain length. With enlarged negative values of S2,trans, the S2,conf term becomes larger in magnitude. The configurational entropy is small above the transition region which implies that S2,trans + S2,conf correlates linearly with S2,trans with a chain length independent slope close to one. A better estimate of the excess entropy can be obtained by taking into account also the third-order term of the entropy expansion, i.e., the three-body entropy S3. It depends on the relative distance between the three pairs of particles in the particle triplets as well as on the angles between two joining vectors. The numerical accuracy in the evaluation of that term, however, is usually quite poor.27,39 It suffers from insufficient sampling over all possible angles between the two joining vectors. This is the reason why it has been neglected in most of the studies so far. A way to circumvent this problem is to adopt an approximate factorization of the three-body correlation function.27,39,40 The simplest scheme was proposed by Kirkwood many decades ago27 on the basis of a “superposition approximation”. In this scheme one assumes that the three-body correlation function g3 does not depend on the absolute position vectors of the particles forming a triplet, r1, r2, and r3, respectively. It is instead a function of the distance between the three possible pairs that can be formed between two particles of the triplet: ||r1 − r2||, ||r1 − r3|| and ||r2 − r3||. Any angular dependence of g3 is neglected in this degree of sophistication. Thus, it can be expressed as the product of the pair correlation function g2,ave(r):

simulations. From the computational results we deduce that all studied chains can be described by a single master curve which reads Sex = 0.806 S2,trans − 1.796. Despite a S2,trans prefactor smaller than unity, the magnitude of Sex exceeds the one of S2,trans. It follows from the large negative constant on the righthand side of the given correlation. The squared regression coefficient R2 amounts to 0.94. When restricting the linear fit to the sample with 70 monomers, we derive Sex = 0.816 S2,trans − 1.780 and R2 = 0.99 which differs only slightly from the master curve sampled over six chain lengths. The theoretical findings from Figure 3 can be summarized as follows: Although S2,trans accounts only for the translational part of the two-body entropy S2, it nevertheless renders possible a qualitative reproduction of Sex. An improvement of the two-body entropy is possible by taking into account inhomogeneities in the pair correlation function. Lazaridis and Karplus37 have proposed a splitting of S2 into two terms; one taking into account the translational contribution, already defined in eq 7, and a second configurational two-body contribution S2,conf = −

ρ 2





π

∫r=0 g2,ave(r) ∫φ=0 ∫θ=0 g2,conf (r , θ , φ)

ln(g 2,conf (r , θ , φ)) sin θ dθ dφ dr

(8)

which accounts also for angular variations in the pair correlation. In line with this philosophy, they have split the pair correlation function g2 appearing in the entropy definition of eq 5, into two terms which appeared already in eq 7 or eq 8: g 2 (r , φ , θ ) = =

V d N (r , φ , θ ) N d V (r , φ , θ ) V dN (r ) 4π d2N (r , φ , θ ) N 4πr 2 dr dN (r) dφ dθ

= g2,ave(r )g2,conf (r , φ , θ )

(9)

g 3(r1, r2, r3) = g2,ave(|| r1 − r2 ||)g2,ave(|| r1 − r3 ||)g2,ave(|| r2 − r3 ||)

More elaborate expressions taking into account the orientational two-body entropy have been also derived.38 They are of particular interest in systems where certain orientations of a molecule or a group of atoms are energetically favored.3 Figure 4 presents the variation of S2,trans + S2,conf, defined by eqs 7 and 8, as a function of S2,trans. Chains with lengths between 20 and 70 monomers have been studied. The densities and temper-

(10)

This scheme implies that the quantity δg3(r1,r2,r3) in eq 5 is equal to 1. The second integral appearing on the right-hand side of the entropy expansion in eq 5, which contains the logarithm of δg3(r1,r2,r3), is thus equal to zero. The third integral of this entropy expansion can be expressed in terms of the pair correlation function g2,ave(r): S3,trans =

8π 2ρ2 3



∫0 ∫0



(g2,ave(x)g2,ave(y)g2,ave(|x − y|)

− 3g2,ave(x)g2,ave(y) + 3g2,ave(x) − 1)x 2y 2 dx dy

(11)

In the present study, the two-dimensional integral is evaluated numerically using the trapezoidal rule. The main advantage of this approach is that only knowledge of the pair correlation function is required. Although there are more elaborate expressions for the δg3(r1,r2,r3) term, such as the ones proposed in,39,40 the Kirkwood approximation is accurate enough if the temperature of the system is not too low and the density high enough. In Figure 5 the correlation between S3,trans + S2,trans and S2,trans is presented for all systems studied. A linear fit allows an accurate correlation between S3,trans and S2,trans; it reads S3,trans = 0.217 S2,trans + 0.062 with R2 = 0.99. Thus, S3,trans is strongly correlated to the more dominant S2,trans term. Again,

Figure 4. Correlation of the sum of the translational and configurational pair correlation entropy S2,trans + S2,conf with the translational pair correlation entropy S2,trans for chain lengths between 20 and 70 monomers. The data have been derived for the temperature and density values considered in Figure 1 E

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

Figure 5. Correlation of the sum of the translational pair correlation entropy S2,trans and the translational three-body correlation entropy S3,trans with the translational pair correlation entropy S2,trans. The data have been derived for temperature and density values considered in Figure 1

Figure 6. Correlation between the excess entropy as predicted by thermodynamic integration, Sex, with the excess entropy predicted by the geometric progression of the translational terms. The data have been derived for temperature and density values considered in Figure 1

this might be caused by the adopted pair approximation for S3,trans. A theoretical justification for this behavior in the present degree of sophistication follows from a reformulation of eq 11 in terms of S2,trans. For instance, when adding and subtracting simultaneously the factor g2,ave(x)g2,ave(y) ln g2,ave(y) + 1 to the second term of eq 11, we derive 8π 2ρ2 3



∫0 ∫0



= −4πρS2,trans − 8π 2ρ2



3.3. Self-Associating Fluid Theory. A different route to estimate the excess entropy is to employ an appropriate EOS description, such as the one developed in SAFT.25 Such an approach has the advantage of estimating the whole excess entropy rather than approximating it through the first term of a series expansion, as done in the previous procedure. On the downside, the entropy cannot be better than the model assumptions of the EOS. A variety of SAFT EOS have been developed in order to describe a wide range of systems, ranging from simple Lennard-Jones chains to very complicated ones, such as ionic liquids41 and polymers.42,43 The SAFT EOS attempts to model the Helmholtz free energy of the system of interest. The excess entropy is estimated by differentiation of this free energy with respect to the temperature. The SAFT EOS has several adjustable parameters which are fine-tuned against experimental data.41,43 Its parametrization is most sensitive to experimental vapor pressure and liquid−liquid data. The force field parameters used in atomistic simulations cannot be used directly in the fine-tuning process. As a result of the difficulties in computing liquid−vapor and liquid−liquid equilibrium quantities for atomistic systems and the large number of required PVT data in the parametrization process, it is not possible to guarantee consistency between the employed parameters of an atomistic simulation and the resulting SAFT EOS. The excess Helmholtz free energy Aex,SAFT of a LennardJones chain system in the SAFT model is written as25,26 Aex,SAFT = A ref + Ach (13)

−3g2,ave(x)g2,ave(y)x 2y 2 dx dy

∫0

∫0 ∫0





g2,ave(x)x 2 dx

(g2,ave(x)g2,ave(y) ln g2,ave(y) + 1)x 2y 2 dx dy

(12)

An interesting feature of the S3,trans−S2,trans correlation is that neither its slope nor its constant term depend on the chain length, n. Thus, the predominance of S2,trans is conserved also for long entangled chains. The correlation between S3,trans and S2,trans indicates that the former term is roughly equal to one-fifth of the latter. Inspired by the fact that a geometric progression relationship holds for the translational entropy terms of the ideal gas,36 we assume a similar behavior in the case of the Lennard-Jones chains. All translational entropy terms are members of a geometric progression with a common ratio of 1/5. However, the geometric progression for the ideal gas has a nonconstant common ratio which decreases with increasing order of the translational terms. This means that the “total” translational entropy estimated by employing the concept of a geometric progression with a constant ratio overestimates the translational contribution to the excess entropy. Nevertheless, it can be expected that such an effect is absent in the studied systems with their reduced degrees of freedom. If a geometric progression approximation is valid, the total translational contribution, Strans, to Sex is equal to (1−(1/5)n)/(1 − 1/5) S2,trans ≅ 1.25S2,trans, where n is the chain length. The correlation between the excess entropy as calculated by thermodynamic integration and Strans as estimated by the proposed geometric progression is shown in Figure 6. The predicted Strans is still smaller than the excess entropy computed by thermodynamic integration. The neglected configurational and other nontranslational entropy terms, however, are partially compensated by the employed geometric approximation for the translational entropy. In line with the findings from the S2,trans − Sex correlation, all studied systems are described by a single master curve.

where Aref is the excess Helmholtz free energy of an atomic Lennard-Jones fluid at the same thermodynamic conditions as the chain system, and Ach is the excess Helmholtz free energy which takes into account the existence of bonds in the system. The excess Helmholtz free energy of a Lennard-Jones atomic fluid can be computed accurately by a modified Benedict− Webb−Rubin EOS.44 The temperatures and densities of the systems studied here lie within those used to estimate the parameters of this EOS. The Ach term in eq 13 is derived on the basis of first-order thermodynamic perturbation theory (TPT1).25 Any information on bond angles and any kind of intramolecular distribution beyond the bond distribution is ignored. The theory does not F

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

with S2,trans. Thus, a master curve can be derived for all simulation data. On the basis of the qualitative agreement between S2,trans and Sex as well as the established linear correlation between them, which was presented in Figure 3, it was expected that, in a similar fashion, S2,trans can be correlated linearly with Sex,SAFT. Similar to the S2,conf − S2,trans case, the correlation of Sex,SAFT and S2,trans can be described accurately by two linear regimes. The transition is rather smooth and takes place in the S2,trans range between −1.7 and −1.6. Above the crossing region we derive Sex,SAFT = 1.519 S2,trans + 0.642 with R2 = 0.98. For the linear regime below the crossing region we observe Sex,SAFT = 0.942 S2,trans − 0.315 with R2 = 0.99. The Helmholtz free energy of a given Lennard-Jones chain system, Aex,SAFT, is written as the sum of the free energy of the reference Lennard-Jones fluid, Aref, and the free energy due to the existence of bonds, Ach, in eq 13. Since the former free energy term is modeled very accurately by the employed modified Benedict−Webb−Rubin EOS,44 any qualitative deviation from a linear correlation between S2,trans and Sex,SAFT should be attributed to the way the second free energy term Ach is estimated. The TPT1 theory neglects the intramolecular interactions due to the attractive part of the reference potential.25 The longer the chains under constant T and ρ, the more important this term is. Moreover, this theory neglects all intramolecular correlations higher than the two body ones. For low temperatures or high densities these correlations become stronger and cause angular inhomogeneities in the system, which are quantified through deviations of the angledependent pair correlation function g2,conf(r,φ,θ) from unity. This element was defined in eq 9 by the formula g2,conf(r,φ,θ) = 4π d2N(r,φ,θ)/(dN(r) dφ dθ) and triggers increased negative values of the S2,conf term. This might be the reason why the crossover region in the Sex,SAFT − S2,trans scheme occurs in the same region as encountered in the S2,conf − S2,trans correlation.

take into account intramolecular interactions described by the attractive part of the reference potential. The larger the chain length, the more important these interactions are. The Ach term is estimated via Ach =

1−n T ln(g2,ave(1)) n

(14)

The excess entropy is computed by differentiating the excess Helmholtz free energy with respect to the temperature. The final equation for the SAFT excess entropy Sex,SAFT is Sex,SAFT(ρn ) = Sex,ref (ρ) +

⎡ ∂[ln(g 2,ave(1))] ⎤ n − 1⎢ ⎥ ln(g2,ave(1)) + T ⎥⎦ ∂T n ⎢⎣

(15)

where ρn = ρ/n is the chain number density; Sex,ref(ρ) is the excess entropy of a Lennard-Jones monatomic fluid at the same thermodynamic conditions, i.e., temperature and density, as the Lennard-Jones chain system. The variation of the SAFT-EOS-based excess entropy, Sex,SAFT, as a function of the chain length and density has been presented already in parts c and d of Figure 2 for different temperatures. In Figure 2c the density ρ is equal to 1.0 while in Figure 2d the chain length n is equal to 50. In analogy to S2,trans in Figure 2a, the density dependence of Sex,SAFT exceeds the chain length dependence. The temperature dependence is similar to the one of Figure 2a. The magnitude of Sex,SAFT is enhanced with decreasing temperature. A direct comparison of the chain length, density and temperature dependence between S2,trans and Sex,SAFT is possible by correlating parts a and b of Figure 2 with parts c and d of Figure 2. Inspection of parts a and c of Figure 2a reveals that both entropy estimators are more sensitive to temperature than to length variations. When comparing the low-density regimes of Figure 2, parts b and d, one sees that Sex,SAFT is smaller in magnitude than S2,trans. These data indicate certain simplifications of the SAFT approach, since the EOS is expected to take into account all entropy contributions in contrast to S2,trans; it is not restricted to the translational two-body part. Nevertheless, the qualitative response of Sex,SAFT and S2,trans to density and chain length variations is similar. Figure 7 presents the correlation of the SAFT excess entropy Sex,SAFT with S2,trans. The SAFT excess entropy correlates well

4. TRANSPORT COEFFICIENTS Prior to the study of the excess entropy scaling of the diffusion coefficient and viscosity, we test whether the system under consideration can be classified as a “strongly correlating liquid”.49,50 To this end we measure whether the fluctuations of the potential energy from its thermodynamic average are strongly correlated to fluctuations of the virial. A fluid is classified as “strongly correlating” if this correlation coefficient is greater than 0.9. It has been shown that phenomenological scalings, such as the excess entropy ones, are valid for strongly correlated systems. Nevertheless, our preliminary results revealed that the correlation coefficient ranges between 0.5 and 0.6 in all cases. Thus, although a statistically significant correlation between the virial and the potential energy exists, our systems cannot be classified as “strongly correlating” ones. The scaling behavior of the diffusion coefficient and viscosity in terms of the different entropy estimators are now examined in detail. These transport coefficients are made dimensionless according to the two main scaling laws considered here. The scaling factors have been derived on the grounds of either continuum-based28 or atomic-based6 theories. We start with an analysis of the scaled diffusion coefficient as determined by the Rosenfeld and Dzugutov laws. Then, we proceed with the discussion of the scaled viscosity by both laws in the second part of this section. 4.1. Diffusion Coefficient Scaling. One of the most important steps in establishing excess entropy scaling laws was the development of scaling factors for the diffusion coefficient

Figure 7. Correlation of the SAFT excess entropy Sex,SAFT with the translational pair correlation entropy S2,trans. The temperatures and densities considered in this diagram coincide with those adopted in the foregoing figures. The continuous line allows a fitting of the data in the upper right region while the dashed curve is valid in the lower left region. The two curves have been extended over the whole entropy interval to show the deviations between the two lines with increasing separation from the crossing region. G

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

Figure 8. Correlation between the natural logarithm of the dimensionless Rosenfeld diffusion coefficient D̅ COM and (a) the excess entropy Sex predicted by thermodynamic integration, (b) the translational two-body entropy S2,trans, (c) the sum S2,trans + S2,conf, and (d) the excess entropy Sex,SAFT as predicted by the SAFT EOS. The diffusion coefficient has been scaled according to the Rosenfeld law for Lennard-Jones chains from 20 to 70 monomers. The temperatures and densities chosen correspond to the same setup as employed in the other correlation schemes.

Table 2. Rosenfeld Excess Entropy Scaling Parameters for the Reduced Diffusion Coefficient D̅ COM = AR exp(αRSex) low entropy regime Sex S2,trans S2,trans + S2,conf

high entropy regime

αR

ln(AR)

αR

ln(AR)

3.2 2.26/n + 2.53 8.704/n + 0.314

42.52/n − 4.611 −48.17/n − 9.941 58.95/n − 13.85

1.8 −9.95/n + 1.53 −3.640/n + 2.266

91.12/n − 10.32 −27.57/n − 8.415 57.47/n − 10.94

by Rosenfeld.28 This method makes use of two microscopic factors which have their origin in kinetic gas theory: the mean interparticle distance, l = ρ−1/3, and the thermal velocity, uthermal = (kBT/m)1/2. Both microscopic factors can be calculated from macroscopic observables. The normalization of the diffusion coefficient by these characteristic length and velocity scales leads to

D̅ = D

ρ1/3 kBT /m

The variation of ln(D̅ COM) with Sex from thermodynamic integration is given in Figure 8a. Although there are deviations from a linear behavior, ln(D̅ COM) and Sex are well-correlated for a given chain length. Similar deviations in the scaling of the reduced diffusion coefficient with either S2,trans or S2,trans + S2,conf are observed in Figure 8b and 8c, respectively. An attempt to fit the data to two nonlinear functional forms has been made by employing either a two-term entropy series expansion of the form D̅ COM = AR exp(αR,1Sex + AR exp(αR,2S2ex) or a power law expression which reads D̅ COM = AR exp(αR(−Sex)v) where ν is a fitting exponent. Nevertheless, none of them could describe accurately the scaling of the diffusion coefficient with all studied excess entropy estimators. The series expansion fails to reproduce sharp transitions such as the ones encountered in the scaling of the diffusion coefficient or the viscosity when S2,conf is employed. The fits obtained from a power-law are plagued by a high statistical uncertainty. A more detailed discussion is provided in the Supporting Information. Despite the scatter of the data, we suggest that two regimes can be identified where (near-) linear relations between ln(D̅ COM) and Sex are observed (Figure 8a). The transition between these regimes takes place in the Sex range between −3.0 and −3.2. The excess entropy scaling parameters, αR and ln(AR), are given in Table 2. In contrast to the chain length dependence of pre-exponential factor AR, the slope αR in the

(16)

The main advantage of this law is that it can be applied directly to real materials. In the case of polymers, the interparticle distance is normalized by the chain density, ρn, and the thermal velocity by the total mass of the chain (i.e., m is the monomer mass).28 These changes lead to DCOM = DCOM ̅

ρn1/3 kBT /nm

(17)

Rosenfeld proposed to relate the reduced diffusion coefficient to an exponential of the excess entropy of the system: D̅ COM = AR exp(αRSex) where AR and αR are coefficients. Both of them depend on the nature of the system: in our case, only on the chain length and the way to approximate the excess entropy. H

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

ln(D̅ COM)−Sex correlation seems to be independent of the chain length. This is valid in both linear regimes. The correlation between the logarithmic ratio ln(D̅ COM/AR), where AR is the pre-exponential factor in the low entropy regime, and the excess entropy Sex predicted by thermodynamic integration is given in Figure 9. With the exception of the

phases. Dzugutov’s scaling reflects the idea that for short time scales the atomic dynamics is dominated by a cage effect and local density fluctuations. The rate of such a “collision” process within the cage, ΓE, can be described by a generalization of the Enskog theory45 2 ΓE = 4rmax g2,ave(rmax )ρ πk BT /m

(18)

where rmax is the position of the first maximum in the pair correlation function g2,ave(r). For longer time scales, the dynamics is dominated by the rate at which the local surroundings of each atom, i.e. the “cage”, are modified and the cage structure is relaxed. This rate is proportional to the number of accessible configurations. The difference between an ideal gas state and a condensed (either liquid or solid) one is assumed to be given by the factor exp(S2/kB).6 On the basis of this assumption, Dzugutov proposed a proportionality between D/(r2maxΓE) and exp(S2/kB) 2 D/(rmax ΓE) = 0.049 exp(S2/kB)

(19)

The Dzugutov normalization can be expressed in a form more suitable for polymers by normalizing DCOM and using the polymer mass (nm) instead of the monomer mass in the collision frequency formula: ΓE = 4r2maxg2,ave(rmax)ρn(πkBT/ nm)1/2. The Dzugutov scaling relates the dimensionless diffusion coefficient, D̃ COM, only to the two-body entropy, S2, and not to the total excess entropy as done in the Rosenfeld approach. In contrast to the elaborate calculation of Sex, the S2 entropy can be estimated easily and directly from a molecular simulation. Thus, the Dzugutov scaling appears to be more convenient than Rosenfeld’s. Although the theory predicts that D̃ COM should be proportional to exp(S2/kB), we focus on the development of more general linear correlations between ln(D̃ COM) and Sex of the form ln(D̃ COM) = ln(ADZ) + αDZSex. It is of general interest to test such correlations and to explore their region of validity as well as to cast light on the dependence of ln(ADZ) and αDZ on the chain length. The variation of ln(D̃ COM) with Sex from thermodynamic integration is given in Figure 10a. This correlation cannot be described accurately by a single linear expression as postulated by Dzugutov’s theory. However, the data can be fitted linearly when allowing a fragmentation into two domains. They are smoothly interchanged via a crossover domain that spans between entropy values of −2.9 to −3.0. The width of the crossover domain is nearly independent of the chain length. When comparing with the Rosenfeld scaling, it is found that, in the Dzugutov case, the transition takes place slightly earlier and the width of the domain is narrower. The slope, αDZ, and the pre-exponential parameter, ln(ADZ), correlate very well with the chain length, n. The correlations are given in Table 3. An important feature of these correlations is that αDZ and ln(ADZ)dependent on n−1. This kind of dependence signals that a terminal value for both quantities in the low and high entropy domains exist in the limit of infinitely long chain systems. Figure 10b presents the variation of ln(D̃ COM) as a function of S2,trans. Similarly to the previous case, the whole domain is considered to consist of two linear regions. The transition region between the two linear regimes occurs between −1.5 and −1.7; it depends weakly on the chain length. The αDZ terms in both entropy regimes are expressed in the form of two weakly n-dependent linear equations, where the slopes have a very small, close to zero, value. The correlations indicate also

Figure 9. Correlation between ln(D̅ COM/AR), where AR is the preexponential factor in the low entropy regime, and the excess entropy Sex predicted by thermodynamic integration. The diffusion coefficient for Lennard-Jones chains from 20 to 70 monomers has been scaled according to the Rosenfeld law. The temperatures and densities chosen correspond to the same setup as employed in the other correlation schemes. The symbols for the different chain lengths have been explained in Figure 8.

shortest chain length, i.e., n = 20, which corresponds to the only unentangled system, all data fall into a master curve within the limits of the statistical uncertainties. This fact elevates AR from a merely mathematical quantity to a diffusion coefficient at Sex = 0 under the hypothetical condition that the system is still in the same phase. Figure 8b presents the variation of ln D̅ COM as a function of the S2,trans entropy; see also the second line of Table 2. Similar to the previous case, there are two (near-) linear regimes with a smooth transition between them at S2,trans ∼ −1.6. Because of the linear dependence between S3,trans and S2,trans, it is expected that the ln(D̅ COM)−(S2,trans + S3,trans) correlation will be similar to the ln(D̅ COM)−S2,trans one. The variation of ln(D̅ COM) with S2,trans + S2,conf is shown in Figure 8c. In analogy to the ln(D̅ COM)−S2,trans correlation, there are two different linear regimes in this scaling. However, now the two gradients exhibit a larger difference than the ones encountered in Figure 8b. The transition point near −1.8 appears to be independent of the chain length. The variation of the ln(D̅ COM) as a function of Sex,SAFT is given in Figure 8d. In contrast to the exact Sex scaling by thermodynamic integration, an almost linear dependence is observed for the whole Sex,SAFT range. The calculated slope is found between 1.51 and 1.58. Given the uncertainties when determining αR and its weak n-dependence, we assume that αR is constant. With the exception of the 20mer, the ln(AR) term can be described by the correlation ln(AR) = −56.46/n − 9.093. An alternative approach to establish a correlation between the diffusion coefficient and the excess entropy was proposed by Dzugutov.6 His scaling was derived for simple atomic and molecular systems. It was tested for both liquid and crystal I

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

dependence and the statistical uncertainties in determining its numerical value from regression analysis. In contrast to the Rosenfeld case, ln(AR) displays also a weak n-dependence with the exception of the unentangled system, i.e. the 20mer chains. Here it is approximately equal to −10.4. The qualitative behavior of the ln(DCOM)−S2,trans correlation in both scaling approaches is identical; both predict a coexistence of two linear regimes. When the contribution from radial inhomogeneities in the pair correlation function to the two-body entropy is taken into account, the transition between the two linear ln(DCOM)−S2 regimes is significantly sharper in the Dzugutov than in the Rosenfeld method. The values of D̅ COM and D̃ COM lie approximately in the same range and there is no pronounced difference between the numerical values of the two scaling factors. A comparison of our findings with the ones presented in ref 23 for unentangled chains is possible, when considering only the Rosenfeld scaling of the diffusion coefficient and the excess entropy estimated by SAFT. The qualitative feature to have only one linear regime is valid in both cases. However, for the entangled systems under consideration, this is not consistent with the Rosenfeld scaling using the exact excess entropy from thermodynamic integration. The parameter αR is found to be independent of the chain length, n, while in ref 23, αR correlates with n−1 for various short chain lengths. It might be the case that the chain lengths considered in this study are long enough so that αR reaches the plateau value predicted by the correlation valid for unentangled systems. Nevertheless, the difference between the two actual αR values is greater than the uncertainty due to statistical errors. A much more pronounced difference of the pre-exponential parameter ln(AR) as a function of n is found; while a power-law dependence is reported for the unentangled systems, an exponential inverse n-dependence seems to be valid for longer chain lengths. 4.2. Viscosity Scaling. The original Rosenfeld normalization for the viscosity, η, was developed for simple atomic systems.1,2 The viscosity is scaled by employing the same characteristic mean interparticle distance and thermal particle velocity as already adopted for the diffusion coefficient scaling. Thus, the scaled viscosity, η̅, is defined as

Figure 10. Correlation between the natural logarithm of the dimensionless Dzugutov diffusion coefficient D̃ COM and (a) the excess entropy Sex predicted by thermodynamic integration, (b) the translational two-body entropy S2,trans, (c) the sum S2,trans + S2,conf, and (d) the excess entropy Sex,SAFT as predicted by the SAFT EOS. The diffusion coefficient has been scaled according to the Dzugutov law for Lennard-Jones chains from 20 to 70 monomers. The temperatures and densities chosen correspond to the same setup as employed in the other correlation schemes. The symbols refer to the same chain lengths as, for example, in Figure 8

that αDZ increases with the chain length n in the high entropy region while it is reduced as a function of n in the low entropy regime. Similar to the previous case, the contribution of the constant term in ln(ADZ) is more important than the one of the n-dependent slope. Increasing chain length implies a more negative constant term in both regimes. Figure 10c presents the variation of ln(D̃ COM) as a function of the two-body entropies S2,trans+ S2,conf. Differences between the two regimes are larger than for S2,trans alone (Figure 10b) indicating that the contribution of S2,conf is rather important and that the inhomogeneities in the radial distribution function probably become larger in the high-entropy regime than the low one. The crossover takes place for entropy values between −1.7 to −1.8. The n-dependent kink in the scaled diffusion coefficient − entropy plot has been explained in the literature by arguments from mode coupling theory.12,46 The variation of ln(D̃ COM) as a function of the SAFT estimate Sex,SAFT is given in Figure 10d. The ln(D̃ COM)−Sex,SAFT correlation shows the same qualitative behavior as encountered in the Rosenfeld approach. A single linear dependence is sufficient to describe accurately the scaled diffusion−entropy data for a given chain length. The slope of the proposed correlation, αR, ranges in the small interval between 2.18 and 2.26. Thus, it can be assumed as constant, due to its weak n-

η̅ = η

ρ−2/3 mkBT

(23)

It is possible to express the original reduction formula in a form which is more suitable for molecular systems. The chosen reduction parameters should reflect characteristic intermolecular distances and molecular thermal velocities. Taking m as the monomer mass, the molecular version of the Rosenfeld scaling for the dimensionless viscosity can be written as η̅ = η

ρn−2/3 nmkBT

(24)

Table 3. Dzugutov Excess Entropy Scaling Parameters for the Reduced Diffusion Coefficient D̃ COM = ADZ exp(αDZSex) low entropy regime Sex S2,trans S2,trans + S2,conf

high entropy regime

αDz

ln(ADz)

αDz

ln(ADz)

−18.82/n + 2.716 −11.96/n + 1.762 −13.01/n + 0.979

42.08/n − 0.934 61.15/n − 9.558 26.62/n − 8.125

15.44/n + 4.406 9.645/n + 3.676 18.33/n + 4.213

−44.324/n − 6.315 −49.33/n − 11.18 −49.96/n − 12.04

J

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

encountered in the ln(D̅ COM)−S2,trans scaling. It is found that βR ≅ αR, thus the Stokes−Einstein law is approximately valid. Figure 11c presents the variation of ln(η)̅ with S2,trans + S2,conf. Similar to the two former correlations, there are two linear regimes. The transition is sharper than in the case of the diffusion coefficient. It takes place at an entropy value of approximately −1.8. The detected correlations, however, are not valid for the 20mer chain. A theoretical disadvantage of using S2,trans + S2,conf as an approximation to the total excess entropy is the observed deviation from the Stokes−Einstein law; now we observe βR ≠ −αR. The correlation between ln(η̅) and Sex,SAFT is given in Figure 11d. The ln(η)−S ex,SAFT relationship can be described by a ̅ single linear expression as proposed by Rosenfeld.1,28 BR is correlated to n according to ln(BR) = 65.71/n+4.880. Similar to the case where Sex was estimated by thermodynamic integration, the parameter βR, in the ln(η)−S ex,SAFT correlations ̅ for different chain lengths, cannot be expressed as a function of n. There is no clear trend between βR and n. When comparing this behavior with the corresponding relations for the diffusion coefficient, it seems that the Stokes−Einstein law is valid and that the identity βR = −αR does hold. In contrast to ref 23, where only unentangled systems were investigated, the dependence of BR is inversely proportional to n and not proportional to a fractional power form. Moreover, the exponential parameter βR displays a significantly different ndependence as the one reported in ref 23. The parameter βR is found to be almost constant; it does not correlate with n−1. Although not explicitly formulated in the original publication, it is possible to derive an appropriate viscosity normalization scheme for simple atomic systems when using Dzugutov’s arguments.48 The characteristic short time- and length-scales of the viscosity normalization are still defined by the collision frequency ΓE, which is given in eq 18, and the position of the first maximum, rmax, in the pair correlation function g2,ave(r). The scaled viscosity of atomic systems, η̃, is defined as r η ̃ = η max ΓEm (25)

which differs from eq 23 by the chain length n and the adoption of the chain density ρn. The Rosenfeld theory predicts that the scaled viscosity is related to the exponential of the excess entropy of the system: η̅ = BR exp(βRSex), where BR and βR are coefficients. Thus, it can be used to develop linear relations of the type ln(η)̅ = ln(BR) + βRSex and to explain how BR and βR are related to the chain length, n, for each way of estimating the excess entropy of the system. Figure 11a presents ln(η)̅ as a function of Sex computed by thermodynamic integration. Despite the uncertainties in the Sex

Figure 11. Correlation between the natural logarithm of the dimensionless viscosity η̅ and (a) the excess entropy Sex predicted by thermodynamic integration, (b) the translational two-body entropy S2,trans, (c) the sum S2,trans + S2,conf, and (d) the excess entropy Sex,SAFT as predicted by the SAFT EOS. The viscosity has been scaled according to the Rosenfeld law for Lennard-Jones chains with 20 to 70 monomers. The temperatures and densities chosen correspond to the same setup as employed in the other correlation schemes. The symbols refer to the same chain lengths as, for example, in Figure 8

In analogy to the diffusion coefficient, a suitable adaptation of this scheme is possible so as to be applicable to molecular systems. In such a case, the mass of the whole molecule should be employed as well as the molecular expression of the collision frequency: ΓE = 4r2maxg2,ave(rmax)ρn(πKBT/nm)1/2. Finally, the molecular Dzugutov definition of the dimensionless viscosity, η̃, can be written as r η ̃ = η max ΓEnm (26)

calculation, two linear regimes can be identified. It is possible to establish linear correlations between Sex and ln(η̅) for all chain lengths studied. The transition occurs when Sex is approximately between −3.0 and −3.2. The scaling of BR and βR with n is shown in Table 4. Similar to the correlation between the normalized Rosenfeld diffusion coefficient and Sex, the slope βR seems to be independent of the chain length. When comparing βR with the parameter αR in the scaled Rosenfeld diffusion coefficient, we find, within the limits of the statistical uncertainties, βR = −αR. This observation signifies the validity of the Stokes−Einstein law for the studied systems. The variation of ln(η)̅ with S2,trans is given in Figure 11b. There are two linear regimes with a smooth transition between them. The transition takes place at the same entropy as

The discussion focuses on the development of linear relations of the type ln(η̃) = ln(BDZ) + βDZSex and on a meaningful correlation of the two constant terms BDZ and βDZ with the chain length, n, for each way of estimating or approximating the excess entropy of the system. The present functional form of

Table 4. Rosenfeld Excess Entropy Scaling Parameters for the Reduced Viscosity η̅ = BR exp(βRSex) low entropy regime Sex S2,trans S2,trans + S2,conf

high entropy regime

βR

ln(BR)

βR

ln(BR)

−3.4 −46.58/n − 0.423 −27.31/n + 0.012

36.45/n − 1.713 −183.3/n + 11.58 −160.5/n + 12.39

−2.0 64.21/n − 3.874 47.66/n − 3.217

87.07/n + 1.910 −76.76/n + 7.367 −99.85/n + 8.247

K

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

in the ln(D̃ COM)−S2,trans + S2,conf diagram, i.e., between −1.7 and −1.8. The crossover is sharper than the one observed in the ln(η̃)−S2,trans case. The Stokes−Einstein law is approximately valid in the high entropy regime. The variation of the natural logarithm of the scaled viscosity η̃ as a function of Sex,SAFT is given in Figure 12d. Similarly to all previously studied correlations using Sex,SAFT as an entropy estimator, a linear relationship is found to describe adequately the ln(η̃)−excess entropy correlation in the whole entropy domain. Although the dependence of the slope, βDz, and the constant term, ln(BDZ), on n is weak, two correlations can be established which read βDZ = −22.46/n − 0.402 and ln(BDZ) = 3.277/n − 1.247. They suggest that for sufficiently long chains both βDz and ln(BDZ) have an almost constant value independent of n; the two parameters are roughly equal to −0.402 and 1.247. Significant deviations from the Stokes− Einstein law are found when comparing βDz with αDz values obtained from the Dzugutov scaling of the diffusion coefficient.

the pursued linear relations is more general than the original one.6,48 The variation of ln(η̃) with Sex from thermodynamic integration is given in Figure 12a. The variation of the scaled

5. CONCLUSIONS In the present simulation study of polymers modeled as Lennard-Jones chains we have analyzed two different problems. On the one hand, we have explored whether dynamic properties of entangled linear chains scale with the excess entropy computed via thermodynamic integration. On the other hand we have analyzed possibilities to estimate the excess entropy of these systems by simpler methods. We have compared two different techniques for reducing dynamic properties: (i) the Rosenfeld approach which is based on macroscopic arguments and (ii) the Dzugutov formalism which has microscopic foundations. Two approximation schemes for the excess entropy have been examined; in the first one, configurational information is taken into account while the second is based on SAFT. The capability of the two schemes has been tested. A number of correlations between the studied entropy approximations have been detected. The simplest approximation for the excess entropy is provided by the configurational route and by restricting the entropy expansion to translational elements. The contribution of S2,trans to the “exact” Sex derived by thermodynamic integration is much smaller for Lennard-Jones chains than for atomic Lennard-Jones fluids. Nevertheless, S2,trans is in qualitative agreement with Sex. Furthermore, we have identified a linear relation between S2,trans and S3,trans. Both translational entropies differ by roughly a factor of 5. In addition, Sex,SAFT also correlates with S2,trans. An empirical relation between the two quantities has been established. When working in the Rosenfeld formalism, we have verified a correlation between dynamic properties and the excess entropy computed by thermodynamic integration. Again, we have observed two linear regimes. Similar results have been detected when approximating the excess entropy by a configurational method. Only when the excess entropy is derived by SAFT EOS, a single linear correlation has been observed over almost

Figure 12. Correlation between the natural logarithm of the dimensionless viscosity η̃ and (a) the excess entropy Sex predicted by thermodynamic integration, (b) the translational two-body entropy S2,trans, (c) the sum S2,trans + S2,conf, and (d) the excess entropy Sex,SAFT as predicted by the SAFT EOS. The viscosity has been scaled according to the Dzugutov law for Lennard-Jones chains with 20 to 70 monomers. The temperatures and densities chosen correspond to the same setup as employed in the other correlation schemes. The symbols refer to the same chain lengths as, for example, in Figure 8

viscosity for constant n is considerably smaller than the one experienced in the Rosenfeld scaling. Two linear regimes can be identified despite the noise. The crossover between them takes place in the same entropy range where the crossover in the scaled diffusion coefficient was observed. The correlations of ln(BDZ) and βDZ with the chain length are shown in Table 5. The n-dependence of ln(BDZ) and βDZ in both regimes is smaller than in the Rosenfeld case. The above relations are not valid for the 20mer case which thus should be excluded. Figure 12b presents the variation of ln(η̃) as a function of S2,trans. The ln(η̃)−S2,trans correlation agrees qualitatively with the one based on the exact Sex. The fragmentation of the excess entropy into two domains allows the derivation of accurate excess entropy relationships. Nevertheless, the transition between the two domains is rather smooth in this case. The Stokes−Einstein law, i.e βDZ = −αDZ, is found to be valid. The variation of ln(η̃) with S2,trans + S2,conf is shown in Figure 12c. In line with the ln(D̃ COM)−S2,trans + S2,conf correlation, two linear subdomains can be identified. The interchange between them occurs in the same excess entropy interval as encountered

Table 5. Dzugutov Excess Entropy Scaling Parameters for the Reduced Viscosity η̃ = BDZ exp(βDZSex) low entropy regime Sex S2,trans S2,trans+ S2,conf

high entropy regime

βDz

ln(BDz)

βDz

ln(BDz)

70.83/n − 0.59 12.96/n − 1.472 29.35/n − 2.590

169.5/n − 2.010 115.7/n − 2.912 20.80/n − 2.397

87.76/n − 3.628 10.145/n − 4.676 27.30/n − 1.032

310.2/n − 11.21 67.76/n − 3.094 87.76/n − 1.610

L

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

Jayant Singh and Nico van der Vegt are gratefully acknowledged.

all temperature and density values considered in our simulations. This trend, however, differs from the findings in ref 47, which suggested the existence of two linear regimes in this excess entropy scaling. The parameters of the excess entropy correlations for both diffusivity and viscosity depend on the chain length n. They differ strongly from those established for short unentangled chains.23 For the LennardJones chains studied here, we have seen that the preexponential term (AR or BR) in the scaling relation depends stronger on n than the exponential coefficient (αR or βR). The latter obey βR = −αR, at least approximately. This is in accordance with the Stokes−Einstein law. When scaling the diffusion coefficient by the Dzugutov law, two linear regimes have been detected. Such a behavior is found in all schemes to approximate the pair entropy. This diversity has not been detected in recent studies of simple fluids. Despite the statistical fluctuations in the data, the transition between the two linear regimes becomes somewhat sharper when the configurational two-body term S2,conf is included. Viscosity correlations derived in the context of the Dzugutov law have two major advantages over the Rosenfeld one. The variation of the scaled viscosity with temperature and density under constant chain length are much smaller than those found in the case of Rosenfeld. Moreover, the dependence of the constants of the linear correlations is also much weaker than in the Rosenfeld approach. An approximate master curve can be obtained by excluding the 20mer, which is the only system without entanglements. The translational two-body entropy, S2,trans, provides the best excess entropy approximation for the Lennard-Jones chain systems under consideration. It is strongly correlated with the “exact” excess entropy obtained by thermodynamic integration. Moreover, it is very easy to evaluate. The excess entropy based on SAFT is found to be inconsistent with the “exact” excess entropy. Although the obtained scaling relationships of the diffusion coefficient and the viscosity are simpler when Sex,SAFT is used, it is questionable whether the same behavior will be encountered in other polymeric systems with long chains. The present excess entropy scaling which have been derived for Lennard-Jones chains should serve as a guideline for real polymer systems provided that the Coulomb interactions encountered in more polar polymers are not too important.





(1) Rosenfeld, Y. Chem. Phys. Lett. 1977, 48, 467−468. (2) Rosenfeld, Y.; Nardi, E.; Zinamon, Z. Phys. Rev. Lett. 1995, 75, 2490−2493. (3) Malvaldi, M.; Chiappe, C. J. Chem. Phys. 2010, 132, 244502−1− 244502−6. (4) Jabes, B. S.; Chakravarty, C. J. Chem. Phys. 2012, 136, 144507− 1−144507−6. (5) Chopra, R.; Truskett, T. M.; Errington, J. R. J. Phys. Chem. B 2010, 114, 16487−16493. (6) Dzugutov, M. Nature 1996, 381, 137−139. (7) Bretonnet, J.-L. J. Chem. Phys. 2002, 117, 9370−9373. (8) Chopra, R.; Truskett, T. M.; Errington, J. R. J. Chem. Phys. 2010, 133, 104506−1−104506−11. (9) Samanta, A.; Musharaf Ali, S.; Ghosh, S. K. Phys. Rev. Lett. 2001, 87, 245901−1−245901−4. (10) Santos, A. J. Chem. Phys. 2012, 136, 136102−1−136102−2. (11) Cao, Q.-L.; Kong, X.-S.; Li, Y. D.; Wu, X.; Liu, C. S. Phys. B. Condens. Matter 2011, 406, 3114−3119. (12) Kaur, C.; Harbola, U.; Das, S. P. J. Chem. Phys. 2005, 123, 034501−1−034501−6. (13) Das, S. P. Phys. Rev. E 1996, 54, 1715−1719. (14) Galliero, G.; Boned, C.; Fernández, J. J. Chem. Phys. 2011, 134, 064505−1−064505−8. (15) Vaz, R. V.; Magalhãe, A. L.; Fernandes, D. L. A.; Silva, C. M. Chem. Eng. Sci. 2012, 79, 153−162. (16) Harbola, U.; Kaur, C.; Das, S. P. Phys. Rev. Lett. 2003, 91, 229601−229601. (17) Fragiadakis, D.; Roland, C. M. J. Chem. Phys. 2011, 134, 044504−1−044504−3. (18) Lopez, E. R.; Pensado, A. S.; Comunas, M. J. P.; Padua, A. A. H.; Fernandez, J.; Harris, K. R. J. Chem. Phys. 2011, 134, 144507−1− 144507−11. (19) Mittal, J.; Errington, J. R.; Truskett, T. M. Phys. Rev. Lett. 2006, 96, 177804−1−177804−4. (20) Mittal, J.; Errington, J. R.; Truskett, T. M. J. Phys. Chem. B 2007, 111, 10054−10063. (21) Chopra, R.; Truskett, T. M.; Errington, J. R. Phys. Rev. E. 2010, 82, 041201−1−041201−10. (22) He, P.; Liu, H.; Zhu, J.; Li, Y.; Huang, S.; Wang, P.; Tian, H. Chem. Phys. Lett. 2012, 535, 84−90. (23) Goel, T.; Nath Patra, C.; Mukherjee, T.; Chakravarty, C. J. Chem. Phys. 2008, 129, 164904−1−164904−9. (24) Doi,M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press: Oxford, U.K., 1988. (25) Johnson, J. K.; Muller, E. A.; Gubbins, K. E. J. Phys. Chem. 1994, 98, 6413−6419. (26) Kolafa, J.; Nezbeda, I. Fluid Phase Equilib. 1994, 100, 1−34. (27) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300−313. (28) Rosenfeld, Y. Phys. Rev. A 1997, 15, 2545−2549. (29) Kremer, K.; Grest, G. S. J. Chem. Phys. 1990, 92, 5057−5086. (30) Frenkel, D.; Smit, B. Understanding Molecular Simulation. From Algorithms to Applications, Second ed.; Academic Press: San Diego, CA, 2001. (31) Dennis, J. E.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; SIAM Press: Philadelphia, PA, 1987. (32) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19. (33) Sharma, R.; Agarwal, M.; Chakravartya, C. Mol. Phys. 2008, 106, 1925−1938. (34) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Chem. Phys. Lett. 1994, 222, 529−539. (35) Kahaner, D.; Moler, C.; Nash, S. Numerical Methods and Software; Prentice-Hall: Englewood Cliffs, NJ, 1989. (36) Baranyai, A.; Evans, D. J. Phys. Rev. A 1989, 40, 3817−3822. (37) Lazaridis, T.; Karplus, M. J. Chem. Phys. 1996, 105, 4294−4316.

ASSOCIATED CONTENT

S Supporting Information *

Additional data regarding the equilibration of the LennardJones chains and an attempt to fit nonlinear expressions in the excess entropy correlations. This material is available free of charge via the Internet at http://pubs.acs.org



AUTHOR INFORMATION

Corresponding Author

*E-mail: (E.V.) [email protected]. Notes

The authors declare no competing financial interest.



REFERENCES

ACKNOWLEDGMENTS

The present work has been supported by the EU project Nanomodel (No. 211778) and the Deutsche Forschungsgemeinschaft by the Priority Program 1369 Polymer-Solid Contacts: Interfaces and Interphases. Fruitful discussions with M

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

(38) Zielkiewicz, J. J. Phys. Chem. B 2008, 112, 7810−7815. (39) Fushiki, M. Mol. Phys. 1991, 74, 307−319. (40) Barrat, J. L.; Hansen, J. P.; Pastore, G. Phys. Rev. Lett. 1987, 58, 2075−2078. (41) Karakatsani, E. K.; Economou, I. G.; Kroon, M. C.; Peters, C. J.; Witkamp, G. J. J. Phys. Chem. C 2007, 111, 15487−15492. (42) Gross, J.; Sadowski, G. Ind. End. Chem. Res. 2002, 41, 1084− 1093. (43) Gross, J.; Sadowski, G. Ind. End. Chem. Res. 2001, 40, 1244− 1260. (44) Johnson, K. J.; Zollweg, J. A.; Gubbins, K. E. Mol. Phys. 1993, 78, 591−618. (45) Chapman, S.; Cowling, T. G. The Mathematical Theory of NonUniform Gases; Cambridge University Press: Cambridge, U.K., 1970. (46) Das, S. P.; Dufty, J. W. Phys. Rev. A 1992, 46, 6371−6385. (47) Rosenfeld, Y. J. Phys.: Condens. Matter 1999, 11, 5415−5427. (48) Li, G. X.; Liu, C. S.; Zhu, Z. G. J. Non-Cryst. Solids 2005, 351, 946−950. (49) Gnan, N.; Schrøder, T. B.; Pedersen, U. R.; Bailey, N. P.; Dyre, J. C. J. Chem. Phys. 2009, 131, 234505−1−234505−17. (50) Ingebrigtsen, T. S.; Schrøder, T. B.; Dyre, J. C. Phys. Rev. X 2012, 2, 011011−1−011011−20.

N

dx.doi.org/10.1021/ma401617z | Macromolecules XXXX, XXX, XXX−XXX