Bridging scales with thermodynamics: from nano to macro

Classical thermodynamics applies to systems on micrometer- or larger scales, as in the laboratory or in a chemical plant. Thermodynamic properties of nano-scale systems can be very far from the thermodynamic limit. It is well known that small systems are no longer extensive due to e.g. surface effects. Typical for a small system is that surface energies can be large. Even line and corner energies may contribute to the system energy.

Systems that no longer obey the extensivity we find on the large scale, can be understood as being small. Hill [1] proposed that the effect of such smallness can be handled by introducing an extra term in the classical thermodynamic equations. It is also necessary to specify the control variables, because unlike a large system, a small system has thermodynamic properties which depend on the environment. Using Hill's formalism one can derive equations for the effect of smallness on thermodynamic properties in an ensemble of interest. The small system that we have studied first [2–4] is controlled by the environment in such a way that it has constant chemical potential μ, volume V and temperature T. The thermodynamic limit is usually approached in computer simulations by using periodic boundary conditions. The effect of such boundary conditions has been well studied, see e.g. Landau et al [5] and Siepmann et al [6]. We can take advantage of properties of periodic boundary conditions to create a well-defined environment for the small system. The small system of interest is put inside a reservoir that has constant μ and T, and the effect of smallness is studied. At equilibrium, the number of particles and the energy of the small system will fluctuate. If the small system is sufficiently small compared to the size of the reservoir, however, it will be in the grand-canonical ensemble.

The important difference between the small and large scale system is illustrated in figure 1. Configurations illustrated in figure 1(a) are characteristic for systems with periodic boundary conditions, while the configurations in figure 1(b) are allowed for the small system. This setup creates surface energy for the small system. This presentation will focus on thermodynamic properties that follow from energy- and density fluctuations: the thermodynamic correction factor Γ, Kirkwood–Buff integrals and the molar enthalpy H.

Zoom In Zoom Out Reset image size

The thermodynamic correction factor is a measure of the deviation from ideal behaviour and can be used to convert Maxwell–Stefan diffusivities into Fick diffusivities. For an open single-component system, the thermodynamic correction factor follows from density fluctuations at a constant volume

$$\begin{equation} \frac{1}{\Gamma } = \frac{{\langle {N^2 } \rangle - \langle N \rangle ^2 }}{{\langle N \rangle }}, \end{equation} \tag{ 1 }$$

where brackets denote an ensemble average, and N is the number of particles. The Kirkwood–Buff integral [7] provides derivatives of activity coefficients in mixtures from

$$\begin{equation} G_{ij} = V\frac{{\langle {N_i N_j }\rangle - \langle {N_i }\rangle \langle {N_j }\rangle }}{{\langle {N_i } \rangle \langle {N_j } \rangle }} - \frac{{\delta _{ij} }}{{c_i }}, \end{equation} \tag{ 2 }$$

where V is the volume, δ_ij the Kroenecker delta and c_i is the concentration of component i. The brackets denote averages in the grand-canonical ensemble. The molar enthalpy H can likewise be computed from energy and density fluctuations

$$\begin{equation} \left( {\frac{{\partial H}}{{\partial N}}} \right)_{p,T} = \frac{{\langle {UN} \rangle - \langle N \rangle \langle U \rangle }}{{\langle {N^2 } \rangle - \langle N \rangle ^2 }} + RT, \end{equation} \tag{ 3 }$$

where U is the potential energy of the system, and R is the universal gas constant. Using molecular dynamics simulations, we shall use these formulae to determine thermodynamic properties at different length scales. The characteristic length of the system, L = V^1/d, is varied, see figure 2, where d is the system's dimension, here d = 3.

Zoom In Zoom Out Reset image size

Depending on L, the thermodynamic properties differ from those of a large system, but in a predictable way. We find that the surface contributions scale linearly with 1/L, when surface energies are important [3, 7]:

$$\begin{equation} \frac{1}{\Gamma } = \frac{1}{{\Gamma _\infty }} + \frac{C}{L},\quad G = G_\infty + \frac{D}{L},\quad H = H_\infty + \frac{B}{L} \cdot \end{equation} \tag{ 4 }$$

We shall see how this scaling can be used in practice to determine thermodynamic properties of bulk systems from fluctuations in a system as small as a few particles. Line and nook energies will add terms of second and third order to these equations.

The method can also be used to determine surface thermodynamic data. The first step will be to determine the surface thickness of the adsorbed layer. This is done according to Gibbs [8]. The sampling area is varied, while the thickness of the surface is fixed. In this manner a two-dimensional thermodynamic surface system is created.

For more information regarding the details of the method in general, we refer to Schnell et al [4].

Diffusion is often a rate-limiting process in chemistry or biology, and reliable data of mixtures are scarce. This is so in particular for ternary systems, where no simulations have been done, and only a few experimental results are reported.

The generalized Fick equations for an n-component mixture can be written for the flux of each component i by

$$\begin{equation} J_i = - c_{\rm{t}} \sum\limits_j^{n - 1} {D_{ij} \nabla x_j } . \end{equation} \tag{ 5 }$$

The total concentration is c_t in mol m⁻³, x_j is the mole fraction and D_ij are Fick diffusivities. These equations are frequently the basis for reduction of experimental data. On the other hand, the Maxwell–Stefan (MS) equations are suitable for computations:

$$\begin{equation} - \frac{1}{{RT}}\nabla \mu _i = \sum\limits_{j = 1,j \ne 1}^n {\frac{{x_j (u_i - u_j )}}{\rm{D_{ij} }}} , \end{equation} \tag{ 6 }$$

where the parenthesis gives the difference in the average velocities of components. Upright D_ij are MS diffusivities. As both sets of equations must describe the same physical phenomena, the entropy production is invariant, and we obtain relations between variables. In vector notation we have [9]

$$\begin{equation} [D] = [B]^{ - 1} [\Gamma ], \end{equation} \tag{ 7 }$$

where

$$\begin{equation} B_{ii} = \frac{{x_i }}{{D_{in} }} + \sum\limits_{j = 1,j \ne 1}^n {\frac{{x_j }}{\rm{D_{ij} }}} , \end{equation} \tag{ 8 }$$

$$\begin{equation} B_{ij} = - x_i \left[ {\frac{1}{\rm{D_{ij} }} - \frac{1}{{D_{in} }}} \right] \cdot \end{equation} \tag{ 9 }$$

The thermodynamic factors of the mixture are

$$\begin{equation} \Gamma _{ij} = \delta _{ij} + x_i \left[ {\frac{{\partial \ln \gamma _i }}{{\partial x_j }}} \right]_{T,p,\Sigma } , \end{equation} \tag{ 10 }$$

where γ_i is the activity coefficient of component i, and the symbol Σ is used to indicate that the differentiation is carried out, while keeping constant the mole fractions of all other components, except the nth, so that the mole fractions sum to unity, see Liu et al [9] for further references. From knowledge of one set of diffusivities and corresponding thermodynamic factors, both as functions of concentration, one can therefore compute another set of diffusivities. This has been done using the method described in the previous section, and coefficients from experiments and simulations have been compared and found to agree remarkably well [9, 10], see below.

Knowledge of thermodynamic factors may in this manner facilitate the computational determination of Fick's diffusion coefficients which are difficult to measure.

We tested first the performance of equations (1) and (2) at the small scale using molecular simulations for systems consisting of particles interacting either with a Lennard Jones (LJ) potential (truncated and shifted at 2.5σ, where σ is the molecular diameter) or a Weeks, Chander and Andersen (WCA)-potential (a shifted LJ potential with the attractive tail cut off, see Weeks et al [11]).

For convenience, the simulations were performed in reduced units, see Frenkel and Smit [12] for details. As the unit of length we use the particle diameter σ, and the unit of energy is the LJ-energy parameter, . We simulated the large reservoir of particles in the micro-canonical ensemble, in which the small system of a varying size was embedded. The side of the small system (L = V^1/d, where d is the system's dimension) was smaller than that of the half reservoir (L_t = V_t^1/d), in order to ensure that the small system is unaffected by the periodic boundary conditions of the reservoir. With the WCA and LJ particles, the box side was L_t = 20.

The method was used to calculate the thermodynamic factor of mixtures, for which the diffusion coefficients were known experimentally. The Stefan–Maxwell diffusivities of the mixture components were computed along well defined schemes [9, 10]. The results of the calculations of the Stefan–Maxwell diffusivities were next combined with the results for the thermodynamic factor, using equations (5)–(10), to determine the Fick type diffusion coefficients. The outcome was compared with experimental results [9, 10].

Results for the thermodynamic factor and for Kirkwood–Buff integrals are shown in figure 3. Both quantities are plotted as a function of 1/L. The dependence on 1/L is significant, meaning that the deviation from the thermodynamic limit is large for small systems. Similar results were obtained for the molar enthalpy (not shown).

Zoom In Zoom Out Reset image size

The interesting fact is that the thermodynamic limit can be determined by extrapolation of small systems' properties. We have thus used the observed scaling property to predict properties of bulk systems by investigating properties of small systems! The lines in figure 3 are straight lines fitted to data for small systems with sizes L = 6–16 (0.06 < 1/L < 0.18). The value L_t = 20 prevented the small system from being affected by the periodic boundary conditions of the reservoir. We selected the value L = 6 as the smallest length of a small system to be investigated. With these values for L, one makes sure that we operate in the linear regime. The scaling law is the same for two- or three-dimensional systems (figure 3 top).

The method was also used to study surface adsorption, CO₂ on a graphite surface [8]. We selected then a cylindrical volume element with thickness equal to the surface thickness, as measured by the CO₂ layer adsorbed on the surface. The value of L was then the radius of the cylinder. We extended the size of the total surface in order to investigate a possible effect of L_t on the thermodynamic factor, see figure 4. The surface areas, L_t², are labelled X1, X4 and X9, with X1 ∼ 50 × 50 Å², X4 ∼ 100 × 100 Å² and X9 ∼ 150 × 150 Å². The figure shows results for 50 < N < 700, where N is the number of CO₂ molecules in the system. The areas X4 and X9 give similar results, while X1 give results far from these. The area X4 has L_t equal to 20 (the length of CO₂ ∼ 5 Å). This value corresponds well with the value referred to above. To find the thermodynamic factor, we therefore fitted results for 0.10 < 1/L < 0.40 to the linear equation (4) for N50 (for the other samples the region was slightly smaller), and made the extrapolation to the thermodynamic limit. From the thermodynamic factor, one can further obtain the chemical potential of the surface or the activity coefficient of the surface.

Zoom In Zoom Out Reset image size

The sign of the slope in figure 3 is related to attractive and repulsive forces at the interface. It is not yet properly understood. The scaling is in sharp contrast to the finite-size scaling observed earlier for systems with periodic boundary conditions [6]. In equation (3), the constants C, D and B, are system-specific. Figure 3 reflects the linear term, except that nook- and corner effects are present for very small L, see [4]. Clearly, such effects can become important for very small system sizes.

The above results show that Monte Carlo simulations can be used to verify that Hill's approach to thermodynamics for small systems successfully describes the size dependence of thermodynamic properties of small systems embedded in a reservoir. We find a 1/L finite size scaling behaviour for the inverse thermodynamic correction factor Γ, the Kirkwood–Buff integral and the molar enthalpy H. The values of the slopes show that finite size effects are very important. The scaling behaviour can be used to connect thermodynamic properties obtained at different length scales. This gives a new way to determine thermodynamic properties of large and small systems.

The thermodynamic factor was found for binary mixtures and used to calculate Fick's diffusivities [9, 10], according to the description above, see figures 5 and 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The agreement between computed values (simulations) and experimental values is remarkable. Diffusion coefficients for binary and ternary systems can be obtained by linking the Green–Kubo formulae to the Onsager coefficients. The values can be further used to predict new and evaluate common formulae for diffusion coefficients of components in a mixture. The method was used favourably for ternary mixtures, where experimental data are few, and where computational data are non-existent up to this point.

A new scaling method [3, 4] was used to compute thermodynamic data for the macroscopic limit from knowledge of fluctuations in small systems. The rapid and precise method offers an alternative to more difficult computations of thermodynamic factors, say from Kirkwood–Buff integrals [7]. When multiplied with computed Maxwell–Stefan coefficients, agreement was found between computed predictions and experiments for several binary systems. Diffusion coefficients for binary and ternary systems can be obtained by linking the Green–Kubo formulae to the Onsager coefficients. The outcome formulae can improve/disprove empirical formulae for diffusion coefficients [9, 10].

The success of the procedure has interesting implications. In the first place, the procedure helps to define smallness, and precise thermodynamic relations on the nano-metre scale. In the second place, it helps define local equilibrium or validity of thermodynamic relations, and when this can be expected. Taking example from, say, simulations of zeolites, one can show that a volume element must be an appropriate fraction of a unit cell in order to speak of thermodynamic properties in a meaningful way.

SS, TJHV and SK acknowledge financial support from NWO-CW through an ECHO grant number 700.58.042 and computational resources through NCF grant numbers MP-213-11, MP-213-12 and MP-213-13.