Molecular theory of solvation: Methodology summary and illustrations

Integral equation theory of molecular liquids based on statistical mechanics is quite promising as an essential part of multiscale methodology for chemical and biomolecular nanosystems in solution. Beginning with a molecular interaction potential force field, it uses diagrammatic analysis of the solvation free energy to derive integral equations for correlation functions between molecules in solution in the statistical-mechanical ensemble. The infinite chain of coupled integral equations for many-body correlation functions is reduced to a tractable form for 2- or 3-body correlations by applying the so-called closure relations. Solving these equations produces the solvation structure with accuracy comparable to molecular simulations that have converged but has a critical advantage of readily treating the effects and processes spanning over a large space and slow time scales, by far not feasible for explicit solvent molecular simulations. One of the versions of this formalism, the three-dimensional reference interaction site model (3D-RISM) integral equation complemented with the Kovalenko-Hirata (KH) closure approximation, yields the solvation structure in terms of 3D maps of correlation functions, including density distributions, of solvent interaction sites around a solute (supra)molecule with full consistent account for the effects of chemical functionalities of all species in the solution. The solvation free energy and the subsequent thermodynamics are then obtained at once as a simple integral of the 3D correlation functions by performing thermodynamic integration analytically.


Introduction
Nanoscale properties, phenomena, and processes are profoundly different from macroscopic laws governing the behavior of continuous media and materials. All functional features of nanostructures stem from the microscopic properties of constituting atoms but manifest at length scale from one to hundreds nanometers and time scale up to microseconds and more. By changing size, composition, and fabrication protocol, nanostructure properties and processes can be tuned in the widest range [1]. Therefore, predictive modelling of nanosystems should operate at length scales from an Angström to hundreds nanometers and microns and time scales to milliseconds and seconds, and yet derive their properties from the chemical functionalities of the constituents. (An example is biological cellular systems of protein nanomachines operating in crowded environment.) Explicit molecular modelling of such nanosystems involves millions and billions of molecules and is by far not feasible in a "brute force" approach employing just ab initio quantum chemical methods and/or molecular simulations. Addressing this challenge requires the development and use of multiscale methods coupling several levels of description, from electronic structure methods for building blocks and classical molecular simulations for critical aggregates in the system, to statistical-mechanical theories for their large assemblies and mean properties in a statistical ensemble over characteristic size and time scales, to eventually come up with macroscopic scale properties of the nanostructures and related processes showing up in the "real observable world". A true, genuine challenge of multiscale modelling is the development of a theoretical framework that couples methods at different scales, so that observables at lower-level scales are analytically linked in a self-consistent field description to force fields of more coarse-grained models at higher-level scales [2]. Statistical mechanics itself, in particular, integral equation theory of liquids [3], is an example of such a theoretical coupling between microscopic molecular variables and thermodynamic, macroscopic properties.

Statistical-mechanical, molecular theory of solvation
Molecular theory of solvation, also known as reference interaction site model (RISM) [3,4], is based on the first principles foundation of statistical mechanics and Ornstein-Zernike (OZ) type integral equation theory of molecular liquids [3]. It provides a firm platform to handle complex chemical and biomolecular systems in solution. As distinct from molecular simulations which explore the phase space of a molecular system by direct sampling, RISM theory operates with spatial distributions rather than trajectories of molecules and is based on analytical summation of the free energy diagrams which yields the solvation structure and thermodynamics in the statistical-mechanical ensemble. It yields a solvation structure by solving the RISM integral equations for correlation functions and then the solvation thermodynamics analytically as a single integral of a closed form in terms of the correlation functions obtained. Its three-dimensional (3D) version, 3D-RISM theory gives 3D maps of solvent distributions around a solute macromolecule of arbitrary shape [5][6][7][8][9][10][11][12][13][14]. An important component enabling a utility of 3D-RISM theory for complex systems in solution has been the closure relation proposed by Kovalenko and Hirata (KH approximation) [9,12,14]. For simple and complex solvents and solutions of a given composition, including buffers, salts, polymers, ligands and other cofactors at a finite concentration, the 3D-RISM-KH molecular theory of solvation properly accounts for chemical functionalities by representing in a single formalism both electrostatic and non-polar features of solvation, such as hydrogen bonding, structural solvent molecules, salt bridges, solvophobicity, and other electrochemical, associative and steric effects. For real systems, to solve the 3D-RISM-KH integral equations is far less computationally expensive than to run molecular simulations which should be long enough to sample relevant exchange and binding events. This enables handling complex systems and processes occurring on large space and long time scales, which is problematic and frequently not feasible for molecular simulations. The 3D-RISM-KH theory has been validated on both simple and complex associating liquids and solutes with different chemical functionalities [14][15][16][17][18][19][20][21][22][23][24][25], including ionic liquids [19], and polyelectrolyte gels [20], in a range of fluid thermodynamic conditions [4,[21][22][23][24], in various environments such as interfaces with metal [10,12], metal oxide [26], zeolite [13,27], clay [28,29], and in confinement of carbon nanotubes [16], synthetic organic rosette nanotubes [14,[30][31][32], nanocellulose based bionanomaterials [33][34][35], and biomolecular systems [14,. The latter include case studies ranging from the structure of hydrated alanine dipeptide [36][37][38][39], miniprotein 1L2Y and protein G [39]; structural water, xenon, and ions bound to lysozyme protein [40,41]; selective binding and permeation of water, ions and protons in channels [14,20,[41][42][43][44][45]; salt-induced conformational transitions of DNA [45][46][47]; partial molar volume and pressure induced conformational transitions of biomolecules [48][49][50][51][52]; formation and conformational stability of β-sheet Amyloid-β (Aβ) oligomers [53], HET-s Prion and Aβ fibrillous aggregates [54,55]; ligand-binding affinities 32601-2 of seven biotin analogues to avidin in aqueous solution [56]; binding modes of inhibitors of pathologic conversion and aggregation of prion proteins [20,55], binding modes of thiamine against extracytoplasmic thiamine binding lipoprotein MG289 [57,58], binding modes and ligand efflux pathway in multidrug transporter AcrB [59], and maltotriose ligand binding to maltose-binding protein in which structural water triggers binding modes with protein domain motion resulting in the holo-closed structure that binds the ligand [60]; to aqueous electrolyte solution at physiological concentrations in biomolecular systems as large as the Gloebacter violaceus pentameric ligand-gated ion channel (GLIC) homologue in a lipid bilayer [14,20,57] and structural water promoting folding in the GroEL/ES chaperonin complex [61]. The 3D-RISM-KH theory provided an insight into such soft matter phenomena as the structure and stability of gels formed by oligomeric polyelectrolyte gelators in different solvents [20].

3D-RISM integral equations for molecular liquid structure
The solvation structure is represented by the probability density ρ γ g γ (r) of finding interaction site γ of solvent molecules at 3D space position r around the solute molecule (which can be both a macromolecule and supramolecule), as given by the average number density ρ γ in the solution bulk times the normalized density distribution, or 3D distribution function, g γ (r). The values of g γ (r) > 1 and g γ (r) < 1 indicate the areas of density enhancement or depletion, respectively, relative to the average density at a distance from the solute in the solution bulk where g γ → 1. The 3D distribution functions of solvent interaction sites around the solute molecule are determined from the 3D-RISM integral equation [5][6][7][8][9][10][11][12][13][14] h γ (r) = α dr c α (r − r )χ αγ (r ) , (2.1) where h γ (r) and c γ (r) are, respectively, the 3D total and direct correlation functions of solvent site γ around the solute molecule; χ αγ (r ) is the site-site susceptibility of solvent (comprising both the intraand intermolecular terms) which is an input to the 3D-RISM theory; and indices α and γ enumerate all interaction sites on all sorts of solvent species. The diagrammatic analysis relates the density distribution function to the total correlation function as [3] g γ (r) = h γ (r) + 1 , (2.2) and so the latter has the meaning of the normalized probability density of 3D spatial correlations between the solute and solvent molecules, or normalized deviations of solvent site density around the solute molecule from its average value in the solution bulk. As follows from the diagrammatic expansion of the direct correlation function [3], the leading term of its asymptotics beyond the short-range region D sr of the solute-solvent repulsive core and attractive well (first solvation shell maximum), is given by the 3D interaction potential u γ (r) between the whole solute molecule and solvent interaction site γ, scaled by the Boltzmann factor k B times temperature T . Inside the repulsive core, the direct correlation function strongly deviates from the asymptotics (2.3) and assumes the values related to the solvation free energy of the solute molecule immersed in the solvent.

KH closure approximation
To obtain correlation functions, the 3D-RISM integral equation (1) should be complemented with another relation between the two functions h γ (r) and c γ (r) which is called a closure and which also involves the 3D solute-solvent site interaction potential u γ (r) specified with a molecular force field. The exact closure has a nonlocal functional form which can be expressed as an infinite diagrammatic series in terms of multiple integrals of the total correlation function [3]. However, these diagrams are cumbersome and the series suffers from a usual drawback of poor convergence, which makes useless to truncate the series at reasonably low-order diagrams and thus renders the exact closure computationally intractable. Therefore, it is replaced in practice with amenable approximations having analytical features that properly 32601-3 represent physical characteristics of the system, such as the long-range asymptotics of the correlation functions related to electrostatic forces and their short-range features related to the solvation structure and thermodynamics. Examples of closure relations to Ornstein-Zernike type integral equations (including RISM) suitable for liquids of polar and charged species are the so-called hypernetted chain (HNC) closure and the mean spherical approximation (MSA) [3], since they both enforce the long-range asymptotics like (2.3), which is critical for systems with charges. Such other well-known approximations as the Percus-Yevick, Modified Verlet, Martynov-Sarkisov, and Ballone-Pastore-Galli-Gazzillo closures reproduce solvation features for species with short-range repulsion but do not properly account for electrostatics in the interaction potential. By generalization, the 3D-HNC and 3D-MSA closures to the 3D-RISM integral equation (2.1) are constructed [5][6][7][8][9]12]. While enforcing the asymptotics (2.3), the 3D-HNC closure strongly overestimates the association effects and, therefore, the 3D-RISM-HNC equations diverge for macromolecules with considerable site charges solvated in polar solvents or electrolyte solutions, which is almost always the case for biomolecules in aqueous solution. The 3D-MSA closure is free from that shortcoming but suffers from another serious deficiency of producing nonphysical negative values of the distribution function in the areas adjacent to the associative peaks. Those drawbacks are overcome with the approximation proposed by Kovalenko and Hirata (KH closure), the 3D version of which reads [9,12,14] The 3D-KH closure relation (2.4) couples in a nontrivial way the HNC approximation and the MSA linearization, the latter automatically applied to spatial regions of solvent density enhancement g γ (r) > 1, including the repulsive core, and the latter to spatial regions of solvent density enrichment g α (r) > 1 such as strong peaks of association and long-range tails of near-critical fluid phases, and the former to spatial regions of density depletion g α (r) < 1 (including the repulsive core), while keeping the right asymptotics (2.3) peculiar in both HNC and MSA. The distribution function and its first derivative are continuous at the joint boundary g α (r) = 1 by construct. The KH approximation consistently accounts for both electrostatic and non-polar (associative and steric) effects of solvation in simple and complex liquids, non-electrolyte and electrolyte solutions, compounds and supramolecular solutes in various chemical , soft matter [20], synthetic organic supramolecular [14,[30][31][32], biopolymeric [33][34][35], and biomolecular [14, systems. The 3D-KH closure underestimates the height of strong associative peaks of the 3D site distribution functions due to the MSA linearization applied to them [9,62]. However, it somewhat widens the peaks and so 3D-RISM-KH quite accurately reproduces the coordination numbers of the solvation structure in different systems, including micromicelles in water-alcohol solutions [14,[23][24][25], solvation shells of metal-water [9,12], metal oxide-water [26] and mixed organic solvent-clay [28,29] interfaces, and structural water solvent localized in biomolecular confinement [41,60,61]. For instance, the coordination numbers of water strongly bound to the MgO surface are calculated from the 3D-RISM-KH theory with a 90% accuracy and the peak positions within a 0.5 Å deviation, compared to molecular dynamics (MD) simulations [26]. The 3D solvation map g γ (r) of function-related structural water in the GroEL chaperon complex (shown to be strongly correlated to the rate of protein folding inside the chaperon cavity) obtained in an expensive MD simulation with explicit solvent involving 1 million atoms is reproduced from the 3D-RISM-KH theory in a relatively short calculation on a workstation with an accuracy of over 90% correlation for the 3D density map and about 98% correlation for the 3D density maxima [61].
Note that different approximate 3D bridge functions, or bridge corrections, to the 3D-HNC approximation can be constructed, such as the 3D-HNCB closure which reproduced MD simulation results for the height of the 3D distribution peaks of water solvent around the bovine pancreatic trypsin inhibitor (BPTI) protein [62]. A related promising approach is the partial series expansions of order n (PSE-n) of the HNC closure [63]. The PSE-n closures interpolate between the KH and HNC approximations and, therefore, combine numerical stability with enhanced accuracy, as demonstrated for aqueous solutions of alkali halide ions [63,64]. However, a particular appeal of the 3D-KH closure is that it provides an adequate accuracy and the existence of solutions if expected from physical considerations for a wide class of solution systems "from the wild", including various chemical and biomolecular solutes in solution with different solvents, co-solvents, and electrolytes , and moreover, with ligand molecules or fragments consid-ered as part of solvent for efficient 3D mapping of binding affinities in such problems as molecular recognition and fragment based drug design [41,58,59]. Construction of bridge functions for such systems with different solvent components other than water constitutes a significant challenge not addressed so far, and the 3D-KH closure offers a reliable choice for a given system.

Dielectrically consistent RISM theory for site-site correlations of solvent
The radially dependent site-site susceptibility of bulk solvent breaks up into the intra-and intermolecular terms, χ αγ (r ) = ω αγ (r ) + ρ α h αγ (r ) . (2.5) The intramolecular correlation function ω αγ (r ) normalized as 4π r 2 dr ω αγ (r ) = 1 represents the geometry of solvent molecules. [Note that ω αγ (r ) = 0 for sites α and γ on different species.] For rigid species with site separations l αγ , it is represented in the direct space in terms of δ-functions as ω αγ (r ) = δ αγ (r − l αγ )/(4πl 2 αγ ) and for numerical calculations it is specified in reciprocal k-space as where j 0 (x) is the zeroth-order spherical Bessel function. The intermolecular term in (2.5) to be input into (2.1) is given by the site-site radial total correlation function h αγ (r ) between all sites α and γ on all sorts of molecules in bulk solvent. The latter are obtained in advance to the 3D-RISM-KH calculation from the dielectrically consistent RISM theory [65,66] coupled with the KH closure (DRISM-KH approach) [12,14] which can be applied to a bulk solution of a given composition, including polar solvent, co-solvent, electrolyte, and ligands at a given concentration. The DRISM integral equation reads [65,66] h αγ (r ) =ω αµ (r ) * c µν (r ) * ω νγ (r ) + ρ νhνγ (r ) , where c αγ (r ) is the site-site direct correlation function of pure bulk solvent, and both the intramolecular correlation functionω αγ (r ) and the total correlation functionh αγ (r ) are renormalized due to a dielectric bridge correction in a particular analytical form [65,66] that ensures the given phenomenological value as well as consistency of the dielectric constant determined through the three different routes related to the solvent-solvent, solvent-ion, and ion-ion effective interactions in electrolyte solution, ω αγ (r ) = ω αγ (r ) + ρ α ζ αγ (r ) , (2.7b) h αγ (r ) = h αγ (r ) − ζ αγ (r ) . (2.7c) The renormalized dielectric correction enforcing the given phenomenological value of the dielectric constant and the proper orientational behavior and consistency of the dielectric response in the electrolyte solution is obtained in the analytical form specified in the reciprocal k-space as follows [65,66], where j 0 (x) and j 1 (x) are the zeroth-and first-order spherical Bessel functions, r α = (x α , y α , z α ) are the Cartesian coordinates of partial charge q α of site α on species s with respect to its molecular origin; both sites α and γ are on the same species s, and its dipole moment d s = α∈s q α r α is oriented along the zaxis: d s = (0, 0, d s ). Note that the renormalized dielectric correction (2.8) is nonzero only for polar solvent species of electrolyte solution which possess a dipole moment and thus are responsible for the dielectric response in the DRISM approach. The envelope function h c (k) has the value at k = 0 determining the dielectric constant of the solution and is assumed in a smooth non-oscillatory form quickly falling off at wavevectors k larger than the inverse characteristic size l of liquid molecules [65,66], polar is the amplitude related to the dielectric constant of the electrolyte solution which is a phenomenological parameter specified at input. The form (2.8) has also been extended to solvent mixtures [67], where the total number density of polar species ρ polar = s∈polar ρ s and the dielectric 32601-5 susceptibility of the mixture is given in the general case by (2.10) The parameter l in the envelope function (2.9) specifies the characteristic separation from a molecule in solution below which the dielectric correction (2.8) is switched off so as not to distort the short-range solvation structure. It can be chosen to about l = 1 Å for water solvent; however, in solvent of larger molecules or in the presence of such co-solvent, it should be increased so as to avoid "ghost" associative peaks appearing in the radial distributions if the dielectric correction (2.8) interferes with the intramolecular structure of the large solvent species. For example, for octanol solvent, it should be set to about l = 10 Å.
The DRISM integral equation (2.7a) is complemented with the KH closure [12,14] which in the site-site radial version reads where the site-site interaction potential u αγ (r ) (all-atom force field) is typically modelled with Coulomb and Lennard-Jones terms. The DRISM-KH equations (2.7a) and (2.11) keep the same dielectrically consistent asymptotics (2.8) as the original DRISM-HNC theory [65,66], but extend the description to solutions with strong associative species in a wide range of compositions and thermodynamic conditions not amenable to the HNC closure. (The latter overestimates the associative attraction and density inhomogeneities and usually diverges for such systems.)

Analytical expressions for solvation thermodynamics
Much as the HNC approximation, the 3D-KH closure (2.4) to the 3D-RISM integral equation (2.1) has an exact differential of the solvation free energy, and thus allows one to analytically perform Kirkwood's thermodynamic integration gradually switching on the solute-solvent interaction. This gives the solvation free energy of the solute molecule (or supramolecule) in multicomponent solvent in a closed analytical form in terms of the 3D site total and direct correlation functions h γ (r) and c γ (r) [9,12,14]: where the sum goes over all the sites of all solvent species, and Θ(x) is the Heaviside step function.
This offers several substantial advantages over molecular simulations using free energy perturbation techniques with multiple productive runs to perform the thermodynamic integration or solute mutation, which is extremely time consuming [68]. Once the 3D-RISM-KH equations (2.1), (2.4) are converged for the correlation functions c α (r) and h α (r), all the solvation thermodynamics is readily available from the expression (2.12) and analysis of its components. Note that accurate calculation of the solvation free energy from the functional (2.12) requires proper analytical account of electrostatic asymptotics in all the correlation functions, both 1D and 3D, in the integral equations (2.1)-(2.11) [10][11][12][13][14].
The integrand Φ KH γ (r) in (2.12a) is interpreted as the solvation free energy density (3D-SFED) coming from interaction site γ of solvent molecules around the solute. The solvation free energy of the solute macromolecule ∆µ is obtained by summation of the 3D-SFED partial contributions from all solvent species and spatial integration over the whole space. Note that the integration volume V in the form (2.12) comprises both the solvation shells and the solute-solvent molecular repulsive cores, since the direct correlation functions c α (r) inside the latter are related to the free energy of creation of a cavity to accommodate the solute excluded volume. In this way, Φ KH γ (r) characterizes the intensity of effective solvation forces in different 3D spatial regions of the solvation shells and indicates where they contribute the most/least to the entire solvation free energy.
The solvation free energy in the form (2.12) can be split up into components from each solvent species by grouping the terms for the corresponding sites. Further, integrating the 3D-SFED values Φ KH γ ( r) over space regions geometrically attributed to the constituting chemical groups of the solute by a Voronoi decomposition (including the repulsive core region) resolves partial contributions of functional groups of the solute molecule (or supramolecule) to the solvation free energy, which is called spatial decomposition analysis [69,70].
With the solvation free energy obtained in the form (2.12), the potential of mean force w(1, n) between i = 1, n solute molecules (or equally, supramolecule parts or nanoparticles) in solution at given distances and orientations in an arrangement (1, n) is calculated directly by definition as their interaction potential u(1, n) plus the change in the solvation free energy from the sum of ∆µ(i) of separate particles to nonadditive ∆µ(1, n) of the particles brought together [10][11][12][13][14], Similarly to the solvation free energy (2.12), the PMF (2.13) can be split up by spatial decomposition analysis into contributions of solute chemical functionalities and into terms mediated by each solvent species. This route provides an efficient and accurate statistically representative way to calculate effective interactions of various chemical species [10-14, 27, 28], nanomaterials [30][31][32][33][34][35], and biomolecular nanosystems [53][54][55][56][57].
Furthermore, PMFs between a solute molecule / supramolecule / nanoparticle and solvent molecules averaged over molecular orientations centered at solvent interaction sites can be conveniently obtained in the 3D site representation from the definition of the 3D site density distribution functions, (2.14) The solute-solvent PMFs so obtained in a single 3D-RISM-KH calculation provide 3D spatial maps of site binding affinity of solvent species to the solute molecule (positive in attractive wells and negative in repulsive barriers of PMFs). This constitutes a new approach toward molecular recognition and computational fragment-based drug design [41]. With fragmental decomposition of flexible ligands treated as distinct species in solvent mixture of arbitrary complexity, the computed density functions and the corresponding PMFs (2.14) of solvent, ions, and ligand fragments around a biomolecule obtained from the 3D-RISM-KH theory on discrete 3D spatial grids uniquely define the scoring function which can be interfaced to a docking protocol [41,42,57,58], e.g., in the 3D-RISM-Dock approach implemented in the AutoDock suite for automated ranking of docked conformations [58]. As a validation, the 3D-RISM-Dock tool predicted the location and residency times of the modes of binding of a flexible thiamine molecule to the prion protein at near-physiological conditions in an excellent agreement with experiment [57,58]. The capabilities of 3D maps of ligand binding affinity calculated from (2.14) include the prediction of functions in biomolecular systems, such as drug efflux pathway in multidrug transporter AcrB [59].
Further, the solvation free energy (2.12) can be split up into the energetic and entropic contributions [3,71,72], by calculating the solvation entropy at constant volume as and the internal energy of the solute ("u") -solvent ("v") interaction as where the remaining term ∆ε vv gives the energy of solvent reorganization around the solute molecule.
In the same way, the PMF (2.13) or (2.14) can be decomposed into energetic and entropic terms.

32601-7
Further to the solvation free energy, the 3D-RISM-KH theory also provides analytical expressions for other thermodynamic functions, including the partial molar volume (PMV) of the solute molecule obtained in the Kirkwood-Buff theory in terms of the 3D direct correlation functions from the 3D-RISM integral equation [48][49][50][51][52] (2.18) where χ T is the isothermal compressibility of bulk solvent which is analytically obtained in the so-called compressibility route from the RISM integral equation as [23] in terms of the site-site direct correlation functions of bulk solvent c αγ (r ) calculated along with the solvent susceptibility χ αγ (r ) from the DRISM-KH equations, and ρ = s ρ s is the total number density of species s in the solution bulk.
As noted above, the solvation free energy expression (2.12) follows from the functional form of the KH closure (2.4) to the 3D-RISM integral equation (2.1) by carrying out the thermodynamic integration analytically. In practice, the 3D-RISM-KH theory can be used with a different functional of the solvation free energy. In particular, Chandler and co-workers derived the so-called Gaussian Fluctuation (GF) free energy functional based on the assumption of Gaussian fluctuations for the solvent, which is similar to the HNC functional form but without the h 2 γ term [73,74]. The solvation free energy calculated from the GF functional using site-site correlation functions obtained from the RISM-HNC integral equations were generally in reasonable agreement with experiment and in better agreement than those obtained with the HNC functional [75,76]. In the context of 3D-RISM theory, the GF functional reads [36,56] The GF functional using the 3D correlation functions obtained from the 3D-RISM-KH theory provides the hydration free energy in better agreement with experiment in some cases [56,77].
Recently, it was found that the 3D-RISM-KH theory supplemented with the solvation free energy correction based on the partial molar volume of the compound (also obtained from the 3D-RISM-KH theory) provides the hydration free energy for a large diverse set of compounds with an accuracy comparable to much more computationally demanding molecular dynamics simulations with explicit solvent [77].
The correction to the solvation free energy ∆µ calculated from either the KH functional (2.12) or the GF one (2.20) is constructed as a linear function of the PMV in turn calculated in the form (2.18) using the correlation functions from the 3D-RISM-KH theory, ∆µ UC = ∆µ + α ρV + β . The constants α and β in (2.21) are obtained from the linear regression analysis. This partial molar volume correction technique, also called a "Universal Correction" (UC), was parameterized with the GF functional on hydration free energies of a training set of 65 molecules and was tested on a set of 120 molecules [77]. The UC to the KH functional has also been parameterized for the hydration free energy on a set of 504 organic molecules [78]. In the same study, another correction to the KH functional based on the Ng modified free energy functional [79] was introduced and shown to accurately predict hydration free energies of the same set of 504 compounds [78].
With the partial molar volume correction (or the UC term) re-parameterized accordingly, the 3D-RISM-KH theory also provides an excellent agreement with the experimental data for the solvation free energy in non-polar solvent (1-octanol) of a large library of 172 small compounds with diverse functional groups, and so accurately predicts the octanol-water partition coefficients [25]. The best agreement with the experimental data for octanol-water partition coefficients is obtained with the KH-UC solvation free energy functional.

Analytical treatment of electrostatic asymptotics
To properly treat electrostatic forces in electrolyte solution with polar solvent and ionic species in 3D-RISM / DRISM theory requires analytical treatment of the electrostatic asymptotics of the radial sitesite total and direct correlation functions of bulk solvent in the DRISM-KH equations (2.7a), (2.11), as well as of the 3D site total and direct correlation functions in the 3D-RISM-KH equations (2.1), (2.4), the solvation free energy functional (2.12)  (Typical box grid sizes range from 64×64×64 to 256×256×256 nodes.) Analytical forms of the non-periodic electrostatic asymptotics in the direct and reciprocal space are separated out from all the correlation functions before and then added back after the 3D-FFT. Note that although the solvent susceptibility χ αγ (r ) has a long-range electrostatic part, no aliasing occurs in the backward 3D-FFT of the short-range part of h γ (k) on the 3D box supercell since the short-range part of the convolution product h γ (r) contains merely 2-3 oscillations (solvation shells) and vanishes outside the 3D box for physical reasons [13,14,18]. Accordingly, the electrostatic asymptotics terms in the thermodynamic integral (2.12) are handled analytically and reduced to one-dimensional integrals easy to compute [10][11][12][13][14]18].

Accelerated numerical solution
The 3D-RISM-KH integral equations (2.1), (2.4) are converged typically to a root mean square accuracy of 10 −4 −10 −8 and the DRISM-KH equations (2.7a), (2.11) to an accuracy of 10 −6 −10 −12 by using the modified direct inversion in the iterative subspace (MDIIS) accelerated numerical solver [10][11][12][13]80]. MDIIS is an iterative procedure accelerating the convergence of integral equations of liquid state theory by optimizing each iterative guess in a Krylov subspace of typically last 10-20 successive iterations and then making the next iterative guess by mixing the optimized solution with the approximated optimized residual [10,12,13,80]. It is closely related to Pulay's DIIS approach for quantum chemistry equations [81]. and other similar algorithms like the generalized minimal residual (GMRes) solver [82]. The MDIIS solver combines the simplicity and relatively small memory usage of an iterative approach with the efficiency of a direct method. Compared to damped (Picard) iterations, MDIIS provides substantial acceleration with quasiquadratic convergence practically throughout the whole range of root mean square residual, and is robust and stable. Of particular importance is that MDIIS ensures convergence (provided a solution exists) for complex charged systems with strong associative and steric effects, which is usually not achievable with Picard iterations and constitutes a challenging task in the case of 3D integral equations on large 3D grids. Converging the 3D-RISM-KH equations takes minutes to hours on a HPC workstation, depending on the 3D grid size, whereas the DRISM-KH equations converge in seconds for water but much slower, up to hours, for multi-component solvent mixtures with complex species.

Examples of 3D-RISM-KH calculations for the solvation structure and thermodynamics of solution systems and nanoparticles
Sketched below are the capabilities of the 3D-RISM-KH molecular theory of solvation in predictive modelling of some representative molecules and nanosystems in solution: structure and potentials of mean force in aqueous electrolyte solutions, activity coefficients of simple and molecular ions in electrolyte solutions, structural transitions in mixtures of water and amphiphile co-solvent, compound partitioning between water and octanol solvents, and effective nanoscale forces between cellulose nanoparticles in hemicellulose hydrogel maintaining the resilient structure of plant cell walls. The discussion below concerns the main features, capabilities and conclusions, whereas the force fields used and different setup details can be found in the corresponding original work cited here.

Simple ions in aqueous electrolyte solution
One of the basic cases of electrolyte solution systems is aqueous solution of sodium chloride. Figure 1 presents the 3D-RISM results for the solvation structure and potential of mean force for all the pairs of Na + and Cl − ions in aqueous electrolyte solution at infinite dilution and at a relatively high concentration of 1.069 mol/l [10,11]. This generic system presents an illustrative test for a solvation theory to reproduce most of the essential features of a variety of chemical and biomolecular effects in solution, and the 3D-RISM theory succeeds in that.
The 3D site distributions of the water oxygen (O) and hydrogen (H) site distributions around the pair of the Na + and Cl − ions in aqueous solution at infinite dilution are shown in figure 1 (a). Both the O and H distributions form high crowns of the first solvation shell around the ions, with the high O and lower H peaks near the contacts of the crowns corresponding to water molecules bridging the ions. This structure is then followed by the shallow second solvation shells. The corresponding potential of mean force (PMF) calculated using equation (2.13) is shown in figure 1 (c). The H distribution has two crowns around the Cl − ion; their separation from the ion shows that the inner H crown gives the water hydrogens in contact with the Cl − ion while the outer H crown corresponds to the other water hydrogens looking outwards but at the angle determined by the tetrahedral hydrogen bonding structure of the water. On the other hand, there is a single H crown around the Na + ion and the separations of the O and H crowns from Na + show that both the water hydrogens are looking outward Na + , tilted at the same angle. These are typical arrangements of water around a cation oriented like a dipole and water bonded with one hydrogen to an anion. The arrangements are visualized in the cartoons in figure 1 (b) for the 3D water site distributions around both the contact ion pair (CIP) and solvent separated ion pair (SSIP) of the Na + and Cl − ions. Water molecules in contact with both the cation and anion form the dipole like association with the former and the hydrogen bonding with the latter, thus creating a water bridge of strongly associated water molecules located in a ring between the ions. This bridge strongly deepens the solvation contribution to the PMF at the CIP arrangement; interplaying with the ion-ion core repulsion, this results in a significant shift of the first minimum on the PMF to a shorter ion-ion separation, compared to the primitive solvation model just uniformly reducing the Coulomb attraction between the ions by the water dielectric constant The decomposition of the PMFs for the pairs of Na + and Cl − ions is also shown in figure 1 (e), (f). The non-electrostatic (or non-polar) part of the PMFs can be defined either as a difference of the full PMF and its Born electrostatic part (solid curve without symbols in the middle right-hand part for the Na + -Cl − ion pair) or as a PMF obtained from the 3D-RISM calculations for the ion pair solute with the same Lennard-Jones short-range potentials but with the charges of the ions entirely switched off while keeping the charges of water molecules (solid curve with symbols in the middle right-hand part). The difference between these two curves illustrates a strong coupling and interplay of the electrostatic and non-electrostatic forces in the solvation structure. Note that in an aqueous solution or other highly polar solvent, the electrostatic part of the mean solvation forces responsible for dielectric screening always opposes and almost cancels out the Coulomb interaction in the ionic pair and so the resulting PMF comes out effectively scaled down by the high dielectric constant of the polar solvent.
The PMFs in the Na + -Na + and Cl − -Cl − pairs of like ions have the same features of the first and second minima and the barrier between them due to the interplay of the associative forces and molecular structure of the ions and solvent molecules. However, at infinite dilution, the strength of the solvent bridges is not sufficient to overcome the electrostatic repulsion and to stabilize the like ion pairs, and both the first and second minima are local [ figure 1 (c)]. (Note that this can be very different for large molecular ions with weaker electrostatic attraction at the separation determined by their size.) The picture changes for the electrolyte solution at a high concentration of 1 mol/l; numerous salt bridges form in addition to water hydrogen bridges, and the like ion pairs get stabilized in both CIP and SSIP arrangements [ figure 1 (d)].

32601-10
For example, the Cl − -Cl − ion pair in the aqueous solution at this concentration is bridged by several Na + ions and water molecules, forming a cluster depicted in figure 1 (g). The structure of such a cluster follows from the analysis of the 3D-RISM results for the 3D distribution functions of solution species and the corresponding coordination numbers of Na + , Cl − , and water around each ion pair (Na + -Cl − , Cl − -Cl − , and Cl − -Cl − ). Shown for comparison is also the potential of mean force of the Na + -Cl − ion pair in the structureless dielectric continuum (primitive model) of water solvent. Decomposition of the potential of mean force of the ion pairs at infinite dilution into: (d) electrostatic term for the Na + -Na + , Cl − -Cl − , and Na + -Cl − ion pairs (two upper solid curves and lower dashed curve, respectively); (e) non-electrostatic (non-polar) term for the Na + -Cl − ion pair at full ion charges and with ion charges switched off (solid curves with and without symbols, respectively). Part (g): Visualization of a cluster hydrated ions, based on the peaks positions on the 3D maps of density distributions: Cl − ions (shown in gray) in a CIP arrangement stabilized due to bridging by associated Na + ions (in blue) as well as by hydrogen bonded water molecules (water O in red, water H in white) in aqueous electrolyte solution of concentration 1M.

Activities of ions in electrolyte solution
The predictive capability of the 3D-RISM-KH theory (DRISM-KH for simple ions as solutes) in reproducing the solvation thermochemistry of electrolyte solutions in a wide range of concentrations has been validated and shown to be in good agreement with experiment, in particular, for aqueous electrolyte solutions of molecular as well as simple ions in ambient conditions at concentrations from infinite dilution

32601-11
to high ionic strength [64,83]. In particular, the theory yields the activity coefficient of the ion γ i in terms of its solvation free energy ∆µ(c) in solution with solvent density ρ s (c) at electrolyte concentration c with respect to those at infinite dilution c = 0 [64,83], For illustration, figure 2 presents a comparison of the activity coefficient in ambient aqueous solution of sodium chloride at concentrations up to 1 mol/kg calculated from the DRISM-KH theory against experimental data. The calculated activity coefficient is in a good agreement with experiment in the whole range of concentrations.

Structure of water-alcohol mixtures
A further important example of crucial importance in chemistry and biomolecular nanosystems is the formation of nanostructures in solution driven by hydrophobic attraction. Figure 3 illustrates the RISM theory predictions for the structure of the ambient mixtures of water and tert-butyl alcohol (TBA) in the whole range of concentrations [23,24]. TBA is a generic example of primitive surfactant with a hydrophobic head of four carbons and a hydrophilic "tale" represented by the hydroxyl group. The solvation structure of this system successively goes through several stages with TBA concentration in water changing from infinite dilution pure TBA: a separate TBA molecule embedded in a water tetrahedral hydrogen bonding cage at infinite dilution changes to micromicelles of four to six TBA molecules in the head-to-head arrangement incorporated in a water hydrogen bonding cage at about 4% TBA molar fraction; then, the tetrahedral hydrogen bonding structure of water gets disrupted at about 40% TBA molar fraction to be replaced by the zigzag hydrogen bonding structure of alcohol, with separate water molecules embedded in it at infinite dilution of water in TBA. The RISM theory predicted both the structure and thermodynamics of these mixtures, in particular, the concentration and temperature dependence of the compressibility, including the isosbestic point and minimum corresponding to the formation of micromicelles [23,24], in an excellent agreement with the findings from neutron diffraction experiment [84,85] and compressibility measurements (see references in [24]).

Solvation free energy of small compounds in water and octanol
The octanol-water partition coefficient characterizes hydrophobic (lipophilic) and hydrophilic properties of chemical compounds [86,87]. For dilute solutions it is defined as the ratio of molar concentrations of a compound in octanol and water, P O/W = [C O ]/[C W ]. In practice, it is more common to use logarithm of the molar concentration ratio due to the large range of changes of the partition coefficient.  [23,24] versus experimental findings [84,85] (and references in [24]). Part (a): tertbutyl alcohol molecule at infinite dilution incorporated in a cage of water tetrahedral hydrogen bonding; part (b): micromicelle of four tert-butyl alcohol molecules in a cage of the water hydrogen bonding. Part (c): water molecules incorporated in the alcohol zigzag hydrogen bonding structure at infinite dilution of water; part (d): the same as in part (c), but at low water concentration. Part (e): solvation structure features corresponding to the water tetrahedral hydrogen bonding at low alcohol -high water and the alcohol zigzag hydrogen bonding at high alcohol -low water concentrations. Part (f): isosbestic point and minimum of the compressibility corresponding to the formation of TBA micromicelles in water.

32601-12
As an equilibrium constant, the partition coefficient is directly related to the transfer free energy of the compound from water (polar solvent) to octanol (non-polar solvent), The partition coefficient is one of the most important physico-chemical characteristics used in pharmacology, environmental studies, food industry, etc. In particular, in pharmacology, the partition coefficient is used to predict distribution of drugs within the body. It is a crucial parameter that defines drug-likeness of a compound. The partition coefficient, along with other descriptors, defines the efficiency of a drug crossing the blood brain barrier and reaching the central nervous system [88].
Due to the practical importance, many theoretical methods have been developed to predict the partition coefficient, ranging from the phenomenological approaches based on continuum solvation models or different descriptors such as quantitative structure-activity relationship to the sophisticated approaches based on first principle calculations of the transfer free energy between different solvents, including explicit solvent MD simulations and calculations based on different solvation models with account for quantum mechanical effects (see references in [25]). A major drawback of the continuous solvation methods is their inability to treat specific solute-solvent and solvent-solvent interactions such as hydrogen bonding, and their non-transferability (re-parametrization is required for a new solvent composition and possibly thermodynamic conditions). In principle, these difficulties can be avoided by using all-atom explicit sol-

32601-13
vent MD simulations. However, such simulations are computationally demanding and cannot be used in most practical applications such as virtual screening in rational drug design. An alternative to explicit solvent MD simulations, the 3D-RISM-KH molecular theory of solvation with the partial molar volume correction provides an excellent agreement with the experimental data for the solvation free energy in both water and 1-octanol, and so accurately predicts the octanol-water partition coefficient [25].
The 3D correlation functions of solvent around the solute compound are obtained by converging the 3D-RISM-KH integral equations (2.1), (2.4), and then used to calculate the solvation free energy from the KH functional (2.12) supplemented with the UC term (2.21) in turn calculated from the partial molar volume (2.18) using the 3D-RISM-KH correlation functions. Figures 4 and 5 illustrate the performance of the 3D-RISM-KH theory using the KH solvation free energy functional (2.12) without and with the UC term (2.21) on solvation free energies in water and octanol [25] for a large set of small compounds with diverse chemical groups against experimental data [89,90]. The linear coefficients in the UC expression (19) parameterized for the solvation free energy in water and octanol calculated with the KH solvation free energy functional are displayed in Table 3.4. The use of the UC to the KH free energy functional significantly improves the accuracy of the solvation free energy obtained from the 3D-RISM-KH theory (root mean square error (RMSE) drops from 22.84 to 1.975 kcal/mol), and thus makes it comparable to the results of explicit solvent MD simulations [77]. A considerably better performance of the 3D-RISM-KH theory with the KH free energy functional alone is observed in the case of solvation in octanol compared to hydration. Even in this case, the use of the UC further improves the accuracy of the results compared to experimental data (the RMSE from 1.37 down to 1.37 kcal/mol, see Table 3.4). Table 1. Linear coefficients in the universal correction (UC) (2.21) to the Kovalenko-Hirata solvation free energy functional (2.12) for the two solvents [25]. Shown is the root mean square error (RMSE) between the calculated results (with and without UC) and experimental data for the library compounds [89,90] Finally, figure 6 makes a comparison of the octanol-water partition coefficient predicted from the 3D-RISM-KH theory with the KH free energy functional supplemented with the UC term for the library of compounds76 against experimental data [89,90]. Shown for comparison are the results obtained for this library of small compounds by using the generalized Born solvent accessible surface area (GBSA) continuum solvation model widely popular in biomolecular calculations [25]. Compared to experimental data, the prediction accuracy level of 3D-RISM-KH with UC reaches the root mean square error as low as 1.28 which is almost twice better than RMSE of 2.32 for the GBSA results.

Effective interactions between cellulose nanoparticles in hydrogel
One representative illustration of the predictive capabilities of the 3D-RISM-KH molecular theory of solvation comes from the quest for fundamental understanding of the chemically-driven effective nanoscale forces that maintain plant cell wall structure [33][34][35]. Efficient conversion of lignocellulosic biomass to second-generation biofuels and valuable chemicals requires decomposition of resilient cell wall structure of plants. Overcoming biomass recalcitrance constitutes the most fundamental unsolved problem of plant-based green technologies [91,92]. Plants naturally evolved to withstand harsh external mechanical, thermal, chemical, and biological factors. Their secondary cell walls are composed of cellulose microfibrils embedded in a complex non-cellulosic matrix comprised mainly of hemicellulose and lignin which are responsible for the cohesive forces within the cell wall that entail structural support to plants [93]. Current technological applications demand decomposition of this resilient structure in order to extract cell wall components for production of second-generation biofuels and other valuable chemical commodities. It is well known that cell wall recalcitrance varies among plant species and even within different phenotypes of the same plant. Close relationship between recalcitrance and chemical composition of a non-cellulosic matrix suggests that cell wall strength could be tuned by carefully controlling the and with (bottom) the UC term using the correlation functions from the 3D-RISM-KH theory [25]. A comparison against experimental data for the library compounds [89,90]. matrix composition [91,[93][94][95]. Full understanding of the chemical interactions within the cell walls is thus fundamental to gain control of the lignocellulosic biomass recalcitrance [92].
MD simulations have been used to investigate decrystallization of cellulose [96,97], the interactions between cellulose and non-cellulosic components of plant cell walls [98], and the structure and dynamics of lignin [99]. However, extremely time-consuming and costly simulations are required to obtain adequate statistical sampling, addressing both solvation structure and thermodynamics of effective interactions in cell walls based on molecular forces. The 3D-RISM-KH molecular theory of solvation bridges the gap between molecular structure and effective forces on multiple length scales, and is uniquely capable of predicting the chemistry-driven effective interactions in plant cell walls [33][34][35].
Glucuronoarabinoxylan -hemicellulose type most abundant in important lignocellulosic grasses used for biofuel production, such as sugar cane and corn -consists of a xylan backbone decorated with branches of mainly arabinose and glucuronic acid whose amount and ratio substantially varies with the plant genotype [100]. Genetic manipulation of glucuronic acid branching has been shown to significantly improve xylan extractability from cell walls without impairing plant growth [93]. All the chemical changes in the xylan structure regard the branches, mainly arabinose and glucuronic acid ( figure 7, top).
The modelling of structure and stability of plant cell walls consisted of the following stages [33]. (i) Calculation began with constructing a model of a primary cell wall represented with two 4-chain 8-glucoselong cellulose fragments, or cellulose nanocrystallites (CNs), immersed in aqueous solution of arabinose, glucuronic acid, and glucuronate monomers at different concentrations ( figure 7, bottom). The Iβ CNs were built by using the Cellulose-Builder [101] in such a way as to have both the hydrophobic and hydrophilic faces exposed to the solvent. (ii) With an input of the interaction potentials between the species in this cell wall model, the 3D-RISM-KH integral equations were solved for the 3D site correlation func-32601-15 Figure 5. Solvation free energy in octanol obtained from the KH solvation free energy functional without (top) and with (bottom) the UC term using the correlation functions from the 3D-RISM-KH theory [25]. A comparison against experimental data for the library compounds [89,90].
tions of hemicellulose monomers and water solvent ("hydrogel") in the solvation structure of two CNs at different separations and arrangements. (iii) The ensemble-averaged 3D spatial maps of the solvation free energy density (SFED) of the hemicellulose moieties around the CNs as well as the solvation free energies were obtained for each arrangement of the CNs. (iv) Finally, the PMFs between the CNs immersed in the hemicellulose hydrogel were calculated from (2.13).
Figures 8 (a) and 8 (b) depict the pathways of aggregation of CNs approaching each other with their hydrophilic faces (top row) and with their hydrophobic ones (bottom row), respectively. The corresponding PMFs against the separation between the CNs along these pathways are shown in figures 8 (c) and 8 (d), respectively. The PMFs for disaggregation of the CNs in both the hydrophilic (part a) and hydrophobic (part b) face arrangements exhibit two well defined local minima. In both cases, the first minimum at the separation marked as r fc ≈ 0 Å corresponds to two aggregated CNs in direct face contact, whereas the PMF second minimum for CNs 3 Å apart refers to CNs separated by a solvent layer and reaches the well depth of −7 kcal/mol in pure water. These results support the suggestion that CNs aggregation through hydrophilic faces is preferable to the hydrophobic faces in primary cell walls [94].
The stability of aggregated CNs indicates how difficult it is to disrupt their arrangement within a primary cell wall. With increase of the glucuronate concentration, the face contact aggregation free energy ∆G agg = PMF(r fc ) decreases similarly to the first maximum PMF(r bar ), while the disaggregation barrier ∆G dis = PMF(r bar ) − PMF(r fc ) remains the same [figures 8 (c) and 8 (d)]. (The PMFs of CNs in glucuronic acid and arabinose hydrogel have similar shapes.) The larger the amount of hemicellulose, the stronger the primary plant cell wall microstructure. This is consistent with the fact that drastic structure change or absence of hemicellulose prevents plant growth due to structural collapse [100].  Finally, figure 10 shows the isosurfaces of the 3D spatial maps of the solvation free energy density (3D-SFED) coming from glucuronate at the cellulose surface. The solvation free energy (2.12a) of the cellulose nanofibril immersed in hemicellulose hydrogel is obtained by integrating the 3D-SFED contributions (2.12b) from all hydrogel species over the solvation shells. The glucuronate contribution in 3D-SFED varies from large negative values for the thermodynamically most favorable arrangements of glucuronate at the cellulose surface to small negative values for less favorable arrangements. The larger negative value isosurface [figure 10 (a)] is highly localized around polar sites on the cellulose surface and indicates hydrogen bonding of hemicellulose to cellulose. The smaller negative value isosurface [figure 10 (b)] indicates a diffuse second layer of hemicellulose monomers stacking over the cellulose surface. The stacking interactions in the second layer are weaker and less specific than hydrogen bonding and are due to hydrophobic as well as enhanced inter-molecular C-H. . . O interactions found in the crystal structure of cellulose [100]. Although the stacking interactions play a considerable role, the hemicellulose-   cellulose binding is controlled mainly by the site-specific H-bonds. The 3D-SFEDs for arabinose and glucuronic acid differ quantitatively, following the trends for the PMF and ∆G agg in figures 8 and 9.

Conclusion
Based on the first principles foundation of statistical mechanics and Ornstein-Zernike type integral equation theory of molecular liquids, the three-dimensional reference interaction site model with the Kovalenko-Hirata closure relation (3D-RISM-KH molecular theory of solvation) [9][10][11][12][13][14] constitutes an essential component of multiscale modelling of chemical and biomolecular nanosystems in solution. As distinct from molecular simulations which explore the phase space of a molecular system by direct sampling, the molecular theory of solvation operates with spatial distributions rather than trajectories of molecules and is based on analytical summation of the free energy diagrams which yields the solvation 32601-19 structure and thermodynamics in the statistical-mechanical ensemble. It provides the solvation structure by solving the integral equations for the correlation functions and then the solvation thermodynamics analytically as a single integral of a closed form in terms of the correlation functions obtained.
The 3D-RISM-KH theory yields the solvation structure in terms of 3D maps of density distribution functions of solvent interaction sites around a solute molecule with full and consistent account for the effects of chemical functionalities of all solution species. The solvation free energy and subsequent thermodynamics are then obtained at once as a simple integral of the correlation functions by performing the thermodynamic integration analytically. Moreover, the possibility of analytical differentiation of the free energy functional enables: (i) self-consistent field coupling of 3D-RISM-KH with both ab initio type and density functional theory (DFT) quantum chemistry methods for multiscale description of electronic structure in solution [9,12,[14][15][16], (ii) using 3D maps of potentials of mean force as scoring functions for molecular recognition [41,56] and protein-ligand binding in docking protocols for fragment based drug design [41,42,57,58], and (iii) hybrid molecular dynamics simulations performing quasidynamics of a biomolecule steered with 3D-RISM-KH mean solvation forces [17,[36][37][38][39].
Based on the first principles of statistical mechanics, the 3D-RISM-KH theory consistently reproduces, at the level of fully converged molecular simulation, the solvation structure and mean solvation forces in complex chemical and biomolecular nanosystems: both electrostatic forces (hydrogen bonding, other association, salt bridges, dielectric and Debye screening, ion localization) and nonpolar solvation effects (desolvation, hydrophobic hydration, hydrophobic interaction), as well as subtle interplays of these, such as preferential solvation, molecular recognition and ligand binding. This is very distinct from the continuum solvation schemes such as the Poisson-Boltzmann (PB) and Generalized Born (GB) models combined with the solvent accessible surface area (SASA) empirical nonpolar terms and additional volume and dispersion integral corrections, which are parameterized for hydration free energy of biomolecules but are neither really applicable to solvation structure effects in complex confined geometries nor transferable to solvent systems with cosolvent or electrolyte solutions at physiological concentrations.