Extension of Kirkwood-Buff Theory to the Canonical Ensemble
David M. Rogers
cond-mat.stat-mech




 University of South Florida, 4202 E. Fowler Ave., CHE 205, Tampa,
FL 33620, USA
 davidrogers@usf.edu




   Kirkwood-Buff (KB) integrals are notoriously difficult to converge from a canonical
simulation because they require estimating the grand-canonical radial distribution.
The same essential difficulty is encountered when attempting to estimate the direct
correlation function of Ornstein-Zernike theory by inverting the pair correlation functions.
We present a new theory that applies to the entire, finite, simulation volume,
so that no cutoff issues arise at all.
The theory gives the direct correlation function for closed systems, while
smoothness of the direct correlation function in reciprocal
space allows calculating canonical KB integrals via a well-posed extrapolation to the origin.
The present analysis method represents an improvement over previous work because
it makes use of the entire simulation volume and its convergence can be accelerated
using known properties of the direct correlation function.
Using known interaction energy functions can make this extrapolation
near perfect accuracy in the low-density case.
Because finite size effects are stronger in the canonical than the grand-canonical
ensemble, we state ensemble correction formulas for the chemical potential
and the KB coefficients.
The new theory is illustrated with both analytical and simulation results on the
1D Ising model and a supercritical Lennard-Jones fluid.  For the latter, the
finite-size corrections are shown to be small.




 Extension of Kirkwood-Buff Theory to the Canonical Ensemble
     David M. Rogers
    December 30, 2023
============================================================





§  INTRODUCTION


  Kirkwood and Buff's theory of solutions exploits a Gibbs relation to
determine derivatives of activities from fluctuations in particle number at constant volume,

    ( βμ_αn_γ)_n_η,V,T
      = [n_α' n_γ' | μ,V,T - V^2 ρ_α'ρ_γ']_αγ^-1,

with the notation ρ_α = n_α | μ, V, T/V, the average number of molecules
of type α in the grand-canonical ensemble (where μ, V, T are given),
and using μ∈ℝ^ν to
denote a vector of chemical potential values and n ∈ℕ^ν, a
vector of molecule numbers for all ν solution components.
The inverse indicated above is the Moore-Penrose pseudoinverse of the ν×ν matrix of
number fluctuations.
KB theory is exact even when applied to periodic simulation volumes with finite unit cells,
but requires an open system (i.e. data from the μ,V,T ensemble).
In this work, we consider only such finite unit cells and use the additional
assumption that no fixed external potential is present that breaks the
translational symmetry.  We introduce a theory that applies equally to open and closed systems
and then investigate its zero-frequency limit (which is only rigorous for the open system).
Also, we maintain constant temperature throughout
and hence omit it in further notation.

  KB theory is commonly cast in terms of the grand-canonical radial distribution function,

    g_αγ(r | μ, L) ≡ρ̂_α(r) ρ̂_γ(0) | μ, L - δ_αγδ(r) ρ_α/ρ_αρ_γ,

whose integral reports on number fluctuations in the grand-canonical ensemble.
Here, r ∈ℝ^3 / L is a coordinate vector in the periodic unit cell
with lattice vectors given as rows of a 3×3 matrix, L,
and

    ρ̂_α(r) ≡∑_i=1^n_αδ(r - r_α,i)

 is a sum of delta-functions located at the centers
of mass for all particles of type α (so that
∫_V ρ̂_α(r) dr = n_α and ρ_α = ρ̂_α(0) | μ, L).
Directly integrating Eq. <ref> gives Kirkwood and Buff's Eq. 7,<cit.>
confirming that it applies to finite volumes with the entire volume as the domain of integration.

  A recurring problem is that KB integrals estimated from closed systems (i.e. those under
NVT or NPT conditions) contain substantial finite-size effects.<cit.>
Proposed corrections fall into two categories: those aimed at computing number
density fluctuations directly in small sub-volumes,<cit.>
and those addressing the cutoff error in truncating the KB integral
at a finite maximum radius.<cit.>
Both categories recognize that two independent corrections are necessary.
First, the sub-volume or truncation radius is typically extrapolated
to infinite size.  Second, the radial distribution (or equivalently the number correlations)
themselves must be corrected for the effect of using a closed system.
These corrections were defined as implicit and explicit by
Salacuse et. al. <cit.>, who derived and compared
calculations of the closed system correction.
However, the finite truncation radius remained as a source of ambiguity,
which more recent works have addressed with renewed interest.<cit.>

  In this work, we re-introduce the structure factor as the preferred method for computing both the
direct correlation function and the KB integral as its zero-frequency limit.
This corrects a minor defect in Ref.  by using the unique vectors in reciprocal space
that correspond to a finite simulation volume, and yields exact direct correlation functions
corresponding to the canonical ensemble simulated.
For the problem of estimating the thermodynamic limit, our method improves on earlier works as it
remains invariant to shifts of g(r|n, L) by a constant,
and makes use of the entire available simulation volume.
In the process, we find two major conclusions.
First, extrapolation of the direct correlation function for closed systems
to zero frequency gives good estimates
of the second derivatives of the Helmholtz free energy.
Second, both closed (canonical) and open (grand-canonical) KB integrals
can be estimated under near-critical conditions from box lengths on the order
of 5 to 10 times the correlation length.



§  THEORY




 §.§  Definitions

  We introduce the structure factors using Fourier transforms of the density operators,

    Ŝ_α(m)    ≡ℱ[ρ̂_α](m) ≡∫_V ρ̂_α(r) e^-2π i m· r  dr
    .
    Explicit formulas for these transforms are,Ŝ_α(m)    = ∑_j=1^n_α e^-2π i m· r_α,j,  ρ̂_α(r) = 1/V∑_m Ŝ_α(m) e^2π i m· r
    ,

where the sum over m runs over the infinite reciprocal
lattice with vectors m = L^-1 u (u∈ℤ^3),
and the volume of the periodic unit cell is V = |L|.
These sums can be computed efficiently using the fast Fourier
transform in molecular simulations.<cit.>
The results in this work were computed with our own
implementation of the method,
which is available with documentation in Ref. .

  For any fixed-volume ensemble with partition function Z(Γ,L)
(e.g. Γ = n,T for canonical or Γ = μ,T for grand-canonical),
and probability distribution P(τ | Γ,L) (where τ
represents the combined positions and momenta of all molecules present),
define an extended ensemble that characterizes its response
to a set of externally applied fields, ϕ = {ϕ_α}, using the partition function,

    Θ(ϕ, Γ, L)    = Z(Γ, L) ∫  
         e^ρ̂, ϕ  P(τ | Γ, L) dτ
    X,Y   ≡∑_α∫_V X_α(r)^† Y_α(r) dr 
       = 1/V∑_α∑_m ℱ[X_α](m)^†ℱ[Y_α](m)
    .

The last two lines define a convenient notation for inner products
between two functions of r in the unit cell and state the appropriate
Poisson summation formula.<cit.>
We also define the function,

    ϕ̃_α(m)    ≡ℱ[ϕ_α](m) / V 
    to absorb the normalization constant in Eq. <ref>.  This way,ρ̂, ϕ   = ∑_α∑_m Ŝ_α^†(m) ϕ̃_α(m), 
    and ϕ̃_α(0) / β   = μ_α,

the chemical potential of species α.
The first and second derivatives of lnΘ generate the densities and correlations,

    lnΘϕ̃_α(m)^†   = Ŝ_α(m)≡ S_α(m)
    
    ^2 lnΘϕ̃_α(m)^†ϕ̃_γ(m')   = S_α(m)ϕ̃_γ(m')
       = V δ_m,m' Q_αγ(m) 
    
    Q_αγ(m)    ≡1/VΔŜ_α(m) ΔŜ_γ(m)^†
       = ℱ[ρ̂_α(r) ρ̂_γ(0)](m) - δ_m,0Vρ_αρ_γ,

where ΔŜ_α(m) ≡Ŝ_α(m) - S_α(m).
All quantities in Eqns. <ref>-<ref> are conditional on the macrostate (ϕ, Γ, L)
and evaluated at points where ϕ is spatially uniform.
The vanishing of m  m' cross-derivatives in Eq. <ref> can be proven
from translational symmetry in this case.

  The matrix Q/ρ is the structure factor of scattering theory.<cit.>
It is connected to the Fourier transform of the radial distribution function
via Eqns. <ref> and <ref>.
It is also related to the integral theory of solutions via the definition,<cit.>

    Q^-1_αγ(r) = δ(r) βμ^id_αρ_γ - c_αγ(r)
    ,

of the Ornstein-Zernike direct correlation function, c_αγ(r).
Here, μ_α^id is the chemical potential of an ideal gas,
whose derivative is ∂βμ_α^id/∂ρ_α = 1/ρ_α.



 §.§  Canonical Kirkwood-Buff Coefficients


  It is well-known that the matrix inverse of Eq. <ref>
gives the shift in the potential required to maintain a small
change in the density profile appropriate for an ensemble with fixed (S, Γ, L).
Eqns. <ref>-<ref> apply exactly to both canonical (Γ = n,T)
and grand-canonical (Γ = μ,T) ensembles.

  From this connection, the most relevant quantities to canonical KB theory
are the fit coefficients, c_0 and c_2, for

    lim_m→ 0 Q_αγ^-1(m | n,L) = βμ^id_αρ_γ
    + c_0,αγ(0) + c_2,αγ (2π m)^2 + O(m^4)
    .

There is no linear term since Q must be an even function of m by symmetry.
For a sufficiently large number of molecules of type α, the first
coefficient provides the density dependence of their excess chemical potential,

    lim_L /L_c →∞ c_0,αγ(n, V) = β V (
    μ^ex_α(n,V)
    - μ^ex_α(n - δ_γ,V)
    )
    ,

where we define

    βμ^ex_α(n,V) ≡ -ln e^-βΔ U_α | n, V
    ,

to be the test particle insertion free energy in the canonical ensemble.<cit.>
The coefficient c_0,αγ
also converges to ∂βμ_α^ex / ∂ρ_γ
in the thermodynamic limit.
Convergence to the particle insertion free energy
is faster though, because the Helmholtz free energy
can be rigorously stated as the sum of n successive test particle insertions,

    β A(n,V,T) = β A^id(n,V,T) + ∑_j=0^n-1βμ^ex(j,V,T)
    ,

while the ideal term is indistinguishable from n lnρ after n ≥ 100.

  As further justification of the limit, we invoke the well-known
result that density fluctuations must become uncorrelated as the
wavelength, 1/m grows much larger than the correlation length.<cit.>

Thus, under these conditions Q^-1(m) is continuous as m approaches zero
and its extrapolation becomes easier as the
cell size increases and the range of interactions becomes shorter.

  Away from m=0, the matrix Q reports on structural features that cannot be calculated
from bulk thermodynamic data.  The second coefficient, c_2, in Eq. <ref>
provides the correlation length for the αγ radial distribution,<cit.>

    L_c^(αγ) = √( |
    Q_αγ(0) c_2,αγ |)
    .

The coefficient c_2 is generally negative for oscillatory decay
and positive for monotonic decay of the radial distribution, g(r).


  The identification of Eq. <ref> as the proper limiting form
for the canonical ensemble allows us to carry through the rest of the
Kirkwood-Buff analysis unchanged.
All the other derived quantities of Kirkwood-Buff theory can be stated by noting that
the matrix Q_αγ(0|μ, L) is identical to
B_αγ defined in Eq. 9 of Ref. .
The supplementary information provides explicit expressions.
In the thermodynamic limit, the density-derivatives of the chemical potential
are the same for the canonical and grand-canonical ensembles.
However, at finite volume, the grand-canonical ensemble provides
observables closer to the thermodynamic limit than the canonical.

  From this identification, we conclude that the most effective method
for estimating KB integrals from canonical simulations is to use density
fluctuations over the entire simulation cell to estimate Q(m) and then
to extrapolate the value of Q(0).  Our numerical results below
indeed show this works, but caution that there are two sources
of error.  The first is due to the difficulty of extrapolating to m=0.
The second is due to finite size corrections, which become increasingly
important for successive density derivatives.
This work shows that the first type of error can be entirely eliminated
within the canonical ensemble when the simulation box size is larger than
about ten correlation lengths (which can be easily estimated using Eq. <ref>).










 §.§  Finite Size Corrections


  Finite size corrections (which extrapolate a result at finite volume to
the thermodynamic limit) are more severe for the canonical
than the grand-canonical ensemble.  This section derives
relations connecting averages in the canonical and grand-canonical
ensembles.  It is usually assumed that the grand-canonical
value is close enough to the thermodynamic limit that this ensemble
correction is a good estimate of the finite size correction.

  The correction to the distribution function was given by Lebowitz <cit.>,

    Q_αγ   (m | μ, L) - Q_αγ(m | n, L) 
       = 1/2V∑_ηζ Q_ηζ(0 | μ, L)
          ^2ρ_η∂ρ_ζ Q_αγ(m | μ, L)
      ,

which is correct to relative order 1/V unless m=0.
At the point m=0, Q_αγ(0 | n, L) = 0, and the entire
contribution to Q_αγ(0 | μ, L) comes from the second term.
Eq. <ref> was shown to match the thermodynamic
limit in Ref. .

  Using the same method, we derive in the appendix the correction
to the excess chemical potential valid to first order in 1/V,

    βμ_α^ex   (μ,V) - βμ_α^ex(n,V) 
       =
       1/2V(
             2/ρ_α - Q^-1_αα(0)
             + ∑_γQ^-1_αγβμ_γ(0)
        ).

For a 1-component system, Q(0) = ρ^2κ / β is proportional to the isothermal
compressibility, κ, and the correction (Eq. <ref>) specializes to
the derivative of the excess compressibility,

    βΔμ^ex(μ,V) = -1/2Vκ( κ - β/ρ )ρ
    .

This result agrees with Ref. , where it was first derived and tested
against the thermodynamic limit for hard-sphere fluids.



§  LATTICE GAS MODEL


  To illustrate the convergence properties of the structure factor approach,
we apply it to a simple one-dimensional, one-component periodic lattice gas
on L sites with nearest-neighbor interaction energy β J
This system is isomorphic to the 1D Ising model.<cit.>
The operators, ρ̂, are defined as indicator functions
on the L sites, and all Fourier transforms
become discrete Fourier transforms over L.  Otherwise, the theory above
goes through as expected except for a change in the
ideal chemical potential, <cit.>

    βμ^idρ = 1/ρ(1-ρ)
    .




  This system, along with several generalizations, were shown to have a strictly nearest-neighbor
direct correlation function in a series of original studies,<cit.> 
yet it still 
poses a nontrivial estimation problem since the correlation length
goes to infinity at strong coupling.
It is ideal for illustrating estimation of the direct correlation function
because the densities and pair distributions can be
computed analytically in the grand-canonical ensemble using the transfer
matrix method.<cit.>
We also calculated the canonical ensemble densities and pair distributions
exactly with the use of the absolutely convergent cluster summation technique
of Ref. .
The most important properties of this model depend on
the parameter y ≡±exp(-1/L_c), where L_c is the correlation length
over which the radial distribution functions decay exponentially.
The complete expression for y in terms of β J and ϕ (= βμ)
is given in the supplementary information.
It can vary from -1 in an antiferromagnetic system when
β J → -∞ to +1 in a ferromagnetic system when β J → +∞.



  Fig. <ref> compares the radial distribution function computed in the canonical
ensemble with the corresponding grand-canonical distributions
(at the chemical potential ϕ = -β J, for which n | μ, L = L/2).
The convergence is slow with increasing L, even though the canonical correlation functions
appear to flatten out to a constant at large distance.
These shifts in g(r | n, L) by a constant do not affect our analysis,
since they merely change the structure factor at m=0 –
which our analysis disregards.  The interesting ensemble effect
shown in this plot is that the r=1 point appears to have a different vertical
shift than the large-distance asymptote.
For the antiferromagnetic case (a), the short-range behavior is relatively insensitive to
the ensemble correction.  In contrast (and disregarding the smallest simulation size),
ensemble dependence shows up most strongly
at short range in the ferromagnetic case (b).

  Fig. <ref> compares the inverse correlation function, Q^-1(m) for the
three systems above in the canonical and grand-canonical ensembles.
The function Q^-1 can be calculated exactly and simplifies
in the limit as L/L_c →∞,



    Q^-1(m|μ) = 1/ρ(1-ρ)(
       1-y/1+y - 2 y/1-y^2( cos(2π m)-1 )
     )
     .


  Fig. <ref> compares two estimates for the KB coefficients
corresponding to extrapolating to m=0 following Eq. <ref>
(panels a,b) and correcting for ensemble dependence, fitting every two successive
lengths to an empirical expression, 

    Q^-1_extrap(1/L) ∼ Q_∞^-1(0) + c/L
    ,

for panels (c,d).
The latter two plots used Eq. <ref> rather than Eq. <ref>
because it gave much faster convergence.
Several coupling values, β J, and several simulation sizes are shown,
and fall on a single curve when scaled by their correlation lengths.
The lines in (a,b) plot Eq. <ref> for second derivatives (times L)
of the Helmholtz free energy for the lattice gas as n is varied.
Its agreement with the extrapolated values
shows that the canonical ensemble Q^-1(0 | n,V)
converges to <ref> first, and
then converges to ∂βμ^ex/∂ρ as V increases.



  The error visible at β J = -1 (panel a) gives another important
message for estimating the KB coefficients.
Rather than fit to Eq. <ref>, we should have subtracted off the
known contribution to the direct correlation function, c_LR = -2β J cos(2π m),
in Eq. <ref> and fit the remainder.
Since the correlation length is very short for this weakly coupled limit,
the resulting estimate would give perfect agreement with the line for β J = -1
(not shown).

  The lower panels (c,d) show extrapolation to 1/L = 0 using the
assumption<cit.> that ^2ρ^2 Q(m | μ, L)
is independent of L in Eq. <ref> and that Q(1/L) ∼ Q(0).
Attempts to combine both types of extrapolation, in either order,
resulted in worse convergence (not shown).

  Comparing the ferromagnetic and anti-ferromagnetic cases shows that
the oscillating radial distribution function of the anti-ferromagnetic
system leads to better overall convergence.  This is associated with
better short-range behavior of g(r) (Fig. <ref>a).
The strong ensemble-dependence
of the ferromagnetic g(r) carries over to Fig. <ref>a, where it appears
as worse agreement at small inverse-distance, m.  It is because
of the slow convergence of canonical and grand-canonical chemical potentials
that Fig. <ref>b shows Q(0 | n,L) is 10% lower than Q(0 | μ, L),
even at very large lengths.
Luckily, this issue appears to be limited to the 1D case.















§  LENNARD JONES MODEL


  This section presents results of applying the formalism to a supercritical Lennard-Jones (LJ) gas.
To mirror the conditions of the lattice model, we fix the density to a near-critical value,
ρ^* σ^3 ≈ 0.36, and set β^-1 = k_B T = 1.5 ϵ.
This is above the critical temperature, T^* = 1.36 ϵ/k_B.
In the lattice model analogy, the coupling constant, |β J|, would be ϵ / k_B (T - T^*)
if the temperature varied while ϵ remained constant.
Fluctuations in this system are interesting because they report on
cluster formation processes.<cit.>
After summarizing target `thermodynamic limit' results computed in the grand-canonical ensemble
and using analytical formulas, we present the scaling behavior in the canonical ensemble.

  The LJ gas has become an indispensable simulation model, for which
analytical approximations many properties are available.
The most important for the present case
are the equation of state, which gives Q(0) via the compressibility, and
the radial distribution function, which gives Q(m) via Eq. <ref>.
We compare our simulation results to the MBWR equation of state,<cit.>
which was parameterized from molecular dynamics and Monte Carlo simulation data,
and includes the effects from tail corrections.  We compare the radial distribution
function with a numerical solution of the Percus-Yevick (PY) closure of the Ornstein-Zernike
equation, computed following the method of Fries and Patey.<cit.>

  Our grand-canonical Monte Carlo (GC) simulations were performed with Towhee,<cit.>
at an excess chemical potential of μ^ex = -1.868 ϵ in
a cubic box of length 10.357σ.  Pair interactions were truncated at 5σ
after which long-range corrections were added for the energy and pressure.
Simulations were run long enough
to collect 10,000 configurations for analysis.  The number of trial moves between
configurations was 40,000, since the potential energy showed an exponential
autocorrelation function with a decay time of 13,200 trial moves.  Moves were selected
randomly with 25% probability for particle insertion/deletion and 75% probability for
translation.


  Because our GC simulation gave a slightly different density
(0.356 σ^-3), the GC properties summarized in Table <ref>
were obtained at ρ = 0.36σ^-3 by differentiating
the partition function estimate,

    β P V    = Θ(βμ) ≈Θ(βμ_0) e^C(βμ - βμ_0)
    
    C(λ | μ, L)    ≡ln∑_n=0^∞  f_n(μ,L) e^λ n
    ,

where f_n is the frequency of observing n particles in the GC simulation.
Eq. <ref> is a maximum likelihood estimate appropriate for a large number of samples,
S, when the probability for observing a histogram, f, is given by Sanov's theorem.
We computed the estimation error using the saddle-point approximation
(taking the second derivative of the relative entropy) to yield the covariance,

    δ C(λ)δ C(λ') S = e^C(λ + λ') - C(λ) - C(λ') - 1
    .




  The MBWR equation of state is fairly accurate near the state point studied here,
but predicts Q^-1(0) that is too low by 0.06 σ^3 and ∂ Q^-1(0)/∂ (βμ)
that is too low by 0.7 σ^-3.  Although its parameterization<cit.> did not include the particular
state point studied here, it can be noted that the MBWR slightly over-estimates
β P / ρ in this region.  That quantity reaches a minimum near the present simulation
conditions, and it is the shape of the minimum which the derivatives in Table <ref> report on.
We estimated the finite-size correction for Q(m) in Eq. <ref>,
from GC, PY and MBWR methods and all three agree on Q(0)/2^2 Q(0)ρ^2 = -20,
although the error in the GC estimate is 50%.





  Canonical simulations were performed with Towhee for L=6.525σ
to L=16.44σ.  The smallest two simulations used pair interaction
cutoffs of 3.25σ and 4.2σ, respectively.  All larger cell sizes
used a cutoff of 5σ.  All results include long-range corrections for the
energy and pressure.  The largest two simulations were run with LAMMPS
using a timestep of 0.003 LJ units and collecting 10^4 samples
– each separated in time by 10^4 dynamics steps.
The Q(0) and L_c values reported in Table <ref>
were fit to the smallest 3 m-points following Eq. <ref>.
For the L=6.525σ simulation, only the smallest 2 points
were used in the fit.



  Figure <ref> shows a running integral of the radial distribution function
from several of the simulations.  The solid lines use the correction
proposed by Krüger et. al.,<cit.>

    Q_est = ρ + ρ^2 ∫_0^L/2 4π r^2 (g(r)-1) (1 - (rL/2)^3 )
    ,

where L is the length of one side of the simulation cell.
The dotted lines were not corrected (setting the 1 - 2r/L
term to 1 in Eq. <ref>).
The figure also shows that the PY closure yields a reasonable estimate of the RDF
up to 4σ, while the remainder is not possible to judge based
on our simulation sizes.  The long correlation length of this near-critical point
system makes for a challenging test case.



  Fig. <ref>a plots the structure factor, Q, using spectroscopy units, k = 2π m.
It should be compared with Figs. 3 and 4 of Ref. ,
where the relationship is ρ^2 h̃(k) = Q(k/2π) - ρ.
That reference shows the shape of Q(m) near m=0 was very flat, and
easy to extrapolate.  The present case is much more difficult because it is near
the critical point and a very large peak has appeared near m=0.

  Fig. <ref>b shows the inverse of Fig. <ref>a near the origin.
It is the analogue of Fig. <ref> for the lattice gas.
For clarity, only the first three points with smallest m from each simulation are shown.
These particular m-points correspond to averages over plane wave shaped density fluctuations
oriented along faces, edges, and corners of the cubic unit cell.
All five of the canonical simulation points shown at m=0 are extrapolations
from these points.  Their numerical values are in Table <ref>.

  Fig. <ref>b also shows two other important points.
First, Q(m) computed from the largest simulation size
is artificially noisy and scaled downward at high wavenumbers.
The scale artifact is due to the inaccuracy of B-spline smoothing
past K_max/2, where K_max = 128/L is set by
our numerical grid, which contained 128^3 points for every L.
A similar artifact also appears for the L=12.65 simulation just
past the right boundary of Fig. <ref>b.  The noise is
caused by the discrete nature of the 128^3 grid points, which appear
at irregularly spaced distances from the origin.  It could be eliminated
by resamping the data to uniform intervals in |m|.

  Second, the long-range portion of the LJ potential can be Fourier-transformed
according to Ref.  to give,

    c_LR(x = m/η) = 4ϵ/3[
    (1-2x^2) e^-x^2 + 2 x^3 √(π)(x)
    ]
    ,

where η is an inverse distance corresponding to a
division between short-range and long-range forces.
Recognizing that this term makes an additive contribution to Eq. <ref>
can accelerate the convergence of the fit in Eq. <ref>.
Indeed, we find that subtracting this term from Q^-1
makes the first few points much more linear.  The result can be visualized
from Fig. <ref>b as the distance between the red and black lines.
However, to keep our discussion simple, we have not included it in our estimates.

  Similar to the antiferromagnetic lattice gas in Fig. <ref>a,
the estimates of Q(0) are larger than the grand-canonical value,
and decrease with increasing cell size.  The former figure showed that
the extrapolated points, Q(0), quickly converged to second differences of
the Helmholtz free energy (Eq. <ref>).
We cannot test this numerically here, since the MD estimates of βμ^ex
at finite volumes have uncertainties on the order of the volume.

  Instead, we have tested the finite-size corrections to βμ^ex
in Eq. <ref>.  Using the value of Q^-1 and its derivative
from Table <ref> gives Δβμ^ex = 4.18σ^3 / V.
This is a good fit to the size-dependence of μ^ex seen in simulations.
In the limit as m → 0, the finite-size corrections to Q
and to βμ^ex are related by a density derivative.
Computing the finite-size correction to Q(0) from Eq. <ref>
gives lim_m→ 0Δ Q(m) = -20 / V, which agrees with
the density derivative of Eq. <ref>.
The simulation points do not allow a single estimation of the
slope near m=0 for the largest simulations,
but the slope of the c_0 values vs. 1/V appears to be +200 rather
than -20 for simulations at intermediate cell-sizes.



  Because finite size corrections are relatively more important for canonical than grand-canonical
simulations, we have plotted the finite size correction to the free energy
predicted by the MBWR equation of state in Fig. <ref>.
The result reassures us that the chemical potential calculated by
test particle insertion has negligible volume dependence in both gas and liquid
states.
Near the solidification line, test particle insertion should be expected to slightly overestimate the
excess chemical potential, while we get an underestimation for high-temperature
liquids near the critical point.



§  CONCLUSIONS













  Kirkwood-Buff integrals are notoriously difficult to calculate from canonical, closed
simulations.  We have shown that this problem
can be helpfully mapped onto the problem of estimating the direct correlation
function of Ornstein-Zernike theory in a canonical ensemble.
All earlier work used truncated real-space estimates, for which the process of extrapolating
to infinite truncation radius, R, is ill-posed.<cit.>
Our exact theory for the canonical ensemble uses the entire simulation volume at once,
is manifestly invariant to shifts of the radial distribution function by a constant, and replaces
the problem of extrapolating to infinite R with the well-posed problem
of extrapolating to zero m within a finite simulation volume.

  The density-dependence of the radial
distribution function becomes important when
estimating the thermodynamic limit.<cit.>
These higher-order derivatives are very difficult to estimate accurately
from simulation data, and it is simpler to test for finite size effects
by scaling up the system.<cit.>
This work showed that 

the principle difficulty in reaching
the thermodynamic limit is the size-dependence of the canonical ensemble itself.
We then provided a simple, general formula for estimating finite-size corrections
of the excess chemical potential from KB coefficients.
These corrections are small for the LJ fluid.

  Our data in Fig. <ref> show excellent agreement for the
entire correlation function at low wave-vectors in reciprocal space.
Higher wave-vectors require smaller grid-spacing to avoid numerical
artifacts, though.  Fits using Eq. <ref> provide an estimation of both
the KB coefficient and the correlation length which are more accurate than
radial-distribution function integrals and fit well into the context of integral theories
of solution.  We have shown that the theory is robust by
treating difficult cases with long-range correlations.  These can be identified
in simulation work from the shape of Q(m) near m=0.

  We have shown that the zero-frequency limit of the canonical KB
theory gives second derivatives of the Helmholtz free energy for simulation sizes
larger than about 5 times the correlation length.
Further, this theory has several nice properties.
In the low density limit, the direct correlation function can be predicted
from the long-range form of the interaction energy function.
This leads to even better extrapolation to zero frequency.
Further information on the applicability of this process is given in the supplementary material.
Alternatively, the direct correlation function can be computed from simulation data
and used to test common assumptions for OZ closure relations.<cit.>
A companion paper uses this process to compute the spatial
dielectric response of water.<cit.>
This use of Q to estimate pair interaction energies
was first proposed by Madden and Rice <cit.>.



§  SUPPLEMENTARY MATERIAL

  The supplementary material includes derivations and extended discussion of several formulas from
the main text.  It provides a definition of the finite-volume radial distribution, free from cutoffs,
and consistent with Eq. <ref>.
It also provides explicit expressions for the other derived quantities of Kirkwood-Buff
theory in terms of the matrix Q_αγ(0|μ).
It derives the ensemble correction to the free energy (Eq. <ref>).
Then it gives computable expressions for the correlation function of the finite-size
Ising model and an approximate expression for the correlation length of the Lennard-Jones fluid.
Both correlation lengths are plotted as a function of state point.
Those results provide additional intuitive insight on the long-range correction proposed above.



§  ACKNOWLEDGMENTS

  This work was supported by the USF Research Foundation and NSF MRI CHE-1531590.



10

kb
John G. Kirkwood and Frank P. Buff.
The statistical mechanical theory of solutions. i.
J. Chem. Phys., 19(6):774–777, 1951.

nmatu96
Nobuyuki Matubayasi and Ronald M. Levy.
Thermodynamics of the hydration shell. 2. excess volume and
  compressibility of a hydrophobic solute.
J. Phys. Chem., 100(7):2681–2688, 1996.

pgang13
Pritam Ganguly and Nico F. A. van der Vegt.
Convergence of sampling Kirkwood-Buff integrals of aqueous
  solutions with molecular dynamics simulations.
J. Chem. Theory Comput., 9(3):1347–1355, 2013.

mrove90
M. Rovere, D. W. Heermann, and K Binder.
The gas-liquid transition of the two-dimensional lennard-jones fluid.
J. Phys. Condens. Matter, 2:7009, 1990.

sschn11
Sondre K. Schnell, Xin Liu, Jean-Marc Simon, André Bardow, Dick Bedeaux,
  Thijs J. H. Vlugt, and Signe Kjelstrup.
Calculating thermodynamic properties from fluctuations at small
  scales.
J. Phys. Chem. B, 115(37):10911–10918, 2011.

rcort16
R. Cortes-Huerto, K. Kremer, and R. Potestio.
Communication: Kirkwood-Buff integrals in the thermodynamic limit
  from small sized molecular dynamics simulations.
J. Chem. Phys., 145:141103, 2016.

pkrug13
Peter Krüger, Sondre K. Schnell, Dick Bedeaux, Signe Kjelstrup, Thijs J. H.
  Vlugt, and Jean-Marc Simon.
Kirkwood-Buff integrals for finite volumes.
J. Phys. Chem. Lett., 4(2):235–238, 2013.

jsala96
J. J. Salacuse, A. R. Denton, and P. A. Egelstaff.
Finite-size effects in molecular dynamics simulations: Static
  structure factor and compressibility. i. theoretical method.
Phys. Rev. E, 53:2382–2390, 1996.

jsala96b
J. J. Salacuse, A. R. Denton, P. A. Egelstaff, M. Tau, and L. Reatto.
Finite-size effects in molecular dynamics simulations: Static
  structure factor and compressibility. ii. application to a model krypton
  fluid.
Phys. Rev. E, 53(3):2390–2401, 1996.

uessm95
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen.
A smooth particle mesh Ewald method.
J. Chem. Phys., 103:8577–8592, 1995.

droge17d
David M. Rogers.
EwaldCorrel.
GitHub, 2017.
https://github.com/frobnitzem/EwaldCorrel.

lax
Peter D. Lax.
The Poisson Summation Formula, chapter 30, page 348.
Wiley-Interscience, 2002.

pdt6
T. L. Beck, M. E. Paulaitis, and L. R. Pratt.
The Potential Distribution Theorem and Models of Molecular
  Solutions, chapter 6, pages 123–141.
Cambridge, New York, 2006.

pdt2
T. L. Beck, M. E. Paulaitis, and L. R. Pratt.
Statistical thermodynamic necessities, chapter 2, pages 23–31.
Cambridge, New York, 2006.

jlebo67
J. L. Lebowitz, J. K. Percus, and L. Verlet.
Ensemble dependence of fluctuations with application to machine
  computations.
Phys. Rev., 153(1):250–254, Jan 1967.

davis9
H. Ted Davis.
Statistical Mechanics of Phases, chapter 9, pages 442–448.
VCH Publishers, New York, 1996.

isiep92
J. I. Siepmann, I. R. McDonald, and D. Frenkel.
Finite-size corrections to the chemical potential.
J. Phys. Cond. Matter, 4(3):679–691, 1992.

sfrie17
Sacha Friedli and Yvan Velenik.
The Ising Model, chapter 3, pages 97–183.
Cambridge University Press, Cambridge, UK, 2017.

jrein00
J. Reinhard, W. Dieterich, P. Maass, and H. L. Frisch.
Density correlations in lattice gases in contact with a confining
  wall.
Phys. Rev. E, 61(1):422–428, 2000.

cborz87
C. Borzi, G. Ord, and J. K. Percus.
The direct correlation function of a one-dimensional Ising model.
J. Stat. Phys., 46(1–2):51–66, 1987.

cteje87
Carlos F Tejero.
One-dimensional inhomogeneous Ising model: A new approach.
J. Stat. Phys., 48(3):531–538, 1987.

davis
H. Ted Davis.
Statistical Mechanics of Phases, chapter 3, pages 128–130.
VCH Publishers, New York, 1996.

jvavr01
Juraj Vavro.
Exact solution for the lattice gas model in one dimension.
Phys. Rev. E, 63:057104, 2001.

hmart96
Hernan L. Martinez, R. Ravi, and Susan C. Tucker.
Characterization of solvent clusters in a supercritical
  Lennard-Jones fluid.
J. Chem. Phys., 104:1067, 1996.

jnico79
J. J. Nicolas, K. E. Gubbins, W. B. Streett, and D. J. Tildesley.
Equation of state for the Lennard-Jones fluid.
Mol. Phys., 37(5):1429–1454, 1979.

jjohn93
J. Karl Johnson, John A. Zollweg, and Keith E. Gubbin s.
The Lennard-Jones equation of state revisited.
Molecular Physics, 78(3):591–618, 1993.

pfrie86
P. H. Fries and G. N. Patey.
The solution of the Percus-Yevick approximation for fluids with
  angle-dependent pair interactions.
J. Chem. Phys., 85(12):7307–7311, 1986.

mmart13
Marcus G. Martin.
MCCCS Towhee: a tool for Monte Carlo molecular simulation.
Mol. Sim., 39(14–15):1212–1222, 2013.

jweek71
John D. Weeks, David Chandler, and Hans C. Andersen.
Role of repulsive forces in determining the equilibrium structure of
  simple liquids.
J. Chem. Phys., 54:5237, 1971.

droge18a
David M. Rogers.
Fluctuation theory of ionic solvation potentials.
2018.
submitted.

wmadd80
William G. Madden and Stuart A. Rice.
The mean spherical approximation and effective pair potentials in
  liquids.
J. Chem. Phys., 72:4208, 1980.













§  ENSEMBLE DEPENDENCE


  Explicit expressions for the radial distribution functions
help to show how ensemble dependence arises in the KB integral.
Eq. 1 of the main text can be written as an average over
canonical simulations,

    ρ̂_α(r) ρ̂_γ(0) | μ, L   = ∑_n P(n | μ, L)
    ρ̂_α(r) ρ̂_γ(0) | n, L
    ρ̂_α(r) ρ̂_γ(0) | n, L   = n_γ/V∑_j=1^n_α P(r_α,j = r | r_γ,1 = 0, n, L)

The conditional probability in the last line gives the probability that the j^th
molecule of type α has its center of mass at r given that
the first molecule of type γ is centered at the origin.
When α = γ, this has the effect of decomposing the sum over
j into j=1, where we know r_j is already
at the origin, and j  1, where j is distinct
from the molecule at the origin.
This gives a delta function at the origin plus n_α - 1 times
the pair distribution function for fixed n.  Explicitly,

    P(r = r_α,j | r_γ,1 = 0, n, L) = δ(r),    α = γ, j = 1 
    
    g_αγ(r | n, L)/V,    o.w.

and

    ρ̂_α(r) ρ̂_γ(0) | n, L
    =  n_γ/V( n_α - δ_αγ/V g_αγ(r | n, L)
    + δ_αγδ(r) )


  The identification of g(r | n, L) in Eq. <ref> agrees with its method of
computation from canonical ensemble simulation data by estimating the conditional probability,
P(r = r_α,j | r_γ,1 = 0, n, L).
Of course, when the two particles are independent, the conditional
probability is 1/V, and g(r | n, L) is 1.
However, even for real molecules with excluded volume, g(r | n, L) always
integrates to the volume, V.  This makes the integral of g(r | n, L)
uninformative for number fluctuations – since there are none in
a canonical simulation.
Explicitly, if we insert Eqns. <ref> and <ref> into Eq. 1
and integrate, then the g(r | n, L) term drops out and the
entire Kirkwood-Buff integral
is due to the ensemble correction.

  For completeness, we re-state the derived quantities of Kirkwood-Buff
theory in terms of the ν×ν matrix Q = [Q_αγ(0|μ,L)] and using
the vector of ν densities, ρ, below.

    κ/β   = 1 / ρ^T  Q^-1ρ
    V̅   = κ/β Q^-1ρ
    β V [ μ_αn_γ]_n,P,T   =  Q^-1 - β/κV̅V̅^T 
    [ Πρ_α]_μ_1,ρ_γ 1,T   =  Q_2^-1ρ_2

These formulas express the isothermal compressibility, vector
of partial molar volumes, matrix of constant-pressure chemical potential derivatives,
and vector of ν-1 osmotic pressure derivatives (respectively).
The derivative of the osmotic pressure requires a special notation, Q_2, for the
ν-1×ν-1 sub-matrix of Q discarding the first row and column, and a
similar notation, ρ_2 for the ν-1 component vector of solute densities.
Note that the entries are discarded from Q before inverting.
The equation for ∂μ/∂ρ at constant temperature and
pressure requires using the relation dn_1 V̅_1 + dn_2 V̅_2 = 0
for a constant pressure variation of n_1 at constant (arbitrary) volume,
which does not generalize beyond 2-components.



§  ENSEMBLE CORRECTION FOR THE EXCESS CHEMICAL POTENTIAL



  The exponent of the excess chemical potential
can be written in the grand-canonical ensemble as,

    e^-α = e^β(μ^id_α - μ_α)≡ Z
    ,

which we call Z for compactness.
The ensemble correction for Z requires alternating
derivatives,[15]

    Z(μ) = Z(n) + V/2 V∑_γ(βμ_γ)n_γ Z(μ)
    ,

evaluating the first at constant chemical potentials (μ_α)
and the second at constant densities (μ_α^id).
The volume is constant in both, and can be used to transform
n_γ / V = ρ_γ.
The first derivative is,

    ρ_γ Z(μ) = Z(μ) (
         αρ_γ
        - Q^-1_αγ(0) )
    .

The next step requires the other derivative,

    (βμ_γ) Z(μ) = Z(μ) (
      αρ_α Q_αγ(0)
      - δ_αγ)
    ,

which uses the fact that μ^id_α does not
depend on any densities other than ρ_α.
Using these in Eq. <ref> yields, (to order 1/V)

    Z(μ)/Z(n) = 1 + 1/2V(
      - 2 αρ_α
      + Q_αα(0) [
                (αρ_α)^2
                + ^2αρ_α^2]
      + Q^-1_αα(0)
      - ∑_γQ^-1_αγ(βμ_γ))
    .


Since ∂α / ∂ρ_α = 1/ρ_α,
the Q_αα(0) term cancels.
The final ensemble correction given in the main text is the first term
in the expansion of -ln Z(μ)/Z(n).
Since Q^-1 is naturally a function of ρ, while its derivative in
Eq. <ref> is taken at constant ρ,
it is useful to change the independent
variables from μ to ρ,

    Q^-1_αγ(βμ_γ)
      = ∑_δ Q_γδ(0) Q^-1_αγρ_δ
    .

Another alternate formula is,

    Q^-1_αγ(βμ_γ)
     = -∑_ij Q^-1_α iQ_ij(βμ_γ) Q^-1_jγ
    .




§  ISING MODEL DISTRIBUTION FUNCTION


  The lattice gas in the main text has potential energy function,

    β H(ρ̂| L) = -x ∑_j=1^L ρ̂_j-1ρ̂_j
    ,

where we have defined β J ≡ x and used
periodic boundary conditions, ρ̂_L = ρ̂_0.
To compute powers of the transfer matrix,

    T ≡[ e^ϕ + x   e^ϕ/2;   e^ϕ/2       1 ]
    ,

we write it as T = α I + r_x σ_x + r_z σ_z,
where σ_x and σ_z are 2× 2 Pauli
spin matrices and

    α≡ (e^ϕ + x + 1)/2,   
    r ≡ [ e^ϕ/2, 0, (e^ϕ + x-1)/2 ]^T
    .

We also define the scale parameter,

    y ≡α - |r|/α + |r|.


  The partition function can then be written as,

    Z(ϕ, L)    =  T^L  = (α + |r|)^L (1 + y^L) 
    T^m/(α + |r|)^m/2   = (1 + y^m) I + (1 - y^m) (r_x σ_x + r_z σ_z)/|r|
    , 
    and it is easy to find expectation values
    of σ̂_z(j),σ̂_z(0)   =  T^L σ_z /Z = ( 1 - y^L/1 + y^L) r_z/|r|
    σ̂_z(0) σ̂_z(j)   = T^L-jσ_z T^jσ_z/Z 
       = r_z^2/|r|^2
    + r_x^2/|r|^2 y^j + y^L-j/1 + y^L
    .

These can be transformed into average densities and radial
distribution functions by substituting,
σ̂_z(j) = 2 ρ̂(j) - I.



  Fig. <ref> shows the prefactor of the cosine-term in
Q^-1(m), 2y/(1-y^2).
Since the prefactor dictates the correlation length estimated
by Eq. 16 of the main text, it is important to compare
to the analytical L_c in this model.
In the high temperature limit, |x| → 0, and
the prefactor simplifies to 2x ρ (1-ρ).
At large coupling, the prefactor approaches
the analytical correlation length, L_c.
Both limits are shown in the figure.

  The finite-size effect on L_c was calculated by estimating
the prefactor for canonical, closed systems.
This estimation is very nearly exact, since Q^-1 for closed
systems maintains its cosine shape, even at high coupling.
Those estimates are shown by crosses in the figure.
All these estimates show a larger effective correlation
at finite sizes.  Our simulation lengths were spaced
logarithmically in L, so the figure suggests
that the prefactor scales linearly with L/L_c.



§  CORRELATION LENGTH OF THE LENNARD-JONES SYSTEM


  The main text suggests estimating Q(m) near m=0
by subtracting a long-range estimate for
the direct correlation function,

    -c̃_LR(m) = ∫ e^-2π i m· r u_LR(r) dr
    .

In a companion paper,[30]
we show that this estimate is near quantitative accuracy for
ionic systems away from any phase transitions.




  To apply this ansatz here, we modify the WCA approximation
by using a fictitious reference system.
WCA showed that the hard-sphere reference system gets
the high-wavenumber behavior of c̃(m) right.
Since c̃_LR(m) is small at high wavenumbers,
we simply add it to the hard-sphere solution in the
mean spherical approximation (MSA),

    c_A(r) = c_d(r; ρσ^3) - βϵ u_LR(r)
    .


  In variance with the main text, we let u_LR
be the long-range part of the WCA potential,

    u_LR(r) = 
    -1,     r < 2^1/6
    
    4(r^-12 - r^-6),     r ≥ 2^1/6 ,

and r is in units of σ.
The function c_d is the cubic polynomial on [0,d)
given by the MSA solution
of the hard sphere system for a diameter chosen
by a consistency condition on Q(0).








We then multiply the solution by e^-βϵ u_0 as described
in Ref. [29]
to smooth the first peak around r = d.
The solution has better agreement with c̃(m)
at small wavenumber, but still
predicts a first peak that is lower than simulation
(and lower than the PY solution used in the main text).

  However, it is a good approximation for predicting the correlation
length.  It can be found from Eq. 20 in the main text, using

    1 - ρ c_0    ≡ρ Q^-1(0) = β/ρκ = 1 - 4πρ∫_0^∞ r^2 c(r) dr 
       = (2η + 1)^2/(1 - η)^4 - 329πρβϵ√(2)
    ρ c_2    = 4πρ/6∫_0^∞ r^4 c(r) dr
    = 4835πρβϵ 2^5/6 - d^2/(1-η)^4(  
          (2η + 1)^2 (η20 - 316) + 3(η + 18)/2
    ) 
    .

The second part in Eq. <ref> is the zeroth-order perturbation, ∫ u_LR(r) dr,
which matches the corresponding term
in the WCA perturbation energy.[29]



  Fig. <ref> plots Eq. <ref>-<ref>
for the correlation length along a few isotherms.
The numerical value at β^-1 = 1.5ϵ
and ρ = 0.36σ^-3 is L_c = 2 σ,
which is higher than simulation, but smaller than the PY
approximation.
The correlation length shows the basic features of the phase
diagram, which is remarkable because the hard-sphere
system does not have a liquid-vapor transition.

  Fig. <ref> is based on an approximation
to the WCA model, so its phase behavior is only qualitative.
The qualitative agreement, however, allows us to conclude
that Eq. <ref> is a good approximation to
the small-wavenumber behavior of c̃(m)
and its inclusion will aid extraction of c̃(m)
from simulation.

  It also shows that L_c is only few molecular diameters
away from phase transitions.  This short-ranged behavior translates
to a small second derivative of Q, which makes estimation of Q(0)
possible from the method in the main text
without corrections under ordinary circumstances.

