Communication: Generalizing Rosenfeld's excess-entropy scaling to predict long-time diffusivity in dense fluids of Brownian particles: From hard to ultrasoft interactions

Computer simulations are used to test whether a recently introduced generalization of Rosenfeld's excess-entropy scaling method for estimating transport coefficients in systems obeying molecular dynamics can be extended to predict long-time diffusivities in fluids of particles undergoing Brownian dynamics in the absence of interparticle hydrodynamic forces. Model fluids with inverse-power-law, Gaussian-core, and Hertzian pair interactions are considered. Within the generalized Rosenfeld scaling method, long-time diffusivities of ultrasoft Gaussian-core and Hertzian particle fluids, which display anomalous trends with increasing density, are predicted (to within 20%) based on knowledge of interparticle interactions, excess entropy, and scaling behavior of simpler inverse-power-law fluids.

Dense suspensions of hard-sphere colloids have long served as useful experimental models for exploring how constraints of particle packing affect the properties of condensed phases.^1–3 Nonetheless, many technologically relevant suspensions contain high concentrations of “particles” (e.g., micelles, microgel colloids, dendrimers, or star polymers), which interact with one another through considerably softer effective potentials.^4–6 Such systems are of fundamental interest because they exhibit novel properties that cannot be readily understood in terms of concepts derived from idealized hard-particle systems such as “excluded volume,” “crowding,” and “collision rate,” etc. Examples highlighted in theoretical and computer simulation studies include re-entrant and cluster crystalline phases,^7–17 as well as anomalous dependencies of fluid-phase thermodynamic,^18–22 structural,^20,22–29 and dynamic^{14,18,20–27,29} properties on particle concentration.

To date, much of the theoretical work on these systems has been devoted to developing and testing statistical mechanical approaches for predicting their structural and thermodynamic behavior.^4,30–32 However, some of the focus is beginning to shift toward predicting dynamics. For example, it has recently been shown that generalized Langevin theories^24,33,34 can qualitatively (but not yet quantitatively) predict the anomalous density dependence of the long-time self-diffusivity observed in Brownian dynamics simulations of the Gaussian-core model. Likewise, mode-coupling theory has been able to capture some of the nonmonotonic dynamic trends displayed by fluids of particles that interact via star-polymer-like,²³ harmonic,³⁵ Hertzian,²⁹ or Gaussian-core^27,36 pair potentials. Kinetic theory has also been used to gain insights into the nontrivial temperature- and density-dependent trends in the long-time molecular dynamics of particles with bounded, penetrable-sphere interactions.³⁷ 

Another promising, albeit more heuristic, approach for predicting the dynamic properties of soft-particle fluids is a recently proposed generalization²¹ of the excess-entropy scaling method of Rosenfeld.^38,39 Here, excess entropy s^ex refers to the difference between the entropy per particle of the fluid and that of an ideal gas of particles with the same number density ρ. Excess entropy is a negative quantity, and its magnitude characterizes the extent to which static interparticle correlations—present due to interparticle interactions— reduce the number of microstates accessible to the fluid. Although a rigorous statistical mechanical derivation of excess entropy scalings for dynamics has proven to be elusive,^21,39–41 the qualitative expectation—corroborated by data from both computer simulations and experiments (see, e.g., Refs. 39 and 42–48)—is that changes to a fluid that strengthen its interparticle correlations also result in slower dynamic relaxation processes.

In the generalized Rosenfeld (GR) scaling approach,²¹ information about interparticle interactions is used to recast the species’ long-time tracer diffusivities as dimensionless (GR-reduced) combinations that are, by construction, single-valued functions of s^ex in the low ρ limit. For fluids of particles that interact solely via an inverse-power-law pair (IPL) repulsion, one can easily show that GR-reduced diffusivities remain single-valued functions of s^ex for all values of ρ and temperature T. Interestingly, Krekelberg et al. have further discovered that GR-reduced tracer diffusivities obtained from molecular dynamics simulations of more complex model systems—equilibrium fluid mixtures of (additive or nonadditive) hard-sphere or Gaussian-core particles at different values of ρ, T, and composition—approximately follow a similar excess-entropy based scaling relation.²¹ From a practical viewpoint, this means that knowledge of the interparticle interactions and the thermodynamic excess entropy for these model fluids allows one to predict key aspects of their long-time dynamic behavior.

Can the GR scaling method be extended to treat systems with different types of dynamics (e.g., with equations of motion that incorporate dynamic effects of “solvent” surrounding the particles)? In this work, we take a step toward addressing this question. Specifically, we present simulation data that tests whether the GR scaling method can be recast in a manner useful for predicting long-time tracer diffusivities of suspended particles undergoing Brownian (i.e., overdamped Langevin) dynamics in the absence of interparticle hydrodynamic forces. The model fluids considered here comprise particles that interact via pair potentials

$V(r)$ of the IPL [

$\mathcal {V}(r)=\epsilon (\sigma /r)^{\mu }$

$V(r)=ε(σ/r)μ$], Gaussian-core [

$\mathcal {V}(r)=\epsilon \exp \lbrace -(r/\sigma )^2\rbrace ]$

$V(r)=εexp{−(r/σ)2}]$⁠, and Hertzian [

$\mathcal {V}(r)=\epsilon (1-r/\sigma )^{5/2}$

$V(r)=ε(1−r/σ)5/2$] forms. Here ε and σ represent characteristic energy and length scales of these interactions, respectively, and μ determines the “softness” of the IPL repulsion. As we discuss below, a key result from our study is that GR-reduced tracer diffusivities, sampled across a wide range of T and ρ for these models, approximately collapse onto a single curve when plotted versus s^ex. A consequence is that the tracer diffusivities of the ultrasoft Gaussian-core and Hertzian particle fluids, which display anomalous trends with increasing ρ,^24,26 can be estimated (to within 20% relative error)—with no adjustable parameters—from knowledge of (1) the interparticle interactions, (2) the value of s^ex at the state point of interest, and (3) the GR scaling behavior of simpler IPL fluids.

For the case of Brownian dynamics considered here, it is convenient to express the GR-reduced long-time diffusivity as²¹ 

where B is the second virial coefficient, β ≡ 1/k_BT, k_B is Boltzmann's constant, and D₂ = −(∂[D/D₀]/∂ρ)_T|_{ρ = 0} quantifies how pair interactions at low particle concentration modify the tracer diffusivity D relative to the infinite dilution value D₀. As has been shown elsewhere,²¹ (1 − D/D₀)^GR = −s^ex/k_B for densities low enough that

$O(ρ2)$ quantities can be neglected.

Based on Eq. (1), it is evident that four quantities must be computed to test the GR excess-entropy scaling relation for state points across the T-ρ plane: B, D₂, D/D₀, and s^ex. The first two depend solely on β and

$V(r)$ and can be calculated from simple theoretical considerations. Specifically, B can be obtained by integrating the following expression numerically:

$B=-(1/2) \int d {\bf r} \lbrace \exp [-\beta \mathcal {V}(r)] -1\rbrace$

$B=−(1/2)∫dr{exp[−βV(r)]−1}$⁠. Likewise, D₂ can be computed via a one-dimensional numerical integration,

Equation (2) is derived by the relaxation method, an approach outlined in detail elsewhere.^49–52 The function Q(r) is calculable from knowledge of

$V(r)$ via⁶¹ 

The physical significance of Q(r) is that it characterizes the extent to which application of a small external force f to a particle in the low-ρ fluid would distort its pair correlation function with its neighbors,

$g(r)=exp[−βV(r)]{1+βQ(r)r·f/2r}$⁠, and hence modify its effective mobility. One boundary condition on Q(r)—valid for all

$\mathcal {V}(r)$

$V(r)$ considered here—is that the force-induced distortion of the pair correlation function must decay at large distances from the central particle, Q(∞) = 0. The second boundary condition depends on the form of the short-range interparticle repulsion. For interactions studied previously with a hard-sphere repulsion at a separation r = σ,^49–52 the following condition holds: (dQ/dr)|_{r = σ} = −1. We have verified that this condition also accurately describes the behavior of particles with softer IPL repulsions if it is applied at an effective hard-sphere “Boltzmann”⁵³ diameter

$r=\sigma _{\text{B}}$

$r=σB$ defined by

$\mathcal {V}(\sigma _{\text{B}})= 10\, k_{\mathrm{B}} T$

$V(σB)=10kBT$⁠. To our knowledge, this work is the first to consider D₂ for ultrasoft potentials like the Gaussian-core or Hertzian models that are bounded at zero interparticle separation. From Eq. (3), it is apparent that the condition Q(0) = 0 applies for such models.

To compute D/D₀, we use Brownian dynamics simulations, where particle translations are governed by the Langevin equation, solved in the overdamped limit. The algorithm we employ, wherein the position r_i of each particle i ∈ [1, N] updates at each time step according to

$ri(t+Δt)=ri(t)−D0Δt∇∑j≠iNβV(rij(t))+Δri$⁠, is detailed elsewhere.^54,55 The quantity Δr_i denotes a random displacement due to solvent collisions; its magnitude is taken to be Gaussian-distributed with a mean of zero and a variance of 6D₀Δt.

The Brownian dynamics simulations track N particles in a periodically-replicated cubic simulation cell whose volume V sets ρ = N/V. We use N = 1000, 2000, and 3000 particles for the IPL, Hertzian, and Gaussian-core fluids, respectively. Interparticle potentials were truncated at r = 1.4 σ– 3.9 σ,⁶²σ, and 3.71 σ, respectively, for these systems. For simulations of the IPL fluids, we set D₀ = 0.001(σ²ε/m)^1/2 and the Δt = 0.025τ_B, where τ_B ≡ mD₀/k_BT is a characteristic Brownian time scale. For simulations of Gaussian-core and Hertzian particle systems, we set D₀ = 0.0001(σ²ε/m)^1/2 and Δt = 0.1τ_B. Long-time tracer diffusivities are computed from the mean-squared displacements via the relation, D = 〈Δr²〉/6t as t → ∞.

We employ free-energy-based simulation techniques to evaluate s^ex. In the first step, we obtain the ρ dependence of the Helmholtz free energy at high T using grand canonical transition matrix Monte Carlo simulation.⁵⁶ In the second step, we perform a canonical temperature expanded ensemble simulation⁵⁷ with a transition matrix Monte Carlo algorithm⁵⁸ to evaluate the change in Helmholtz free energy with T at constant ρ. Collectively, these simulations provide values for excess Helmholtz free energy and excess energy, and hence s^ex, at a state point of interest. Interested readers can find additional details in our earlier papers.^59,60

We now examine the simulated behavior of six IPL fluids with exponents of μ = 6, 8, 10, 12, 18, and 36. Specifically, Fig. 1(a) shows how (1 − D/D₀) and its GR-reduced counterpart (1 − D/D₀)^GR vary with density ρσ³ and −s^ex/k_B, respectively, for these fluids at T = ε/k_B. As should be expected, systems with softer interparticle repulsions (i.e., lower μ) display considerably faster single-particle dynamics at a given ρσ³. For example, the μ = 6 fluid must realize a density twice that of the μ = 36 fluid for both systems to exhibit D/D₀ = 0.2. In contrast, (1 − D/D₀)^GR for these fluids shows only minor systematic variation when plotted versus −s^ex/k_B, with less than 20% difference between the μ = 6 and μ = 36 values.

To put to the GR scaling relationship for Brownian dynamics to a more stringent test, we further analyze the behaviors of model fluids with ultrasoft Hertzian and Gaussian-core interactions [see Figs. 1(b) and 1(c), respectively]. In particular, we explore a broad range of ρσ³ for these two systems along isotherms spanning from the moderately supercooled liquid (T ⩽ 0.004ε/k_B) to temperatures high enough (T ⩾ 0.01ε/k_B) for the models to exhibit mean-field-like behavior.^4,5,7,26 Note that this data set includes state points for which D/D₀ increases with increasing ρσ³ in both systems, behavior that is anomalous when compared to the trends exhibited by simpler atomic and colloidal fluids. The main conclusion to be drawn from Figs. 1(b) and 1(c) is that, despite the anomalous relationship between single-particle dynamics and ρ in these systems, the GR excess entropy scaling relation looks strikingly similar to the quasi-universal form of the IPL fluids. In fact, the Hertzian and Gaussian-core data deviate by less than 20% from the least-squares fit to the (1 − D/D₀)^GR versus s^ex data from Fig. 1(a).

One consequence of this quasi-universal scaling behavior is that one can estimate Brownian dynamics D/D₀ data, even for nontrivial systems with ultrasoft potentials, simply from knowledge of their interparticle interactions, the excess entropy, and, e.g., a least-squares fit of the IPL data of Fig. 1(a). Figure 2 illustrates the quality of the predictions made by this approach for the Hertzian and Gaussian-core model fluids. As can be seen—with no adjustable parameters—the estimated values of diffusivity are accurate for low-to-moderate ρ and remain qualitatively reliable for higher ρ where the anomalous trends emerge. Collectively, these results illustrate a strong apparent link between static interparticle correlations and long-time dynamics for model fluids with either hard or ultrasoft interactions.

Finally, whether the approach outlined here can provide an accurate estimate to the “thermodynamic” component of diffusivity for systems which also have strong hydrodynamic interparticle contributions to dynamics is an open question, which deserves attention in future studies. Furthermore, it will be interesting to test whether the GR scaling for diffusivity in Brownian dynamics also holds for systems with more complex interparticle interactions (e.g., with more than one characteristic length scale).

We thank Professor David S. Dean for suggesting the relaxation method for computing D₂. T.M.T. acknowledges support of the Welch Foundation (F-1696) and the David and Lucile Packard Foundation. J.R.E. acknowledges financial support of the National Science Foundation, Grant No. CBET-0828979. M.J.P. acknowledges the support of the Thrust 2000—Harry P. Whitworth Endowed Graduate Fellowship in Engineering. The Texas Advanced Computing Center (TACC), the University at Buffalo Center for Computational Research, and the Rensselaer Polytechnic Institute Computational Center for Nanotechnology Innovations provided computational resources for this study.