Fluctuation Theory of Ionic Solvation Potentials
David M. Rogers
cond-mat.stat-mech cond-mat.soft








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
















  This work presents a rigorous statistical mechanical theory of solvation free energies,
specifically useful for describing the long-range nature of ions in an electrolyte solution.
The theory avoids common issues with field theories by writing
the excess chemical potential directly as a maximum-entropy
variational problem in the space of solvent 1-particle density functions.
The theory was developed to provide a simple physical picture of
the relationship between the solution's spatial dielectric function,
ion screening, and the chemical potential.
The key idea is to view the direct correlation function of molecular
Ornstein-Zernike theory as a Green's function for both longitudinal
and transverse electrostatic dipole relaxation of the solvent.
Molecular simulation data is used to calculate these direct correlation functions,
and suggests that the most important solvation effects can be captured with only a
screened random phase approximation.
Using that approximation predicts both the Born solvation free energy
and a Debye-Hückel law in close agreement with
the mean spherical approximation result.
These limiting cases establish the simplicity and generality
of the theory, and serve as a guide to replacing local dielectric and
Poisson-Boltzmann approximations.


 Fluctuation Theory of Ionic Solvation Potentials
     David M. Rogers
    December 30, 2023
=================================================






§  INTRODUCTION


  The importance and many applications of implicit solvent methods hardly requires an introduction.
A complete equilibrium answer to this problem would be given by a computable expression
for the free energy of transferring a set of molecules with a fixed configuration into solution.
However, it is more common (and intuitive) to present the solvation free energy in terms
of a coarse-grained solvent density.  This recognizes both that the solvent has special, `relevant,'
degrees of freedom and that they respond to the solute in a self-consistent way.

  Most heavily used solvation methods focus on electrostatics, carrying
over the Maxwell theory to atomistic volumes by assuming the solvent
dipole density responds to a local electric
field with a slope that is essentially the solvent dielectric coefficient,
ϵ_r ∼ 80 in water.<cit.>
For ionic solutions, the theory is extended with a mean-field assumption for the ion
densities that leads to the Poisson-Boltzmann equation and the Debye-Hückel
(DH) limiting law for ionic solvation free energies.  If these theories are used as a starting point,
the goal is to find transferrable models for predicting the spatial behavior of ϵ_r(k),
and corrections to ion densities, from which all other quantities are derived.
However, comparing theory and experiment within this paradigm requires
contorting direct measurements from both sides to find “the dielectric,” which can
depend sensitively on boundary conditions.

  Experimentally, the procedure for measuring dielectric response is well-defined
on a macroscopic scale,<cit.> and corresponding molecular dynamics
calculations of spatial dielectrics, ϵ_r(k),
have appeared for many bulk liquids.<cit.>
The latter show that the manipulations required to transform
dipole-dipole or current-current correlation functions, Σ(k), into ϵ_r(k)
have the form Σ / (β^-1 - Σ/ϵ_0) (where β = 1/k_B T
is the thermal energy scale), and can exacerbate numerical errors
– even though solvation energies scale as βΣ.
Calculations for geometries including surfaces have also appeared,<cit.>
as well as methods aimed at computing ϵ_r in real space as a function of
distance from a solute.<cit.>

  These latter calculations highlight the issues faced by turning the Maxwell theory,
meant for describing macroscopic scales, into an atomistically
detailed calculation device.  First and foremost, the picture of electromagnetic
wave reflection at a macroscopic dielectric interface does not scale down to an
atomistic theory.<cit.>  It is well-known that solvent dipoles are
oscillatory at an interface due to orientational saturation,<cit.> and that water has
heightened in-plane correlations near hydrophobic solutes.<cit.>
These boundary layer effects motivate treatment of interfacial water
as a chemically distinct species.<cit.>
Describing the electric field experienced by molecular solutes can be
done within the Maxwell theory only by altering to the usual macroscopic
boundary conditions to account for these effects.<cit.>

  This paper explores a theory of solvent density response that is distinct
from the Maxwell dielectric picture.  It eliminates the solvent dielectric in
favor of a direct prediction of the local force (including electric and molecular fields)
to which solvent and solute dipoles respond.  The approach mirrors Ref. 
and the results are closely tied to integral theories of solution.<cit.>
We deviate from those theories by only
describing solvent response, with a focus on electrostatics.
The major results are summarized in section <ref> with a sketch
presenting solvation free energy components for an exactly solvable model.
In sec. <ref>, we present the general theory by
expressing the free energy in terms of a density functional,
explaining our method for calculating correlation functions of real fluids (<ref>)
and then demonstrating their use in determining the potential distribution (<ref>).
Section <ref> shows two results.  First, the effective local field (molecular potential
or direct correlation function) is computed for electrolyte solutions of varying concentrations.
These show that the contribution of nearby charges is screened out almost
exactly following an error-function like splitting.
Second, sec. <ref> shows that, with the error-function screening, Born and Debye-Hückle
solvation laws are recovered even using a random phase-like approximation.
The discussion in sec. <ref> compares these results with recent literature,
and we conclude by listing some of the unexplored consequences of the present development.

  Our final results support the molecular Ornstein-Zernike 
perspective that solvent position and orientational distributions should be described by
effective energy functions in close connection to pair correlations.



This is the same conclusion as an earlier nonlocal response function theory
developed for describing electron transfer.<cit.>




The present work was developed independently,
and is more simply motivated by the structure of the inverse pair correlations.
Our main result constitutes a novel and rigorous foundation for using density functional
methods to replace dielectric solvation theories for calculating free energies.
The focus on the dielectric is appropriate because, even using the mean spherical approximation
(MSA) solution for the primitive model of electrolytes,
the state-dependence of the solvent dielectric constant
is a major difficulty when comparing to experiment.<cit.>
Recent progress in integral theories has been made in predicting nontrivial
density-dependence for solvent dielectric response.<cit.>
Our approach of analyzing correlation functions has close connections
to both dipolar fluctuation theories<cit.>
and inference on bridge functions for RISM models.<cit.>
However, it differs from both in its focus on predicting solvent
density distributions directly through a simple maximum entropy procedure, as
opposed to parameterizing any particular theory.



§  LATTICE MODEL


  We will sketch our main results by recalling the simple, exactly solvable
system made up of a lattice of polarizable point dipoles, p_i,
at fixed positions, r_i, in the presence of fixed ions.<cit.>
Our purpose is to ground the discussion by demonstrating that this system contains
all electrostatic contributions to the solvation free energy –
dielectric self-energy, screening and solvent dispersion energies.
The result also reveals the root cause of some issues with “local dielectric” models.
Formulas which apply only to this system will have the superscript “ref.”
As a matter of convenience, statistical mechanical averages are written
using double-angle brackets, as in p for the average of the dipole vector, p.
Single angle brackets are reserved for a bra-ket notation for matrix inner products
(defined in the Appendix <ref>).

  The long-range part of the electrostatic energy of the system can be written as,<cit.>

    E_int^ref = 1/2∑_ij (q_i + p_i ·∇_i) (q_j - p_j·∇_i) G(r_i - r_j) 
    
    - η/4πϵ_0 √(π)( 2η^2/3∑_i |p_i^2| + ∑_i q_i^2 )
    
    .

In order to include the self (i=j) term, we only consider the long-range part of the
Coulomb potential, so that

    G(r) = ∑_n (η |r + n|)/4πϵ_0 |r + n|
    ,

where the sum runs over 3D lattice vectors, n, of the system's unit cell.
For a uniform distribution of dipoles, the screening effectively
ignores solvent near each molecule – creating a molecular cavity around each one.
If we assume the dipoles can take any magnitude with internal energy given by,

    E^ref_0 = 1/2∑_i |p_i|^2 / α_i
    ,

then the potential energy function for the collection of
dipole moments, p, is equivalent to that of a forced Harmonic oscillator.

  Since this system is Gaussian, the dipole fluctuations (Δ p_i ≡ p_i - p_i)
exactly satisfy the matrix equality,

    [Σ^-1]_ij/β = [G_p]_ij + δ_ij/α_i
      ,

for any position of the ions.
Here Σ_ij = Δ p_i Δ p_j^T is the
dipole-dipole correlation function, ρ_w is the density of dipoles, and
[G_p]_i,j≡ -(1-δ_ij)∇_i ∇_i^T G(|r_i-r_j|) is
the interaction energy between dipoles at r_i and r_j
(The Kronecker delta function omits the i=j terms).
Eq. <ref> gives a trivial derivation of the formula for the spatial dielectric response of a medium.
Typically, one would infer the dielectric function, ϵ≡ϵ_0 ϵ_r,
by comparing the linear response predicted by
statistical mechanics, Δ p = βΣE⃗^ext,
to the phenomenological equation expressing polarization
by a local electric field, ρ_w Δ p = ϵ (E⃗^ext - G_pΔ p).
Then the comparison reads,

    ρ_w Δ p = ϵ (E⃗^ext - G_p Δ p) = βρ_w ΣE⃗^ext

and is solved when ρ_w ϵ^-1 = (βΣ)^-1 - G_p.
Substituting the exact correlation function from Eq. <ref> shows that the phenomenological equation
gives the diagonal matrix, ϵ_ij^-1 = δ_ij/ρ_w α_i.




If, in addition to electrostatic interactions, there are other local interactions between dipoles,
then (βΣ)^-1 - G_p will still be localized.  However, its inverse will not.
This non-locality is the central problem with adapting the Maxwell theory of dielectric response
to atomistic scales.

  This simple, extensively studied,<cit.>
picture of a lattice of dipoles provides a new perspective on the role of dielectric
in integral theories of electrolyte solutions.
The direct correlation functions between water dipoles
provide the effective interaction energy between dipoles at every separation in a fluid,
and are more transferrable than their inverse “dielectric.”



The solvent response part of dielectric theory is replaced with a rigorous
maximum entropy theory for predicting the solvation free energy.
The density functional to be maximized is the natural target for
developing approximate physical theories.
If the interactions are restricted to long-range,
avoiding large high-wavenumber perturbations,
then a Gaussian approximation is shown to be very accurate.

  To finish the simple example, note that the classical free energy for any configuration of ions,
A^ref(β, r_α), can be found by integrating the partition function
over the vector of solvent dipoles,

    A^ref(β, r_α)    = 12βln |Σ^-1 / 2π| + E_eff^ref
    In the limit of a uniform distribution of N dipole locations and letting
    η→∞ while 1/α̅= 1/α - η^3/3ϵ_0π^3/2,ln |Σ^-1/2π|    = 3N ln (β / 2πα̅) 
          + (N-1) ln (ϵ_r)
      + ln|I_3 + α̅ρ_w J | 
    
    E_eff^ref   = 1/2ϵ_r∑_ij q_i q_j G^ref(r_i - r_j) 
         - η/4πϵ_0√(π)∑_i q_i^2 
    ϵ_r    = 1 + α̅ρ_w/ϵ_0 
    .

Here J is the depolarization tensor that expresses the surface energy,
E_surf = M^T J M/2V, due to a net system dipole, M,
surrounded by a bulk medium of dielectric ϵ_e.<cit.>
For a spherical boundary at infinity, J = I_3/ϵ_0 (2ϵ_e + 1),
with I_3 the 3× 3 identity matrix.

  Eq. <ref> contains electrostatic screening, dielectric self-energy and solvent
dispersion energies.
The screening appears directly in E_eff via ϵ_r > 1.
The dielectric self-energy is the free energy of solvation for a single ion.
It appears as the difference between the i=j term,
where G_ref is η/2ϵ_0 π^3/2,
and the right-hand side.
Identifying this with the Born solvation free energy lets us make the definition
R_B≡√(π)/2 η.
The solvent dispersion energy appears in the
normalization constant, |Σ^-1|.<cit.>
We have left the mass contribution out of  <ref>
and present only the classical free energy.
Adding ions does not remove polarizable centers in this picture, so there is no dispersion energy
of solvation unless the ion polarizability differs from bulk water.



  The result is exact for the reference system, but cannot be generalized to solution
density functions because integration over all possible densities is ill-defined.
Physically, the space of all possible density functions is much larger
than the configuration space of the system and introduces non-physical degrees of freedom.<cit.>
The Hubbard-Stratonvich transformation provides one alternative route,
but also meets difficulties with integration over infinite field variables.
Even when those integrations can be defined with a convergent limit,
the field functional cannot be given in a closed form and must be approximated
in a way directly paralleling density functional theories.<cit.>
It should be noted that this procedure has been carried out rather clearly
for fluids with soft pairwise interactions in Ref. ,
and both it and an earlier work<cit.>
contain parallels to many of the results from our Gaussian approximation
in Sec. <ref>.  Implementations of that theory have
also helpfully cast the solution process as a maximum entropy problem.<cit.>

  This work avoids issues with field variable integration by applying large deviation theory to
find a representative density that yields the exact distribution of solute interaction energies.
This approach appears to be a new alternative to the random phase approximation route
to electrolytes solution structure recently shown by Frydel and co-workers.<cit.>
Comparing the results of this procedure with the exact Eq. <ref> is interesting
because it helps eliminate misconceptions and provides an intuitive
context for all of the quantities that will appear below.







§  GENERAL MODEL




  The molecular picture of a uniform dipolar fluid above works well because
local interactions were screened out by the replacement 1/r →(η r)/r.
Locally, the system behaves like a fluid with only short-range order while we focus on treating
long-range correlations correctly.  This leads us to the thermodynamic cycle of Fig. <ref>.
Solvation of a charged species is divided into a long-range (LR) and a local (SR)
contribution.  The long-range step is performed first.
According to the potential distribution theorem, the excess chemical potential of
a solvated molecule, α, is then expressed exactly as,

    βμ^ex_α   = -ln e^-βΔÛ^LR_α_0 - ln e^-βΔÛ^SR_α_LR
    ΔÛ^LR_α   ≡⟨ρ̂_q(r') |(η_α |r' - r|)/4πϵ_0 |r' - r|⟩ρ̂_α(r) 
    ΔÛ^SR_α   ≡ΔÛ_α - ΔÛ^LR_α .

The first term on the right of Eq. <ref> is defined as βμ^ex,LR_α,
the focus of this work.
Here, ρ̂_q(r) is the instantaneous charge density at point r
corresponding to a single solvent
microstate, while ρ_q(r) is the average charge density which is zero by symmetry.
Eq. <ref> for the long-range interaction energy uses a bra-ket notation
for the double-integral over dr and dr'.
Hats are used used on quantities that depend on the molecular coordinates.
Although this long-ranged interaction just contains the screened Coulomb energy,
we note that it is simple to generalize η_α to η_α(r), so that
different atoms can have different effective radii.  The total solvation free energy,
μ^ex_α, will not be affected by the choice of ΔÛ_α^LR,
but our results confirm that choosing √(π)/2η = R_B
minimizes the importance of the second, short-range step.
The theory in Sec. <ref> is general enough to handle any
choice for ΔÛ^LR_α, including the full
ΔÛ_α itself.



 §.§  Correlation Functions


  We have recently presented a simple Fourier-space method
for inferring the direct correlation function from molecular dynamics data.<cit.>
This method of inverting the dipole correlation function to find an “effective
dipole-dipole interaction energy” is synonymous with computing the direct correlation
function of  Ornstein-Zernike theory and the polarization structure factor of nonlocal
response function theory.<cit.>
Formally, define a vector, F_α, characterizing the orientation of a molecule of type α.
For example, the vector characterizing
a molecule containing a point dipole could be F = [1, p_x, p_y, p_z]^T.
Correlations between the density operators,

    F̂_α(r) = ∑_j=1^N_α F_α,jδ(r - r_j)
    ,

where j indexes the N_α molecules of type α,
then report on both scalar and vector interactions.




  The correlation function for these vectors is organized into a matrix where
molecule indices are combined with vector indices and Fourier transformed,

    Q_αγ(k) ≡∫_V e^-i r· kΔF̂_α(r) ΔF̂_γ(0)^T  dr
    .

The integral ranges over the unit cell with volume V.
To minimize error in the low-frequency components and to compute Q efficiently,
Eq. <ref> is estimated from simulation data by
averaging squared Fourier transforms (ψ̂_α(k) ≡ℱ[F̂_α(r)](k)),
so that Q = Δψ̂Δψ̂^†/V.
Also, we define a potential after dividing by the volume
so ϕ_α(k) ≡ℱ[ Φ_α(r)](k) / V.
The Appendix gives helpful relations for the Fourier transform used here.





  The correlation function is related to the Ornstein-Zernike definition,<cit.>

    Q_αγ(k)^-1 = βμ^id_αρ_γ - c_αγ(k)
    .

Here,

    βμ^id_α = ln (ρ_αΛ^3_α / q^int_α)

is the ideal gas expression for the chemical potential of species α.<cit.>
Approximations to Q are available both from analysis of experimental data<cit.>,
and from solutions to integral equations for charged or dipolar hard-sphere liquids.<cit.> 
The results section presents the correlation functions and inverses computed
from MD simulations of a series of 1:1 electrolyte solutions.





 §.§  Solvation Potential Distribution


  The basic quantity of the present theory is an exponential average,

    Z[Φ ; V,β] = expρ̂Φ
    ,

which is taken over the grand-canonical ensemble with μ
describing the system before coupling to the solute (when Φ = 0).
A constant shift of Φ has the same effect as changing the
chemical potential in Eq. <ref>.
The inner product in the exponent is defined in Appendix <ref> and carries units of volume
if both sides are in real space and inverse volume if both sides are in Fourier space.
The long-range excess chemical potential can be
written in terms of Eq. <ref> as
βμ^ex,LR_α = -ln Z[-βΦ_α ; n,V,β],
which requires Legendre transformation of ln Z into the nVT ensemble,
and defines Φ_α to be the pair interaction field produced by the solute
(by comparison to Eq. <ref>) as,

    ΔÛ^LR_α = Φ_αρ̂ = ϕ_αψ̂
    .






  Equation <ref> is a moment generating function for the density field, ρ̂,
and we can define its corresponding density functional as the Legendre transform,

    ℛ[ρ ; V,β] ≡inf_Φ[ ln Z[Φ ; V,β] - ρΦ]
    .

Eq. <ref> is a generalized entropy, and is also called a (negative) rate function
for the empirical distribution, ρ̂_V,β,Λ
in large deviation theory.<cit.>
It is also the negative of the traditional “density functional.”<cit.>
The functional ℛ is maximized under extra constraints Λ
at constant volume and temperature when ρ = ρ̂_V,β,Λ.
Stated in terms of Fourier-transformed densities and potentials,
ℛ satisfies

    ℛ[ψ]ψ^†_α(k)   = -ϕ_α(k;Λ), 
    
    V ϕ_α(k; Λ)ψ_γ(k')    = δ_k,k' Q^-1_αγ(k;Λ)
    .

The second equality requires translational invariance of the constrained ensemble (V,T,Λ)
and explicitly acknowledges that the fluctuations in Eq. <ref> depend on that ensemble.

  Because Z and ℛ form a Legendre transform pair, βμ^ex,LR
is simple to state in terms of ℛ,

    -βμ^ex,LR_α = sup_ψ[ ℛ[ψ] - ψβϕ_α]
    .

This equation may prove to be the most useful result in this paper.

It shows that an approximation for Eq. <ref> can be used to create a density
functional that avoids the “charging”
integration process of density-functional theory.<cit.>
The maximization over densities in Eq. <ref> can be carried out under
fixed n = ψ(0), but might also include other conditions.
Those other constraints can be used to change the ensemble or
even to fix the density to zero near the origin (as for the mean-spherical approximation).
Because ℛ was defined as a minimization problem,
constraints on ρ will be reflected physically by deviations of
the pair potential, Φ, which solves ℛ away from the
default (ρ-unconstrained) solution at Φ_α.
Rigorously, constrained maximization is justified by the
Gibbs conditioning principle, which states more formal conditions on
the constraint and density spaces.<cit.>




  This development makes an important departure from traditional Ornstein-Zernike theory.
Rather than seeking to provide an extra closure between direct and indirect correlation
functions, our primary target is to model the rate function (Eq. <ref>).
The Ornstein-Zernike relation is embedded in the structure of Eq. <ref>,
so that the density fluctuations (radial distribution functions) come out as a consequence
of a proper model for ℛ.







  So far, the development from Eq. <ref> through Eq. <ref> has been formal and exact.
To present analytical results in Sec. <ref>, we use the two-moment approximation,

    ℛ[ψ] ≃ - 1/2⟨Δψ|Q^-1⟩Δψ
    ,

where Δψ≡ψ - ψ_0 and ψ_0
is the (Fourier) density of the uncoupled system.
Technically, the domain of ψ is limited to positive densities,
which could also be included either as a better approximation to ℛ[ψ]
or as additional constraints.
The results in this work do not include such a constraint.

  Earlier work showed the utility of a Gaussian perturbation theory for the long-range part
of the solvation free energy after forming a sufficiently large
cavity.<cit.>

    μ^ex,LR_α = ΔÛ^LR_α_0
      -  βσ_LR^2 / 2 + O(β^2)
     .

Before forming a cavity, the first term in the expansion is zero by symmetry.
We recover this result by inserting Eq. <ref> into Eq. <ref>, adding ϕ(0)^T ψ(0)
at k=0 to constrain the number of solvent molecules (with undetermined
Lagrange multiplier ϕ(0)), and maximizing.
This leads to the identifications,

    Δψ(k)    = -β Q(k) ϕ_α(k) + δ_k,0 Q(0) ϕ(0)
    , 
    and
    -2 μ^ex,LR_α   = βσ_LR^2
      = δ_k,0ϕ(0) β  - ϕ_αΔψ
       = - ϕ_αΔψ

and shows that when Eq. <ref> applies, the linear response
approximation actually gives the ensemble average of the 1-particle density functions
in the coupled state, and proves that the distribution of interaction energies,
ΔÛ^LR_α, is also Gaussian.
The last simplification step in Eq. <ref> happens because
each component of the multiplier, ϕ_j(0), is chosen so that the corresponding
molecule number does not change (Δψ_j(0) = Δ n_j = 0).

  For describing ionic solvation, specialize the solution components
to a 1:1 electrolyte (e.g. NaCl) in water (modeled as a dipole so F_H_2O is a 3-component
vector).  Each density vector, ψ(k), then has 5 components – cation and anion
densities, ψ_q, plus water dipole vectors, ψ_p.
In this basis, the screened interaction energy of an additional ion of charge z_α
with solution (Eq. <ref> and Eq. <ref>)
has the form βϕ_α(k) = z_α f(k) [1,-1,-ik]^T
where f(k) ≡β e^-k^2 / 4η^2 / k^2 ϵ_0 is the Fourier-transform
of the screened Coulomb operator (Eq. <ref>).



  Equation <ref> is the linear response of both water dipoles and ion densities
to charged solutes.  It is given, in the Gaussian approximation,
by the convolution of the indirect correlation functions, Q,
with the potential, ϕ_α.  Since the potential is smooth and long-ranged,
the small k behavior of the correlation functions
is the most important aspect of solvation.  Second most is improving ℛ
beyond the Gaussian approximation to describe excluded
volume effects.



§  RESULTS


  We show two key results.  First, for sodium chloride solutions, the direct correlation
functions exhibit error-function like screened forms with little dependence on ionic strength.
Second, the long-range limits of these screened Coulomb forms motivate a random-phase
like approximation whose analytical predictions for solution dielectric and ionic response re-derive
the Born law and extend the Debye-Hückel theory with an ionic size-dependence that
closely mirrors the full mean-spherical approximation result.
Together, these results justify abandoning the Maxwell theory and its associated extensions of the
Poisson equation.
Instead, direct charge correlation functions, available from both experiment and simulations,
form the basis for computable methods of solution response and solvation free energies.



 §.§  Simulation Correlation Functions


  We simulated sodium chloride solutions in SPC/E using the Kirkwood-Buff forcefield model
for ions<cit.> replicating the simulation
conditions in Ref. <cit.> but extending all simulations to at least 10 ns.
Figures <ref> and <ref>
show the water position and dipole indirect correlation functions (F = [1,p_x,p_y,p_z])
for pure SPC/E<cit.> and TIP5P<cit.> water models under NVT conditions at 300 K.
These were computed by inverse Fourier-transform of the correlations in reciprocal space.
Deviations of h(r) below the exact value of -1 near the origin indicate the amount of numerical error
due to the band-limiting inherent in our method.

  For all vector-scalar interactions plotted in real-space, the geometry should be pictured as
with the vector fixed at the origin, pointing along the +z axis.
Positive interaction energies or low densities then describe unfavorable interactions
of the dipole with molecules on the +z hemisphere of radius r.
Vector-vector interactions (e.g. between dipoles p_1 and p_2)
along a separation direction r̂
in real space are broken into parallel
and isotropic components scaling with
p_1^T r̂r̂^T p_2 or p_1^T p_2, respectively.



  Figure <ref> plots the components of water's the inverse correlation function, Q^-1(k).
Panels b and d should be compared with Figs. 7 and  15 of Ref. ,
after noting that each point, Q^-1(k), plotted here is the result of a 4× 4
matrix inverse.
The two symmetry components of the dipole-dipole direct
correlation function, χ≡ Q^-1_pp,
in Fig. <ref>b and <ref>d are,

    χ(k) = χ^0(|k|) k̂k̂^T + χ^1(|k|) (I_3 - k̂k̂^T)
    .

Since the field produced by a dipole at the origin is proportional to k̂k̂^T,
χ^0 represents the effective energy governing the
`longitudinal' fluctuations of the second dipole in the direction
parallel to the field lines created by the first water dipole, while χ^1 represents
`free' or `tranverse' fluctuations in the perpendicular direction.
Without screening, these would be χ^0ϵ_0/β = 1 + ϵ_0 / ρ_w α
and χ^1 ϵ_0/β = ϵ_0 / ρ_w α in the lattice model.
With screening, the 1 in χ^0 is
replaced by the Fourier-transform of the screening charge distribution
while χ^1 is unchanged.




This is the usual justification for using χ^0(0) = χ^1(0) = V/M_x^2 (when J = 0)
to compute ϵ_r(0) = β / χ(0) ϵ_0 + 1.
However, away from k=0, only χ^1 varies continuously.
Fig. <ref>b and <ref>d clearly show that χ^0
starts at 1, but decreases exponentially as k increases.
The fitted line uses 1/η = 1.8Å for SPC/E, which is equivalent to R_B = 1.6Å.
Interestingly, the χ^1 component starts near 1/79, but increases to a larger constant
value around 0.0540 (SPC/E) or 0.0556 (TIP5P).  Using this in Eq. <ref>
gives only 19.5 or 19.0 for the dielectric, respectively.
This behavior also appears in the MSA solution, and is more pronounced
for dipolar-only solvents.<cit.>



  Figures <ref>-<ref> show the ion-ion, ion-water, and water-water
radial distribution functions and their transformations as a function of salt concentration.
The direct interaction functions in real-space (Figs. <ref>, <ref>, <ref>, and <ref>)
follow screened Coulomb electrostatic laws so closely that we have subtracted out the screened forms
and plotted the differences instead.
Subtracted parts are shown so that the two curves can be visually added to get
the total Q^-1.
Ion-water (Figs. <ref> and <ref>) and ion-ion interactions (Fig. <ref>) show negligible concentration dependence.



  Figure <ref> shows the ion-ion interactions in reciprocal
space to have the expected, ideal, 1/ρ
contribution and the k^-2 divergence near the origin.<cit.>
The k → 0 limit in Fig. <ref>a and <ref>b can be used to find
the Kirkwood-Buff coefficients for single-ions.<cit.>
A log-log plot (Fig. <ref>c,d) shows the k^-2
divergence transitions to an exponential-like screening around k=1Å^-1 before
flattening out due to the ideal (1/ρ) contribution.  These profiles were used to identify
appropriate screening lengths for each ion, resulting in
1/η_++ = 2.2Å and 1/η_– = 2.65Å. 
Note that Born radii are expected to be temperature and pressure-dependent.
However, our results indicate that they are not sensitive to
salt concentration in the 0-4M range.



  The ion center of mass direct correlation functions in Fig. <ref>
are surprisingly featureless apart from the 1/r interaction.  An unexplainable, small linear trend
shows up in Fig. <ref>b at long range with a positive slope for like-charged
ion interaction and a negative slope for cation-anion interactions.
At short-range, the like-charged interaction shows numerical noise before a steep,
repulsive increase.  This does not appear in the cation-anion interaction.
This short-range behavior is not directly linked to interaction energies.
In the mean-spherical approximation, it is treated as a fitting parameter to
achieve excluded volume in g(r).





  Fig. <ref> returns our focus to the number and dipole
density response behavior of water.  The dipole-dipole interaction energies
for co-linear aligned and in-plane parallel orientations (Fig. <ref>c,d, respectively)
both follow the expected trends.  Within this dipolar response function description,
local effects only appear within the first 6Å.  Interesting salt-concentration
effects show up in the the shape of the water solvation
shell around a water fixed 
at 
the origin (Fig. <ref>b).  This may be linked to electrostriction
and specific-ion effects.<cit.>



 §.§  Approximate long-range limiting laws


  We show the utility of the theory above by deriving simple
analytical expressions for long-range
ionic chemical potentials (Eq. <ref>) using a “screened MSA”
(for a 1:1 electrolyte with equal cation and anion concentrations, ρ),

    Q^-1(k) = [  ρ^-1 + f        -f   i k^T f;       - f  ρ^-1 + f - i k^T f;    -i k f     i k f         χ ]
    ,

where χ = I_3 β/ρ_w α + f k k^T is the lattice result,
containing contributions from water polarizability
and the screened dipole-dipole interaction energy.
To properly express the energy contribution from polarizing the surface,
the dipole-dipole energy at k=0, should be replaced by the
depolarization tensor (χ(0) = I_3β/ρ_wα + β J).<cit.>
Eq. <ref> is essentially a mean-spherical approximation, and

is a caricature of the low-k behavior found in Sec. <ref> with
only one η value. 

  This ansatz is simple, supported by our simulation results 
and yields both the Born theory and a mean-spherical approximation-like
Debye-Hückle theory as limits.
However, it is not intended to be used in practice, since Ref. 
contains analytical expressions for Q that match our simulation results
for water much more closely.



 §.§  Born Limit


  The developments above make re-derivation of Eq. <ref> a simple
matter of removing charge-ion interactions.  To accomplish this, delete
the top two rows and columns from Q^-1 in Eq. <ref>
and in Eq. <ref>.  The result is,

    β^2 σ^2_LR,Born   = z_α^2 ⟨k| f^2(k) 
     χ^-1⟩k
       = z_α^2/V∑_k f(k) - f̅(k)
     ,

which introduces a definition for the scaled interaction energy,

    f̅(k) ≡ f / (1 + e^-k^2/4η^2αρ_w /ϵ_0) ≡ f / ϵ_r(k)
    .

The last equality also defines the “dielectric function,” ϵ_r, of the lattice model.
It reduces to Eq. <ref> when η→∞.
The k^-2 divergence at k=0 integrates
to a contribution that scales as 1/V.

  Because the solvent density response is purely along the longitudinal (k̂)
direction, it can be written in two forms,

    Δψ_p = i z_α f χ^-1 k = i z_α k/k^2(1 - 1/ϵ_r )
    .

  It is extremely important to notice the difference between 1/ϵ_r(k),
defined for mathematical convenience in Eq. <ref>,
and the indirect dipole-dipole correlation function, χ^-1,
which gives the both transverse and longitudinal contributions to the orientation
response.  Both these correlations are long-ranged, while their inverses appear
short-ranged.  However, ϵ_r is only short-ranged in real-space when there are
no other interactions besides electrostatic ones (hence the k^-1 prefactor),
while χ(r) is not short-ranged until the electrostatic 1/r^3 term is subtracted.



 §.§  Debye-Hückel Limit


  Maximizing only the dipoles, Δψ_p, in Eq. <ref>
leads to the linear response equation,

    Q^-1_pp | Δψ_p ⟩ + Q^-1_pq | Δψ_q ⟩ + | βϕ_p ⟩ = 0
    ,

with obvious notation for 2- and 3-component sub-blocks of Eq. <ref>.
Replacing this back in Eq. <ref> and simplifying
leads to a charge-only system with Q^-1_qq
scaled by the effective dielectric (Eq. <ref>),

    Q̅^-1_qq(k)    = 
    [ ρ^-1 + f̅       -f̅;      - f̅ ρ^-1 + f̅;           ]
    ,
    and  βϕ̅_q    = [f̅, -f̅]^T .


  Inverting Eq. <ref> yields the ionic charge density
induced by an external potential from a charge at the origin,

    Δψ_q(k) = - z_α(
      ρf̅/1 + 2ρf̅ - δ_k )
    [  1; -1 ]
    .

Ref.  also found Eq. <ref> by applying the
random phase approximation to a fluid of screened, soft-core ions.
They obtain a closed form for its Fourier transform (h(r))
and show that it predicts a cross-over from exponential to damped oscillatory
behavior at high ionic strength.

  Collecting all the terms in Eq. <ref> created by our procedure,
the fluctuations in the solvation potential
of an ion with charge z_α is given by the sum,

    σ^2_LR   = σ^2_LR,Born + σ^2_LR,DH
    β^2 σ^2_LR,DH   = z_α^2/V∑_k 02 ρf̅^2(k)/1 + 2ρf̅
     .

Using the relations for the inner product developed in Appendix <ref>,
the summation in Eq. <ref> goes over to the inverse Fourier transform
as the cell volume increases,

    β^2 σ^2_LR,DH   →β z_α^2/(2π)^3 ϵ_rϵ_0∫_0^+^∞4πκ^2 e^-k^2/2η^2/k^2 + κ^2 e^-k^2/4η^2  dk 
       = β z_α^2 κ/4 π^2 ϵ_rϵ_0∫_-∞^∞e^-k^2 u^2/k^2 e^k^2 u^2 + 1  dk
    .

Here κ^2 ≡ 2βρ/ϵ_r ϵ_0 and u^2 ≡κ^2/4η^2 = κ^2 R_B^2/π.  When u goes to zero (zero ion radius),
the integral becomes π, and we recover the classical Debye-Hückel result for the
solvation free energy using Eq. <ref> with Eq. <ref>.

  For finite η, ϵ_r is k-dependent, and
the integral is strictly smaller than π.  It decreases with increasing
κ R – making ionic solvation less favorable.
We calculated the integral numerically and verified that it compares well
with the mean spherical approximation result for the electrostatic free energy component with
ion radius R_B.<cit.>





§  DISCUSSION


  The assumption in Eq. <ref> used to derive Born and modified DH laws is essentially
the random phase approximation (RPA) of fluid density functional theory.
A similar derivation has recently been presented (starting from the RPA)
by Frydel and Ma.<cit.>  They do not include screening and thus arrive at the linearized
Poisson-Boltzmann equation.  Earlier work showed the utility of this approximation
for describing charge reversal near a charged surface.<cit.>
In later work, they caution that the RPA, being equivalent
to a variational Gaussian approximation, has the same weaknesses as a mean-field approximation
and loses applicability at strong coupling.<cit.>
The primary difference in this work is that we do not use the adiabatic connection to
charge up the whole system at once, but only calculate the effect of
adding a single ion in a self-consistent way.  This frees approximations like
Eq. <ref> from having to represent the density of the entire system.
The RPA works well because the solute potential, Φ_α,
is concentrated at long range, near k=0.

  We view the focus on addition of a single molecule to a pre-existing fluid as the
major reason for success of the present theory.  An entire fluid constructed from
soft ions has very different properties.  It has been investigated as the
“ultrasoft restricted primitive model” (URPM) of penetrable electrolytes
such as charged polymers.<cit.>
The screening prevents Coulomb collapse in the URPM, but allowing ion
overlap removes the competition between long-range attraction and short-range excluded
volume underpinning most density functional closures,
explaining their difficulty in predicting its phase diagram.








  It is hard to compare our simulation results with the body of literature
on the Hubbard-Stratonovich transformation, which is a Fourier instead of a Legendre
transform.
The major difficulty is technical.  Since the transformation introduces an imaginary,
auxiliary field with an infinite number of degrees of freedom, most of its expressions
are related to distributions of the imaginary potential and cannot be easily
compared to molecular simulations.  Analytically, the theory is commonly used in combination
with a mean-field type expression for the potential or density functional
which can result in various modified Poisson-Boltzmann theories depending
on this approximation.<cit.>
When employing the variational Gaussian approximation within the statistical
field theory, the results of Sec. <ref> appear almost exactly
with the same conclusions.<cit.>  Nevertheless, the picture there is of the URPM model
and is more difficult to extend to non-pairwise interactions.
In addition, entropic terms similar to the
primitive available volume approximation (-ln(1 - V_0 ρ))
have been successfully used to account for the free energy of cavity formation
in related works.<cit.>






  The physical picture of solvent response gained from Fig. <ref>
is that the “local electric field,” which determines the force felt by
a solvent dipole has special contributions from the first two or three solvation layers.
Beyond that, it has the expected form of an integral over dipole-dipole electrostatic
1/r^3 terms.  Two maxima surround a sharp minimum in χ^0
near the first peak of the water-water radial distribution (Fig. <ref>c)
– so that the dipole-dipole interaction is slightly more repulsive for
waters deviating from the optimal distance.
This behavior likely reflects quadrupolar forces causing correlations between the position and
orientation angle of first-shell waters.  The isotropic part of the local electric field
(that scales with the cosine of the water-water angle, μ·μ) in
Fig. <ref>d mainly follows a typical screened form,
but shows an extra tendency for
waters closer than 3Å to take antiparallel orientations.<cit.>

  For ion-water association, Fig. <ref> shows a single peak where
the usual ion-dipole interaction is strengthened, while Fig. <ref>
shows a local minimum for the center-to-center interaction.  This strengthening occurs
for both ions and both interactions right at the maximum of the ion-water
radial distribution function.  The extra dipolar response is direct evidence
for the hypothesis that the first solvent dipole layer at the solute/solvent interface
should be scaled from the Maxwell theory.<cit.>  The extra
center-to-center interaction provides an energetic basis for electrostriction.
Both are associated with specific-ion effects.<cit.>



§  CONCLUSION


  The limiting cases studied throughout Sec. <ref> establish the utility
of this theory, and serve as a guide for translating the language of
dielectric polarization into structural, energetic terms.
Future work should explore the key features of the correlation functions
responsible for experimental evidence of specific ion effects.

  We have taken one step, showing the
direct correlation function found from simulations has a surprisingly good agreement
with the random phase and mean spherical approximations.
They are indistinguishable past the second solvation
shell of water.  Even within the first and second solvation shells, the effective energy only oscillates
about the error-function screened form.  One way to rationalize the far-field result is to note
that long-range electrostatic ordering (and charge density fluctuations) does not depend sensitively
on local charge ordering.  Rigorously, higher-order multipoles have a rapidly decreasing interaction
radius.






  We reiterate here the importance of distinguishing between the ideal dielectric continuum
theory, for which ϵ_r can be usefully defined, and real solutions containing additional
local interactions besides electrostatics.  As Eq. <ref> and <ref> illustrate,
the inverse correlation function is more informative than an inferred dielectric function.

  This observation helps give a simple mathematical structure to a nonlocal response theory
for solution density response to the solute's molecular field.  It is well-known that a majority
of the solvation free energy in polar liquids is due to long-range solvent electrostatic response.
Most dielectric theories spend a lot of effort finding a solute-solvent
boundary that makes bulk dielectric polarization theories give accurate free energies.
The density response paradigm parallels density functional theories by focusing
on reproducing average solvent dipole (and ionic) densities conditional on interaction with the solute.
Eq. <ref> then provides the excess chemical potential exactly.
Using a screened electrostatic interaction removes non-physical issues
with singularities of the Coulomb potential at the origin, while capturing long-range
contributions.

  The same division into short and perturbative long-range interactions is a crucial step in
local molecular field (LMF) theory,<cit.>.
While LMF theory takes an additional step to find a long-range
Coulomb field self-consistently,
the present work has the more modest
goal of describing only single-solute energetics.
Our connection with solution spatial response functions allows for detailed tests of the theory,
and suggests simple generalizations to self-consistent
spatial and time-dependent problems, for which response data
are already available.<cit.>

  The technical foundation for the theory applied the Gibbs conditioning principle
to transform the grand potential into an exact expression for the solvation free energy.<cit.>
This appears to be a novel route for eliminating
the Hubbard-Stratonovich transform that proceeds directly to a computable density functional,
and for expressing the free energy in conventional density functional theories.<cit.>
Although few density functionals treat molecular orientation in an angular expansion,<cit.>
many ideas from both integral equations and density functional theory
can still be carried over into the present formalism.
Another obvious addition is to include the contribution of energetic degeneracy
of solvent dipole orientations belonging to a given dipole density in Eq. <ref>.
This would begin to address the dispersion contribution identified
for the lattice model in Eq. <ref>.

  There are many important questions which might be addressed with this theory.
We have not attempted to calculate the second step in Fig. <ref>,
which involves forming a cavity at the center of a screened field.
It has, however, been explicitly computed by others<cit.>
who have shown that cavity formation in water creates a positive relative
potential at the center.<cit.>
That asymmetry lowers the free energy of anionic solvation and,
by path equivalence, makes cavity formation in step 2 of Fig. <ref> easier
after polarizing water with a negative potential.

Including both polarization and cavity formation in the first step
(by adding an occupancy constraint)<cit.>
provides a possible route to combining and comparing these ideas.
The limit, Q^-1(0), provides ∂βμ/∂ρ via
the Kirkwood-Buff  theory,<cit.> and, given a predictive theory for Q,
gives a more traditional route to solvation free energies.<cit.>
Further, cluster expansion of Eq. <ref> could be used to include explicit solvent molecules
near the solute as a means to systematically improve approximations to Eq. <ref>
in a quasi-chemical way.<cit.>



§ ACKNOWLEDGMENTS

This work was supported by the USF Research Foundation.










87
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Ren et al.(2012)Ren,
  Chun, Thomas, Schnieders,
  Marucho, Zhang, and Baker]pren12
  author author P. Ren, author J. Chun, author D. G. Thomas, author
  M. J. Schnieders, author
  M. Marucho, author J. Zhang,  and author N. A.  Baker, 10.1017/S003358351200011X
  journal journal Quart. Rev. Biophys. volume 45, pages 427 (year
  2012)NoStop
[Mennucci et al.(1997)Mennucci, Cancès, and Tomasi]bmenn97
  author author B. Mennucci, author E. Cancès,  and author J. Tomasi, 10.1021/jp971959k journal
  journal J. Phys. Chem. B volume 101,  pages 10506 (year 1997)NoStop
[Schnieders and Ponder(2007)]mschn07
  author author M. J. Schnieders and author J. W. Ponder, 10.1021/ct7001336 journal
  journal J. Chem. Theory Comput. volume
  3, pages 2083 (year 2007)NoStop
[Duignan et al.(2013)Duignan, Parsons, and Ninham]tduig13
  author author T. T. Duignan, author D. F. Parsons,  and author B. W. Ninham, 10.1021/jp403595x journal
  journal J. Phys. Chem. B volume 117,  pages 9412 (year 2013)NoStop
[Buchner and Barthel(2001)]rbuch01
  author author R. Buchner and author J. Barthel, 10.1039/B101629F journal
  journal Annu. Rep. Prog. Chem. C volume
  97, pages 349 (year 2001)NoStop
[Skaf and Ladanyi(1995)]mskaf95
  author author M. S. Skaf and author B. M. Ladanyi, 10.1063/1.469368 journal
  journal J. Chem. Phys. volume 102,  pages 6542 (year 1995)NoStop
[Omelyan(1997)]iomel97
  author author I. P. Omelyan, 10.1016/S0378-4371(97)00404-4 journal journal Physica A volume
  247, pages 121 (year 1997)NoStop
[Bopp et al.(1996)Bopp,
  Kornyshev, and Sutmann]pbopp96
  author author P. A. Bopp, author A. A. Kornyshev,
   and author G. Sutmann, 10.1103/PhysRevLett.76.1280 journal journal Phys. Rev. Lett. volume 76, pages 1280 (year 1996)NoStop
[Bopp et al.(1998)Bopp,
  Kornyshev, and Sutmann]pbopp98
  author author P. A. Bopp, author A. A. Kornyshev,
   and author G. Sutmann, 10.1063/1.476884 journal journal J.
  Chem. Phys. volume 109, pages 1939
  (year 1998)NoStop
[Stern and Feller(2003)]hster03
  author author H. A. Stern and author S. E. Feller, 10.1063/1.1537244 journal
  journal J. Chem. Phys. volume 118,  pages 3401 (year 2003)NoStop
[Ballenegger and Hansen(2005)]vball05
  author author V. Ballenegger and author J.-P. Hansen, 10.1063/1.1845431 journal
  journal J. Chem. Phys. volume 122,  pages 114711 (year 2005)NoStop
[Zasetsky and Svishchev(2001)]azase01
  author author A. Y. Zasetsky and author I. M. Svishchev, 10.1063/1.1381055 journal
  journal J. Chem. Phys. volume 115,  pages 1448 (year 2001)NoStop
[Friesen and V.Matyushov(2011)]afrie11
  author author A. D. Friesen and author D. V.Matyushov, 10.1016/j.cplett.2011.06.031 journal journal Chem. Phys. Lett. volume 511, pages 256 (year
  2011)NoStop
[Schaaf and Gekle(2015)]cscha15
  author author C. Schaaf and author S. Gekle,  10.1103/PhysRevE.92.032718 journal journal Phys. Rev. E volume 92, pages 032718 (year 2015)NoStop
[Schaaf and Gekle(2016)]cscha16
  author author C. Schaaf and author S. Gekle,  10.1063/1.4960775 journal journal J. Chem. Phys. volume 145, pages 084901 (year 2016)NoStop
[Matyushov(2014)]dmaty14
  author author D. V. Matyushov, 10.1063/1.4882284 journal
  journal J. Chem. Phys. volume 140,  pages 224506 (year 2014)NoStop
[Bagchi(2001)]bbagc01
  author author B. Bagchi, 10.1021/cr020661+ journal
  journal Chem. Rev. volume 105,  pages 3197 (year 2001)NoStop
[Martin et al.(2014)Martin,
  Fioretto, and Matyushov]dmart14
  author author D. R. Martin, author D. Fioretto,   and author D. V. Matyushov,  10.1063/1.4861965 journal journal J. Chem. Phys. volume 140, pages 035101 (year 2014)NoStop
[Martin and Matyushov(2017)]dmart17
  author author D. R. Martin and author D. V. Matyushov, 10.1063/1.4989641 journal
  journal J. Chem. Phys. volume 147,  pages 084502 (year 2017)NoStop
[Heyden et al.(2012)Heyden,
  Tobias, and Matyushov]mheyd12
  author author M. Heyden, author D. J. Tobias,
   and author D. V. Matyushov,  10.1063/1.4772000 journal journal J. Chem. Phys. volume 137, eid 235103 (year 2012), 10.1063/1.4772000NoStop
[Remsing and Weeks(2016)]rrems16
  author author R. C. Remsing and author J. D. Weeks, 10.1021/acs.jpcb.6b02238 journal
  journal J. Phys. Chem. B volume 120,  pages 6238 (year 2016)NoStop
[Beglov and Roux(1997)]dbegl97
  author author D. Beglov and author B. Roux,  10.1021/jp971083h journal journal J. Phys. Chem. B volume 101, pages 7821 (year 1997)NoStop
[Matyushov(2004)]nrft
  author author D. V. Matyushov, 10.1063/1.1676122 journal
  journal J. Chem. Phys. volume 120,  pages 7532 (year 2004)NoStop
[Martin and Matyushov(2008)]dmart08
  author author D. R. Martin and author D. V. Matyushov, 10.1063/1.3006313 journal
  journal J. Chem. Phys. volume 129,  pages 174508 (year 2008)NoStop
[Perry et al.(1988)Perry,
  Massie, and Cummings]rperr88
  author author R. L. Perry, author J. D. Massie,   and author P. T. Cummings,  @noop  journal journal Fluid Phase
  Equil. volume 39, pages 227 (year 1988)NoStop
[Kovalenko and Hirata(2000)]akova00
  author author A. Kovalenko and author F. Hirata, 10.1063/1.481676 journal
  journal J. Chem. Phys. volume 112,  pages 10391 (year 2000)NoStop
[Maribo-Mogensen et al.(2012)Maribo-Mogensen, Kontogeorgis, and Thomsen]bmari12
  author author B. Maribo-Mogensen, author G. M. Kontogeorgis,  and author K. Thomsen, 10.1021/ie2029943 journal
  journal Ind. Eng. Chem. Res. volume
  51, pages 5353 (year 2012)NoStop
[Raineri and Stell(2001)]frain01
  author author F. O. Raineri and author G. Stell, 10.1021/jp0121163 journal
  journal J. Phys. Chem. B volume 105,  pages 11880 (year 2001)NoStop
[Dyer et al.(2008)Dyer,
  Perkyns, Stell, and Pettitt]kdyer08
  author author K. M. Dyer, author J. S. Perkyns,
  author G. Stell,  and author B. M. Pettitt, 10.1063/1.2976580 journal journal J. Chem.
  Phys. volume 129, pages 104512
  (year 2008)NoStop
[Overduin and Patey(2012)]sover12
  author author S. D. Overduin and author G. N. Patey, 10.1021/jp3075749 journal
  journal J. Phys. Chem. B volume 116,  pages 12014 (year 2012)NoStop
[Zhao et al.(2013)Zhao,
  Liu, Liu, and Wu]szhao13
  author author S. Zhao, author Y. Liu, author H. Liu,  and author
  J. Wu, 10.1063/1.4817784
  journal journal J. Chem. Phys. volume 139, pages 064509 (year
  2013)NoStop
[Chuev et al.(2014)Chuev,
  Vyalov, and Georgi]gchue14
  author author G. N. Chuev, author I. Vyalov,   and author N. Georgi, 10.1002/jcc.23586 journal journal J.
  Comput. Chem. volume 35, pages 1010
  (year 2014)NoStop
[Sheng and Wu(2016)]sshen16
  author author S. Sheng and author J. Wu,  10.1080/00268976.2016.1143567 journal
  journal Mol. Phys. volume 114,  pages 2351 (year 2016)NoStop
[Sergey V. Novikov(1995)]snovi95
  author author A. V. V.  Sergey V. Novikov, 10.1021/j100040a001
  journal journal J. Phys. Chem. volume 99, pages 14573 (year
  1995)NoStop
[Ravichandran and Bagchi(1996)]sravi96
  author author S. Ravichandran and author B. Bagchi, 10.1103/PhysRevE.54.3693 journal
  journal Phys. Rev. E volume 54,  pages 3693 (year 1996)NoStop
[Papazyan and Warshel(1997)]apapa97
  author author A. Papazyan and author A. Warshel, 10.1021/jp971632j journal
  journal J. Phys. Chem. B volume 101,  pages 11254 (year 1997)NoStop
[Essmann et al.(1995)Essmann, Perera, Berkowitz, Darden, Lee, and Pedersen]uessm95
  author author U. Essmann, author L. Perera,
  author M. L. Berkowitz, author T. Darden, author
  H. Lee,  and author
  L. G. Pedersen, @noop  journal journal J. Chem. Phys. volume
  103, pages 8577 (year 1995)NoStop
[Toukmaji et al.(2000)Toukmaji, Sagui, Board, and Darden]atouk00
  author author A. Toukmaji, author C. Sagui,
  author J. Board,  and author T. Darden, 10.1063/1.1324708 journal journal J. Chem.
  Phys. volume 113, pages 10913
  (year 2000)NoStop
[Zhou et al.(1992)Zhou,
  Bagchi, Papazyan, and Maroncelli]hzhou92
  author author H.-X. Zhou, author B. Bagchi,
  author A. Papazyan,  and  author M. Maroncelli, 10.1063/1.463307 journal journal J.
  Chem. Phys. volume 97, pages 9311
  (year 1992)NoStop
[Ballenegger(2014)]vball14
  author author V. Ballenegger, 10.1063/1.4872019 journal
  journal J. Chem. Phys. volume 140,  pages 161102 (year 2014)NoStop
[Johnson(2011)]sjohn11
  author author S. G. Johnson, in @noop  booktitle Casimir
  Physics, series Lecture Notes in Physics, Vol. volume 834 (publisher Springer, year 2011)  Chap. chapter 6, pp. pages 175–218NoStop
[Dyson(1972)]fdyso72
  author author F. J. Dyson, @noop  journal journal Bull.
  Amer. Math. Soc. volume 78, pages 635
  (year 1972)NoStop
[Diehl et al.(1997)Diehl,
  Barbosa, and Levin]adieh97
  author author A. Diehl, author M. C. Barbosa,
   and author Y. Levin,  @noop  journal journal Phys. Rev. E  volume 56, pages 619 (year
  1997)NoStop
[Martin et al.(2016)Martin,
  Li, Delaney, and Fredrickson]jmart16
  author author J. M. Martin, author W. Li, author K. T. Delaney,  and author G. H. Fredrickson, 10.1063/1.4964680 journal journal J.
  Chem. Phys. volume 145, pages 154104
  (year 2016)NoStop
[Wang(2010)]zwang10
  author author Z.-G. Wang, 10.1103/PhysRevE.81.021501 journal
  journal Phys. Rev. E volume 81,  pages 021501 (year 2010)NoStop
[Pujos and Maggs(2015)]jpujo15
  author author J. S. Pujos and author A. C. Maggs, 10.1021/acs.jctc.5b00023 journal
  journal J. Chem. Theory Comput. volume
  11, pages 1419 (year 2015)NoStop
[Frydel and Ma(2016)]dfryd16
  author author D. Frydel and author M. Ma,  10.1103/PhysRevE.93.062112 journal journal Phys. Rev. E volume 93, pages 062112 (year 2016)NoStop
[Xiang and Frydel(2017)]yxian17
  author author Y. Xiang and author D. Frydel,  10.1063/1.4983129 journal journal J. Chem. Phys. volume 146, pages 194901 (year 2017)NoStop
[Rogers(2018)]droge17e
  author author D. M. Rogers, @noop  title Extension of
  Kirkwood-Buff theory to the canonical ensemble,  (year
  2018), note submittedNoStop
[Beck et al.(2006a)Beck, Paulaitis,  and Pratt]pdt6
  author author T. L. Beck, author M. E. Paulaitis,
   and author L. R. Pratt,  title Statistical tentacles, in @noop 
  booktitle The Potential Distribution Theorem and Models
  of Molecular Solutions (publisher Cambridge, New York,  year 2006) Chap. chapter 6, pp. pages
  123–141NoStop
[Beck et al.(2006b)Beck, Paulaitis,  and Pratt]pdt2
  author author T. L. Beck, author M. E. Paulaitis,
   and author L. R. Pratt,  title Statistical thermodynamic necessities, in  @noop  booktitle The Potential Distribution
  Theorem and Models of Molecular Solutions (publisher
  Cambridge, New York, year 2006) Chap. chapter
  2, pp. pages 23–31NoStop
[Wertheim(1971)]mwert71
  author author M. S. Wertheim, 10.1063/1.1676751 journal
  journal J. Chem. Phys. volume 55,  pages 4291 (year 1971)NoStop
[Wertheim(1984)]mwert84
  author author M. S. Wertheim, @noop  journal journal J.
  Stat. Phys. volume 35, pages 19
  (year 1984)NoStop
[Varadhan(2008)]svara08
  author author S. R. S.  Varadhan, 10.1214/07-AOP348 journal journal Ann. Prob. volume
  36, pages 397 (year 2008)NoStop
[Touchette(2009)]htouc09
  author author H. Touchette, 10.1016/j.physrep.2009.05.002 journal journal Phys. Rep. volume
  478, pages 1 (year 2009)NoStop
[Chandler et al.(1986)Chandler, McCoy, and Singer]dchan86
  author author D. Chandler, author J. D. McCoy,
   and author S. J. Singer,  10.1063/1.451510 journal journal
  J. Chem. Phys. volume 85, pages 5971
  (year 1986)NoStop
[Singer and Chandler(1985)]ssing85
  author author S. J. Singer and author D. Chandler, 10.1080/00268978500101591 journal journal Molecular Physics volume 55, pages 621 (year 1985)NoStop
[Léonard and Najim(2002)]cleon02
  author author C. Léonard and author J. Najim, @noop  journal journal
  Bernoulli volume 8, pages 721
  (year 2002)NoStop
[Rogers and Beck(2008)]droge08
  author author D. M. Rogers and author T. L. Beck, 10.1063/1.2985613 journal journal J. Chem. Phys. volume 129, eid 134505 (year 2008)NoStop
[Rogers and Beck(2010)]droge10
  author author D. M. Rogers and author T. L. Beck, 10.1063/1.3280816 journal journal J. Chem. Phys. volume 132, eid 014505 (year 2010)NoStop
[Weerasinghe and Smith(2003)]sweer03
  author author S. Weerasinghe and author P. E. Smith, 10.1063/1.1622372 journal
  journal J. Chem. Phys. volume 119,  pages 11342 (year 2003)NoStop
[Gee et al.(2011)Gee,
  Cox, Jiao, Bentenitis,
  Weeerasinghe, and Smith]mgee11
  author author M. B. Gee, author N. R. Cox,
  author Y. Jiao, author
  N. Bentenitis, author
  S. Weeerasinghe,  and author
  P. E. Smith, 10.1021/ct100517z journal journal J. Chem.
  Theory Comput. volume 7, pages 1369
  (year 2011)NoStop
[Berendsen et al.(1987)Berendsen, Grigera, and Straatsma]hbere87
  author author H. J. C.  Berendsen, author J. R.  Grigera,  and author T. P.  Straatsma, @noop  journal journal J. Phys. Chem. volume 91, pages 6269 (year 1987)NoStop
[Mahoney and Jorgensen(2000)]mmaho00
  author author M. W. Mahoney and author W. L. Jorgensen, 10.1063/1.481505 journal
  journal J. Chem. Phys. volume 112,  pages 8910 (year 2000)NoStop
[Attard(1993)]patta93
  author author P. Attard, 10.1103/PhysRevE.48.3604 journal
  journal Phys. Rev. E volume 48,  pages 3604 (year 1993)NoStop
[Schnell et al.(2013)Schnell, Englebienne, Simon, Krüger, Balaji, Kjelstrup,
  Bedeaux, Bardow, and Vlugt]sschn13
  author author S. K. Schnell, author P. Englebienne, author J.-M. Simon, author P. Krüger,
  author S. P. Balaji, author S. Kjelstrup, author
  D. Bedeaux, author A. Bardow,  and author T. J.  Vlugt, 10.1016/j.cplett.2013.07.043
  journal journal Chem. Phys. Lett. volume 582, pages 154 (year
  2013)NoStop
[Mazzini and Craig(2017)]vmazz17
  author author V. Mazzini and author V. S. J. Craig, 10.1039/C7SC02691A journal
  journal Chem. Sci. volume 8,  pages 7052 (year 2017)NoStop
[Nikoubashman et al.(2012)Nikoubashman, Hansen, and Kahl]aniko12
  author author A. Nikoubashman, author J.-P. Hansen,  and author G. Kahl, 10.1063/1.4748378 journal journal J. Chem. Phys. volume 137, pages 094905 (year 2012)NoStop
[Blum(1975)]lblum75
  author author L. Blum, @noop  journal journal Mol.
  Phys. volume 30 (year
  1975)NoStop
[Frydel and Levin(2013)]dfryd13
  author author D. Frydel and author Y. Levin,  10.1063/1.4802994 journal journal J. Chem. Phys. volume 138, pages 174901 (year 2013)NoStop
[Buyukdagli and Ala-Nissila(2013)]sbuyu13
  author author S. Buyukdagli and author T. Ala-Nissila, 10.1103/PhysRevE.87.063201 journal journal Phys. Rev. E volume
  87, pages 063201 (year 2013)NoStop
[Levy et al.(2013)Levy,
  Andelman, and Orland]alevy13
  author author A. Levy, author D. Andelman,   and author H. Orland, 10.1063/1.4826103 journal journal J.
  Chem. Phys. volume 139, pages 164909
  (year 2013)NoStop
[Borukhov et al.(1997)Borukhov, Andelman, and Orland]iboru97
  author author I. Borukhov, author D. Andelman,
   and author H. Orland, 10.1103/PhysRevLett.79.435 journal journal Phys. Rev. Lett. volume 79, pages 435 (year 1997)NoStop
[Koehl et al.(2009)Koehl,
  Orland, and Delarue]pkoeh09
  author author P. Koehl, author H. Orland,   and author M. Delarue, 10.1021/jp9010907 journal journal J.
  Phys. Chem. B volume 113, pages 5694
  (year 2009)NoStop
[Wachter et al.(2005)Wachter, Kunz, Buchner, and Hefter]wwach05
  author author W. Wachter, author W. Kunz,
  author R. Buchner,  and  author G. Hefter, 10.1021/jp053299m journal journal J.
  Phys. Chem. A volume 109, pages 8675
  (year 2005)NoStop
[Beck(2011)]tbeck13
  author author T. L. Beck, 10.1021/jp204883h journal journal J. Phys. Chem. B volume 115, pages 9776 (year 2011)NoStop
[Chen and Weeks(2006)]ychen06
  author author Y.-G. Chen and author J. D. Weeks, 10.1073/pnas.0600282103 journal
  journal Proc. Nat. Acad. Sci. USA volume
  103, pages 7560 (year 2006)NoStop
[Rodgers and Weeks(2008)]jrodg08
  author author J. M. Rodgers and author J. D. Weeks, 10.1073/pnas.0807623105 journal
  journal Proc. Nat. Acad. Sci. USA volume
  105, pages 19136 (year 2008)NoStop
[Bagchi and Chandra(1989)]bbagc89
  author author B. Bagchi and author A. Chandra, http://dx.doi.org/10.1063/1.456213 journal journal J. Chem. Phys. volume
  90, pages 7338 (year 1989)NoStop
[Salacuse and Egelstaff(2001)]jsala01
  author author J. J. Salacuse and author P. A. Egelstaff, @noop  journal journal
  Phys. Rev. E volume 64, pages 051201
  (year 2001)NoStop
[Rogers(2015)]droge15
  author author D. M. Rogers, @noop  journal journal J.
  Chem. Phys. volume 142, pages 074101
  (year 2015)NoStop
[Ashbaugh and Pratt(2006)]hashb06
  author author H. S. Ashbaugh and author L. R. Pratt, 10.1103/RevModPhys.78.159 journal
  journal Rev. Mod. Phys. volume 78,  pages 159 (year 2006)NoStop
[Molavi Tabrizi et al.(2017)Molavi Tabrizi, Goossens, Mehdizadeh Rahimi, Cooper, Knepley, and  Bardhan]mtabr17
  author author A. Molavi Tabrizi, author S. Goossens, author A. Mehdizadeh Rahimi, author C. D. Cooper, author M. G. Knepley,
   and author J. P. Bardhan,  10.1021/acs.jctc.6b00832 journal journal J. Chem. Theory Comput. volume 13,  pages 2897 (year 2017), note pMID:
  28379697, http://arxiv.org/abs/http://dx.doi.org/10.1021/acs.jctc.6b00832
  http://dx.doi.org/10.1021/acs.jctc.6b00832 NoStop
[Duignan et al.(2017)Duignan, Baer, Schenter, and  Mundy]tduig17
  author author T. T. Duignan, author M. D. Baer,
  author G. K. Schenter,  and  author C. J. Mundy, 10.1063/1.4994912 journal journal J.
  Chem. Phys. volume 147, pages 161716
  (year 2017)NoStop
[Hummer et al.(1996)Hummer,
  Garde, García, Pohorille,  and Pratt]ghumm96
  author author G. Hummer, author S. Garde,
  author A. E. García, author A. Pohorille,  and author L. R. Pratt, http://www.pnas.org/content/93/17/8951.abstract journal
  journal Proc. Nat. Acad. Sci. USA volume
  93, pages 8951 (year 1996), http://arxiv.org/abs/http://www.pnas.org/content/93/17/8951.full.pdf+html
  http://www.pnas.org/content/93/17/8951.full.pdf+html NoStop
[Pratt and Ashbaugh(2003)]lprat03
  author author L. R. Pratt and author H. S. Ashbaugh, @noop  journal journal
  Phys. Rev. E volume 68, pages 021505
  (year 2003)NoStop
[Rogers et al.(2012)Rogers,
  Jiao, Pratt, and Rempe]droge12
  author author D. M. Rogers, author D. Jiao,
  author L. R. Pratt,  and  author S. B. Rempe, in 10.1016/B978-0-444-59440-2.00004-1 booktitle
  Annu. Rep. Comp. Chem., Vol. volume 8, editor
  edited by editor R. A. Wheeler (year 2012) Chap. chapter 4, pp.  pages 71–127NoStop




§  INNER PRODUCT DEFINITION


  The definition of the inner product in Eq. <ref> is straightforward, but extremely
useful both for simplifying our notation and performing computations.
Denote a vector of ν functions each in ℝ^3 →ℂ^d, as f
and another vector as g.  The inner product is the sum-integral,

    fg = ∑_α = 1^ν∫_V f_α(r)^† g_α(r)   dr
    ,

where the integration ranges over the 3D unit cell with volume, V.
The dagger denotes the complex-transpose and is used so that f_α(r)^† g_α(r)
is a scalar.  Note that fg is the complex conjugate of gf,
motivating our adoption of Dirac brackets.

  Defining Fourier transforms of each function as,

    ℱ[f_α](k) ≡∫_V e^-ik· r f_α(r)   dr,

the Fourier-Plancherel theorem proves the equivalence of the inner product with the infinite sum,

    fg = 1/V∑_α = 1^ν∑_k ℱ[f_α](k)^†ℱ[g_α](k)
    .

In Eq. <ref>, the summation is over the reciprocal lattice, k ∈{2π L^-1 z | z ∈ℤ^3},
where L is a 3× 3 matrix whose rows are unit cell translation vectors.
In the text we use the same inner product symbol for both representations, so
ℱ[f]ℱ[g]≡fg.

  In the limit of infinite volume, the reciprocal lattice sum goes over to,

    lim_V→∞1/V∑_k r(k) = ∫_ℝ^3r(k)   dk/(2π)^3
    .


  This work uses hats (as in ΔÛ_α)
to denote functions which depend on the system microstate.
This is distinct from their usual quantum-mechanical interpretation
as operators.  Mathematically, the hat just denotes a random variable here.

  Operator notation is defined similarly to the inner-product notation,

    ⟨f|M⟩g   = ∑_αγ∫_V ∫_V f_α(r) M_αγ(r - r') g_γ(r') dr dr' 
       = 1/V∑_αγ∑_k ℱ[f_α](k)^†ℱ[M_αγ](k) ℱ[g_γ](k)
    .


