Open AccessReview

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

Maziar Heidari

¹,

Kurt Kremer

¹,

Raffaello Potestio

^1,2,3 and

Robinson Cortes-Huerto

^1,*

Max Planck Institute for Polymer Research, Ackermannweg 10, 55128 Mainz, Germany

Physics Department, University of Trento, via Sommarive 14, I-38123 Trento, Italy

INFN-TIFPA, Trento Institute for Fundamental Physics and Applications, I-38123 Trento, Italy

Author to whom correspondence should be addressed.

Entropy 2018, 20(4), 222; https://doi.org/10.3390/e20040222

Submission received: 1 March 2018 / Revised: 21 March 2018 / Accepted: 22 March 2018 / Published: 24 March 2018

(This article belongs to the Special Issue Thermodynamics and Statistical Mechanics of Small Systems)

Download

Browse Figures

Versions Notes

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.

Keywords:

computer simulations; finite-size effects; calculation of free energies; thermodynamic limit

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

κ_{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]:

χ_{T}^{\infty} = \frac{〈 N^{2} 〉 - {〈 N 〉}^{2}}{〈 N 〉},

(1)

with

〈 N 〉

the average number of particles contained in a volume V of the fluid. The reduced isothermal compressibility

χ_{T}^{\infty} = ρ k_{B} T κ_{T}

is the ratio between the bulk isothermal compressibility of the system,

κ_{T}

, and the isothermal compressibility of the ideal gas

{(ρ k_{B} T)}^{- 1}

with

ρ = 〈 N 〉 / 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]:

{〈 N^{k} 〉}_{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

χ_{T} (L, L_{0})

[14,16,17,25]:

χ_{T} (L, L_{0}) = \frac{{〈 N^{2} 〉}_{L, L_{0}} - {〈 N 〉}_{L, L_{0}}^{2}}{{〈 N 〉}_{L, L_{0}}} .

(3)

The finite-size reduced isothermal compressibility,

χ_{T} (L, L_{0})

, can be extrapolated to the reduced isothermal compressibility in the TL,

χ_{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 ≫ ζ

with

ζ

the correlation length of the system, it has been proposed that the difference between

χ_{T} (L, L_{0})

and

χ_{T}^{\infty}

is related to boundary effects associated with the finite-size of the subdomains. This difference takes the form [16,17]:

χ_{T} (L, L_{0} \to \infty) = χ_{T}^{\infty} + \frac{c}{L} + O (\frac{1}{L^{2}}),

(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} / σ = 2^{1 / 6}

, and shifted. The parameters

ϵ

σ

and m, define the units of energy, length and mass, respectively. All the results are expressed in LJ units with time

σ {(m / ϵ)}^{1 / 2}

, temperature

ϵ / k_{B}

and pressure

ϵ / σ^{3}

. Various system sizes, namely

N_{0} = 10^{4}, 10^{5}

and

10^{6}

, are considered, and the density is fixed at

ρ σ^{3} = 0.864

, thus defining the linear size of the simulation box

L_{0}

. The systems are equilibrated at

k_{B} T = 1.2 ϵ

, enforced with a Langevin thermostat with damping coefficient

γ (σ {(m / ϵ)}^{1 / 2}) = 1.0

, for

2 \times 10^{6}

molecular dynamics (MD) steps using a time step of

δ t / (σ {(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 / σ < L_{0} / σ

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

χ_{T} (L, L_{0})

as a function of

σ / L

. The linear behavior predicted in Equation (4) is apparent for

L ≪ 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}

χ_{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.,

σ / 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} / σ = 105

for

N_{0} = 10^{6}

where

ζ / σ \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

ρ = N_{0} / L_{0}^{3}

, should equate to the bulk isothermal compressibility

κ_{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:

χ_{T} (L, L_{0}) = χ_{T}^{\infty} (1 - {(\frac{L}{L_{0}})}^{3}) .

(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_{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

ρ σ^{3} = 0.1, \dots, 1.0

for systems of size

N_{0} = 10^{5}

particles. Results are presented in Figure 3 for the cases

ρ σ^{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 ≪ 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

ρ σ^{3} = 0.3

This is also seen in Figure 4, where for a system with density

ρ σ^{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:

χ_{T} (L, L_{0}) = χ_{T}^{\infty} (1 - {(\frac{L}{L_{0}})}^{3}) + \frac{c}{L} + O (\frac{1}{L^{2}}),

(6)

captures the two finite-size effects, ensemble and boundary [25]. By neglecting the

O (1 / L^{2})

terms, defining

λ = L / L_{0}

and multiplying everything by

λ

, we obtain:

λ χ_{T} (λ) = λ χ_{T}^{\infty} (1 - λ^{3}) + \frac{c}{L_{0}} .

(7)

Equation (7) is more convenient to analyze because in the limit

λ \to 0

, provided that

ζ < L < L_{0}

λ^{3}

is negligible, and this expression can be approximated to a linear function in

λ

with slope

χ_{T}^{\infty}

and y-intercept equal to

c / L_{0}

. In particular, we use a simple linear regression in the interval

0.0 < λ < 0.3

, with the fluctuations data for

N_{0} = 10^{5}

, to find

χ_{T}^{\infty} = 0.0295 (5)

and

c = 0.415 (5) σ

. Results of the scaled fluctuations

λ χ_{T} (λ)

minus

c / L_{0}

are presented in Figure 5, where the intensive character of the constant c becomes clear. By replacing the calculated values

χ_{T}^{\infty}

and c in Equation (7), we obtain the black curve that superimposes on the simulation data in the full range

0 < λ < 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

λ χ_{T} (λ)

for

λ \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 < ρ σ^{3} < 1.0

, we compute the scaled fluctuations

λ χ_{T} (λ)

and determine, as before, the ratio

χ_{T}^{\infty} = κ_{T} / κ_{T}^{I G}

as a function of the density, with

κ_{T}^{I G} = {(ρ k_{B} T)}^{- 1}

the isothermal compressibility of the ideal gas (see Figure 6). As expected for this system at

k_{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:

κ_{T} = \frac{1}{ρ^{2}} {\frac{\partial ρ}{\partial μ}|}_{T},

(8)

which can be rearranged, in terms of the chemical potential

μ

, as:

δ μ = \int_{ρ_{0}}^{ρ} \frac{d ρ^{'}}{ρ^{' 2} κ_{T}}

(9)

with

δ μ = μ - μ_{0}

and

μ_{0}

the chemical potential of the system at the reference density

ρ_{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.):

δ μ^{ex} = δ μ - k_{B} T ln ρ,

(10)

obtained by subtracting from

δ μ

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

μ_{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

ρ σ^{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

ρ_{0} σ = 0.6

has been used as the reference value. Once

δ μ^{e x}

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

ρ

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{Δ^{2} (N)}{〈 N 〉} = 1 + \frac{ρ}{V} \int_{V} \int_{V} [g^{o} (r_{1}, r_{2}) - 1] d r_{1} d r_{2},

(11)

where

Δ^{2} (N) / 〈 N 〉

are the fluctuations of the number of particles,

Δ^{2} (N) = 〈 N^{2} 〉 - {〈 N 〉}^{2}

and

g^{o} (r_{1}, r_{2})

is the pair correlation function of the open system and

r_{1}, 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

〈 N 〉 \to \infty

with

ρ = 〈 N 〉 / V =

constant. An infinite, homogeneous and isotropic system is translationally invariant; therefore, we rewrite Equation (11) as [1]:

χ_{T}^{\infty} = \frac{Δ^{2} (N)}{〈 N 〉} = 1 + 4 π ρ \int_{0}^{\infty} (g^{o} (r) - 1) r^{2} d r,

(12)

with

χ_{T}^{\infty} = ρ k_{B} T κ_{T}

κ_{T}

being the isothermal compressibility of the bulk system. We have replaced

g^{o} (r_{1}, r_{2})

with

g^{o} (r)

the radial distribution function (RDF) of the open system, with

r = | r_{2} - 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:

χ_{T} (V, V_{0}) = \frac{Δ^{2} (N; V, V_{0})}{{〈 N 〉}_{V, V_{0}}} = 1 + \frac{ρ}{V} \int_{V} \int_{V} [g^{c} (r_{12}) - 1] d r_{1} d r_{2},

(13)

where

g^{c} (r_{12})

r_{12} = | r_{2} - r_{1} |

, is the pair correlation function of the closed system with total number of particles

N_{0}

, and

Δ^{2} (N; V, V_{0}) = {〈 N^{2} 〉}_{V, V_{0}} - {〈 N 〉}_{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

ρ

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{matrix} g^{c} (r) = & g^{o} (r) - \frac{χ_{T}^{\infty}}{N_{0}} . \end{matrix}

(14)

As a matter of fact, the expansion includes terms that depend on the partial derivative of

g^{o} (r)

with respect to the density. However, we anticipate here that for the present analysis, their contribution is negligible [6]. By replacing

g^{c} (r)

in the integral on the r.h.s of Equation (13), we obtain:

\frac{ρ}{V} \int_{V} \int_{V} (g^{c} (r_{12}) - 1) d r_{1} d r_{2} = I_{V, V} - \frac{V}{V_{0}} χ_{T}^{\infty},

(15)

where:

I_{V, V} = \frac{ρ}{V} \int_{V} \int_{V} (g^{o} (r_{12}) - 1) d r_{1} d r_{2},

(16)

and we use that

ρ = 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{matrix} I_{V, V_{0}} & = \frac{ρ}{V} \int_{V} \int_{V_{0}} (g^{o} (r_{12}) - 1) d r_{1} d r_{2} \\ I_{V, V_{0} - V} & = \frac{ρ}{V} \int_{V} \int_{V_{0} - V} (g^{o} (r_{12}) - 1) d r_{1} d r_{2} . \end{matrix}

As pointed out by Rovere, Heermann and Binder [17], the two integrals

I_{V, V}

and

I_{V, V_{0}}

are equal when

r_{1}

and

r_{2}

are both within the volume V. When

r_{12} > ζ

, the integrand

(g^{o} (r_{12}) - 1) = 0

, and it does not contribute to the integrals. Close to the boundary of the subdomain V, for

r_{12} < ζ

, and in particular when

r_{1}

lies inside and

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} + {(\frac{c_{2}}{L})}^{2} + O (\frac{1}{L^{3}}),

(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

ζ < L < L_{0}

. Since we assume PBCs, the system is translationally invariant. Hence, upon applying the transformation

r_{12} \to r = r_{2} - r_{1}

, the expression:

I_{V, V_{0}} = ρ \int_{V_{0}} (g^{o} (r) - 1) d r = χ_{T}^{\infty} - 1

(18)

is obtained, where we assume that

g^{o} (r > ζ) = 1

, thus ignoring fluctuations of the RDF beyond the volume V. By combining these two results, we obtain:

I_{V, V} = χ_{T}^{\infty} - 1 + \frac{c_{1}}{L} + {(\frac{c_{2}}{L})}^{2},

(19)

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

\frac{ρ}{V} \int_{V} \int_{V} (g^{c} (r_{12}) - 1) d r_{1} d r_{2} = χ_{T}^{\infty} (1 - {(\frac{L}{L_{0}})}^{3}) - 1 + \frac{c_{1}}{L} + {(\frac{c_{2}}{L})}^{2} .

(20)

Finally, this expression becomes:

χ_{T} (L, L_{0}) = χ_{T}^{\infty} (1 - {(\frac{L}{L_{0}})}^{3}) + \frac{c_{1}}{L} + {(\frac{c_{2}}{L})}^{2},

(21)

and by defining

λ = L / L_{0}

, we write:

λ χ_{T} (λ) = λ χ_{T}^{\infty} (1 - λ^{3}) + \frac{c_{1}}{L_{0}} + {(\frac{c_{2}}{L_{0}})}^{2} \frac{1}{λ} .

(22)

Equations (7) and (22) differ in the

c_{2}^{2} / L_{0}^{2} λ

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} / σ = 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 < ρ σ^{3} < 0.70

. Two temperatures were considered,

k_{B} T = 2.00 ϵ

and

1.15 ϵ

. The critical point of this system has been reported at

ρ_{c} σ^{3} = 0.319

and

k_{B} T_{c} = 1.086 ϵ

[38].

We report the reduced fluctuations

λ χ_{T} (λ)

as a function of

λ

for

ρ σ^{3} = 0.3

in Figure 8. In the case

k_{B} T = 2.00 ϵ

, the effect of the

λ^{- 1}

term in Equation (22) is negligible, and a linear approximation in the region

λ < 0.3

seems to be well justified. However, for the case close to the critical point, i.e.,

k_{B} T = 1.15 ϵ

, the effect of this term is evident and should be included in the extrapolation to

χ_{T}^{\infty}

Finally, upon extrapolating to

χ_{T}^{\infty}

, an interesting behavior is observed for the bulk isothermal compressibility

κ_{T}

as a function of density (Figure 9). In the case

k_{B} T = 2.00 ϵ

, as expected, a monotonically-decreasing behavior with increasing density is observed. More interestingly, in the case

k_{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_{i j}^{o} = V (\frac{〈 N_{i} N_{j} 〉 - 〈 N_{i} 〉 〈 N_{j} 〉}{〈 N_{i} 〉 〈 N_{j} 〉} - \frac{δ_{i j}}{〈 N_{i} 〉}) = \frac{1}{V} \int_{V} \int_{V} [g_{i j}^{o} (r_{12}) - 1] d r_{1} d r_{2},

(23)

with

δ_{i j}

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,

〈 N_{i} 〉

is the average number of i-particles inside V, or

ρ_{i} = 〈 N_{i} 〉 / V

g_{i j}^{o} (r_{12})

is the multicomponent radial distribution function (RDF) of the infinite system, with

r_{12} = r_{2} - 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_{i j} (L, L_{0}) = V (\frac{{〈 N_{i} N_{j} 〉}^{'} - {〈 N_{i} 〉}^{'} {〈 N_{j} 〉}^{'}}{{〈 N_{i} 〉}^{'} {〈 N_{j} 〉}^{'}} - \frac{δ_{i j}}{{〈 N_{i} 〉}^{'}}) = \frac{1}{V} \int_{V} \int_{V} [g_{i j}^{c} (r_{12}) - 1] d r_{1} d r_{2} .

(24)

The finite-size KBI

G_{i j} (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

{〈 N_{i} 〉}^{'} \equiv {〈 N_{i} 〉}_{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_{i j}^{c} (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_{i j}^{c} (r) = g_{i j}^{o} (r) - \frac{1}{V_{0}} (\frac{δ_{i j}}{ρ_{i}} + G_{i j}^{\infty}),

(25)

based on the asymptotic limit

g_{i j}^{c} (r ≫ ζ) = 1 - (δ_{i j} / ρ_{i} + G_{i j}^{\infty}) / V_{0}

discussed in [2]. As expected, when the total volume

V_{0} \to \infty

, we recover

g_{i j}^{c} (r) = g_{i j}^{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:

λ G_{i j} (λ) = λ G_{i j}^{\infty} (1 - λ^{3}) - λ^{4} \frac{δ_{i j}}{ρ_{i}} + \frac{α_{i j}}{L_{0}},

(26)

with

λ \equiv L / L_{0}

and

α_{i j}

an intensive parameter with units of length. In the limit

L_{0} \to \infty

, the following expression is obtained:

G_{i j} (L, L_{0} \to \infty) = G_{i j}^{\infty} + \frac{α_{i j}}{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} σ

. The potential parameters are chosen as

σ_{A A} = σ_{B B} = σ_{A B} = σ

, and

ϵ_{A A} = 1.2 ϵ

ϵ_{B B} = 1.0 ϵ

with

ϵ_{A B} = (ϵ_{A A} + ϵ_{B B}) / 2 = 1.1 ϵ

. All the results are expressed in LJ units with energy

ϵ

, length

σ

, mass

m_{A} = m_{B} = m

, time

σ {(m / ϵ)}^{1 / 2}

, temperature

ϵ / k_{B}

and pressure

ϵ / σ^{3}

. As before, simulations are carried out using ESPResSo++ [32] with a time step of

δ t / (σ {(m / ϵ)}^{1 / 2}) = 10^{- 3}

. Constant temperature

k_{B} T = 1.2 ϵ

is enforced through a Langevin thermostat with damping coefficient

γ (σ {(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, \dots, 1.0

. The pressure is fixed at

P σ^{3} / ϵ = 9.8

by adjusting the number density of the system at values around

ρ σ^{3} \approx 0.86

(or

L_{0} / σ \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_{i j} (λ)

, we select 800 frames per trajectory and for each frame identify 1000 randomly-positioned subdomains with linear sizes ranging from

2 < L / σ < L_{0} / σ

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_{A B}

(green circles) tend to zero when

λ \to 1

, as indicated by the horizontal green lines. By contrast,

G_{A A} \to 1 / ρ_{A}

(indicated by horizontal red lines) when

λ \to 1

. The region

λ < 3

, indicated by vertical black lines, is where simple linear regression is used to find

G_{i j}^{\infty}

and

α_{i j}

. 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 < λ < 1

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

κ_{T} = \frac{1 + ρ_{A} G_{A A} + ρ_{B} G_{B B} + ρ_{A} ρ_{B} (G_{A A} G_{B B} - G_{A B}^{2})}{k_{B} T (ρ_{A} + ρ_{B} + ρ_{A} ρ_{B} (G_{A A} + G_{B B} - 2 G_{A B}))},

(28)

with

ρ_{A, B}

the number density of the corresponding species.

Results for the isothermal compressibility obtained from the

G_{i j}^{\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

κ_{T} = (1 - x_{A}) κ_{T}^{B} + x_{A} κ_{T}^{A}

, with

κ_{T}^{A} ϵ / σ^{3} = 0.012 (1)

and

κ_{T}^{B} ϵ / σ^{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

ρ_{A}

using the expression [39]:

\frac{1}{k_{B} T} {(\frac{\partial μ_{A}}{\partial ρ_{A}})}_{P, T} = \frac{1}{ρ_{A}} + \frac{G_{A B} - G_{A A}}{1 + ρ_{A} (G_{A A} - G_{A B})},

(29)

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

δ μ_{A} = k_{B} T \int_{ρ_{A}^{0}}^{ρ_{A}} [\frac{1}{ρ_{A}^{'}} + \frac{G_{A B} - G_{A A}}{1 + ρ_{A}^{'} (G_{A A} - G_{A B})}] d ρ_{A}^{'} .

(30)

This is the chemical potential shifted by a reference chemical potential computed at density

ρ_{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:

δ μ_{A}^{e x} = δ μ_{A} - k_{B} T ln (x_{A}) - k_{B} T ln (ρ_{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.

Acknowledgments

We thank Debashish Mukherji and Roberto Menichetti for many stimulating discussions. We also thank Claudio Perego and Nancy C. Forero Martinez for a critical reading of the manuscript. Maziar Heidari, Kurt Kremer and Raffaello Potestio gratefully acknowledge funding from SFB-TRR146 of the German Research Foundation (DFG). Robinson Cortes-Huerto gratefully acknowledges the Alexander von Humboldt Foundation for financial support.

Author Contributions

Robinson Cortes-Huerto conceived the study. Kurt Kremer, Raffaello Potestio and Robinson Cortes-Huerto planned the computer simulations. Maziar Heidari and Robinson Cortes-Huerto carried out the numerical work. All authors discussed the results, helped with their interpretation and contributed to the final manuscript. All authors have read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

SBA	Spatial block analysis
TSLJ	Truncated and shifted Lennard–Jones
MD	Molecular dynamics
TL	Thermodynamic limit
PBCs	Periodic boundary conditions
OZ	Ornstein–Zernike
KB	Kirkwood–Buff

References

Hansen, J.P.; McDonald, I.R. Theory of Simple Liquids, 3rd ed.; Academic Press: Cambridge, MA, USA, 2006. [Google Scholar]
Ben-Naim, A. Molecular Theory of Solutions; Oxford University Press: Oxford, UK, 2006. [Google Scholar]
Adams, D. Chemical potential of hard-sphere fluids by Monte Carlo methods. Mol. Phys. 1974, 28, 1241–1252. [Google Scholar] [CrossRef]
Mukherji, D.; Kremer, K. Coil-Globule-Coil Transition of PNIPAm in Aqueous Methanol: Coupling All-Atom Simulations to Semi-Grand Canonical Coarse-Grained Reservoir. Macromolecules 2013, 46, 9158–9163. [Google Scholar] [CrossRef]
Wang, H.; Hartmann, C.; Schütte, C.; Site, L.D. Grand-Canonical-like Molecular-Dynamics Simulations by Using an Adaptive-Resolution Technique. Phys. Rev. X 2013, 3, 011018. [Google Scholar] [CrossRef]
Heidari, M.; Kremer, K.; Potestio, R.; Cortes-Huerto, R. Finite-size integral equations in the theory of liquids and the thermodynamic limit in computer simulations. Submitted.
Galata, A.A.; Anogiannakis, S.D.; Theodorou, D.N. Thermodynamic analysis of Lennard–Jones binary mixtures using Kirkwood–Buff theory. Fluid Phase Equilibria 2017, in press. [Google Scholar] [CrossRef]
Site, L.D.; Ciccotti, G.; Hartmann, C. Partitioning a macroscopic system into independent subsystems. J. Stat. Mech. Theory Exp. 2017, 083201. [Google Scholar] [CrossRef]
Strom, B.A.; Simon, J.M.; Schnell, S.K.; Kjelstrup, S.; He, J.; Bedeaux, D. Size and shape effects on the thermodynamic properties of nanoscale volumes of water. Phys. Chem. Chem. Phys. 2017, 19, 9016–9027. [Google Scholar] [CrossRef] [PubMed]
Dawass, N.; Krüger, P.; Schnell, S.K.; Bedeaux, D.; Kjelstrup, S.; Simon, J.M.; Vlugt, T.J.H. Finite-size effects of Kirkwood–Buff integrals from molecular simulations. Mol. Simul. 2018, 44, 599–612. [Google Scholar] [CrossRef]
Milzetti, J.; Nayar, D.; van der Vegt, N.F.A. Convergence of Kirkwood–Buff Integrals of Ideal and Nonideal Aqueous Solutions Using Molecular Dynamics Simulations. J. Phys. Chem. B 2018. [Google Scholar] [CrossRef] [PubMed]
Rogers, D.M. Extension of Kirkwood–Buff theory to the canonical ensemble. J. Chem. Phys. 2018, 148, 054102. [Google Scholar] [CrossRef] [PubMed]
Dawass, N.; Krüger, P.; Simon, J.M.; Vlugt, T.J.H. Kirkwood–Buff integrals of finite systems: Shape effects. Mol. Phys. 2018. [Google Scholar] [CrossRef]
Binder, K. Finite size scaling analysis of ising model block distribution functions. Z. Phys. B 1981, 43, 119–140. [Google Scholar] [CrossRef]
Binder, K. Critical Properties from Monte Carlo Coarse Graining and Renormalization. Phys. Rev. Lett. 1981, 47, 693–696. [Google Scholar] [CrossRef]
Rovere, M.; Hermann, D.W.; Binder, K. Block Density Distribution Function Analysis of Two-Dimensional Lennard–Jones Fluids. EPL 1988, 6, 585. [Google Scholar] [CrossRef]
Rovere, M.; Heermann, D.W.; Binder, K. The gas-liquid transition of the two-dimensional Lennard–Jones fluid. J. Phys. Condens. Matter 1990, 2, 7009–7032. [Google Scholar] [CrossRef]
Rovere, M.; Nielaba, P.; Binder, K. Simulation studies of gas-liquid transitions in two dimensions via a subsystem-block-density distribution analysis. Z. Phys. B 1993, 90, 215–228. [Google Scholar] [CrossRef]
Weber, H.; Marx, D.; Binder, K. Melting transition in two dimensions: A finite-size scaling analysis of bond-orientational order in hard disks. Phys. Rev. B 1995, 51, 14636–14651. [Google Scholar] [CrossRef]
Román, F.L.; White, J.A.; Velasco, S. Block analysis method in off-lattice fluids. EPL 1998, 42, 371. [Google Scholar] [CrossRef]
Salacuse, J. Particle fluctuations within sub-regions of an N-particle, three-dimensional fluid: Finite-size effects and compressibility. Phys. A 2008, 387, 3073–3083. [Google Scholar] [CrossRef]
Sengupta, S.; Nielaba, P.; Rao, M.; Binder, K. Elastic constants from microscopic strain fluctuations. Phys. Rev. E 2000, 61, 1072–1080. [Google Scholar] [CrossRef]
Lebowitz, J.L.; Percus, J.K. Long-Range Correlations in a Closed System with Applications to Nonuniform Fluids. Phys. Rev. 1961, 122, 1675–1691. [Google Scholar] [CrossRef]
Salacuse, J.J.; Denton, A.R.; Egelstaff, P.A. Finite-size effects in molecular dynamics simulations: Static structure factor and compressibility. I. Theoretical method. Phys. Rev. E 1996, 53, 2382–2389. [Google Scholar] [CrossRef]
Román, F.L.; White, J.A.; Velasco, S. Fluctuations in an equilibrium hard-disk fluid: Explicit size effects. J. Chem. Phys. 1997, 107, 4635. [Google Scholar] [CrossRef]
Villamaina, D.; Trizac, E. Thinking outside the box: Fluctuations and finite size effects. Eur. J. Phys. 2014, 35, 035011. [Google Scholar] [CrossRef]
Rowlinson, J.S. The equation of state of dense systems. Rep. Prog. Phys. 1965, 28, 169. [Google Scholar] [CrossRef]
Schnell, S.K.; Vlugt, T.J.; Simon, J.M.; Bedeaux, D.; Kjelstrup, S. Thermodynamics of a small system in a μT reservoir. Chem. Phys. Lett. 2011, 504, 199–201. [Google Scholar] [CrossRef]
Lebowitz, J.L.; Percus, J.K. Thermodynamic Properties of Small Systems. Phys. Rev. 1961, 124, 1673–1681. [Google Scholar] [CrossRef]
Hill, T.L. Thermodynamics of Small Systems; Courier Corporation: Dover, UK, 1963. [Google Scholar]
Humphrey, W.; Dalke, A.; Schulten, K. VMD—Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38. [Google Scholar] [CrossRef]
Halverson, J.D.; Brandes, T.; Lenz, O.; Arnold, A.; Bevc, S.; Starchenko, V.; Kremer, K.; Stuehn, T.; Reith, D. ESPResSo++: A modern multiscale simulation package for soft matter systems. Comput. Phys. Commun. 2013, 184, 1129–1149. [Google Scholar] [CrossRef]
Román, F.L.; González, A.; White, J.A.; Velasco, S. Fluctuations in the number of particles of the ideal gas: A simple example of explicit finite-size effects. Am. J. Phys. 1999, 67, 1149–1151. [Google Scholar] [CrossRef]
Román, F.L.; White, J.A.; González, A.; Velasco, S. Fluctuations in a small hard-disk system: Implicit finite size effects. J. Chem. Phys. 1999, 110, 9821. [Google Scholar] [CrossRef]
Román, F.; White, J.; González, A.; Velasco, S. Theory and Simulation of Hard-Sphere Fluids and Related Systems; Chapter Ensemble Effects in Small Systems; Springer: Berlin/Heidelberg, Germany, 2008; pp. 343–381. [Google Scholar]
Heidari, M.; Kremer, K.; Cortes-Huerto, R.; Potestio, R. Spatially Resolved Thermodynamic Integration: An Efficient Method to Compute Chemical Potentials of Dense Fluids. ArXiv, 2018; arXiv:cond-mat.soft/1802.08045. [Google Scholar]
Ornstein, L.S.; Zernike, F. Accidental deviations of density and opalescence at the critical point of a single substance. Proc. Akad. Sci. (Amsterdam) 1914, 17, 793–806. [Google Scholar]
Thol, M.; Rutkai, G.; Span, R.; Vrabec, J.; Lustig, R. Equation of State for the Lennard–Jones Truncated and Shifted Model Fluid. Int. J. Thermophys. 2015, 36, 25–43. [Google Scholar] [CrossRef]
Kirkwood, J.G.; Buff, F.P. The Statistical Mechanical Theory of Solutions. I. J. Chem. Phys. 1951, 19, 774–777. [Google Scholar] [CrossRef]
Pierce, V.; Kang, M.; Aburi, M.; Weerasinghe, S.; Smith, P.E. Recent Applications of Kirkwood–Buff Theory to Biological Systems. Cell Biochem. Biophys. 2008, 50, 1–22. [Google Scholar] [CrossRef] [PubMed]
Kjelstrup, S.; Schnell, S.K.; Vlugt, T.J.H.; Simon, J.M.; Bardow, A.; Bedeaux, D.; Trinh, T. Bridging scales with thermodynamics: From nano to macro. Adv. Nat. Sci. Nanosci. Nanotechnol. 2014, 5, 023002. [Google Scholar] [CrossRef]
Ben-Naim, A. Theoretical aspects of self-assembly of proteins: A Kirkwood–Buff-theory approach. J. Chem. Phys. 2013, 138, 224906. [Google Scholar] [CrossRef] [PubMed]
Mukherji, D.; Marques, C.M.; Stuehn, T.; Kremer, K. Depleted depletion drives polymer swelling in poor solvent mixtures. Nat. Commun. 2017, 8, 1374. [Google Scholar] [CrossRef] [PubMed]
Cortes-Huerto, R.; Kremer, K.; Potestio, R. Communication: Kirkwood–Buff integrals in the thermodynamic limit from small-sized molecular dynamics simulations. J. Chem. Phys. 2016, 145, 141103. [Google Scholar] [CrossRef] [PubMed]
Schnell, S.K.; Liu, X.; Simon, J.; Bardow, A.; Bedeaux, D.; Vlugt, T.J.H.; Kjelstrup, S. Calculating Thermodynamic Properties from Fluctuations at Small Scales. J. Phys. Chem. B 2011, 115, 10911. [Google Scholar] [CrossRef] [PubMed]
Ganguly, P.; van der Vegt, N.F.A. Convergence of Sampling Kirkwood–Buff Integrals of Aqueous Solutions with Molecular Dynamics Simulations. J. Chem. Theory Comput. 2013, 9, 1347–1355. [Google Scholar] [CrossRef] [PubMed]

Figure 1. Snapshot of the simulation box for a system of particles interacting via a TSLJ potential at density

ρ σ^{3} = 0.1

and temperature

k_{B} T = 1.2 ϵ

. In this particular example, a box of linear size

L_{0}

has been subdivided into blocks of linear dimension

L = L_{0} / 5

as indicated by the different color shades. The figure has been rendered with the Visual Molecular Dynamics (VMD) program [31].

Figure 1. Snapshot of the simulation box for a system of particles interacting via a TSLJ potential at density

ρ σ^{3} = 0.1

and temperature

k_{B} T = 1.2 ϵ

. In this particular example, a box of linear size

L_{0}

has been subdivided into blocks of linear dimension

L = L_{0} / 5

as indicated by the different color shades. The figure has been rendered with the Visual Molecular Dynamics (VMD) program [31].

Figure 2. Fluctuations of the number of particles

χ_{T} (L, L_{0})

as a function of

σ / L

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Data corresponding to system sizes

N_{0} = 10^{4}, 10^{5}

and

10^{6}

are presented using red squares, blue triangles and green circles, respectively. The vertical lines indicate the limit

σ / L_{0}

at which fluctuations become zero. The black horizontal dashed line indicates the value

χ_{T}^{\infty} = ρ k_{B} T κ_{T} = 0.0295

with

κ_{T}

the bulk compressibility obtained with the method described in [6].

Figure 2. Fluctuations of the number of particles

χ_{T} (L, L_{0})

as a function of

σ / L

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Data corresponding to system sizes

N_{0} = 10^{4}, 10^{5}

and

10^{6}

are presented using red squares, blue triangles and green circles, respectively. The vertical lines indicate the limit

σ / L_{0}

at which fluctuations become zero. The black horizontal dashed line indicates the value

χ_{T}^{\infty} = ρ k_{B} T κ_{T} = 0.0295

with

κ_{T}

the bulk compressibility obtained with the method described in [6].

Figure 3. Fluctuations of the number of particles

χ_{T} (L, L_{0})

as a function of the ratio

L / L_{0}

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Results corresponding to systems of

N_{0} = 10^{5}

particles with densities

ρ σ^{3} = 0.1, 0.2

and 0.3 are presented using red squares, blue triangles and green circles, respectively. The theoretical prediction presented in the text is plotted using the corresponding value for

χ_{T}^{\infty}

, obtained as described in [6], and solid-line curves with the same color code.

Figure 3. Fluctuations of the number of particles

χ_{T} (L, L_{0})

as a function of the ratio

L / L_{0}

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Results corresponding to systems of

N_{0} = 10^{5}

particles with densities

ρ σ^{3} = 0.1, 0.2

and 0.3 are presented using red squares, blue triangles and green circles, respectively. The theoretical prediction presented in the text is plotted using the corresponding value for

χ_{T}^{\infty}

, obtained as described in [6], and solid-line curves with the same color code.

Figure 4. Fluctuations of the number of particles

χ_{T} (L, L_{0})

as a function of the ratio

L / L_{0}

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Results corresponding to sizes

N_{0} = 10^{4}, 10^{5}

and

10^{6}

, with density

ρ σ^{3} = 0.864

, using red squares, blue triangles and green circles, respectively. The theoretical prediction presented in the text is plotted as the black dashed curve using

χ_{T}^{\infty} = 0.0295

Figure 4. Fluctuations of the number of particles

χ_{T} (L, L_{0})

as a function of the ratio

L / L_{0}

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Results corresponding to sizes

N_{0} = 10^{4}, 10^{5}

and

10^{6}

, with density

ρ σ^{3} = 0.864

, using red squares, blue triangles and green circles, respectively. The theoretical prediction presented in the text is plotted as the black dashed curve using

χ_{T}^{\infty} = 0.0295

Figure 5. Scaled fluctuations of the number of particles

λ χ_{T} (L, L_{0})

, minus

c / L_{0}

, versus the ratio

λ = L / L_{0}

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Results corresponding to sizes

N_{0} = 10^{4}, 10^{5}

and

10^{6}

, with density

ρ σ^{3} = 0.864

, using red squares, blue triangles and green circles, respectively. The theoretical prediction Equation (7) presented in the text is plotted as the black solid curve using

χ_{T}^{\infty} = 0.0295

and

c = 0.415 σ

Figure 5. Scaled fluctuations of the number of particles

λ χ_{T} (L, L_{0})

, minus

c / L_{0}

, versus the ratio

λ = L / L_{0}

for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Results corresponding to sizes

N_{0} = 10^{4}, 10^{5}

and

10^{6}

, with density

ρ σ^{3} = 0.864

, using red squares, blue triangles and green circles, respectively. The theoretical prediction Equation (7) presented in the text is plotted as the black solid curve using

χ_{T}^{\infty} = 0.0295

and

c = 0.415 σ

Figure 6. Ratio

χ_{T}^{\infty} = κ_{T} / κ_{T}^{I G}

k_{B} T = 1.2 ϵ

as a function of the density for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

, with

κ_{T}^{I G} = {(ρ k_{B} T)}^{- 1}

the isothermal compressibility of the ideal gas. The red curve is a guide to the eye.

Figure 6. Ratio

χ_{T}^{\infty} = κ_{T} / κ_{T}^{I G}

k_{B} T = 1.2 ϵ

as a function of the density for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

, with

κ_{T}^{I G} = {(ρ k_{B} T)}^{- 1}

the isothermal compressibility of the ideal gas. The red curve is a guide to the eye.

Figure 7. Excess chemical potential

μ^{e x} / ϵ

k_{B} T = 1.2 ϵ

as a function of the density for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. Red squares indicate the data obtained with the spatially-resolved thermodynamic integration (SPARTIAN) method [36], and the blue triangles are the data points obtained with the method outlined in the text.

Figure 7. Excess chemical potential

μ^{e x} / ϵ

k_{B} T = 1.2 ϵ

as a function of the density for systems described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

Figure 8. Reduced fluctuations as a function of

λ

for systems described by a TSLJ potential with

r_{c} / σ = 2.5

with density

ρ σ^{3} = 0.3

at temperatures

k_{B} T = 2.00 ϵ

and

1.15 ϵ

. For the latter case, it is apparent that the contribution proportional to

λ^{- 1}

is not negligible. The inset shows the full range

0 < λ < 1

. The black curves are the result of fitting the data to Equation (22).

Figure 8. Reduced fluctuations as a function of

λ

for systems described by a TSLJ potential with

r_{c} / σ = 2.5

with density

ρ σ^{3} = 0.3

at temperatures

k_{B} T = 2.00 ϵ

and

1.15 ϵ

. For the latter case, it is apparent that the contribution proportional to

λ^{- 1}

is not negligible. The inset shows the full range

0 < λ < 1

. The black curves are the result of fitting the data to Equation (22).

Figure 9. Bulk isothermal compressibility

κ_{T}

as a function of the density

ρ

k_{B} T = 1.15 ϵ

(red circles) and

k_{B} T = 2.00 ϵ

(green squares) for systems described by a TSLJ potential with

r_{c} / σ = 2.5

. The vertical black line indicates the location of the critical density

ρ σ^{3} = 0.319

[38].

Figure 9. Bulk isothermal compressibility

κ_{T}

as a function of the density

ρ

k_{B} T = 1.15 ϵ

(red circles) and

k_{B} T = 2.00 ϵ

(green squares) for systems described by a TSLJ potential with

r_{c} / σ = 2.5

. The vertical black line indicates the location of the critical density

ρ σ^{3} = 0.319

[38].

Figure 10. Scaled finite-size Kirkwood–Buff integrals

λ G_{i j} (λ)

as a function of

λ

for different mole fractions: (a)

x_{A} = 0.20

; (b)

x_{A} = 0.30

; (c)

x_{A} = 0.50

and (d)

x_{A} = 0.80

, for mixtures described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. For clarity, only the cases

G_{A A}

(red squares) and

G_{A B}

(green circles) are plotted. In the asymptotic case

λ \to 1

G_{A B} \to 0

and

G_{A A} \to 1 / ρ_{A}

, as indicated by the horizontal green and red lines, respectively. The black curves correspond to Equation (26) with

G_{i j}^{\infty}

and

α_{i j}

obtained from a simple regression analysis in the interval

λ < 0.3

Figure 10. Scaled finite-size Kirkwood–Buff integrals

λ G_{i j} (λ)

as a function of

λ

for different mole fractions: (a)

x_{A} = 0.20

; (b)

x_{A} = 0.30

; (c)

x_{A} = 0.50

and (d)

x_{A} = 0.80

, for mixtures described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. For clarity, only the cases

G_{A A}

(red squares) and

G_{A B}

(green circles) are plotted. In the asymptotic case

λ \to 1

G_{A B} \to 0

and

G_{A A} \to 1 / ρ_{A}

, as indicated by the horizontal green and red lines, respectively. The black curves correspond to Equation (26) with

G_{i j}^{\infty}

and

α_{i j}

obtained from a simple regression analysis in the interval

λ < 0.3

Figure 11. Isothermal compressibility at

k_{B} T = 1.20 ϵ

and

P σ^{3} / ϵ = 9.8

as a function of the mole fraction of type-A particles

x_{A}

for mixtures described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. The horizontal black lines indicate the compressibility for a pure system of type-A particles

κ_{T}^{A} ϵ / σ^{3} = 0.012 (1)

and for a pure system of type-B particles

κ_{T}^{B} ϵ / σ^{3} = 0.0281 (8)

. The red line is a guide to the eye. The ideal case corresponds to

κ_{T} = (1 - x_{A}) κ_{T}^{B} + x_{A} κ_{T}^{A}

Figure 11. Isothermal compressibility at

k_{B} T = 1.20 ϵ

and

P σ^{3} / ϵ = 9.8

as a function of the mole fraction of type-A particles

x_{A}

for mixtures described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

. The horizontal black lines indicate the compressibility for a pure system of type-A particles

κ_{T}^{A} ϵ / σ^{3} = 0.012 (1)

and for a pure system of type-B particles

κ_{T}^{B} ϵ / σ^{3} = 0.0281 (8)

. The red line is a guide to the eye. The ideal case corresponds to

κ_{T} = (1 - x_{A}) κ_{T}^{B} + x_{A} κ_{T}^{A}

Figure 12. Excess chemical potential of type-A particles as a function of the mole fraction

x_{A}

for mixtures described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

k_{B} T = 1.2 ϵ

and

P σ^{3} / ϵ = 9.8

. Data points obtained with the method in [36], in particular for

x_{A} = 0.3

, are used as a reference for the data points obtained with Equations (30) and (31).

Figure 12. Excess chemical potential of type-A particles as a function of the mole fraction

x_{A}

for mixtures described by a TSLJ potential with

r_{c} / σ = 2^{1 / 6}

k_{B} T = 1.2 ϵ

and

P σ^{3} / ϵ = 9.8

. Data points obtained with the method in [36], in particular for

x_{A} = 0.3

, are used as a reference for the data points obtained with Equations (30) and (31).

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Heidari, M.; Kremer, K.; Potestio, R.; Cortes-Huerto, R. Fluctuations, Finite-Size Effects and the Thermodynamic Limit in Computer Simulations: Revisiting the Spatial Block Analysis Method. Entropy 2018, 20, 222. https://doi.org/10.3390/e20040222

AMA Style

Heidari M, Kremer K, Potestio R, Cortes-Huerto R. Fluctuations, Finite-Size Effects and the Thermodynamic Limit in Computer Simulations: Revisiting the Spatial Block Analysis Method. Entropy. 2018; 20(4):222. https://doi.org/10.3390/e20040222

Chicago/Turabian Style

Heidari, Maziar, Kurt Kremer, Raffaello Potestio, and Robinson Cortes-Huerto. 2018. "Fluctuations, Finite-Size Effects and the Thermodynamic Limit in Computer Simulations: Revisiting the Spatial Block Analysis Method" Entropy 20, no. 4: 222. https://doi.org/10.3390/e20040222

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. Boundary and Ensemble Finite-Size Effects

3. Finite-Size Ornstein–Zernike Integral Equation

4. Mixtures

5. Summary and Concluding Remarks

Acknowledgments

Author Contributions

Conflicts of Interest

Abbreviations

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI