Fluctuations, Finite-Size Effects and the Thermodynamic Limit in Computer Simulations: Revisiting the Spatial Block Analysis Method

Abstract

The spatial block analysis (SBA) method has been introduced to efficiently extrapolate thermodynamic quantities from finite-size computer simulations of a large variety of physical systems. In the particular case of simple liquids and liquid mixtures, by subdividing the simulation box into blocks of increasing size and calculating volume-dependent fluctuations of the number of particles, it is possible to extrapolate the bulk isothermal compressibility and Kirkwood–Buff integrals in the thermodynamic limit. Only by explicitly including finite-size effects, ubiquitous in computer simulations, into the SBA method, the extrapolation to the thermodynamic limit can be achieved. In this review, we discuss two of these finite-size effects in the context of the SBA method due to (i) the statistical ensemble and (ii) the finite integration domains used in computer simulations. To illustrate the method, we consider prototypical liquids and liquid mixtures described by truncated and shifted Lennard–Jones (TSLJ) potentials. Furthermore, we show some of the most recent developments of the SBA method, in particular its use to calculate chemical potentials of liquids in a wide range of density/concentration conditions.

1. Introduction

In the last few decades, computational studies of soft matter have gained ground in the no-man’s land between purely theoretical studies and experimental investigations. Arguably, this success is due to the use of statistical mechanics relations between macroscopic thermodynamic properties and microscopic components and interactions of a physical system in the thermodynamic limit (TL) [1,2]. However, and apart from a few examples [3,4,5], computer simulations are mainly constrained to consider closed systems with a finite and usually small number of particles $N_0$. These limitations introduce spurious finite-size effects, apparent in the simulation results, that in spite of the current computing capabilities are still the subject of intense investigations [6,7,8,9,10,11,12,13].

A meaningful comparison between computer simulations of finite systems and experimental results has been always a difficult task. In principle, it is possible to extrapolate the simulation data to the quantities of interest in the thermodynamic limit by considering systems of increasing size and performing simulations for each of them. The SBA method has been proposed as a more efficient alternative where only one system is examined and then subdivided into blocks of different size from which the data are extracted. The method is rather general since it was originally proposed to study the critical behavior of Ising systems [14,15] and then extended to study liquids [16,17,18,19,20,21] and even the elastic constants of model solids [22].

In this paper, we examine the SBA method focusing on the extrapolation of bulk thermodynamic properties of simple liquids. We use prototypical liquids and mixtures described by truncated and shifted Lennard–Jones (TSLJ) potentials to discuss the original ideas [16,17] and explore the background [20,21,23,24,25,26] for the most recent developments [6,7,9] of the method. The simple examples presented here, in addition to the results available in the literature [6], suggest that the method is suitable for the calculation of trends in the chemical potential of complex liquids in a wide range of density/concentration conditions.

The paper is organized as follows: In Section 2, we introduce the relevant finite-size effects present in standard computer simulations. In Section 3, we introduce the finite-size integral equations for liquids and illustrate the procedure to extrapolate thermodynamic quantities. In Section 4, we discuss the extension of the block analysis method to liquid mixtures. We conclude the paper in Section 5.

2. Boundary and Ensemble Finite-Size Effects

Statistical mechanics establishes the connection between macroscopic thermodynamic properties and the microscopic components and interactions of a physical system. An interesting example of this relation is provided by the compressibility equation that identifies the density fluctuations of a system in the grand canonical ensemble with the bulk isothermal compressibility ${\kappa }_T$ [27]. In the thermodynamic limit (TL), the isothermal compressibility of a homogeneous system is related to the fluctuations of the number of particles via the expression [1]:

$${\chi }_T^{\infty }=\frac{\langle N^2\rangle -{\langle N\rangle }^2}{\langle N\rangle },$$

(1)

with $\langle N\rangle$ the average number of particles contained in a volume V of the fluid. The reduced isothermal compressibility ${\chi }_T^{\infty }=\rho k_{\mathrm{B}}T{\kappa }_T$ is the ratio between the bulk isothermal compressibility of the system, ${\kappa }_T$, and the isothermal compressibility of the ideal gas ${(\rho k_{\mathrm{B}}T)}^{-1}$ with $\rho =\langle N\rangle /V$.

Various finite-size effects can be included in the block analysis aiming at extrapolating interesting thermodynamic quantities. In practice, let us consider a system of $N_0$ particles where the simulation box of volume $V_0=L_0^3$ is divided into subdomains of volume $V=L^3$, as illustrated in Figure 1. By evaluating the fluctuations of the number of particles in these subdomains, it is possible to obtain the distribution $P_{L,L_0}(N)$ of the number of particles, with k-moments given by [25]:

$${\langle N^k\rangle }_{L,L_0}=\sum _{N=0}^{N_0}{N}^k{P}_{L,L_0}(N).$$

(2)

The second moment of the distribution is related to the reduced isothermal compressibility of the finite system ${\chi }_T(L,L_0)$ [14,16,17,25]:

$${\chi }_T(L,L_0)=\frac{{\langle N^2\rangle }_{L,L_0}-{\langle N\rangle }_{L,L_0}^2}{{\langle N\rangle }_{L,L_0}}.$$

(3)

The finite-size reduced isothermal compressibility, ${\chi }_T(L,L_0)$, can be extrapolated to the reduced isothermal compressibility in the TL, ${\chi }_T^{\infty }$, taking the limits $L,L_0\to \infty$. Originally [17,25], by applying periodic boundary conditions (PBCs) to the total linear size $L_0$ and taking into account volumes such that $L\gg \zeta$ with $\zeta$ the correlation length of the system, it has been proposed that the difference between ${\chi }_T(L,L_0)$ and ${\chi }_T^{\infty }$ is related to boundary effects associated with the finite-size of the subdomains. This difference takes the form [16,17]:

$${\chi }_T(L,L_0\to \infty )={\chi }_T^{\infty }+\frac{c}{L}+O\left(\frac{1}{L^2}\right),$$

(4)

with c a constant. Recently, Equation (4) has been obtained [28] using arguments based on the thermodynamics of small systems [29,30], underpinning the consistency of the result.

To investigate this expression, we consider a liquid system whose potential energy is described by a 12–6Lennard–Jones potential truncated, with cutoff radius $r_c/\sigma =2^{1/6}$, and shifted. The parameters $ϵ$, $\sigma$ and m, define the units of energy, length and mass, respectively. All the results are expressed in LJ units with time $\sigma {(m/ϵ)}^{1/2}$, temperature $ϵ/k_{\mathrm{B}}$ and pressure $ϵ/{\sigma }^3$. Various system sizes, namely $N_0={10}^4,{10}^5$ and ${10}^6$, are considered, and the density is fixed at $\rho {\sigma }^3=0.864$, thus defining the linear size of the simulation box $L_0$. The systems are equilibrated at $k_{\mathrm{B}}T=1.2ϵ$, enforced with a Langevin thermostat with damping coefficient $\gamma (\sigma {(m/ϵ)}^{1/2})=1.0$, for $2\times {10}^6$ molecular dynamics (MD) steps using a time step of $\delta t/(\sigma {(m/ϵ)}^{1/2})={10}^{-3}$. Production runs span ${10}^6$ MD steps. All the simulations have been performed with the ESPResSo++ [32] simulation package.

To use the block analysis method, we compute the fluctuations of the number of particles. In particular, we choose domains of size $1<L/\sigma <L_0/\sigma$ to scan continuously the fluctuations as a function of domain size. To increase the amount of statistics, we use 100 randomly-positioned subdomains per simulation frame.

In Figure 2, we report ${\chi }_T(L,L_0)$ as a function of $\sigma /L$. The linear behavior predicted in Equation (4) is apparent for $L\ll L_0$. There are evident deviations from the linear behavior, which are not included in Equation (4), since this equation has been obtained for a system in the grand canonical ensemble. As a matter of fact, the deviations from linearity are mainly related to the fixed size of the system because when $L\to L_0$, ${\chi }_T(L_0,L_0)=0$, that is, the fluctuations of the number of particles for a closed system are equal to zero. In principle, the isothermal compressibility in the TL can be extracted by extrapolating a line to the y-axis, i.e., $\sigma /L\to 0$, and determining the y-intercept. This procedure, however, might lead to ambiguous and strongly-size-dependent results as suggested by the same plot.

From the previous discussion, Equation (4) satisfactorily describes the boundary size effects present in a system described in the grand canonical ensemble. However, ensemble size effects, i.e., the fact that we are computing quantities defined in the grand canonical ensemble using information obtained from a system in a canonical ensemble, are important even in cases where the size of the system might appear to be enormous ($L_0/\sigma =105$ for $N_0={10}^6$ where $\zeta /\sigma \approx 10$).

It is thus clear that the isothermal compressibility of a finite-size system in the TL, i.e., $L,L_0\to \infty$ with $\rho =N_0/L_0^3$, should equate to the bulk isothermal compressibility ${\kappa }_T$. An elegant analysis using probabilistic arguments for the ideal gas case [26,33] shows that the finite-size reduced isothermal compressibility can be written as:

$${\chi }_T(L,L_0)={\chi }_T^{\infty }\left(1-{\left(\frac{L}{L_0}\right)}^3\right).$$

(5)

In spite of the simplicity of the system chosen in this study, it cannot be identified with the ideal gas. However, at very low densities and temperature $k_{\mathrm{B}}T=1.2ϵ$, the system behaves more like a real gas, and a meaningful trend could be identified. Therefore, to investigate Equation (5), we consider the density range $\rho {\sigma }^3=0.1,\cdots ,1.0$ for systems of size $N_0={10}^5$ particles. Results are presented in Figure 3 for the cases $\rho {\sigma }^3=0.1,0.2$ and 0.3. The three datasets follow the theoretical prediction in Equation (5) with deviations from this behavior for $L\ll L_0$, thus indicating the signature of boundary finite-size effects. As expected, the data presented also suggest that upon increasing density, the deviations from the ideal gas behavior become more evident, as can be seen in the case $\rho {\sigma }^3=0.3$.

This is also seen in Figure 4, where for a system with density $\rho {\sigma }^3=0.864$, the deviations from the ideal gas case are much more evident. As a matter of fact, even for the largest size considered ($N_0$ = ${10}^6$), it is not possible to convincingly reproduce the ideal gas behavior.

Nonetheless, one intuitively could imagine that the following expression:

$${\chi }_T(L,L_0)={\chi }_T^{\infty }\left(1-{\left(\frac{L}{L_0}\right)}^3\right)+\frac{c}{L}+O\left(\frac{1}{L^2}\right),$$

(6)

captures the two finite-size effects, ensemble and boundary [25]. By neglecting the $O(1/L^2)$ terms, defining $\lambda =L/L_0$ and multiplying everything by $\lambda$, we obtain:

$$\lambda {\chi }_T(\lambda )=\lambda {\chi }_T^{\infty }\left(1-{\lambda }^3\right)+\frac{c}{L_0}.$$

(7)

Equation (7) is more convenient to analyze because in the limit $\lambda \to 0$, provided that $\zeta <L<L_0$, ${\lambda }^3$ is negligible, and this expression can be approximated to a linear function in $\lambda$ with slope ${\chi }_T^{\infty }$ and y-intercept equal to $c/L_0$. In particular, we use a simple linear regression in the interval $0.0<\lambda <0.3$, with the fluctuations data for $N_0={10}^5$, to find ${\chi }_T^{\infty }=0.0295(5)$ and $c=0.415(5)\sigma$. Results of the scaled fluctuations $\lambda {\chi }_T(\lambda )$ minus $c/L_0$ are presented in Figure 5, where the intensive character of the constant c becomes clear. By replacing the calculated values ${\chi }_T^{\infty }$ and c in Equation (7), we obtain the black curve that superimposes on the simulation data in the full range $0<\lambda <1$.

In addition to the explicit finite-size effects discussed above, there is another type of effect related to the periodicity of the simulation box. This is the case of implicit finite-size effects that appear due to anisotropies in the pair correlation function of the system, generated by the use of PBCs [34,35]. These effects, extremely important for small simulation setups, appear as oscillations in $\lambda {\chi }_T(\lambda )$ for $\lambda \approx 1$ caused by short range interactions between the system and its nearest neighbor images. However, given the large sizes of the systems considered here, implicit finite-size effects can be safely ignored in the present discussion.

With the trajectories of the system with $N_0={10}^5$ particles in the density interval $0.1<\rho {\sigma }^3<1.0$, we compute the scaled fluctuations $\lambda {\chi }_T(\lambda )$ and determine, as before, the ratio ${\chi }_T^{\infty }={\kappa }_T/{\kappa }_T^{IG}$ as a function of the density, with ${\kappa }_T^{IG}={(\rho k_{\mathrm{B}}T)}^{-1}$ the isothermal compressibility of the ideal gas (see Figure 6). As expected for this system at $k_{\mathrm{B}}T=1.2ϵ$, a monotonically-decreasing behavior is observed since the system becomes less compressible as the density increases.

The isothermal compressibility as a function of the density allows one to investigate more interesting thermodynamic properties, as has been recently demonstrated [6,7]. For example, the isothermal compressibility can be written as:

$${\kappa }_T=\frac{1}{{\rho }^2}{\left.\frac{\partial \rho }{\partial \mu }\right|}_T,$$

(8)

which can be rearranged, in terms of the chemical potential $\mu$, as:

$$\delta \mu ={\int }_{{\rho }_0}^{\rho }\frac{d{\rho }^{\prime }}{{\rho }^{\prime 2}{\kappa }_T}$$

(9)

with $\delta \mu =\mu -{\mu }_0$ and ${\mu }_0$ the chemical potential of the system at the reference density ${\rho }_0$. In practice, one usually is interested in the excess chemical potential (In this context, the word excess should be replaced with residual. The residual chemical potential is the difference between the chemical potential of the target system and that of an ideal gas at the same density, temperature and composition. We misuse the expression excess chemical potential to match the modern literature.):

$$\delta {\mu }^{\mathrm{ex}}=\delta \mu -k_{\mathrm{B}}Tln\rho ,$$

(10)

obtained by subtracting from $\delta \mu$ the density-dependent part of the chemical potential of the ideal gas.

To validate the results obtained using Equation (10), it is necessary to use a different computational method to evaluate ${\mu }_0$. For that purpose, any computational method aiming at calculating chemical potentials could be used. In particular, we use the spatially-resolved thermodynamic integration (SPARTIAN) method [36], recently implemented by us. In SPARTIAN, the target system, described with atomistic resolution, is embedded in a reservoir of ideal gas particles. An interface between the two subdomains is defined such that molecules are free to diffuse, adapting their resolution on the fly. A uniform density across the simulation box is guaranteed by applying a single-molecule external potential that is identified with the difference in chemical potential between the two resolutions, i.e., the excess chemical potential of the target system. This method has been validated by calculating excess chemical potentials for Lennard–Jones liquids, mixtures, as well as for simple point-charge (SPC) and extended simple point-charge (SPC/E) water models and aqueous sodium chloride solutions, all in good agreement with state-of-the-art computational methods.

For the comparison, we consider the same system at the same temperature with densities $\rho {\sigma }^3$ = 0.2, 0.4, 0.6, 0.8 and 1.0. Results for the excess chemical potential as a function of the density are presented in Figure 7 where the value of ${\rho }_0\sigma =0.6$ has been used as the reference value. Once $\delta {\mu }^{ex}$ is rescaled, it becomes clear that the agreement between the two methods is remarkable. This result suggests that the simple calculation of the fluctuations of the number of particles, used in combination with Equation (7), provides us with an efficient and accurate method to compute the chemical potential of simple liquids, which can be extended to more complex fluids [6].

In this section, Equation (7) has been introduced in a rather intuitive manner. However, the presented results suggest that it encompasses the relevant finite-size effects of the system and allows one to compute bulk thermodynamic quantities. In the following section, we derive Equation (7) more rigorously and explore, using a different example, its range of validity.

3. Finite-Size Ornstein–Zernike Integral Equation

Fluctuations of the number of particles are related to the local structure of a liquid. Let us consider a molecular liquid of average density $\rho$ at temperature T in equilibrium with a reservoir of particles, i.e., an open system. The fluctuations of the number of molecules are related to the local structure of the liquid via the Ornstein–Zernike integral equation [1,37]:

$$\frac{{\Delta }^2(N)}{\langle N\rangle }=1+\frac{\rho }{V}{\int }_V{\int }_V[g^{\mathrm{o}}({\mathbf{r}}_1,{\mathbf{r}}_2)-1]d{\mathbf{r}}_{1}d{\mathbf{r}}_2,$$

(11)

where ${\Delta }^2(N)/\langle N\rangle$ are the fluctuations of the number of particles, ${\Delta }^2(N)=\langle N^2\rangle -{\langle N\rangle }^2$ and $g^{\mathrm{o}}({\mathbf{r}}_1,{\mathbf{r}}_2)$ is the pair correlation function of the open system and ${\mathbf{r}}_1,{\mathbf{r}}_2$ the position vectors of a pair of fluid particles. To solve the integral in Equation (11), one assumes that the fluid is homogeneous, isotropic and that the system is in the thermodynamic limit (TL), i.e., $V\to \infty$, $\langle N\rangle \to \infty$ with $\rho =\langle N\rangle /V=$ constant. An infinite, homogeneous and isotropic system is translationally invariant; therefore, we rewrite Equation (11) as [1]:

$${\chi }_T^{\infty }=\frac{{\Delta }^2(N)}{\langle N\rangle }=1+4\pi \rho {\int }_0^{\infty }(g^{\mathrm{o}}(r)-1)r^{2}dr,$$

(12)

with ${\chi }_T^{\infty }=\rho k_{\mathrm{B}}T{\kappa }_T$, ${\kappa }_T$ being the isothermal compressibility of the bulk system. We have replaced $g^{\mathrm{o}}({\mathbf{r}}_1,{\mathbf{r}}_2)$ with $g^{\mathrm{o}}(r)$ the radial distribution function (RDF) of the open system, with $r=|{\mathbf{r}}_2-{\mathbf{r}}_1|$.

An alternative version of the OZ integral equation for finite systems has been introduced [25]. For a finite system with total volume $V_0$ with PBCs we have:

$${\chi }_T(V,V_0)=\frac{{\Delta }^2(N;V,V_0)}{{\langle N\rangle }_{V,V_0}}=1+\frac{\rho }{V}{\int }_V{\int }_V[g^{\mathrm{c}}(r_{12})-1]d{\mathbf{r}}_{1}d{\mathbf{r}}_2,$$

(13)

where $g^{\mathrm{c}}(r_{12})$, $r_{12}=|{\mathbf{r}}_2-{\mathbf{r}}_1|$, is the pair correlation function of the closed system with total number of particles $N_0$, and ${\Delta }^2(N;V,V_0)={\langle N^2\rangle }_{V,V_0}-{\langle N\rangle }_{V,V_0}^2$. The fluctuations of the number of particles thus depend on both subdomain and simulation box volumes.

For a single component fluid of density $\rho$ at temperature T with fixed number of particles $N_0$ and volume $V_0$, its RDF can be written in terms of an expansion around $N_0$ as [23,24,25,26,33]:

$$\begin{array}{cc}\hfill g^{\mathrm{c}}(r)=& g^{\mathrm{o}}(r)-\frac{{\chi }_T^{\infty }}{N_0}.\hfill \end{array}$$

(14)

As a matter of fact, the expansion includes terms that depend on the partial derivative of $g^{\mathrm{o}}(r)$ with respect to the density. However, we anticipate here that for the present analysis, their contribution is negligible [6]. By replacing $g^{\mathrm{c}}(r)$ in the integral on the r.h.s of Equation (13), we obtain:

$$\frac{\rho }{V}{\int }_V{\int }_V(g^{\mathrm{c}}(r_{12})-1)d{\mathbf{r}}_{1}d{\mathbf{r}}_2=I_{V,V}-\frac{V}{V_0}{\chi }_T^{\infty },$$

(15)

where:

$$I_{V,V}=\frac{\rho }{V}{\int }_V{\int }_V(g^{\mathrm{o}}(r_{12})-1)d{\mathbf{r}}_{1}d{\mathbf{r}}_2,$$

(16)

and we use that $\rho =N_0/V_0$.

Next, we include explicitly the second finite-size effect, i.e., the fact that the volume V is finite and embedded into a finite volume $V_0$ with PBCs. For this, we rewrite $I_{V,V}$ as [17]:

$$I_{V,V_0-V}=I_{V,V_0}-I_{V,V},$$

with:

$$\begin{array}{cc}\hfill I_{V,V_0}& =\frac{\rho }{V}{\int }_V{\int }_{V_0}(g^{\mathrm{o}}(r_{12})-1)d{\mathbf{r}}_{1}d{\mathbf{r}}_2\hfill \\ \hfill I_{V,V_0-V}& =\frac{\rho }{V}{\int }_V{\int }_{V_0-V}(g^{\mathrm{o}}(r_{12})-1)d{\mathbf{r}}_{1}d{\mathbf{r}}_2.\hfill \end{array}$$

As pointed out by Rovere, Heermann and Binder [17], the two integrals $I_{V,V}$ and $I_{V,V_0}$ are equal when ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$ are both within the volume V. When $r_{12}>\zeta$, the integrand $(g^{\mathrm{o}}(r_{12})-1)=0$, and it does not contribute to the integrals. Close to the boundary of the subdomain V, for $r_{12}<\zeta$, and in particular when ${\mathbf{r}}_1$ lies inside and ${\mathbf{r}}_2$ outside the volume V, there are contributions missing in $I_{V,V}$, which are present in $I_{V,V_0}$. Therefore, the difference between the two integrals $I_{V,V_0-V}=I_{V,V_0}-I_{V,V}$ must be proportional to the surface volume ratio of the subdomain V [17], i.e.,

$$I_{V,V_0-V}=\frac{c_1}{L}+{\left(\frac{c_2}{L}\right)}^2+O\left(\frac{1}{L^3}\right),$$

(17)

with $c_1$, $c_2$ proportionality constants with units of length that, at this point, we assume to be intensive.

To compute $I_{V,V_0}$, we require that $\zeta <L<L_0$. Since we assume PBCs, the system is translationally invariant. Hence, upon applying the transformation ${\mathbf{r}}_{12}\to \mathbf{r}={\mathbf{r}}_2-{\mathbf{r}}_1$, the expression:

$$I_{V,V_0}=\rho {\int }_{V_0}(g^{\mathrm{o}}(r)-1)d\mathbf{r}={\chi }_T^{\infty }-1$$

(18)

is obtained, where we assume that $g^{\mathrm{o}}(r>\zeta )=1$, thus ignoring fluctuations of the RDF beyond the volume V. By combining these two results, we obtain:

$$I_{V,V}={\chi }_T^{\infty }-1+\frac{c_1}{L}+{\left(\frac{c_2}{L}\right)}^2,$$

(19)

and by including this result in Equation (15), we arrive at the following expression:

$$\frac{\rho }{V}{\int }_V{\int }_V(g^{\mathrm{c}}(r_{12})-1)d{\mathbf{r}}_{1}d{\mathbf{r}}_2={\chi }_T^{\infty }\left(1-{\left(\frac{L}{L_0}\right)}^3\right)-1+\frac{c_1}{L}+{\left(\frac{c_2}{L}\right)}^2.$$

(20)

Finally, this expression becomes:

$${\chi }_T(L,L_0)={\chi }_T^{\infty }\left(1-{\left(\frac{L}{L_0}\right)}^3\right)+\frac{c_1}{L}+{\left(\frac{c_2}{L}\right)}^2,$$

(21)

and by defining $\lambda =L/L_0$, we write:

$$\lambda {\chi }_T(\lambda )=\lambda {\chi }_T^{\infty }\left(1-{\lambda }^3\right)+\frac{c_1}{L_0}+{\left(\frac{c_2}{L_0}\right)}^2\frac{1}{\lambda }.$$

(22)

Equations (7) and (22) differ in the $c_2^2/L_0^2\lambda$ term that appears from considering the boundary finite-size effects. One possible scenario in which this difference might play a role is in the case of simulations near critical conditions where the correlation length of the system tends to infinity.

To test this expression, we perform simulations of systems with potential energy described by the truncated, at $r_c/\sigma =2.5$, and shifted 12–6 Lennard–Jones potential. We consider systems with $N_0=24,000$ particles, with densities spanning the range $0.05<\rho {\sigma }^3<0.70$. Two temperatures were considered, $k_{\mathrm{B}}T=2.00ϵ$ and $1.15ϵ$. The critical point of this system has been reported at ${\rho }_c{\sigma }^3=0.319$ and $k_{\mathrm{B}}{T}_c=1.086ϵ$ [38].

We report the reduced fluctuations $\lambda {\chi }_T(\lambda )$ as a function of $\lambda$ for $\rho {\sigma }^3=0.3$ in Figure 8. In the case $k_{\mathrm{B}}T=2.00ϵ$, the effect of the ${\lambda }^{-1}$ term in Equation (22) is negligible, and a linear approximation in the region $\lambda <0.3$ seems to be well justified. However, for the case close to the critical point, i.e., $k_{\mathrm{B}}T=1.15ϵ$, the effect of this term is evident and should be included in the extrapolation to ${\chi }_T^{\infty }$.

Finally, upon extrapolating to ${\chi }_T^{\infty }$, an interesting behavior is observed for the bulk isothermal compressibility ${\kappa }_T$ as a function of density (Figure 9). In the case $k_{\mathrm{B}}T=2.00ϵ$, as expected, a monotonically-decreasing behavior with increasing density is observed. More interestingly, in the case $k_{\mathrm{B}}T=1.15ϵ$, the monotonically-decreasing behavior is interrupted by a singularity in the isothermal compressibility in the vicinity of the critical density. This cusp in the curve is expected since the isothermal compressibility of a fluid at the critical point is infinite.

The use of finite-size integral equations is general enough to admit generalizations of other systems of interest. In the next section, we describe one of such possible extensions: the study of binary mixtures.

4. Mixtures

Kirkwood–Buff (KB) theory [39] is arguably the most successful framework to investigate the properties of liquid mixtures that relates the local structure of a system to density fluctuations in the grand canonical ensemble. These quantities are in turn related to equilibrium thermodynamic quantities such as the compressibility, the partial molar volumes and the derivatives of the chemical potentials [2]. Formulated more than sixty years ago, KB enjoys renewed interest in the computational soft-matter and statistical physics communities [6,7,9,10,11,12,13]. Recent works have shown promising applications related to solvation of biomolecules [40] and potential uses to compute multicomponent diffusion in liquids [41] and to study complex phenomena such as self-assembly of proteins [42] and polymer conformation in complex mixtures [4,43].

For a multicomponent fluid of species $i,j$ in equilibrium at temperature T, the Kirkwood–Buff integral (KBI) is defined as:

$$G_{ij}^{\mathrm{o}}=V\left(\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}}{\langle N_i\rangle }\right)=\frac{1}{V}{\int }_V{\int }_V[g_{ij}^{\mathrm{o}}({\mathbf{r}}_{12})-1]d{\mathbf{r}}_{1}d{\mathbf{r}}_2,$$

(23)

with ${\delta }_{ij}$ the Kronecker delta. The superscript (o) indicates that this definition holds for an open system, i.e., a system in the grand canonical ensemble. In practice, we compute fluctuations of the number of particles in a subdomain of volume V embedded in a reservoir whose size goes to infinity. Thus, $\langle N_i\rangle$ is the average number of i-particles inside V, or ${\rho }_i=\langle N_i\rangle /V$. $g_{ij}^{\mathrm{o}}({\mathbf{r}}_{12})$ is the multicomponent radial distribution function (RDF) of the infinite system, with ${\mathbf{r}}_{12}={\mathbf{r}}_2-{\mathbf{r}}_1$.

Let us recall that in computer simulations one considers systems with total fixed number of particles $N_0$ and volume $V_0$ with PBCs. In this case, we have [35]:

$$G_{ij}(L,L_0)=V\left(\frac{{\langle N_i{N}_j\rangle }^{\prime }-{\langle N_i\rangle }^{\prime }{\langle N_j\rangle }^{\prime }}{{\langle N_i\rangle }^{\prime }{\langle N_j\rangle }^{\prime }}-\frac{{\delta }_{ij}}{{\langle N_i\rangle }^{\prime }}\right)=\frac{1}{V}{\int }_V{\int }_V[g_{ij}^{\mathrm{c}}({\mathbf{r}}_{12})-1]d{\mathbf{r}}_{1}d{\mathbf{r}}_2.$$

(24)

The finite-size KBI $G_{ij}(L,L_0)$ is evaluated by computing fluctuations of the number of particles in finite subdomains of volume V inside a simulation box of volume $V_0$. The average number of i-particles ${\langle N_i\rangle }^{\prime }\equiv {\langle N_i\rangle }_{V,V_0}$ depends on both subdomain and simulation box volumes. Moreover, the integral on the r.h.s. of Equation (24) should be evaluated for the RDF of the finite system $g_{ij}^{\mathrm{c}}({\mathbf{r}}_{12})$ with volume $V_0$ by using a finite integration domain V.

As has been done for the single component case, we include in this example both, ensemble and boundary, finite-size effects. For the former, the following correction has been suggested [44]:

$$g_{ij}^{\mathrm{c}}(r)=g_{ij}^{\mathrm{o}}(r)-\frac{1}{V_0}\left(\frac{{\delta }_{ij}}{{\rho }_i}+G_{ij}^{\infty }\right),$$

(25)

based on the asymptotic limit $g_{ij}^{\mathrm{c}}(r\gg \zeta )=1-({\delta }_{ij}/{\rho }_i+G_{ij}^{\infty })/V_0$ discussed in [2]. As expected, when the total volume $V_0\to \infty$, we recover $g_{ij}^{\mathrm{c}}(r)=g_{ij}^{\mathrm{o}}(r)$. By including Equation (25) in the integral on the r.h.s. of Equation (24) and evaluating the finite-size integral as for the single component case, we finally obtain:

$$\lambda G_{ij}(\lambda )=\lambda G_{ij}^{\infty }\left(1-{\lambda }^3\right)-{\lambda }^4\frac{{\delta }_{ij}}{{\rho }_i}+\frac{{\alpha }_{ij}}{L_0},$$

(26)

with $\lambda \equiv L/L_0$ and ${\alpha }_{ij}$ an intensive parameter with units of length. In the limit $L_0\to \infty$, the following expression is obtained:

$$G_{ij}(L,L_0\to \infty )=G_{ij}^{\infty }+\frac{{\alpha }_{ij}}{L},$$

(27)

that describes the finite-size effects on the KBIs for a system in the grand canonical ensemble. Consistent with this limiting case in Equation (26), Equation (27) has been obtained from the thermodynamics of small systems [45,46].

For the investigation of Equation (26), we perform simulations for binary mixtures $(A,B)$ of Lennard–Jones (LJ) fluids. We use a purely repulsive 12-6 LJ potential truncated and shifted with cutoff radius $2^{1/6}\sigma$. The potential parameters are chosen as ${\sigma }_{AA}={\sigma }_{BB}={\sigma }_{AB}=\sigma$, and ${ϵ}_{AA}=1.2ϵ$, ${ϵ}_{BB}=1.0ϵ$ with ${ϵ}_{AB}=({ϵ}_{AA}+{ϵ}_{BB})/2=1.1ϵ$. All the results are expressed in LJ units with energy $ϵ$, length $\sigma$, mass $m_A=m_B=m$, time $\sigma {(m/ϵ)}^{1/2}$, temperature $ϵ/k_{\mathrm{B}}$ and pressure $ϵ/{\sigma }^3$. As before, simulations are carried out using ESPResSo++ [32] with a time step of $\delta t/(\sigma {(m/ϵ)}^{1/2})={10}^{-3}$. Constant temperature $k_{\mathrm{B}}T=1.2ϵ$ is enforced through a Langevin thermostat with damping coefficient $\gamma (\sigma {(m/ϵ)}^{1/2})=1.0$. The size of the system is $N_0=23,328$ in the range of mole fractions of A-molecules $x_A=0.1,\cdots ,1.0$. The pressure is fixed at $P{\sigma }^3/ϵ=9.8$ by adjusting the number density of the system at values around $\rho {\sigma }^3\approx 0.86$ (or $L_0/\sigma \approx 30$). We perform equilibration runs of $64\times {10}^6$ MD steps and production runs of $2\times {10}^6$ MD steps. To compute $G_{ij}(\lambda )$, we select 800 frames per trajectory and for each frame identify 1000 randomly-positioned subdomains with linear sizes ranging from $2<L/\sigma <L_0/\sigma$.

In Figure 10, results for finite-size KBIs are presented for four mole fractions, namely (a) $x_A=0.20$; (b) $x_A=0.30$; (c) $x_A=0.50$ and (d) $x_A=0.80$. Plots of $G_{AB}$ (green circles) tend to zero when $\lambda \to 1$, as indicated by the horizontal green lines. By contrast, $G_{AA}\to 1/{\rho }_A$ (indicated by horizontal red lines) when $\lambda \to 1$. The region $\lambda <3$, indicated by vertical black lines, is where simple linear regression is used to find $G_{ij}^{\infty }$ and ${\alpha }_{ij}$. By replacing such values in Equation (26), we obtained the black curves that, in all cases, superimpose on the simulation data for the full interval $0<\lambda <1$.

The bulk KBIs are related to various thermodynamic quantities. For example, the isothermal compressibility is given by [39]:

$${\kappa }_T=\frac{1+{\rho }_A{G}_{AA}+{\rho }_B{G}_{BB}+{\rho }_A{\rho }_B(G_{AA}{G}_{BB}-G_{AB}^2)}{k_{\mathrm{B}}T({\rho }_A+{\rho }_B+{\rho }_A{\rho }_B(G_{AA}+G_{BB}-2G_{AB}))},$$

(28)

with ${\rho }_{A,B}$ the number density of the corresponding species.

Results for the isothermal compressibility obtained from the $G_{ij}^{\infty }$ values are presented in Figure 11. Single component cases corresponding to systems composed by only type-A and type-B particles are indicated by the horizontal black lines. As expected, the system composed by strongly interacting particles, i.e., the type-A, has a lower compressibility. The behavior of the isothermal compressibility is nearly ideal since it follows closely the relation ${\kappa }_T=(1-x_A){\kappa }_T^B+x_A{\kappa }_T^A$, with ${\kappa }_T^Aϵ/{\sigma }^3=0.012(1)$ and ${\kappa }_T^Bϵ/{\sigma }^3=0.0281(8)$, as indicated by the solid black line.

Finally, the extrapolated KBIs have been used to compute the derivative of the chemical potential of type-A particles with respect to the number density ${\rho }_A$ using the expression [39]:

$$\frac{1}{k_{\mathrm{B}}T}{\left(\frac{\partial {\mu }_A}{\partial {\rho }_A}\right)}_{P,T}=\frac{1}{{\rho }_A}+\frac{G_{AB}-G_{AA}}{1+{\rho }_A(G_{AA}-G_{AB})},$$

(29)

that, as has been done for the single component case, can be integrated to obtain [6]:

$$\delta {\mu }_A=k_{\mathrm{B}}T{\int }_{{\rho }_A^0}^{{\rho }_A}\left[\frac{1}{{\rho }_A^{\prime }}+\frac{G_{AB}-G_{AA}}{1+{\rho }_A^{\prime }(G_{AA}-G_{AB})}\right]d{\rho }_A^{\prime }.$$

(30)

This is the chemical potential shifted by a reference chemical potential computed at density ${\rho }_A^0$ [4,43]. By removing the density and concentration terms of the chemical potential of an ideal mixture, the excess chemical potential can be written as:

$$\delta {\mu }_A^{ex}=\delta {\mu }_A-k_{\mathrm{B}}Tln(x_A)-k_{\mathrm{B}}Tln({\rho }_A).$$

(31)

We compare the results obtained using Equations (30) and (31) with the results obtained with the SPARTIAN method [36] and use the excess chemical potential result from $x_A=0.3$ to find the reference value. We present the results in Figure 12 where a good agreement between the two datasets is apparent. To conclude this section, it has been shown that the block analysis method constitutes a robust strategy to compute chemical potentials of liquids and mixtures in a wide range of density/concentration conditions.

5. Summary and Concluding Remarks

In general, a direct comparison between a real system and a finite-size simulation is prevented by the fixed and relatively small number of particles used in the latter. As has been encoded in the title of the paper, the spatial block analysis method employs a clever combination of finite-size effects, ensemble and boundary, and density fluctuations to extrapolate bulk isothermal compressibilities and Kirkwood–Buff integrals in the thermodynamic limit.

In this work, we have illustrated with prototypical Lennard–Jones liquids and liquid mixtures the working mechanisms of the method. Upon identifying the relevant finite-size effects and assessing their impact on the simulation results, we have intuitively introduced an analytical expression connecting the fluctuations measured in a small subdomain of the simulation box with the bulk isothermal compressibility for a single component fluid.

Subsequently, the same analytical expression has been rigorously obtained from a finite-size version of the Ornstein–Zernike integral equation. Using a challenging system close to critical point conditions, we have tested the range of validity of the method and obtained results in line with theoretical expectations.

Then, for a multicomponent system, we have applied the same protocol to the Kirkwood–Buff integrals. Using the corresponding analytical expression, it is possible to obtain the Kirkwood–Buff integrals in the thermodynamic limit. In both single and multicomponent systems, the method allows one to compute the chemical potential of a liquid/mixture for a wide range of density/concentration conditions, provided a single reference chemical potential has been determined using a different computational method. These results contribute to establishing the spatial block analysis method as a powerful tool to investigate systems where the accurate computation of the chemical potential is of paramount importance.