Note on the physical basis of spatially resolved thermodynamic functions

ABSTRACT

The spatial resolution of extensive thermodynamic functions is discussed. A physical definition of the spatial resolution based on a spatial analogy of partial molar quantities is advocated, which is shown to be consistent with how hydration energies are typically spatially resolved in the molecular simulation literature. It is then shown that, provided the solvent is not at a phase transition, the spatially resolved entropy function calculated by first-order grid inhomogeneous solvation theory [Nguyen CN, Young TK, Gilson MK. Grid inhomogeneous solvation theory: hydration structure and thermodynamics of the miniature receptor cucurbit [7] uril. J Chem Phys. 2012;137(4):044101] satisfies the definition rigorously, whereas that calculated by grid cell theory [Gerogiokas G, Calabro G, Henchman RH, et al. Prediction of small molecule hydration thermodynamics with grid cell theory. J Chem Theory Comput. 2014;10(1):35–48] most likely does not. Moreover, for an ideal gas in an external field, the former theory is shown consistent whereas the latter is not. Finally, consistent with the proposed definition and with the case of an ideal gas in an external field, we derive an approximate expression for the solvent contribution to the free energy of solvation in the limit of infinite dilution from the spatial variation of the density around the solute.

KEYWORDS:

• Solvation thermodynamics
• spatial resolution
• grid cell theory
• grid inhomogeneous solvation theory

1. Introduction

We consider an arbitrary extensive thermodynamic quantity A, say the energy or the entropy. In experiments, typically only changes in A as the system undergoes transformations can be inferred. For model systems, statistical mechanics often allows a calculation of A, either analytically in some special cases or numerically with the aid of computers using the methods of molecular simulation or integral equation theories.

The calculation of A through statistical mechanics naturally involves variables that are hidden in the macroscopic experiment. Mark and van Gunsteren [1] criticised any attempt at decompositions of the entropy (or free energy) in terms of particular atomic interactions for non-ideal systems. Nevertheless, many authors [2–13] have studied decompositions of the form (in this Note, all integrals are to be interpreted as definite integrals over the entire liquid phase), (1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) where $\overline{A}\left(\overrightarrow{x}\right)$ is the local density of A at position $\overrightarrow{x}$ in the system. Since A is an extensive quantity, such a decomposition appears natural. The problem is that the conclusions drawn from an analysis of the spatial behaviour of $\overline{A}\left(\overrightarrow{x}\right)$ depend on the choice of this function, as it is not uniquely defined by Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) ) alone.

To stress the main point of this Note: in theory any spatial decomposition of an extensive thermodynamic function that is consistent with the integrated (macroscopic) value can be obtained and none can be said to be more ‘correct’ than the other from an experimental point of view. In other words, any conclusion can in principle be drawn from an analysis of spatially resolved thermodynamic functions without there being, even in theory, any rigorous way to discern ‘good’ from ‘erroneous’ conclusions.

As a more minor point, we shall appeal to the theory of mixtures and seek to introduce a general axiomatic form for $\overline{A}\left(\overrightarrow{x}\right)$. For simplicity, our starting point will be a solvation process of a monoatomic solute, so that changes in A correspond to changes in the solvent only. Our definition for $\overline{A}\left(\overrightarrow{x}\right)$ will be shown to coincide with how the energy is usually decomposed spatially in the literature. Consistency dictates that the entropy should be spatially resolved in the same manner, not the least to ensure the correct interpretation of any entropy-energy compensation [14].

2. Spatially resolved thermodynamics of solvation

Provided with Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) ), it follows that the total change in A in an arbitrary solvation process is given by, (2) $\mathrm{\Delta }A=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{A}\right(\overrightarrow{x})-{\overline{A}}_{\mathrm{\infty }}]$(2) where ${\overline{A}}_{\mathrm{\infty }}$ is the A-density function evaluated for the bulk solvent (infinitely far removed from the solute). Whereas Equation (2(2) $\mathrm{\Delta }A=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{A}\right(\overrightarrow{x})-{\overline{A}}_{\mathrm{\infty }}]$(2) ) follows from Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) ) in total, we impose the extra condition (‘meanfield approximation’) that the integrand $\overline{A}\left(\overrightarrow{x}\right)-{\overline{A}}_{\mathrm{\infty }}$ equals the local change in A when the solvent is restructured at point $\overrightarrow{x}$ from its initial (unperturbed) bulk state to its (perturbed) solvation state. If, for simplicity, we restrict ourselves to a monoatomic solvent, this local change can only manifest itself through the local fluid density. We therefore seek to express $\overline{A}\left(\overrightarrow{x}\right)$ as a function of the local density, or particle number, only and ignore any dependencies on derivatives, fluctuations or spatial correlations of this quantity. As we shall see in Section 3.1, the definition which we shall present coincides with what is commonly implicitly used for the energy (enthalpy) [2, 4, 7, 10, 13].

Because of correlations, the individual volume elements of the fluid are not independent subsystems. If they had been, the extensivity of A would imply a linear dependence on the number of molecules in the subsystem – at fixed pressure P and temperature T – and $\overline{A}\left(\overrightarrow{x}\right)v\left(\overrightarrow{x}\right)$ could be identified with the change in A upon addition of one molecule in a small volume $v\left(\overrightarrow{x}\right)$ centred at $\overrightarrow{x}$. For interacting subsystems, this relation is not necessarily linear, and we therefore make an infinitesimal change, so that the local A at position $\overrightarrow{x}$ is (3) ${\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P},$(3) and hence $\overline{A}\left(\overrightarrow{x}\right)$ of Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) ) becomes: (3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) Here $\rho \left(\overrightarrow{x}\right)$ and $N\left(\overrightarrow{x}\right)$ denote, respectively, the number density and number of molecules at position $\overrightarrow{x}$, both of which are fluctuating quantities: appropriate ensemble averaging is implicitly understood for both. By the derivative with respect to the number of particles, we actually understand the finite change when adding one particle to the system. The use of the derivative symbol in this case is normally justified by the number of particles in a thermodynamic system being very large. While this argument is not applicable here since the ensemble averaging entails that $N\left(\overrightarrow{x}\right)$ is – in the technical jargon of statistics – ‘almost surely’ zero, it does not deter us and we simply regard the derivative symbol as a notational shorthand for a finite difference.

That $\overline{A}\left(\overrightarrow{x}\right)$ satisfies Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) ) (in the sense that its spatial integral yields A) can be seen the most easily by using a discrete approximation for $\rho \left(\overrightarrow{x}\right)$ in which space is divided into M cells, indexed by i, of a small but finite volume v. Writing (4) ${\overline{A}}_i=\frac{N_i}{v}{\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }{N}_i}\right)}_{T,P}$(4) where $N_i$ is the number of molecules in cell i, we have the differentials (5) $\mathrm{d}A=\sum _i^M{\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }{N}_i}\right)}_{T,P}\mathrm{d}{N}_i$(5) and (6) $\mathrm{d}{\overline{A}}_i=\frac{\mathrm{d}{N}_i}{v}{\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }{N}_i}\right)}_{T,P}$(6) which combined lead to a discrete analog of Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) ), viz. (7) $A=\sum _i^{M}v{\overline{A}}_i.$(7) Two remarks are justified before we continue. First, the quantity $\left(\mathrm{\partial }A/\mathrm{\partial }N\right(\overrightarrow{x}){)}_{T,P}$ is the partial molar A of solvation for a solvent molecule kept fixed at position $\overrightarrow{x}$ (this position is defined with respect to the location of the solute); this is a theoretical device by which diffusion of the solvent in the liquid is taken to be analogous to chemical interconversion, each solvent molecule being defined not only by its chemical identity but also by its physical location. For the bulk fluid, this expression then equals the partial molar A of solvation in itself tout court because of translational symmetry in the fluid. Only in the inhomogeneous fluid (around a solute) will there be a spatial dependence (solvent molecules close to a solute will behave differently from those far in the bulk). Second, in the particular case that A is the Gibbs energy, then the partial derivative in question is the chemical potential which is also constant throughout the fluid, whether homogeneous or not.

3. Spatial decompositions in the literature

Before continuing the development in Section 4, we here digress slightly to see whether some of the spatial decompositions proposed in the literature conform with the axiomatic approach above. We deal with the spatial resolution of the energy and of the entropy separately for reasons that will be clear shortly.

3.1. Energy

In molecular simulations, the total potential energy may be seen as the sum of solvent-solvent $U_{\mathrm{s}\mathrm{s}}$, solute-solvent $U_{\mathrm{x}\mathrm{s}}$ and kinetic contributions. These last ones are entirely local and thus basically ‘resolve themselves’ spatially. The first two contributions (the potential energy) are computed as a sum over individual terms over all of the molecules. This leads to a very intuitive scheme for their spatial decomposition, in that the energy terms involving molecules are projected onto the corresponding positions in which said molecules are located. For the solvent-solvent contributions, this means (8) ${\overline{U}}_{\mathrm{s}\mathrm{s}}\left(\overrightarrow{x}\right)=\frac{1}{2}〈\sum _i\delta ({\overrightarrow{r}}_i-\overrightarrow{x})\sum _{j\ne i}{u}_{\mathrm{s}\mathrm{s}}\left(\right|{\overrightarrow{r}}_j-{\overrightarrow{r}}_i\left|\right)〉$(8) where $\delta \left(\overrightarrow{r}\right)$ is the three-dimensional Dirac delta-function, $u_{\mathrm{s}\mathrm{s}}(\cdot )$ is the pair potential between solvent molecules and the sums run over all molecules positioned at $\left\{{\overrightarrow{r}}_i\right\}$ and the angle brackets denote proper ensemble averaging over all positions $\left\{{\overrightarrow{r}}_i\right\}$. In the averaging procedure, the delta-function picks out only those contributions for which $\overrightarrow{x}$ coincides with a molecular centre. In other words, the energy associated with spatial position $\overrightarrow{x}$ is taken as the average energy required to remove a molecule from position $\overrightarrow{x}$ to infinity, and the factor of 1/2 corrects for double counting.¹ Since only terms with ${\overrightarrow{r}}_i=\overrightarrow{x}$ contribute to the summations, the sum over j becomes independent of the sum over i, and the expression can be rewritten as² (9) ${\overline{U}}_{\mathrm{s}\mathrm{s}}\left(\overrightarrow{x}\right)=\frac{1}{2}〈\sum _i\delta ({\overrightarrow{r}}_i-\overrightarrow{x})〉〈\sum _j{u}_{\mathrm{s}\mathrm{s}}\left(\right|{\overrightarrow{r}}_j-\overrightarrow{x}\left|\right)〉.$(9) Hence, recalling that $\rho \left(\overrightarrow{x}\right)$ is the ensemble-averaged density, we have (10) ${\overline{U}}_{\mathrm{s}\mathrm{s}}\left(\overrightarrow{x}\right)=\frac{1}{2}\rho \left(\overrightarrow{x}\right)〈\sum _i{u}_{\mathrm{s}\mathrm{s}}\left(\right|{\overrightarrow{r}}_i-\overrightarrow{x}\left|\right)〉$(10) and Equation (13(13) $\overline{U}\left(\overrightarrow{x}\right)={\overline{U}}_{\mathrm{s}\mathrm{s}}\left(\overrightarrow{x}\right)+{\overline{U}}_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)$(13) ) is thus seen to satisfy Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ) up to an error, inversely proportional to the total system size, that disappears in the thermodynamic limit.³

The solute-solvent contribution is decomposed in the same way, (11) ${\overline{U}}_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)=〈\sum _i\delta ({\overrightarrow{r}}_i-\overrightarrow{x})u_{\mathrm{x}\mathrm{s}}\left(\right|\overrightarrow{x}\left|\right)〉$(11) where the monoatomic solute is taken to lie at the origin of the coordinate system and $u_{\mathrm{x}\mathrm{s}}$ is the solute-solvent pairwise interaction energy. This leads to (12) ${\overline{U}}_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right)〈u_{\mathrm{x}\mathrm{s}}\left(\right|\overrightarrow{x}\left|\right)〉.$(12) The spatial decomposition of the total potential energy is simply the sum of these two contributions, (13) $\overline{U}\left(\overrightarrow{x}\right)={\overline{U}}_{\mathrm{s}\mathrm{s}}\left(\overrightarrow{x}\right)+{\overline{U}}_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)$(13) and this is essentially the procedure used in Refs. [2, 4, 7, 10, 13] and many others, and it seems to be the only intuitive one. The kinetic energy density, in turn, is simply given by $\frac{3}{2}T\rho \left(\overrightarrow{x}\right)$ for a monoatomic solvent, when taking Boltzmann's constant to be unity, a convention that we shall observe throughout.

3.2. Entropy

There is no intuitive spatial resolution of the entropy as there is for the intermolecular energy. Consequently, there are different suggestions in the literature, of which we shall briefly comment on three. In Section 4, we shall suggest a fourth alternative.

3.2.1. Grid cell theory

We now consider whether the entropy of grid cell theory [4] (GCT; a spatially resolved variant of Henchman's cell theory [15]) conforms to Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ), in order to be consistent with how the energy is resolved in the cited reference. In this theory, the configuration space is divided into discrete cells and the average magnitude of the force of every molecule within a specific cell is computed. The local entropy density of the cell is then computed as,⁴ (14) ${\overline{S}}_i^{\mathrm{G}\mathrm{C}\mathrm{T}}=\frac{N_i}{v}{S}_i^{\mathrm{H}\mathrm{O}}$(14) where $S_i^{\mathrm{H}\mathrm{O}}$ is the entropy of the harmonic oscillator whose average force magnitude $f_i$ equals that of the molecules in cell i: $S_i^{\mathrm{H}\mathrm{O}}\propto \ln f_i$. This choice of matching the molecule of the fluid to the harmonic oscillator is arbitrary and there is no a priori reason to prefer matching the average force magnitude to, say, the average square of the force. However, the choice will impact the calculated spatial entropy distribution. This entails that if the GCT entropy should satisfy definition (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ) for one particular matching choice, it will not for another one.

Since the total entropy is taken as [to satisfy Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) )] (15) $S_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{G}\mathrm{C}\mathrm{T}}=\sum _j^M{N}_j{S}_j^{\mathrm{H}\mathrm{O}}$(15) then $S_i^{\mathrm{H}\mathrm{O}}=(\mathrm{\partial }{S}_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{G}\mathrm{C}\mathrm{T}}/\mathrm{\partial }{N}_i{)}_{T,P}$ only if (16) $\sum _{j\ne i}{N}_j{\left(\frac{\mathrm{\partial }{S}_j^{\mathrm{H}\mathrm{O}}}{\mathrm{\partial }{N}_i}\right)}_{T,P}=-N_i{\left(\frac{\mathrm{\partial }{S}_i^{\mathrm{H}\mathrm{O}}}{\mathrm{\partial }{N}_i}\right)}_{T,P}$(16) which using the chain rule for derivatives, simplifies to (17) $\sum _{j\ne i}\frac{N_j}{f_j}{\left(\frac{\mathrm{\partial }{f}_j}{\mathrm{\partial }{N}_i}\right)}_{T,P}=-\frac{N_i}{f_i}{\left(\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{N}_i}\right)}_{T,P}$(17) where $f_k>0$ is the average magnitude of the molecular forces in cell k and $N_k>0$ is the corresponding number of molecules. This condition expresses a peculiar, and possibly unphysical, condition in that, up to a numerical factor which is unity for the homogeneous fluid, the average change in magnitude of the molecular forces in all other cells should always exactly cancel the corresponding change in the cell where a molecule is added. Of course, for the homogeneous fluid, all derivatives in Equation (17(17) $\sum _{j\ne i}\frac{N_j}{f_j}{\left(\frac{\mathrm{\partial }{f}_j}{\mathrm{\partial }{N}_i}\right)}_{T,P}=-\frac{N_i}{f_i}{\left(\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{N}_i}\right)}_{T,P}$(17) ) vanish and the equation is satisfied trivially. Seeing as any interesting spatial resolution is predicated on the fluid being inhomogeneous around a solute, we shall not pay further attention to this case.

Let us out of simplicity examine this equation in the special case of a system that is divided into two grid cells that together contain all of the molecules in the system. Equation (17(17) $\sum _{j\ne i}\frac{N_j}{f_j}{\left(\frac{\mathrm{\partial }{f}_j}{\mathrm{\partial }{N}_i}\right)}_{T,P}=-\frac{N_i}{f_i}{\left(\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{N}_i}\right)}_{T,P}$(17) ) may then be cast as (18) $\frac{N_2}{f_2}{\left(\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{N}_1}\right)}_{T,P}=-\frac{N_1}{f_1}{\left(\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{N}_1}\right)}_{T,P}$(18) Clearly, this amounts to a requirement that the change in the average magnitude of the molecular force in one cell is opposite in sign to the corresponding change in the other cell when a molecule is added to either. Since the average total force is zero, Equation (18(18) $\frac{N_2}{f_2}{\left(\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{N}_1}\right)}_{T,P}=-\frac{N_1}{f_1}{\left(\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{N}_1}\right)}_{T,P}$(18) ) implies that the distribution of molecular forces should narrow in one cell if it broadens in the other one, when adding a molecule. While this condition is very stringent, it is difficult to prove explicitly that it is not fulfilled by GCT short of performing numerical simulations.

The inconsistency of GCT is more readily proved in a different way. Consider an ideal gas in an external field, so that the gas density is inhomogeneous, an idealised model of the Earth's atmosphere. Now (19) ${\left(\frac{\mathrm{\partial }{f}_k}{\mathrm{\partial }{N}_i}\right)}_{T,P}=0$(19) for arbitrary i and k since for the ideal gas the force depends only on the external field and not on the other molecules. If we for simplicity assume that the field is homogeneous (as it will be over sufficiently short distances), then the local entropy of GCT becomes a spatial constant. However, this means that the ‘gravito-chemical’ potential (in the continuum limit), (20) $\mu \left(\overrightarrow{x}\right)=u_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)+\frac{3}{2}T-T{S}^{\mathrm{H}\mathrm{O}}$(20) is not constant, since the first term depends on $\overrightarrow{x}$ but the other ones do not. This is in clear contradiction with the postulates of equilibrium thermodynamics.

3.2.2. Grid inhomogeneous solvation theory

While not based directly on Equation (2(2) $\mathrm{\Delta }A=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{A}\right(\overrightarrow{x})-{\overline{A}}_{\mathrm{\infty }}]$(2) ), it is nevertheless instructive to investigate to what extent the spatially resolved solvation entropy, computed by first-order grid inhomogeneous solvation theory (GIST) in Ref. [2], satisfies Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ). Note that GIST provides equations for the change of entropy in the solvation process, whereas GCT provides equations for the absolute entropy. This difference is immaterial to our arguments, since the pure solvent entropy is a spatial constant in any case.

In the present notation, the first-order local GIST solvation entropy density is written,⁵ (21) $\mathrm{\Delta }{\overline{S}}_i^{\mathrm{G}\mathrm{I}\mathrm{S}\mathrm{T}}=-\frac{N_i}{v}\ln \left(\frac{N_i}{v{\rho }_{\mathrm{\infty }}}\right)$(21) where ${\rho }_{\mathrm{\infty }}$ is the bulk density, and so (v implicitly depends on $N_i$ due to the condition of constant pressure) (22) ${\left(\frac{\mathrm{\partial }\mathrm{\Delta }{S}_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{G}\mathrm{I}\mathrm{S}\mathrm{T}}}{\mathrm{\partial }{N}_i}\right)}_{T,P}=-[\ln \left(\frac{N_i}{v{\rho }_{\mathrm{\infty }}}\right)+1-\frac{N_i}{v}{\left(\frac{\mathrm{\partial }v}{\mathrm{\partial }{N}_i}\right)}_{T,P}]$(22) Multiplying this result by $N_i/v$ we do not recover Equation (21(21) $\mathrm{\Delta }{\overline{S}}_i^{\mathrm{G}\mathrm{I}\mathrm{S}\mathrm{T}}=-\frac{N_i}{v}\ln \left(\frac{N_i}{v{\rho }_{\mathrm{\infty }}}\right)$(21) ) as we should by Equation (4(4) ${\overline{A}}_i=\frac{N_i}{v}{\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }{N}_i}\right)}_{T,P}$(4) ), and so this proves that Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ) is not satisfied under arbitrary thermodynamic conditions save for an ideal gas (for which $(\mathrm{\partial }v/\mathrm{\partial }{N}_i{)}_{T,P}=v/N_i$ always). However, in the single-phase region (outside the binodal line, i. e. no phase transitions), Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ) is satisfied also for a non-ideal solvent, because the molecular volume is then uniquely defined by the thermodynamic state which does not change upon molecule insertion (no singularity of the compressibility in the thermodynamic limit). This is far from presenting any actual limitation when it comes to biomolecular simulations of aqueous solutions and GIST is thus arguably superior in this respect compared to GCT.

For completeness, we consider also the case of the ideal gas in the external field. In the case of first-order GIST, we have, (23) $\mu \left(\overrightarrow{x}\right)=u_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)+\frac{3}{2}T+T\ln \left(\frac{\rho \left(\overrightarrow{x}\right)}{{\rho }_{\mathrm{\infty }}}\right)$(23) in the continuum limit. But clearly, (24) $\rho \left(\overrightarrow{x}\right)={\rho }_{\mathrm{\infty }}\exp (-\frac{u_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)}{T}).$(24) This means that first-order GIST is consistent with equilibrium thermodynamics since the last term in Equation (23(23) $\mu \left(\overrightarrow{x}\right)=u_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)+\frac{3}{2}T+T\ln \left(\frac{\rho \left(\overrightarrow{x}\right)}{{\rho }_{\mathrm{\infty }}}\right)$(23) ) may be written, (25) $T\ln [\exp (-\frac{u_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)}{T})]=-u_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right).$(25) Unlike GCT, the GIST entropy depends on the fluid density explicitly.

3.2.3. Short comment on 3D-2PT theory

As for Ref. [10], the local entropy is computed using the 2PT theory of Lin and coworkers [16], and its extension to a discrete spatial grid is referred to as ‘3D-2PT’. In this theory, the entropy is related to the dynamical properties of the fluid. Their sampling requires that the trajectories of the molecules be followed for a certain time, which however leads to a problem of precisely ascertaining their position in space and of ascribing the computed entropy solely to that position, however determined. Moreover, because the entropy depends on dynamic properties, the derivative with respect $N_i$ is especially intractable analytically. The question as to whether the local entropy function of Ref. [10] satisfies Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ) thus remains open, albeit finding that it does seems a priori highly unlikely, given the very stringent requirements.

4. Consistent spatial decomposition of the entropy

As shown in the previous section (with the exception of the spatially resolved 2PT theory, although the same conclusion is strongly suspected), the spatial decompositions of the energy and of the entropy found in the literature are in general rigorously inconsistent with each other in the sense that one and the same general definition for spatial resolution is not applicable to both thermodynamic functions, although the actual discrepancies might be numerically small and, indeed, for the first-order GIST, the definition is satisfied if one avoids phase transitions in the solvent phase diagram. Let us now turn to a simple alternative for the spatially resolved entropy function which is also consistent with Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ) as well as with the ideal gas in an external field. Like this we prove constructively that even when the entropy decomposition is consistent with the energy one, there is still arbitrariness remaining in how to define it.

From the definitions of the Gibbs energy, G, and the chemical potential, we may write the local solvent entropy density as, (26) $\overline{S}\left(\overrightarrow{x}\right)=\frac{\overline{H}\left(\overrightarrow{x}\right)}{T}-\frac{\overline{\mu }\left(\overrightarrow{x}\right)}{T}$(26) where $\overline{\mu }\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right)\left(\mathrm{\partial }G/\mathrm{\partial }N\right(\overrightarrow{x}){)}_{T,P}$ is the chemical potential density of the solvent, and $\overline{H}\left(\overrightarrow{x}\right)$ is the solvent enthalpy density. This definition ensures that $\overline{S}\left(\overrightarrow{x}\right)$ satisfy Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ), it being equivalent to (27) $\overline{S}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }(H/T-G/T)}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}.$(27) Neglecting the difference between the energy and enthalpy (which in condensed phases typically amounts to a negligibly small difference of a few joules per mole), we may write (28) $\overline{S}\left(\overrightarrow{x}\right)\approx \frac{\overline{U}\left(\overrightarrow{x}\right)}{T}+\frac{3}{2}\rho \left(\overrightarrow{x}\right)-\frac{\overline{\mu }}{T}$(28) where $\overline{U}\left(\overrightarrow{x}\right)$ is computed as in Section 3.1. The approximation, although very good for liquids and solids under ambient conditions, of substituting the energy for the enthalpy, is mandated by the problem of spatially resolving the pressure-volume term, a problem that remains without suggestions so far in the literature.

Provided the solute lacks internal degrees of freedom, Equation (26(26) $\overline{S}\left(\overrightarrow{x}\right)=\frac{\overline{H}\left(\overrightarrow{x}\right)}{T}-\frac{\overline{\mu }\left(\overrightarrow{x}\right)}{T}$(26) ) is simple to apply directly. In this way, we may write expressions for the entropy and energy of solvation consistent with each other, through Equation (2(2) $\mathrm{\Delta }A=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{A}\right(\overrightarrow{x})-{\overline{A}}_{\mathrm{\infty }}]$(2) ): (29) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}S=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{S}\right(\overrightarrow{x})-{\overline{S}}_{\mathrm{\infty }}]$(29) (30) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}U=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{U}\right(\overrightarrow{x})-{\overline{U}}_{\mathrm{\infty }}]$(30) where $\overline{U}\left(\overrightarrow{x}\right)$ is computed using Equation (13(13) $\overline{U}\left(\overrightarrow{x}\right)={\overline{U}}_{\mathrm{s}\mathrm{s}}\left(\overrightarrow{x}\right)+{\overline{U}}_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)$(13) ), and $\overline{S}\left(\overrightarrow{x}\right)$ is computed using Equation (26(26) $\overline{S}\left(\overrightarrow{x}\right)=\frac{\overline{H}\left(\overrightarrow{x}\right)}{T}-\frac{\overline{\mu }\left(\overrightarrow{x}\right)}{T}$(26) ). The bulk values (denoted by subscript ∞) are computed for spatial points far away from the solute within the liquid phase. Now, when inserting Equation (28(28) $\overline{S}\left(\overrightarrow{x}\right)\approx \frac{\overline{U}\left(\overrightarrow{x}\right)}{T}+\frac{3}{2}\rho \left(\overrightarrow{x}\right)-\frac{\overline{\mu }}{T}$(28) ) into Equation (29(29) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}S=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{S}\right(\overrightarrow{x})-{\overline{S}}_{\mathrm{\infty }}]$(29) ), we obtain (31) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}S=\int {\mathrm{d}}^3\overrightarrow{x}[\frac{\overline{U}\left(\overrightarrow{x}\right)}{T}-\frac{\overline{\mu }\left(\overrightarrow{x}\right)}{T}-\frac{{\overline{U}}_{\mathrm{\infty }}}{T}+\frac{{\overline{\mu }}_{\mathrm{\infty }}}{T}+\frac{3}{2}\{\rho \left(\overrightarrow{x}\right)-{\rho }_{\mathrm{\infty }}\left\}\right]=\frac{{\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}U}{T}-\frac{{\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu }{T},$(31) where the identifications made in the last equality should be obvious after taking into account Equation (30(30) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}U=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{U}\right(\overrightarrow{x})-{\overline{U}}_{\mathrm{\infty }}]$(30) ). Seeing as the solute that we have considered lacks all internal degrees of freedom, ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu$ may be interpreted as the solvent contribution to the free energy of solvation, as is evident by simple rearrangement of Equation (31(31) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}S=\int {\mathrm{d}}^3\overrightarrow{x}[\frac{\overline{U}\left(\overrightarrow{x}\right)}{T}-\frac{\overline{\mu }\left(\overrightarrow{x}\right)}{T}-\frac{{\overline{U}}_{\mathrm{\infty }}}{T}+\frac{{\overline{\mu }}_{\mathrm{\infty }}}{T}+\frac{3}{2}\{\rho \left(\overrightarrow{x}\right)-{\rho }_{\mathrm{\infty }}\left\}\right]=\frac{{\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}U}{T}-\frac{{\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu }{T},$(31) ). Now $\overline{\mu }\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right)\mu$ and ${\overline{\mu }}_{\mathrm{\infty }}={\rho }_{\mathrm{\infty }}\mu$ where μ is the constant chemical potential. The solvent contribution to the free energy of solvation may hence be written, (32) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu =(\mu -\frac{3}{2}T)\int {\mathrm{d}}^3\overrightarrow{x}\left[\rho \right(\overrightarrow{x})-{\rho }_{\mathrm{\infty }}]$(32) in compliance with Equation (2(2) $\mathrm{\Delta }A=\int {\mathrm{d}}^3\overrightarrow{x}\left[\overline{A}\right(\overrightarrow{x})-{\overline{A}}_{\mathrm{\infty }}]$(2) ).

That Equation (32(32) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu =(\mu -\frac{3}{2}T)\int {\mathrm{d}}^3\overrightarrow{x}\left[\rho \right(\overrightarrow{x})-{\rho }_{\mathrm{\infty }}]$(32) ) expresses the free energy of solvation as a functional of the liquid density, is the result of our assumptions of a monoatomic solute and solvent and is the explicit end result of the goal introduced in the opening paragraph of Section 2. However, Equation (32(32) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu =(\mu -\frac{3}{2}T)\int {\mathrm{d}}^3\overrightarrow{x}\left[\rho \right(\overrightarrow{x})-{\rho }_{\mathrm{\infty }}]$(32) ) is not exact, even for the structureless monoatomic solute within the meanfield approximation. The most obvious deficiency is our neglect of the pressure-volume term in Equation (28(28) $\overline{S}\left(\overrightarrow{x}\right)\approx \frac{\overline{U}\left(\overrightarrow{x}\right)}{T}+\frac{3}{2}\rho \left(\overrightarrow{x}\right)-\frac{\overline{\mu }}{T}$(28) ) which means that Equation (32(32) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu =(\mu -\frac{3}{2}T)\int {\mathrm{d}}^3\overrightarrow{x}\left[\rho \right(\overrightarrow{x})-{\rho }_{\mathrm{\infty }}]$(32) ) is at best valid only at low pressures. Moreover, our neglect of the Gibbs–Duhem equation is a second source of error. Whereas the chemical potential of the solvent is spatially constant at equilibrium, it is not conserved in the solvation process, invalidating Equation (32(32) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu =(\mu -\frac{3}{2}T)\int {\mathrm{d}}^3\overrightarrow{x}\left[\rho \right(\overrightarrow{x})-{\rho }_{\mathrm{\infty }}]$(32) ) except at high dilution where this change is negligible. Finally, Equation (32(32) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu =(\mu -\frac{3}{2}T)\int {\mathrm{d}}^3\overrightarrow{x}\left[\rho \right(\overrightarrow{x})-{\rho }_{\mathrm{\infty }}]$(32) ) only gives the solvent contributions to the free energy of solvation and the numerical value of μ (for a given solvent) is fixed by choice of reference point for the solvent potential energy U only without any regard for solute or solute-solvent energy terms.

With this in mind, we note that Equation (32(32) ${\mathrm{\Delta }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}\mu =(\mu -\frac{3}{2}T)\int {\mathrm{d}}^3\overrightarrow{x}\left[\rho \right(\overrightarrow{x})-{\rho }_{\mathrm{\infty }}]$(32) ) could, for instance, be used to interpret relative contributions to the free energy of solvation of a biomolecule from different parts of the molecule, from the fluid density in the ‘solvation layer’. Provided an unambiguous definition of the ‘solvation layer’ can be established, different analogs could in this way be compared. But great care must be taken in any such analysis to make sure the analysis is physically meaningful. Exactly how such an analysis could proceed is beyond the scope of this note.

5. Conclusion: what is gained by spatial resolution?

The experimental thermodynamic observables correspond to the definite integral in Equation (1(1) $A=\int {\mathrm{d}}^3\overrightarrow{x}\overline{A}\left(\overrightarrow{x}\right)$(1) ), that is the total A, and are indifferent to the precise choice of the decomposition function $\overline{A}\left(\overrightarrow{x}\right)$ of which there is an infinite number. Ideally, spatial decomposition offers possible computational advantages in that a single molecular simulation – to be more precise, a single and physical $6N$-dimensional Hamiltonian corresponding only to the intermolecular potential and kinetic energies – may capture both bulk and solvation properties from different spatial regions of one and the same physical simulation, thus negating the need to run separate simulations (or use interpolating Hamiltonians with non-physical degrees of freedom) to compute the difference in A.

However, if the spatial variations of $\overline{A}\left(\overrightarrow{x}\right)$ themselves are to be interpreted in a physically meaningful way, it ought to be clear to what they correspond. But this is not the only concern: whereas GIST leads to a computed chemical potential that is constant at equilibrium, that obtained by GCT is not, in contradiction with the postulates of thermodynamics, unless a different way of spatially resolving the energy is chosen. However, taking the spatially resolved energy as anything other than the local binding energy would seem to lead to highly contrived analyses if one extends the argument to electronic degrees of freedom.

Nevertheless, although I have argued for maintaining consistency in the way entropy and energy are spatially resolved, and showed the nice – and consistent – properties exhibited by the definition of Equation (3(3) $\overline{A}\left(\overrightarrow{x}\right)=\rho \left(\overrightarrow{x}\right){\left(\frac{\mathrm{\partial }A}{\mathrm{\partial }N\left(\overrightarrow{x}\right)}\right)}_{T,P}$(3) ), it is clear from the derivation in Section 4 that the self-consistency requirement for the chemical potential amounts at most to a necessary, and not a sufficient, physical condition to define the spatial resolution. Therefore, like Mark and van Gunsteren [1] before me, I still strongly caution against the overinterpretation of results from any molecular simulations that report spatially resolved thermodynamic functions, regardless of the way that the spatial resolution is obtained. This will remain the state of affairs until we have enough physical conditions to impose to make such an analysis unique and hence edifying.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 Including this factor leads to the spatial decomposition of the average energy per molecule; otherwise it leads to that of the total binding energy: the energy difference of the system with respect to its isolated constituent molecules. In the literature, it is more common to see the former spatially resolved than the latter, but it is a subtle question which one is actually more appropriate for local thermodynamics.

2 Whether one excludes from the average, terms for which ${\overrightarrow{r}}_j=\overrightarrow{x}$, or simply takes $u_{\mathrm{s}\mathrm{s}}\left(0\right)$ finite by definition, does not affect the converged ensemble average.

3 Since the relative fluctuations of the energy are inversely proportional to the square-root of the total system size, the error in Equation (13(13) $\overline{U}\left(\overrightarrow{x}\right)={\overline{U}}_{\mathrm{s}\mathrm{s}}\left(\overrightarrow{x}\right)+{\overline{U}}_{\mathrm{x}\mathrm{s}}\left(\overrightarrow{x}\right)$(13) ) is negligible next to the thermal fluctuations already for systems of a few thousand molecules such as in molecular simulations.

4 For notational simplicity, we only consider the translational contribution, as we have only dealt with monoatomic molecules elsewhere.

5 Only the translational contribution is considered.