Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Phys Chem B. Author manuscript; available in PMC 2010 October 22.
Published in final edited form as:
PMCID: PMC2766175

Using compressibility factor as a predictor of confined hard-sphere fluid dynamics


We study the correlations between the diffusivity (or viscosity) and the compressibility factor of bulk hard-sphere fluid as predicted by the ultralocal limit of the barrier hopping theory. Our specific aim is to determine if these correlations observed in the bulk equilibrium hard-sphere fluid can be used to predict the self-diffusivity of fluid confined between a slit-pore or a rectangular channel. In this work, we consider a single-component and a binary mixture of hard spheres. To represent confining walls, we use purely reflecting hard walls and interacting square-well walls. Our results clearly show that the correspondence between the diffusivity and the compressibility factor can be used along with the knowledge of the confined fluid's compressibility factor to predict its diffusivity with quantitative accuracy. Our analysis also suggests that a simple measure, the average fluid density, can be an accurate predictor of confined fluid diffusivity for very tight confinements (≈ 2-3 particle diameters wide) at low to intermediate density conditions. Together, these results provide further support for the idea that one can use robust connections between thermodynamic and dynamic quantities to predict dynamics of confined fluids from their thermodynamics.


The behavior of fluids in small spaces, with one of the characteristic dimensions of the order of few fluid particle diameters, can be significantly different from the bulk. Some of the well-studied features of such confined fluids are the formation of layering structure normal to the confining walls,1 inhomogenous pressure equation of state,1 appearance of surface induced phase transitions,2,3 and shift in the glass transition temperature.4,5 The scientific progress on some of the issues related to the behavior of confined fluids has been biased towards understanding the structural changes and associated thermodynamic properties due to confinement. In fact, now we have theories at our disposal which can predict the implications of confinement on the thermodynamic properties quantitatively over a wide range of conditions.6 On the contrary, the situation is not so satisfactory for predicting the dynamic behavior (such as the average self-diffusivity) of confined fluids.79 One of the reasons for this slow advancement on understanding fluid dynamics under confinement is the lack of tractable kinetic theories even for the bulk fluid behavior over the whole equilibrium fluid range.

Recent studies have utilized a different strategy to predict the transport properties of fluids under confinement; their aim is to find empirical correlations between thermodynamic and dynamic quantities which remain approximately unchanged due to confinement.1012 The thermodynamic quantity used to find such correlations can be as simple as the average fluid density,10 although significantly better results are obtained by using quantities such as the excess entropy10,11,13,14 of the fluid with respect to the ideal gas state and the available volume.12,15 The idea is to utilize these correlations obtained from the bulk fluid behavior along with the information about the thermodynamic property of confined fluids (which can be easily obtained via theory, molecular simulation, or laboratory experiment) to calculate the dynamic quantity with quantitative accuracy. This line of enquiry seems rather promising from the studies so far on equilibrium fluids as one can predict the effect of confinement on the self-diffusivity (or alternatively shear viscosity which follows trivially from the Stokes-Einstein relationship for equilibrium fluids) for a wide range of fluid models such as hard-sphere, Lennard-Jones (LJ), square well and others by utilizing “excess entropy-diffusivity”10,11 and “available volume-diffusivity”12 bulk relationships.

The possible explanation for the observed correlations between these quantities of seemingly different origin (thermodynamic versus dynamic) even in the bulk has mostly been qualitative and there is no first principle based justification available. It has been argued that excess entropy should naturally be a quantity which will track the fluid dynamics as it characterizes the reduced number of states due to interparticle correlations.16 But it is not clear what the functional relationship should be17,18 and if there is a universal connection or something more like a parameteric fit with fitting coefficients dependent on the state of the fluid (density, temperature) and the type of fluid (strong versus fragile). Whether these correlations persist in the supercooled state is also somewhat an open question19 with previous work pointing towards the need for an additional adjustable parameter (as a function of density in case of binary Lennard-Jones20 and SPC/E water21).

In this regard, recent results from the single particle barrier hopping theory of Schweizer and co-workers22 may be an important step towards bridging the gap between observations and fundamental understanding. One specific result out of many predictions made by this theory is that the relaxation time (or inverse diffusivity D) for the bulk hard-sphere fluid should scale with the compressibility factor Z = βP/ρ as,23


where β = [kBT]−1, kB is the Boltzmann constant, T is temperature, P is pressure, ρ is the fluid density, D0 is an elementary diffusion scale, a is a numerical prefactor, [var phi] = πρσ3/6 is the fluid packing fraction, and σ is the fluid diameter. A qualitative agreement with the predictions of eq 1 was already found in a previous work on the dynamics of polydisperse hard-sphere fluid which showed that the relaxation time as a function of Z for several different polydispersity values can be approximately collapsed onto a single curve.24 A quantitative comparison of eq 1 with the actual diffusion data will provide a useful starting point for understanding some of the outstanding issues related to the dynamics of supercooled fluids. As the focus of this paper is the equilibrium fluid dynamics under confinement, we will only touch upon this issue very briefly in the beginning of results and discussion section. Detailed results and the discussion on the supercooled fluid behavior from the perspective of eq 1 will be provided in a future publication.

Here, instead we focus on the possiblity of DZ relationship as observed in the bulk for an equilibrium hard-sphere fluid as a predictor of diffusion under confinement by utilizing the knowledge about its Z. For confined systems, pressure is a tensorial quantity (and therefore Z) and for a slit-pore system (fluid confined between walls in one direction) pressure has two components, parallel and normal to the walls. For average diffusion parallel to the walls (which will be the main dynamic quantity of interest in this work), we expect that the compressibility factor based on the average pressure component in the parallel direction will be relevant. As our fluid models, we consider both a single-component and a binary mixture of hard-sphere particles. For confinement models, we consider slit-pore geometry with hard walls and attractive square-well walls and rectangular channels with hard walls. We find that the bulk DZ relationship remains essentially unchanged due to these different types of confinement and can be successfully used to predict the diffusivity of confined fluids with quantitative accuracy for most of the state points. We further notice that any expected deviations from the exact diffusion values are generally negative (under-prediction) for low to intermediate density states and positive (overprediction) for high density states. Data from previous studies show that qualitatively similar deviations are expected from the excess entropy - diffusivity and available volume - diffusivity correlations.12 An attractive aspect of our observed correlation between the compressibility factor and diffusivity, as compared to previous work on other thermodynamic measures, is that one can calculate Z fairly easily using colloidal hard sphere experiments in bulk (and possibly under confinement with appropriate modifications in the analysis)25 and in molecular simulations.

Simulation details

We have used discontinuous molecular dynamics (DMD) simulations26 to track various thermodynamic and dynamic properties of the hard-sphere (HS) fluid model. The DMD simulations each involved N = 10000 HS particles contained within a simulation cell of volume V = Hx × Hy × Hz, where Hi is the box length in a given direction (i = x,y,z). Periodic boundary conditions were applied in all three directions for bulk and in two directions (x and y) for the slit-pore, and in one direction (x) for the rectangular channel. The non-periodic directions in case of confined fluids can be either smooth “hard” walls which are perfectly reflecting or “attractive” walls with a square-well potential between the fluid particles and the confining walls,


where z represents the shortest distance between a given particle center and the wall of interest, and εw is the strength of the effective particle-wall interaction. We only consider hard walls in case of rectangular channels.

The diffusivity D of the fluid was obtained by fitting the long-time behavior of the average mean-squared displacement of the particles to the Einstein relation Δrd2=2dDt, where Δrd2 corresponds to the mean-square displacement per particle in the d periodic directions after time t. To calculate the compressibility factor Z, we use the following form of the virial theorem,27,28


where ri j is the position vector between particles i and j, vi is the change in velocity of particle i after the collision, N is the number of particles, and the summation is over all the collisions during the elapsed time te. In case of confined systems, the vector product includes only the periodic components of ri j and vi for Z parallel to the walls and non-periodic components for Z normal to the walls. We test our Z calculations against the GC-TMMC simulation results29,30 and find favorable comparison for both bulk and confined systems (less than 0.5%) difference.13,14 To simplify the notation in this paper, we have implicitly non-dimensionalized all quantities by appropriate combinations of a characteristic length scale (which we take to be the HS particle diameter σ) and time scale (which we choose to be σmβ, where m is particle mass). As a result, all quantities with dimensions of energy are understood to be “per kBT”, the only energy scale in the problem.

Results and discussion

Bulk hard-sphere fluid

First we check the predictions of eq 1 for the bulk HS fluid. Although the barrier hopping theory in its ultralocal limit is expected to be valid for strongly supercooled liquids,23 we want to see if this can also provide some reasonable predictions for weakly supercooled and equilibrium fluid states. As we are using a monatomic fluid which is hard to supercool and will crystallize readily, it is not possible to access strongly supercooled regime. Figure 1 shows the D data plotted against [var phi] (top panel) as well as (Z−1)2/[var phi] (bottom panel) as expected from eq 1 for [var phi] ≤ 0.53. We are able to fit the data for most of the dense fluid range (D < 0.1 and [var phi] > 0.40) to eq 1. Equation 1 is only valid for fluid states above the ideal glass transition packing fraction ([var phi]i = 0.432) at which the activated barrier hopping dynamics becomes important22 and the deviations seen in Figure 1 at low density are therefore expected. We also show the inverse viscosity 1/η on the same plot and the expected fit from the Stokes-Einstein relationship, i.e., = 1/2π in the bottom panel. To verify that the fluid systems between [var phi] = 0.49 to [var phi] = 0.53 are not phase separated (as coexistence between the fluid-solid phase for the HS model is between [var phi] = 0.494−0.545),31 we check the time evolution of Z as well as several other order parameters, e.g., bond-orientationl order Q6,32 pair-correlation function g(r), translational order parameter s233 and do not find any evidence of phase separation at the time scale of the simulations. Of course, if we run the simulations for long enough, these state points will show fluid and solid phases in coexistence. We use a similar approach to identify stable and metastable (at the simulation timescales) state points for the confined HS system for which the location of the fluid-solid phase boundary may not be known apriori. In future we will investigate whether the trends observed here also hold in the deeply supercooled state by utilizing a polydisperse HS fluid to avoid crystallization.

Figure 1
Bulk data. (Top panel) Diffusivity D and inverse viscosity 1/η data are plotted as a function of [var phi]. (Bottom panel) Diffusivity D and inverse viscosity 1/η data are plotted as a function of (Z−1)2/[var phi] along with a ...

Next, we want to examine if the correspondence between D and Z observed in the bulk is unchanged due to confinement and can therefore be used as a predictive tool for confined fluid diffusivity. To facilitate this process, we fit our bulk data to an analytic functional form34 D = f (Z) as shown in the inset of Figure 1. To simplify the analysis we have used DZ relationship but D−(Z−1)2/[var phi] can also be used for rest of the analysis in this paper. Further, we plot the ratio of the actual DMD data and f (Z) to assess the reliability of our fit, and find overall good agreement with less than 5% error for Z > 5. From now onwards we will use f (Z) as a functional representation of our bulk D data.

Hard-sphere fluid confined between hard walls

First, we test the effect of hard-wall confinement on the DZ relationship. Figure 2 shows the data for bulk HS (solid line), HS confined between slit-pores of size H = Hz, and rectangular channels with dimensions Hy × Hz. Overall, the collapse of confined fluid data on the bulk curve is good suggesting that f (Z) can be used as a predictor of diffusivity if Z for the confined fluid is known. For more restrictive confinements, e.g., in rectangular channels and for a slit-pore of size H = 4, DZ curves have a different slope from the bulk curve which gives rise to underprediction in D from f (Z) at intermediate Z and overprediction in D at high Z as shown in the inset of Figure 2. The Z-range (or alternatively [var phi]-range) considered here covers most of the equilibrium state points as known from the exact fluid-solid phase diagram.35

Figure 2
Hard wall confinement. Diffusivity D versus compressibility factor Z for hard-sphere fluid in bulk (line) and confined between hard-walls (symbols). The ratio D/f (Z) as a function of Z is also shown in the inset.

Hard-sphere fluid confined between attractive walls

To test if DZ correspondence is still valid when we introduce fluid-wall interactions, we calculate these quantities for systems with varying attraction strength εw (eq 2). Figure 3 shows the bulk HS data (line) along with the confined fluid data (symbols) which demonstrates that indeed the DZ relationship remains unchanged. To assess quantitative differences, we plot D/f (Z) in the inset of Figure 3 which shows that difference between the bulk and confined diffusivity based on Z is less than 10% for 75% of the state points and always less than 15%.

Figure 3
Attractive wall confinement. Diffusivity D versus compressibility factor Z for hard-sphere fluid in bulk (line) and confined between attractive walls represented by a square-well potential eq. 2 (symbols). The ratio D/f (Z) as a function of Z is also ...

Binary mixture of hard-sphere fluid confined between hard walls

As most of the fluid systems of practical interest are not single-component but rather a mixture of different sizes, it will be instructive to investigate if the DZ relationship can also be useful for predicting the diffusivity of equilibrium confined fluid mixtures. We have simulated a binary mixture of hard spheres with particle size ratio given by σ21 = 0.75 and equal mass (m2 = m1). We present data for two mixture compositions, equal particle numbers (N1 = N2 = 5000) and equal packing fraction (N1 = 2967, N2 = 7033; [var phi]1 = [var phi]2). Figure 4 shows the diffusivity Di for both the components (i = 1 (black lines and symbols) and i = 2 (red lines and symbols)) as a function of Z in bulk (equal particle numbers: solid lines, equal packing fraction: dashed lines) and confined in a slit-pore of size H = 5 (equal particle numbers: filled symbols, equal packing fraction: empty symbols). Again, we find a one-to-one correspondence between the bulk and confined DZ relationship thereby suggesting that Z can be a reliable predictor of confined fluid diffusivity for fluid mixtures. We fit individual component data for the equimolar mixture in the bulk to an analytic functional form34 Di = fi(Z) to facilitate comparison with the confinement data. This will also help in assesing the dependence of observed correlations on the mixture composition. To quantify the differences between the actual DMD data and predictions based on the bulk fi(Z) and the confined fluid Z, we plot Di/fi(z) in the inset of Figure 4. It is clear that the confinement data for an equimolar mixture is quantitatively described by the bulk relationship based on the same mixture composition but the prediction for a different composition (equal packing fraction) is only semi-quantitatively described at high Z (or alternatively high [var phi]).

Figure 4
Binary hard-sphere mixture. Diffusivity Di versus compressibility factor Z for a binary mixture of hard spheres in bulk (equal particle numbers: solid lines, equal packing fraction: dashed lines) and confined between hard walls (equal particle numbers: ...

Frustrated dynamics due to imperfect layering

It was observed in a previous work that the formation of layering structure normal to the walls can actually facilitate the motion in both parallel14 and normal15 directions. Similarly, frustrated layering formation due to pore-sizes which cannot accomodate an integer number of layers can impede the particle motion or its diffusivity. These observations were utilized to propose a possible route to control mobility in small channels by modifying the wall-particle interactions to impose a predetermined density profile and thereby achieve the required D.36 For a HS fluid confined between hard walls, it was found that the fluctuations in D as a function of H for a given packing fraction can be explained by a similar change in the excess entropy of the fluid.14 But the excess entropy in this case was only able to provide a semi-quantitaive measure of dynamics as seen by an imperfect collapse of confinement data on the bulk diffusivity-excess entropy curve at high [var phi]. Later, it was shown that the available volume can provide a better description of this data and make improved predictions.12

Here, we test if DZ correlations can describe the layering induced fluctuations in confined fluid D qualitatively and further quantitatively. Figure 5 shows the diffusivity data for HS between hard walls as a function of varying H along various isochores. Symbols are the actual simulation data and solid lines are the predictions based on bulk data (f (Z)) and the confined fluid Z. Overall, we find that the predictions match well with the actual data, especially for high [var phi] the agreement is remarkable as opposed to the predictions based on excess entropy.14 Moreover, low to intermediate [var phi] data can actually be described simply based on the average packing fraction as it matches the bulk diffusivity at the same [var phi] (shown by horizontal dashed lines) very well.

Figure 5
Layering induced dynamics. Diffusivity D versus slit-pore size H for hard-sphere fluid confined between hard walls. The molecular dynamics data is shown by symbols and the predictions based on bulk diffusivity-compressibility factor relationship is shown ...


We have shown that the relationship between a thermodynamic and a dynamic quantity, i.e., between “diffusivity and compressibility factor” can be used to predict the effect of confinement on fluid diffusivity. Along with the knowledge of bulk DZ relationship, the information about the confined fluid's compressiblity factor is needed which can be easily estimated from experiments or approximated using existing theoretical approaches such as the density functional theory. This provides an indirect but quantitative route to predict confined fluid diffusivity based on its thermodynamics. One of the interesting features about the observed DZ correlations is that the barrier hopping theory for supercooled fluid dynamics in its ultralocal limit predicts such an expectation for a bulk hard-sphere fluid and a precise functional form for these correlations. We find that this expected functional dependency is satisfied for a single-component hard-sphere fluid in a weakly supercooled state and even in the equilibrium fluid state. In future work, we will focus on understanding the DZ correlations for supercooled fluids in bulk and under confinement by utilizing polydisperse fluids to thwart crystallization and size segregation. We are currently also investigating if DZ correlations can be used to predict the effect of confinement on dynamics for attractive fluids with square-well or Lennard-Jones interactions.


It is my great pleasure to contribute this paper for an issue honoring Prof. Ted Davis as his work on confined fluids has been a great source of inspiration for me. I also acknowledge the kind hospitality of Dr. Robert Best and the Department of Chemistry during my stay at the University of Cambridge where a part of this paper was written. I thank Dr. Artur Adib for supporting a postdoctoral fellowship. This research was supported by the Intramural Research Program of the NIH, NIDDK. This study utilized the high-performance computational capabilities of the Biowulf PC / Linux cluster at the National Institutes of Health, Bethesda, MD ( and the NSF TeraGrid resources provided by TACC.


1. Davis HT. Statistical Mechanics of Phases, Interfaces, and Thin Films. VCH; 1996.
2. Gelb LD, Gubbins KE, Radhakrishnan R, Sliwinska-Bartkowiak M. Rep Prog Phys. 1999;62:1573.
3. Sellers MS, Errington JR. J Phys Chem C. 2008;112:12905.
4. Alcoutlabi M, McKenna GB. J Phys: Condens Matter. 2005;17:R461.
5. Mittal J, Shah P, Truskett TM. J Phys Chem B. 2004;108:19769.
6. Hansen JP, McDonald IR. Theory of Simple Liquids. 3rd. Academic Press; 2006.
7. Magda JJ, Tirrell MV, Davis HT. J Chem Phys. 1985;83:1888.
8. Vanderlick TK, Davis HT. J Chem Phys. 1987;87:1791.
9. Krishnan SH, Ayappa KG. J Phys Chem B. 2005;109:23237. [PubMed]
10. Mittal J, Errington JR, Truskett TM. Phys Rev Lett. 2006;96:177804. [PubMed]
11. Mittal J, Errington JR, Truskett TM. J Phys Chem B. 2007;111:10054. [PubMed]
12. Goel G, Krekelberg WP, Pond MJ, Mittal J, Shen VK, Errington JR, Truskett TM. J Stat Mech. 2009:P04006.
13. Mittal J, Shen VK, Errington JR, Truskett TM. J Chem Phys. 2007;127:154513. [PubMed]
14. Mittal J, Errington JR, Truskett TM. J Chem Phys. 2007;126:244708. [PubMed]
15. Mittal J, Errington JR, Truskett TM, Hummer G. Phys Rev Lett. 2008;100:145901. [PubMed]
16. Dzugutov M. Nature. 1996;381:137.
17. Rosenfeld Y. Phys Rev A. 1977;15:2545.
18. Rosenfeld Y. J Phys: Condens Matter. 1999;11:5415.
19. Dzugutov M. Phys Rev E. 2002;65:032501. [PubMed]
20. Mittal J, Errington JR, Truskett TM. J Chem Phys. 2006;125:076102. [PubMed]
21. Mittal J, Errington JR, Truskett TM. J Phys Chem B. 2006;110:18147. [PubMed]
22. Schweizer KS, Yatsenko G. J Chem Phys. 2007;127:164505. [PubMed]
23. Schweizer KS. J Chem Phys. 2007;127:164506. [PubMed]
24. Sear RP. J Chem Phys. 2000;113:4732.
25. Dullens RPA, Arts DGAL, Kegel WK. Proc Natl Acad Sc. 2006;103:529. [PubMed]
26. Rapaport DC. The Art of Molecular Dynamics Simulation. 2nd. Cambridge University Press; 2004.
27. Alder BJ, Wainwright TE. J Chem Phys. 1959;31:459.
28. Denlinger MA, Hall CK. Mol Phys. 1990;71:541.
29. Errington JR. J Chem Phys. 2003;118:9915–9925.
30. Errington JR, Shen VK. J Chem Phys. 2005;123:164103. [PubMed]
31. Anderson VJ, Lekkerkerker HNW. Nature. 2002;410:811. [PubMed]
32. Steinhardt P, Nelson DR, Ronchetti M. Phys Rev B. 1983;28:784.
33. Truskett TM, Torquato S, Debenedetti PG. Phys Rev E. 2000;62:993. [PubMed]
34. f (Z) = (a + bZ + cZ2 + dZ3)/(1+eZ + fZ2+gZ3), where a, b, c, d, e, f, g are constants.
35. Fortini A, Dijkstra M. J Phys: Condens Matter. 2006;18:L371. [PubMed]
36. Goel G, Krekelberg WP, Errington JR, Truskett TM. Phys Rev Lett. 2008;100:106001. [PubMed]