An EQT-cDFT approach to determine thermodynamic properties of confined fluids

Fluids confined in a narrow pore or channel, especially at length scales ranging from a few angstroms to micrometers, exhibit interesting properties¹ and have many practical applications, such as nanofiltration,^2–6 drug delivery,^7,8 enhanced oil recovery,^9,10 heat management,^11,12 and material synthesis.¹³ In confinement, due to the spatial constraints and competition between the surface-fluid and fluid-fluid atomic interactions, the fluid properties can be significantly different from the bulk fluid. For example, the pressure experienced by a fluid confined in a nanopore can be orders of magnitude higher than the pressure in the bulk fluid.¹⁴ Such unusual behavior of confined fluids can have practical applications, e.g., a carbon nanotube (CNT) can be used as a super-compressor for the synthesis of valuable high-pressure materials, such as KI nanocrystals.¹³ Therefore, study of confined fluids is important to get atomic-level insights into their unusual properties and to enable the design of novel nanofluidic applications.

Study of confined fluids at atomic-level by experimental methods can be difficult and expensive, therefore, we need theoretical and computational methods. Such methods also help explain experimental results and establish the relations between various properties of the fluid. Molecular simulation techniques, such as molecular dynamics (MD) and Monte Carlo (MC) simulations, can be used to study the microscopic properties of confined fluids.^15,16 However, these techniques are computationally expensive, and accessible time and length scales are limited. On the other hand, the classical continuum methods, such as Navier-Stokes equations, are computationally efficient, but fail to accurately describe atomic-level structure and properties of confined fluids.¹⁷ As a result, molecular simulations and classical continuum methods are not applicable for applications that involve multiple time and length scales ranging from the quantum to atomic to continuum scales.¹⁸ To address these issues, we need a multiscale method which is not only as accurate as molecular simulations, but also as fast as classical continuum methods.

An empirical potential-based quasi-continuum theory (EQT)^19–23 is a multiscale theory, that seamlessly integrates interatomic potentials describing various atomic interactions into the continuum framework, to accurately predict the equilibrium density and potential profiles of confined fluids. On the one hand, EQT captures the atomic details, and hence, it is more rigorous than the classical continuum methods. On the other hand, it is based on a continuum framework, therefore, it is computationally faster than molecular simulations.

In this work, we focus on an application of the EQT framework to determine the equilibrium thermodynamic properties of a confined fluid. For a fluid confined to a slit-like channel, the fundamental thermodynamic relation is²⁴ 

where Ω is the grand potential, S is the entropy, T is the temperature, P is the pressure, V is the volume, N is the number of the fluid particles, μ is the chemical potential, γ is the wall-fluid surface tension, A is the surface area of the wall, L is the channel width, and f_S is the solvation force. Given the thermodynamic state of a confined fluid defined by T, V, and μ, we can determine the thermodynamic properties of the fluid by evaluating the gradients of Ω. For example, from Eq. (1), the wall-fluid surface tension can be computed as $γ= 1 2 ∂ Ω ∂ A V , T , L , μ$⁠. Therefore, to determine the thermodynamic properties of a confined fluid, we need an expression for Ω. Classical density functional theory (cDFT)^24–26 provides a framework to determine Ω in terms of the equilibrium density and interaction energy obtained from EQT.

In this work, first, we use EQT to model an expression for Ω in terms of the equilibrium density and potentials. Then, we can use the EQT-cDFT-based grand potential functional to compute various thermodynamic properties of confined fluids. The remainder of the paper is organized as follows. In Sec. II, we describe the theory of the EQT-cDFT approach and methods to compute the thermodynamic properties, such as adsorption, local pressure tensor, surface tension, and solvation force. In Secs. III and IV, we demonstrate the EQT-cDFT approach by simulating Lennard-Jones (LJ) fluids, namely, methane and argon, confined inside slit-like channels of graphene, and compare the results from the EQT-cDFT with the MD simulations. Finally, we draw conclusions in Sec. V.

In EQT, for a slit-channel system, the 1-D steady-state Nernst-Planck (NP) equation,

with boundary conditions

is solved to obtain the self-consistent density and potential profiles of a confined fluid. In Eqs. (2) and (3), ρ is the fluid density, U is the total potential, T is the fluid temperature, R is the ideal gas constant, L is the channel width, ρ_avg is the average density of the confined fluid, and z is the direction normal to the slit-channel walls. The principal idea in EQT is to compute U by a continuum approximation. In the continuum approximation, wall and fluid particles are represented by their local densities, and U is determined by summing the contributions from the wall and fluid as a density weighted integration of the effective quasi-continuum interaction potentials over the surrounding medium, i.e.,

where r and r′ are the position vectors, r = |r − r′|, ρ_w is the wall atom density, u^wf(r) and u^ff(r) are the effective wall-fluid and fluid-fluid pair potentials, respectively. In Eq. (4), first term is the wall-fluid potential, U^wf(r), i.e.,

and second term is the fluid-fluid potential, U^ff(r), i.e.,

Using structurally consistent u^wf(r) and u^ff(r), EQT can accurately predict the molecular density profiles of confined fluids, such as simple LJ fluids,^19,20 CO₂,²¹ and water.^22,23

We note that, in Eq. (6), to determine the fluid-fluid potential energy, U^ff(r), a mean field approximation (MFA) is invoked, i.e., the fluid-fluid pair correlations are neglected.²⁷ In earlier work on EQT,^19,20 an effective fluid-fluid pair potential, u^ff(r), was modeled using a soft-core repulsion model, which was optimized to account for the neglected pair correlations in the MFA. In this work, we use an alternative approach to account for the fluid-fluid pair correlations. In this approach, the fluid-fluid pair correlation is approximated by a hard-sphere radial distribution function (RDF), g^hs(r). Similar approach to approximate the fluid-fluid correlations has been used by Tang and Wu.²⁸ However, the hard sphere RDF approximation may not accurately reproduce the properties of a real fluid. Therefore, to account for the effects of real fluid pair correlations in the hard sphere approximation, we add a correlation correction potential, $u ccp ff ( r )$⁠, and reformulate Eq. (6) as

In Eq. (7), for u^ff(r) we use an atomistic model-based fluid-fluid pair potential, and for g^hs(r), we use the Percus-Yevick (PY) equation for a hard sphere fluid.^29–31 To model $u ccp ff ( r )$⁠, we use uniform cubic B-splines as

where the separation interval from 0 to the cutoff, $R cut ff$⁠, is discretized into n − 1 segments, {r₀, r₁, r₂, …, r_n−1}, of equal size $Δr= R cut ff / ( n − 1 )$ such that r_i = i × Δr (i ∈ (0, …, n − 1)), {c₀, c₁, c₂, …, c_n+1} are the spline knots, the index j satisfies the condition r_j ≤ r < r_j+1, and $t= r − r j Δ r$⁠.

cDFT is a continuum-based technique that describes the properties of inhomogeneous fluids from a microscopic level.^24–26 It is based on the theorem that, for a fluid in an external field, the Helmholtz free energy, F, is a unique functional of the average molecular density profile, ρ(r), independent of the external potential, V_ext(r).^24–26 Therefore, in cDFT, the grand potential is defined as a functional of ρ(r),

To determine Ω from Eq. (9), we require an expression for $F ρ ( r )$⁠. The Helmholtz energy has two parts: (i) the ideal part, $F id ρ ( r )$⁠, and (ii) the excess part, $F ex ρ ( r )$⁠, i.e.,

The ideal part of the Helmholtz energy accounts for the ideal gas free energy,

where k_B is the Boltzmann constant, $Λ= 2 π ħ 2 m k B T 1 2$ is the thermal wavelength, ħ is the reduced Planck’s constant, and m is the mass of an atom. The excess part of the intrinsic Helmholtz energy accounts for the non-bonded interactions between molecules. Modeling of $F ex ρ ( r )$ is the most challenging part of the cDFT. The exact expression for $F ex ρ ( r )$ is in general unknown.²⁶ There exist approximate functionals for $F ex ρ ( r )$⁠, such as fundamental-measure theory (FMT) functionals^32–38 and functionals based on the statistical associating fluid theory (SAFT).^39–41 In this work, we use EQT to formulate $F ex ρ ( r )$ as explained in Sec. II C.

In the EQT-cDFT approach, the fluid-fluid EQT potential model (Eq. (7)) is used to construct the excess part of the intrinsic Helmholtz energy functional as

For a confined fluid system, the external potential is due to the wall-fluid interactions, i.e., V_ext(r) = U^wf(r). Therefore, from Eqs. (9)–(11), (7), (12), and (5), we get the EQT-cDFT-based grand potential functional, $Ω EQT ρ ( r )$⁠, as

We note that the chemical potential, μ, in Eq. (13) is same as the bulk fluid chemical potential, μ_b, at temperature T and density ρ_b. The bulk chemical potential can be defined as

where F_b is the Helmholtz energy per unit volume of the bulk fluid. Applying Eqs. (10)–(12) for a bulk fluid, F_b can be computed as

where U_b is the potential energy per molecule in the bulk fluid, given by

where $R cut ff$ is the cutoff distance for u^ff(r). Hence, in the EQT-cDFT approach, the chemical potential of a confined fluid can be computed as

At equilibrium, Ω is minimum. Therefore, the equilibrium ρ(r) can be obtained by minimizing $Ω ρ ( r )$⁠. From Eq. (13), by minimizing the EQT-cDFT-based grand potential, $Ω EQT ρ ( r )$⁠, with respect to ρ(r), we obtain the equilibrium density profile of a confined fluid as

Therefore, in the EQT-cDFT approach, one can obtain the equilibrium density and potential profiles of a confined fluid by self-consistently solving Eqs. (5), (7), and (18). One of the advantages of the EQT-cDFT approach (Eq. (18)) over the Nernst-Planck approach (Eq. (2)) is that a solution of the EQT-cDFT approach does not require ρ_avg as an input; instead, ρ_avg can be computed as an output of the EQT-cDFT simulation from the equilibrium density profile.

We consider a fluid confined in a slit-like channel, which consists of two infinitely long plane parallel walls placed in the xy plane at z = 0 and z = L. Therefore, the system is periodic in the x and y directions, and we focus only on the z-variation of the properties. From the equilibrium values of ρ(z) and $Ω EQT ρ ( z )$⁠, one can determine various thermodynamic properties of a confined fluid.²⁴ In this work, we compute the properties like total adsorption, local pressure tensor, surface tension, and solvation force as described below.

The total adsorption is the difference between the average number of fluid molecules in the confined region with and without the channel walls. The total adsorption per unit surface area, Γ, can be computed as an integral over the confined region,

where ρ_b is the bulk density of the fluid at a given T and μ, and the factor of $1 2$ is multiplied to account for the two channel walls.

The surface tension, γ, according to the thermodynamic definition, is the isothermal work required to increase the interface by unit area, i.e., $γ= 1 2 ∂ Ω ∂ A T , μ , L$ for a slit-channel system. Alternatively, γ can also be determined, according to the mechanical definition, in terms of the stress transmitted across a strip of unit width normal to the interface. In this work, we use a mechanical definition of the surface tension given by⁴² 

where P_n(z) and P_l(z) are the normal and lateral components of the local pressure tensor. In a bulk fluid, pressure is homogeneous and isotropic, however, in a confined fluid, pressure varies with the position and is anisotropic due to the wall-fluid force field and local variations in the fluid density.^43–45 For a slit-channel system, P_l(z) can be computed as a negative of the local grand potential density, Ω(z).^46,47 Therefore, from Eq. (13),

Also, for a slit-channel system in the steady-state, P_n(z) must be uniform across the channel width to satisfy a mechanical equilibrium condition. Therefore, for a given channel of width L, an average normal pressure value, P_n(L), can be computed using the thermodynamic definition as

where Ω^EQT(L) is the total grand potential of the channel of width L, which can be computed from Eq. (13). To compute ∂Ω^EQT(L)/∂L in Eq. (22), we use the central difference scheme as

where ϵ is the infinitesimal change in the channel width. Eqs. (22) and (23) are analogous to the volume perturbation expressions proposed by de Miguel and Jackson⁴⁸ in the context of vapour-liquid interfaces as an extension of the formalism introduced by Eppenga and Frenkel⁴⁹ and Harismiadis et al.⁵⁰ 

The solvation force, f_S, is the difference between the pressure exerted by a confined fluid on the channel walls and the bulk fluid pressure, P_b.^24,51 For a slit-like system in mechanical equilibrium, the pressure exerted by a confined fluid on the channel walls is equal to the average normal pressure, P_n(L). Therefore, the solvation force can be computed as

To demonstrate the EQT-cDFT approach, we simulate two different confined LJ fluid systems, namely, methane-graphene and argon-graphene slit-channel systems. In both the systems, the LJ fluid is confined between two flat graphene walls in equilibrium with the bulk reservoir. The thermodynamic state of the confined fluid is defined by the bulk reservoir temperature, T, and density, ρ_b. We consider supercritical states of the methane and argon fluids given in Table I.

To validate the accuracy of the EQT-cDFT approach, we compare the EQT-cDFT results with the equilibrium MD simulations for various channel widths from 2σ to 20σ, where σ is the length-scale parameter for the LJ interaction between fluid molecules (see Table II). To perform the reference MD simulations, we use the similar procedure and the interaction parameters given in Ref. 20. The MD simulations are performed in the NVT (canonical) ensemble by GROMACS.⁵² Methane, argon, and graphene carbon atoms are modeled as LJ type spherical particles. The LJ interaction parameters used in MD simulations for various pairs of methane, argon, and graphene carbon particles are given in Table II. The two graphene layers are placed along the xy plane, and the lateral dimensions of the layers are 3.834 00 × 3.689 27 nm². The separation distance between the two graphene layers, i.e., the channel width, is varied from 2σ to 20σ. The average densities of fluid molecules, i.e., ρ_avg = no. of molecules/volume of the channel, in various size channels are given in Table III. Spherical cutoff of 1.38 nm is used for the Lennard-Jones interactions. Wall atoms are kept fixed at their original positions. Periodic boundary conditions are specified in the x, y, and z directions. The simulation box is padded with a vacuum layer of 60σ width in the z dimension to avoid the interactions between periodic images in the z dimension. Temperature is maintained using the Nosé-Hoover thermostat⁵³ with 0.5 ps time constant. All systems are equilibrated for 2 ns and production runs of 8 ns are performed with 1 fs time step. The density profiles are computed using 0.05σ bin size along the z direction. To compute the local pressure tensors from MD, we use the method of Schofield and Henderson⁵⁴ in combination with the Gaussian smoothing kernel similar to Ref. 44. The local pressure determination method is not available in the default GROMACS 4.6.1 version. Therefore, for this work, we modified the GROMACS source code and implemented the method for determining the local pressure tensor in a slit-like geometry. Our implementation of the local pressure tensor method in GROMACS is publicly available on GitHub.⁵⁵ Recently, Vanegas et al.⁵⁶ have also implemented a local pressure calculation method in a custom version of GROMACS, which is based on the Hardy–Murdoch procedure. The surface tension and solvation force values from MD are determined by substituting the MD local pressure values in Eqs. (20) and (24), respectively. To estimate the errors in the properties from MD, we perform 5 different MD simulations with different initial conditions and obtain 5 sets of mean values of the properties. The estimate of error in the properties from MD is found to be less than 1.0%.

In the EQT-cDFT simulations of methane-graphene and argon-graphene slit-channel systems, we model u^wf(r) and u^ff(r) as LJ potentials,

where $C 12 wf$⁠, $C 6 wf$⁠, $C 12 ff$⁠, and $C 6 ff$ are the usual LJ parameters for the wall-fluid and fluid-fluid interactions. For the methane-graphene, argon-graphene, methane-methane, and argon-argon LJ interactions, we use the same LJ parameters as in the MD simulations that are given in Table II. The cutoffs for the wall-fluid, $R cut wf$⁠, and fluid-fluid, $R cut ff$⁠, pair interactions are set to 1.4 nm. For the cubic B-splines-based $u ccp ff ( r )$ (Eq. (8)), we use Δr = 0.04 nm and n = 36 and optimize the spline knot values, {c₀, c₁, c₂, …, c_n+1}, using a systematic approach based on the potential of mean force (PMF) matching, which we developed in Ref. 23. PMF-matching approach optimizes structurally consistent potential parameters. Here, to optimize the spline knot values of $u ccp ff ( r )$ for methane and argon, we use the corresponding 20σ channel atomistic trajectories from the MD simulation as a reference. More details on the PMF-matching optimization scheme can be found in Ref. 23. Fig. 1 shows the optimized $u ccp ff ( r )$ for the methane-methane and argon-argon interactions.

First, we self-consistently solve Eqs. (7), (5), and (18) and obtain the equilibrium density profiles of methane and argon inside graphene-slit channels of various widths. Fig. 2 shows that, for both methane-graphene and argon-graphene systems, the equilibrium density profiles from the EQT-cDFT agree well with the MD simulations for various channel widths. We observe that the methane and argon density profiles are similar and oscillatory, because the confined LJ fluid molecules arrange in the layers near the walls due to the competition between the wall-fluid and fluid-fluid interactions. When L = 3σ, fluid molecules arrange in two layers located around 0.95σ distance from each wall. A layer of particles is added mid-way between the walls with each σ increase in the channel width, and the density of the added layer decreases with increase in the distance from the walls. The maximum number of layers occurs when L = 18σ. Further increase in L only adds flat bulk-like region in the middle of a channel.

Fig. 3 shows the variation of ρ_avg and Γ with L. We observe that, for the channels 2σ–20σ, Γ < 0 and ρ_avg < ρ_b. Due to the strong repulsion from the wall atoms, fluid molecules cannot access the volume very close to the walls. Moreover, the layering of particles not only forms the regions of high (>ρ_b) densities, but also the regions of low (<ρ_b) densities inside a channel. The net effect of the excluded volume and layering is that the total number of fluid molecules inside the channel of a volume V is smaller than the number of molecules in the bulk of the same volume. We also observe that, for the smaller channels with no bulk-like region, Γ oscillates with L. The oscillations in Γ follow the formations of adsorbed layers with increasing L. The minimum in Γ occurs when adsorbed particles form an additional layer to arrange in a closely packed structure and reduce the average density.

Next, we compute the local pressure tensor, surface tension, and solvation force. Fig. 4 shows the lateral pressure profiles in the methane-graphene and argon-graphene systems. It can be observed that the lateral pressure predictions from the EQT-cDFT compare well with the MD for various channel widths. We observe that, in a channel, P_l(z) oscillates similar to ρ(z). The lateral pressure values are much higher near the walls than the bulk pressure. The maximum value of P_l(z) occurs near the first density peak, i.e., 0.95σ from the walls, and it is ≈5 times greater than the bulk pressure. Such high pressures in a confined fluid near the channel walls provide explanations for the confined fluid nanophases,^14,57 such as high pressure solid phases^58,59 and chemical reactions.¹³ Away from the walls, the oscillations in P_l(z) decay towards the bulk value.

Fig. 5 shows the variation of the normal pressure, surface tension, and solvation force as a function of L for both methane-graphene and argon-graphene systems. We observe that the predictions for P_n(L), γ(L), and f_S(L) from the EQT-cDFT simulations compare well with the MD simulations. Fig. 5 shows that the normal pressure oscillates with L for channels less than 9σ and it approaches the bulk pressure value for L > 9σ. The oscillations in the normal pressure are well-known and they arise because of the oscillations in the average density values (see Fig. 3).^{14,43,51,57,58} Similar to the normal pressure, the surface tension and solvation force oscillate for the smaller channels and the amplitudes of oscillations decay rapidly with increasing L.

In this work, we presented a multiscale continuum-based method to predict the structure and thermodynamic properties of confined fluids. The multiscale method is a combination of two approaches, namely, EQT and cDFT. We developed a free energy functional for cDFT based on EQT potentials. We demonstrated the EQT-cDFT approach by simulating methane and argon confined in slit-like graphene channels of various widths. The EQT-cDFT predictions for the structure and thermodynamic properties, like the density, adsorption, local pressure tensor, surface tension, and solvation force, compare well with the MD simulations. Therefore, the EQT-cDFT approach is a promising approach to accurately predict the structure and thermodynamic properties of confined fluids. Finally, we note that there is a scope to further test the EQT-cDFT approach for more complex fluid systems, such as confined water, mixtures, and electrolyte solutions.

This work was supported by NSF under Grant Nos. 1264282 and 1420882 and AFOSR under Grant No. FA9550-12-1-0464. The authors acknowledge the use of the Taub cluster provided by the Computational Science and Engineering Program at the University of Illinois.