Finite-size integral equations in the theory of liquids and the thermodynamic limit in computer simulations

ABSTRACT We present an efficient method to obtain bulk isothermal compressibilities () and Kirkwood–Buff (KB) integrals of single- and multicomponent liquids using fluctuations of the number of molecules obtained from small-sized molecular dynamics simulations. We write finite-size versions of the Ornstein–Zernike and the KB integral equations and include there finite size effects related to the statistical ensemble and the finite integration volumes required in computer simulations. Consequently, we obtain analytical expressions connecting and the KB integrals in the thermodynamic limit (TL) with density fluctuations in the simulated system. We validate the method by calculating various thermodynamic quantities, including the chemical potentials of SPC/E water as a function of the density, and of aqueous urea solutions as a function of the mole fraction. The reported results are in excellent agreement with calculations obtained by using the best computational methods available, thus validating the method as a tool to compute the chemical potentials of dense molecular liquids and mixtures. Furthermore, the present method identifies conditions in which computer simulations can be effectively considered in the TL. GRAPHICAL ABSTRACT


I. Introduction
Statistical mechanics establishes the connection between macroscopic thermodynamic quantities and the microscopic interactions and components of a physical system. In particular, integral equations relate the local structure of a fluid with density fluctuations in the grand canonical ensemble that, in the thermodynamic limit (TL), can be identified with equilibrium thermodynamic quantities such as the compressibility and the derivatives of the chemical potential [1,2]. In computer simulations, this relation is routinely employed in spite of the fact that the systems under consideration are constrained to the canonical ensemble and usually far away from TL conditions. Hence, various finite-size effects need to be  [3][4][5] and their impact should be evaluated for an appropriate interpretation of the simulation results [6][7][8][9][10]. Building on our previous work [11], in this study we explicitly write finite-size versions of the Ornstein-Zernike (OZ) [12] and the Kirkwood-Buff (KB) [13] integral equations of the theory of liquids. By including in these equations size effects related to (i) the statistical ensemble used in computer simulations and (ii) the finite integration volumes required when the system is not in the TL, we derive analytical expressions to obtain isothermal compressibilities and KB integrals (KBIs) in the TL from relatively small-sized simulations. The accuracy of the method presented here allows one to derive other thermodynamic quantities of interest such as the chemical potential of molecular liquids, in particular SPC/E water and aqueous urea solutions, that well reproduce available simulation data for these systems.
The paper is organised as follows: In Section II we discuss the finite-size version of the OZ integral equation and the evaluation of the bulk isothermal compressibility of SPC/E water which allows one to compute the chemical potential as a function of the density. Section III extends these results to multicomponent systems, where we study urea-water liquid mixtures and obtain the KBIs in the TL and subsequently the chemical potential of urea as a function of mole fraction. In Section IV we summarise our results.

II. Simple liquids
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 OZ integral equation [1,12] 2 (N) where 2 (N)/ N are the fluctuations of the number of particles, 2 (N) = N 2 − N 2 and g • (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 (1), one assumes that the fluid is homogeneous, isotropic and that the system is in the TL, i.e. V → ∞, N → ∞ with ρ = N /V = constant. An infinite, homogeneous and isotropic system is translationally invariant, therefore we rewrite Equation (1) as [1]: with χ ∞ T = ρk B Tκ T , κ T being the isothermal compressibility of the bulk system. We have replaced g • (r 1 , r 2 ) with g • (r) the radial distribution function (RDF) of the open system, with r = |r 2 − r 1 |. Equation (2) is only valid in the TL. However, it is routinely used in computer simulations where closed systems with finite volume and fixed number of particles are usually considered. The modified version of Equation (2) used in simulations reads: where g c (r) is the RDF of the closed system and the infinity upper limit in the integral on the r.h.s. of Equation (2) is replaced by a radius ζ R < L with ζ the correlation length of the system and L the size of the simulation box. These implicit assumptions introduce two different finite-size effects: ensemble effects appear from assuming g o (r) = g c (r), which only holds in the TL. If R ζ then it is true that g • (r > R) = 1. By contrast even small fluctuations of g c (r > R) might contribute significantly to the integral on the r.h.s. of Equation (3). Boundary effects appear because in computer simulations the volume V in Equation (1) is finite and therefore the step leading to Equation (2) has to be handled with care [17]. As a matter of fact, the errors introduced by using Equation (3) yield significant deviations that become apparent when comparing standard molecular dynamics results with semi-grand canonical simulations [24]. An alternative OZ integral equation for finite-size systems was proposed a long time ago [14,16,22,25]. Let us consider a closed system with fixed number of particles N 0 and volume V 0 , with periodic boundary conditions (PBCs), where fluctuations of the number of particles are computed for subdomains of volume V ≤ V 0 . Therefore, we define [4] 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; The fluctuations of the number of particles thus depend on both subdomain and simulation box volumes. The question now is: is it possible to compute χ ∞ T using Equation (4) [26]?
To answer this question, we consider first the relation between the RDFs of an open system and its closed counterpart. The original idea is to write the grand canonical RDF as a sum of RDFs, in the canonical ensemble, weighted by the probability that the system has a number of particles N when in contact with a reservoir of particles [21,27]. 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 [21,23,[25][26][27] : In the following, we neglect the O(1/N 2 0 ) terms in Equation (5) because it has been demonstrated that, for sufficiently large systems, their contribution can be safely ignored [25]. By inserting Equation (5) into the integral on the right hand side (r.h.s.) of Equation (4) we obtain where and we use that ρ = N 0 /V 0 . At this point, we include explicitly the second finitesize effect, i.e. the fact that the volume V is finite. For this we rewrite I V,V as [17] To solve these equations, we recall that both integrals, I V,V and I V,V 0 , are equal when r 1 and r 2 are both within the volume V . Moreover, the integrand is zero for r 12 > ζ and therefore it does not contribute to the integral. However, for values of r 12 close to the boundary of the subdomain V, and in particular when r 1 lies inside V and r 2 outside with r 12 < ζ , there are contributions missing in I V,V that are present in I V,V 0 . Thus, the difference I V,V 0 −V = I V,V 0 − I V,V must be proportional to the surface-to-volume ratio S/V of the subdomain, as first discussed in [17], i.e.
with α , α > 0 proportionality constants that at this point, we assume to be intensive. S is the surface area of the subdomain with volume V [17].
To compute I V,V 0 , we require that V ζ < V < V 0 with V ζ = 4πζ 3 /3. Since we assumed PBCs, the system is translationally invariant. Hence, upon applying the transformation r 12 → r = r 2 − r 1 , we obtain where we assume that g • (r > ζ ) = 1 thus ignoring fluctuations of the RDF beyond the volume V . By combining these two results, we obtain and by including this result in Equation (6) we obtain The expression inside the square bracket on the r.h.s. is approximated to 1 since the term with the second derivative with respect to ρ is proportional to V −4/3 ≈ 0. Replacing in Equation (4) we obtain This expression is consistent with the two obvious limiting cases.
In the TL, V, V 0 → ∞, the correct value is recovered, i.e. χ T (λ = 0) = χ ∞ T . In the asymptotic limit λ = 1 we obtain χ T (λ = 1) = 0, that is, the fluctuations of the number of particles for a closed system are, by definition, equal to zero. The two size effects considered above can be easily identified in the previous expression. First, the term proportional to λ 3 takes into account the difference in the ensemble, i.e. g c (r) = g • (r). Second, the term proportional to V −1/3 accounts for the fact that the integration domain V is finite. We emphasise here that this last result can be equally obtained by neglecting the derivative with respect to the density of Equation (5), i.e. by directly using the expression: valid for large values of r and with an asymptotic correction independent of r.
Finally, for practical purposes [22], we multiply Equation (12) by λ and obtain: To make use of Equation (14), one runs a standard NVT simulation and selects a system with volume   (4) for subvolumes V such that V < V 0 , as schematically illustrated in Figure 1. Finally, one plots the results as a function of λ and extracts χ ∞ T from a simple linear regression in the range (V ζ /V 0 ) 1/3 < λ 1. We recall here that extrapolations of thermodynamic properties from finite-size computer simulations of single-component systems have been obtained in the past. Examples include the calculation of the compressibility of the Ising lattice gas [14] and hard disk fluids [22,29], and the investigation of the gas-liquid transition in two-dimensional Lennard-Jones fluids [16,18], as well as the calculation of the elastic constants of model solids [20].
We have used this strategy to calculate the isothermal compressibility of SPC/E [30] water. Molecular dynamics simulations have been carried out with GROMACS 4.5.1 [31] for 8000, 16,000 and 32,000 SPC/E [30] water molecules. To equilibrate the system, we started with a simulation box of size such that the initial density was ≈ 26 waters/nm 3 (≈ 776 kg/m 3 ). Then, optimisation of the system was achieved by performing steepest descent minimisation (50,000 steps are sufficient). Finally, an equilibration run of 3.5 ns has been carried out in the NPT ensemble at 1 bar. During this time, the displacement of a randomly chosen water molecule is comparable to half the linear size of the simulation box.
Following this minimisation step, we proceeded by alternating 3.5 ns (time step = 1 fs) constant pressure (NPT) at P = 1 bar and constant volume (NVT) simulations at T = 300 K. For NPT simulations we used a Berendsen barostat [32] whereas for NVT simulations temperature was enforced by a velocity rescaling thermostat [33]. We continued with this protocol until we verified that in the NPT ensemble the density is 33.5 waters/nm 3 (1000 kg/m 3 ) and that in the NVT simulation pressure fluctuates around the 1 bar value. For all the cases considered in the present study, the last NVT trajectory obtained after this sequence of NPT-NVT equilibration runs was used for the subsequent analysis in terms of the fluctuations of the number of molecules.
We evaluated λχ T (λ) as a function of λ for a system of 8000 water molecules. These data suggest a linear regime for λ < 0.3, consistent with Equation (14), where we used simple linear regression to the data to obtain χ ∞ T and α. We report χ ∞ T = 0.062 (3 % error), and by multiplying this value by the molar mass of water and dividing by ρk B T, we obtain 4.48 ×10 −5 bar −1 in excellent agreement with the widely accepted value for the isothermal compressibility of the SPC/E model (4.50 ×10 −5 bar −1 ) [34]. Figure 2 shows λχ T (λ) − ρα/V 1/3 0 as a function of λ for systems of 8000 (circles), 16,000 (triangles) and 32,000 (squares) water molecules. We replaced χ ∞ T and α back in Equation (14) to obtain the solid curve that superposes on the simulation data for the full interval 0 < λ < 1. Nevertheless, one should keep in mind that for very small values of λ, when the linear size of the subdomain is comparable to the size of the molecules, the statistical analysis is questionable. The three sets of data follow the same curve, indicating that α is indeed an intensive quantity but also suggesting that it is possible to accurately estimate χ ∞ T using the results of a relatively small simulation.
By neglecting the term proportional to λ 3 in Equation (14), we obtained the straight line in Figure 2 that indicates that deviations from linearity indeed become important for λ > 0.3. This limit λ 1, where a linear behaviour is apparent, is rather interesting. In this regime, the volume V V 0 , that is, the system is effectively in the grand canonical ensemble. If V 0 V ζ we obtain χ T (λ 1) = χ ∞ T . Namely, provided that V ζ < V < V 0 in the interval λ < 0.3, both V and V 0 are large enough such that the system seems to be in the TL. A way to challenge this argument involves the RDF: if the system is in the TL then g c (r) = g • (r) and Equation (3) gives χ ∞ T for R > ζ . We choose R = 1.5 nm as a reasonable value since the O-O RDF of the SPC/E model goes to 1 at approximately 1 nm. This choice gives V 1/3 = 1.5 × (4π/3) 1/3 nm, and with the sizes V 1/3 0 ≈ 6.2, 7.8 Figure 2. Fluctuations of the number of water molecules as a function of λ for systems of 8000 (circles), 16,000 (triangles) and 32,000 (squares) corresponding to the volumes indicated in the legends. Using Equation (14) and simple linear regression to the data for 8000 water molecules in the range λ < 0.3, we obtain the curve represented by the dashed line and the solid line (neglecting the λ 3 term in Equation (14)). and 9.9 nm we obtain λ ≈ 0.39, 0.31 and 0.24, respectively. For these systems, we compute g c (r) and evaluate χ R T from Equation (3). We obtain χ R T from the integral by averaging in the interval 1.2 < R < 1.5 nm. As expected, upon increasing the size of the simulation box (with λ = 0.39, 0.31 and 0.24 the integral takes the values 0.093 ± 0.004, 0.072 ± 0.003, 0.061 ± 0.003) the integral converges to the value of 0.062 computed with the present method (see Figure 3).
We conclude this part by saying that this method also allows us to identify simulation conditions in which the TL can be reached: in the case of SPC/E water, the volume at which the RDF is computed should be maximum of the order of 1% of the total volume of the simulation box.
To test the accuracy of the method, we have computed the isothermal compressibility of SPC/E water in the temperature interval 260-360 K. As shown in Figure 4 there is reasonable agreement with previous calculations using directly volume fluctuations [35][36][37]. More important, the efficient and accurate calculation of isothermal compressibility for dense liquids opens the possibility to compute other thermodynamic quantities of great interest. For example, the isothermal compressibility is related to the chemical potential μ via which can be rearranged in terms of μ as:  Isothermal compressibility for SPC/E water as a function of temperature at a fixed pressure (1 bar) using the method described in the text (circles) and direct volume fluctuations: triangles [35] and squares [36,37].
with δμ = μ − μ 0 and μ 0 the chemical potential of the system at the reference density ρ 0 . This expression suggests that, upon an accurate calculation of the compressibility for different densities at fixed temperature, it is possible to obtain the chemical potential of the system, shifted by a constant μ 0 . We have computed the compressibility of SPC/E water in the density interval 0.99 -1.05 g/cm 3 and used the expression (16) to obtain δμ. To estimate the shift and to prove that the method outlined here gives correct results, we have computed the excess chemical potential of SPC/E water in the same density interval using a method recently developed by us which can be interpreted as spatially resolved thermodynamic integration (SPARTIAN) [38]. To obtain the excess chemical potential, we subtract the densitydependent part of the chemical potential of the ideal gas from δμ, i.e.
Results are presented in Figure 5, and upon displacing the point corresponding to 0.995 g/cm 3 we find that both data sets overlap remarkably well within error bars. In this way, we have shown that with the strategy presented here it is possible to compute efficiently the chemical potential of dense liquids. When used in combination with another method that allows one to calculate the chemical potential at a reference state point, it also brings results consistent with state-of-the-art chemical potential calculations. Since in this context the chemical potential is obtained via the simple evaluation of fluctuations of the number of particles, we provide here an alternative strategy useful, in particular, in cases where the high density of the liquid becomes problematic for methods based on particle insertion or thermodynamic integration.

III. Multicomponent systems
From the previous analysis, the generalisation to multicomponent systems is straightforward [11]. For a multicomponent fluid of species i,j in equilibrium at temperature T, analogous to the OZ integral equation, the KBI is defined as [13]: with δ ij the Kronecker delta. As before, the superscript (o) indicates that this definition holds for an open system, i.e. a system described by the grand canonical ensemble. g • ij (r 1 , r 2 ) is the multicomponent pair correlation function of the open system. KB integrals connect the local structure of a multicomponent system to density fluctuations calculated in the grand canonical ensemble, that, in theTL, are related to equilibrium quantities such as the isothermal compressibility and derivatives of the chemical potential [2]. Recently, and to illustrate the robustness of this framework, KB integrals have been used to investigate the thermodynamics of Lennard-Jones binary mixtures [39], solvation of biomolecules [40], diffusion in multicomponent liquids [41,42], complex phenomena in molecular mixtures [43] and denaturation [44] and self-assembly [45] of proteins.
As in the single-component case, by assuming that the fluid is homogeneous and isotropic we obtain the KB integrals in the TL: where g • ij (r) is the multicomponent RDF of the infinite system. However, in computer simulations, the expression is widely used. Here, g c ij (r) is the multicomponent RDF of the finite system and R ζ is a cut-off distance. We explicitly include both finite-size effects: ensemble and boundary effects, into an expression, analogous to Equation (18), defined for a finite system. Let us consider systems with a total fixed number of particles N 0 and volume V 0 with PBCs. We define [4] with g c ij (r 12 ), r 12 = r 2 − r 1 , the pair correlation function of the closed system, and N i ≡ N i V,V 0 the average number of i-particles that depends on V and V 0 . Thus, we define a quantity G ij (V, V 0 ) that can be evaluated by computing fluctuations of the number of particles in finite subdomains of volume V inside a simulation box of volume V 0 .
We have proposed an expression, valid for large separations r, relating the open and closed RDFs of the form [11]: based on the asymptotic limit g c ij (r → ∞) = 1 − (δ ij /ρ i + G ∞ ij )/V 0 discussed in [2]. As expected, when the total volume V 0 → ∞ we recover g c ij (r) = g • ij (r). We include Equation (22) in the integral on the r.h.s. of Equation (21) and obtain: In analogy to the single-component case, in Equation (22) we neglect O(1/V 0 ) contributions that include derivatives of g • ij (r) with respect to ρ i , as well as O(1/V 2 0 ) terms, because their contribution to ensemble size effects becomes negligible upon integration.
Second, for the boundary effects, by explicitly taking into account the finite domain of integration V, we obtain [4,11]: where α ij is a constant that depends only on intensive thermodynamic system properties such as density and temperature. Finally, we rewrite Equation (24) as [11] with λ = (V/V 0 ) 1/3 . And like in the single-component case, we can use Equation (25) to extract the KBIs in the TL, G ∞ ij . In the case of multicomponent fluids, the explicit inclusion of finite-size effects to obtain G ∞ ij has been attempted in the literature by either using arguments from the thermodynamics of small systems [46] where only the 1/V 1/3 term appears [47] or by introducing empirical tail corrections to the RDF and evaluating Equation (20) with a modified kernel [48,49].
Here we identify the relevant finite-size effects at play and use them in a simple and physically sound framework that enables an accurate and efficient calculation of G ∞ ij [11]. Recently, the accuracy of this method has been tested by computing several thermodynamic properties of Lennard-Jones mixtures showing good agreement with theoretical predictions [39].
We rewrite Equation (25) as To validate Equation (26), we have performed simulations of aqueous urea solution [50,51] using the KBderived force field [52] and SPC/E water [30] in GRO-MACS 4.5.1 [31] with a relatively small size of the simulation box (L ∼ 8 nm). In addition to the trajectories from [51] we have considered four more molar concentrations for a total of seven molar concentrations: 2.00, 3.06, 3.90, 5.07, 6.03, 7.10 and 8.03. To equilibrate the system, we have alternated 3.5 ns (time step = 1 fs) constant pressure (NPT) at P = 1 bar and constant volume (NVT) simulations at T = 300 K. For NPT simulations, we used a Berendsen barostat [32], whereas for NVT simulations, temperature was enforced by a velocity rescaling thermostat [33]. We continued with this protocol (38 NPT-NVT cycles) until we verified that in the NVT simulation pressure fluctuates around 1 bar. The last NVT trajectory obtained after this NPT-NVT equilibration sequence was used for the analysis in terms of fluctuations of the number of molecules. Figure 6 shows fluctuations of the number of molecules of species i,j = urea (U), water (W) as a function of λ for a molar concentration of 2.00 M. Triangles correspond to λG UW (λ) that goes to zero when λ → 1, as the horizontal solid line indicates. Circles represent λG UU (λ) that in the limit λ → 1 goes to 1/ρ U = 1/1.208 = 0.828 as indicated by the horizontal dashed line. As in the case of single-component systems, here it is also clear that there is a linear region in the interval 0 < λ < 0.3, indicated by the vertical solid line ( see inset in Figure 6). From a simple linear regression in this region, G ∞ UW , G ∞ UU as well as the values of α ij are obtained. By inserting these values in Equation (26) we obtain the solid curves that overlap the simulation data points in all cases.
Since we can compute g c (r) directly from our trajectories, we compare the fluctuations results with the evaluation of the integral Equation (20). Results for G R UU , G R UW with R = 1.5 nm are presented in Figure 7 using solid curves. To obtain a value for the KBIs we have averaged these curves in the interval 1.3 -1.5 nm as indicated by the dotted horizontal lines. In addition, the dashed lines correspond to the values G ∞ UU and G ∞ UW obtained from the simple linear regression of Equation (26). Differences in G ∞ UU and the average from the curve are significant. The reason behind such differences lies in the fact that the concentration of urea is low, C U = 2.00M, and a larger  simulation box is needed to correctly compute the corresponding RDF. Conversely, the differences for G UW are negligible because the size of the system is large enough to accurately obtain the urea -water RDF.
In our previous paper, we have used these tools to investigate solvation thermodynamics [11]. In ureawater mixtures at pressure P, temperature T and number density ρ U the derivative of the urea activity coefficient can be written as with γ U the activity coefficient and k B T ln γ U the chemical potential of urea. We have computed γ UU as a function  [52,53] solid curve and Exp. 2 [56] dashed line) are also presented for comparison. Concerning other simulation methods, the triangles have been obtained by using Equation (20) [50,51]. Recently, the activity coefficient has been computed by using a modified version of Equation (20) with empirical tail corrections to the RDF. The squares correspond to one data set extracted from [54] with the smallest error bars and that best reproduces the experimental results. For the latter, as well as for our data, the error bars are comparable to the symbol size.
of molar concentration and report the results in Figure 8. These results are in good agreement with the experimental measurements [53] used to parameterise the force field [52]. Recently, a modified version of Equation (20) with empirical tail corrections to the RDF have been used to compute γ UU [54]. To compare with our method, we extract [55], among the various examples provided [54], the data set with the smallest error bars and that best reproduces the experimental results. The results presented in Figure 8 suggest that a simple calculation of the fluctuations of the number of particles used in combination with Equation (26) produce results consistently in better agreement with Exp. 1 [53].
More important, we can also use KBIs to compute the chemical potential of urea in urea-water mixtures. We use the expression , (28) and as in the case of single-component fluids, we integrate this equation to obtain a shifted chemical potential: where we have dropped the subindices P,T to lighten the notation. Using a combination of different methods, the chemical potential of urea in urea-water mixtures as a function of mole fraction has been computed in [57]. In Figure 9 we report these results in the interval 0,0.20 using triangles. We have obtained the shifted chemical potential of urea as a function of mole fraction using Equation (29) and plot the results using the data for the mole fraction x U = 0.13 (C U = 6.03) as a reference state point (circles). The two data points are in excellent agreement. This suggests that the method presented here constitutes an efficient and accurate tool to compute the chemical potential of liquid mixtures for a wide range of concentrations.

IV. Conclusions
During the last decades, we have witnessed an exponentially increasing computational capacity that allowed us to simulate enormous systems. In some cases, like the SPC/E example with V 0 = (9.9 nm) 3 discussed above, the system could practically be considered in the thermodynamic limit. Usually, this situation constitutes the exception rather than the rule, and the computation of bulk thermodynamic quantities from small-sized simulation results remains a challenging task. This limitation did not prevent pioneer computational scientists to investigate these properties for various complex systems. On the contrary, researchers have cleverly used the finite-size effects present in their systems to extrapolate interesting quantities in the TL.
Inspired by these works, we used finite-size variants of integral equations in the theory of liquids, namely for the OZ and KBI equations. By explicitly including ensemble and boundary size effects, we have derived analytical expressions that allow one to accurately calculate isothermal compressibilities and KBIs in the TL. In particular, the accurate estimation of these quantities for various density/concentration conditions allows one to determine, upon integration, the chemical potential for molecular liquids and mixtures.
This protocol to compute chemical potentials simply relies on running a standard computer simulation. Provided that the linear dimension of the simulation box is larger than the correlation length of the system, the obtained trajectories can be easily used to compute the fluctuations of the number of molecules for different subdomain sizes. This can be carried out during the postprocessing and without any external sophisticated computational strategies, in contrast with the standard methods available in the literature.
The efficiency and accuracy of the method have been demonstrated by the calculation of shifted chemical potentials of SPC/E water as a function of the density, and aqueous urea solution as a function of the mole fraction. The reported results exhibit trends in excellent agreement with state-of-the-art calculations. The reference chemical potential can be calculated using an alternative, wellestablished method for a single density/concentration value. Thus, the simple estimation of density fluctuations, as described here, provides a powerful tool to determine the chemical potential of molecular fluids for a wide range of density/concentration conditions.