Multiscale Methods Framework with the 3D-RISM-KH Molecular Solvation Theory for Supramolecular Structures, Nanomaterials, and Biomolecules: Where Are We Going?

: 3D-RISM-KH molecular solvation theory based on statistical mechanics has been an engine of the multiscale methods framework, which also includes molecular simulation techniques. Its applications range from the solvation energy of small molecules to the phase behavior of polymers and biomolecules. Molecular solvation theory predicts and explains the molecular mechanisms and functioning of a variety of chemical and biomolecular systems. This includes the self-assembly and conformational stability of synthetic organic rosette nanotubes (RNTs), the aggregation of peptides and proteins related to neurodegeneration, the binding of ligands to proteins, and the solvation properties of biomolecules related to their functions. The replica RISM-KH-VM molecular solvation theory predicts and explains the structure, thermodynamics


Introduction
Molecular properties often differ significantly from continuous media mostly due to the different scalability of these systems.The governing factors in such different behaviors are rooted in the different forms arising due to size, shape, and composition(s), as well as physical state(s).Thus, from the perspective of the building materials of desired properties, understanding all the underlying interactions between different constituent fragments is essential.This is a difficult task from an experimental chemistry point of view, due to the sheer number of possible combinations between structures and resultant reactivity and/or activity.Theoretical prediction is the most useful tool to address this conundrum, as prediction can extend over the entire range of the size scale and also to the larger part of the time scale.This allows predictive modeling capabilities to address issues related to single molecules from the realm of complex biomolecules to the world of nanomaterials.The limitation of predictive modeling is the incompatibility of sophisticated electronic structure methods and/or molecular-dynamics-based simulation techniques for systems of a realworld size constituting millions of atoms, if not billions, as the computing cost becomes astronomical.A true multiscale predictive method should thus possess sufficient ability to couple lower-scale methods to higher-scale ones by propagating the accuracy of the lower-scale methods to a large extent to the higher-scale systems without the impediment of extensive computational requirements.It is important to note that the sophistication level of gas-phase chemical physics is most often impossible to apply to real-world systems due to not only the size issue but also to the presence of the solvent medium interacting directly.
For a system embedded in a continuum medium, quantum mechanical calculations using continuum solvation models generally provide reliable outcomes, provided the systems are of reasonable size to be treated with at least a triple-ζ basis set.However, the presence of heavy elements renders such a methodology useless.The applicability of the Onsager equation in predicting solvation free energy has been questioned recently [1].An extensive search of the available literature clearly indicates a clear classification of the possible schemes in handing (bio)chemical processes in solution, either by using continuum solvation models, or the so-called cluster continuum models and quasi-chemical models, all with system size limitations [2][3][4][5][6][7][8][9][10].The other class of methods based on statistical mechanics, viz., the reference interaction site model (RISM), has been gaining popularity in addressing molecular systems in liquid media due to its reasonable accuracy, speed of calculation, and extendibility over the entire scale of the molecular size and time scales [11,12].The molecular solvation theory with integral equation formalism is based on the modified Ornstein-Zernike (OZ) theory for the dimensional reduction of molecular liquids [13].The spatial distributions and statistical mechanical ensembles are key to the RISM theory in order to predict solvation thermodynamics and behavior via integration over an infinite number of interaction diagrams.The three-dimensional version, 3D-RISM theory [14][15][16][17][18][19][20][21][22], gives 3D maps of the distributions of a solvent around a solute macromolecule of arbitrary shape by integrating a single integral of correlation functions.A successful implementation of the MOZ equation for the 3D-RISM application needs a mathematical function, a closure relation, which can be generalized as a functional to impose a consistency condition of the path-independent chemical potential µ.Despite over forty years of developments in RISM theory, only a handful of closure relations exist, varying substantially in accuracy and the scope of application [23][24][25].The most promising closure relation, the Kovalenko-Hirata (KH) approximation, has been shown to work for small molecules, biomolecules, nanomaterials, aggregates, and ligand-protein interactions with high accuracy in the whole range of thermodynamic conditions [26][27][28][29][30][31].The replica RISM-KH-VM molecular solvation theory predicts the properties of electrolyte solutions sorbed in disordered nanoporous electrodes [32][33][34][35][36][37].The 3D-RISM-KH theory is also applied to protein-folding problems, the drug efflux mechanism, clay minerals, ionic solvents, and electrolytes in solutions [38][39][40][41][42].It has been shown that the theoretical framework of the 3D-RISM-KH closure significantly accelerates sampling over slow and rare solvation events in complex biomolecular systems, e.g., solvent and ion exchanges and localization in biomolecular simulations, and molecular recognition, including protein-ligand binding, DNA association, and protein-protein (or peptide) interactions [43,44].One can also use electronic structure calculations with the 3D-RISM-KH molecular theory of solvation with various levels of sophistication [45].

Theoretical Background
Molecular solvation theory, in particular, the 3D-RISM-KH integral equation for molecular liquids, relies on the probability density ρ γ g γ (r) of finding interaction site γ of solvent molecules at 3D space position r around a solute molecule or arbitrary shape, with the average solvent site number density ρ γ in the solution bulk and the normalized solvent site density distribution g γ (r).This 3D solvent site distribution function g γ (r) assumes the values g γ (r) > 1 or g γ (r) < 1 in the areas of density enhancement or depletion, respectively.At a distance from the solute, it levels out to g γ →1.The 3D solvent site distribution functions around the solute molecule are obtained from the 3D-RISM integral equation where h γ (r) = g γ (r) − 1 is the 3D total correlation function of solvent interaction site γ, c γ (r) ~−u γ (r)/(k B T) is the 3D direct correlation function which has the asymptotics of the solute-solvent site interaction potential u γ (r), the site-site susceptibility of pure solvent χ αγ (r) is an input to the 3D-RISM theory, indices α and γ enumerate all interaction sites on all sorts of solvent species, T is temperature, and k B is the Boltzmann constant [14][15][16][17][18][19][20][21][22].
A closure relation connecting the 3D total and direct correlation functions complements the 3D-RISM integral of Equation ( 1), providing a computational handle to integrate the infinite chain of diagrams.The exact closure can be formally expressed as a series in terms of multiple integrals of the combinations of the total correlation functions, which is cumbersome, and in practice, is replaced with tenable approximations.For instance, the KH closure approximation accounts, in a consistent manner, for both electrostatic and non-polar features (i.e., associative and steric effects) of solvation in simple and complex liquids, and has the form: where u γ (r) is the 3D interaction potential between the whole solute and solvent site γ specified by the molecular force field.The 3D-KH closure in Equation ( 2) couples, in a nontrivial way, the so-called mean spherical approximation (MSA) applied to spatial regions of solvent density enrichment g γ (r) > 1 and the hypernetted chain (HNC) one applied to the regions of density depletion g γ (r) < 1.The 3D solvent site distribution function and its first derivative are continuous at the joint boundary g γ (r) = 1 by construct.The 3D-RISM-HNC theory is known to overestimate solvation structures in strongly associated systems, thus imparting numerical instability and also diverging for strongly charged systems like electrolytes in solution.On the other hand, the 3D-RISM-KH theory can handle all such systems with ease.There are several other closure relations which have been developed over the years but they are aimed at specific applications, lacking the generality of the KH closure [46][47][48][49].As a critical drawback, KH closure underestimates the height of strong associative peaks.On the other side, it widens the peak to some extent, and so provides the correct thermodynamics and solvation structure [50,51].The site-site susceptibility of solvent breaks up into the intra-and intermolecular terms, where the normalized intramolecular correlation function ω αγ (r) represents the geometry of solvent molecules.For rigid species with site separations l αγ , it has the form, ω αγ (r) = δ(r − l αγ )/(4πl 2 αγ ) specified in the reciprocal k-space in terms of the zeroth-order spherical Bessel function j 0 (x) as and h αγ (r) is the intermolecular radial total correlation function between solvent interaction sites α and γ.The solvent site-site total correlation functions h αγ (r) in Equations ( 1) and ( 3) are obtained in advance from the dielectrically consistent RISM theory [52] coupled with the KH closure (DRISM-KH approach).It is applied to the bulk solution of a given composition, including polar solvent, co-solvent, electrolyte, and ligands at a given concentration.The DRISM integral equation has the form [52] where c αγ (r) is the site-site direct correlation function of a bulk solvent, and both the intramolecular site-site correlation functions ω αγ (r) and the total site-site correlation functions h αγ (r) are renormalized via an analytical dielectric bridge correction, ensuring all inter-and intra-species interactions in liquid state, The renormalized dielectric correction has the following form in the reciprocal k-space [52]: where j 0 (x) and j 1 (x) are the zeroth-and first-order spherical Bessel functions over the positions of each atom r α = (x α , y α , z α ) with a partial site charge q α of site α on species s with respect to its molecular origin, where both sites α and γ are on the same species s, and its dipole moment d s = ∑ α∈s q α r α is oriented along the z-axis, d s = (0, 0, d s ).
The renormalized dielectric correction (8) is nonzero only for polar solvent species of sorbed electrolyte solution which possess a dipole moment and are responsible for the dielectric response in the DRISM approach.The value of the envelope function h c (k) at k = 0 determines the dielectric constant of the solution, and is assumed in the ae smooth non-oscillatory form, quickly falling off at wavevectors k larger than those corresponding to the characteristic size l of liquid molecules and hence, to the dielectric constant (ε) of the solvent, and A = 1 A being the amplitude.The form (8)- (10) applies to mixed solvents with the total number density of solution polar species and the dielectric susceptibility of the solution The parameter l in the dielectric correction (9) specifies the characteristic separation from a liquid molecule, below which this correction is switched off so as not to distort the short-range solvation structure.It can be chosen as l = 1 Å for a water solvent; however, the value of l should be carefully chosen for solvents with a large molecular diameter.
The KH closure relation to the DRISM integral equations has the following form: Notably, the 3D-RISM integral equation has an exact differential of the solvation free energy for both the HNC and KH closures, allowing an analytical expression of Kirkwood's Thermo 2023, 3 thermodynamic integration that gradually switches the solute-solvent interaction on.The solvation free energy of a molecule in multicomponent solvent that follows from the 3D-RISM-KH integral Equations ( 1) and ( 2) is obtained in a closed analytical form as a single integral of the 3D solvent site correlation functions: where Θ(x) is the Heaviside step function.The integrand Φ γ (r) in Equation (13a) can be taken as the 3D-solvation free energy density arising due to all solvent-solute interactions.
The solvation free energy of the solute molecule ∆µ is obtained by summation of the partial contributions over all solvent sites integrated over the whole volume.Other thermodynamic quantities are derived from the solvation free energy (13) via differentiation.This includes the solvation chemical potential which is decomposed into the energetic and entropic components at a constant volume, where entropy at constant volume is the internal energy of the solute-solvent ("uv") interaction is and the remaining term ∆ε vv gives the energy of solvent reorganization around the solute.
The partial molar volume of the solute macromolecule is obtained from the Kirkwood-Buff theory [53] extended to the 3D-RISM formalism as [54,55] where χ T is the isothermal compressibility of bulk solvent obtained in terms of the site-site direct correlation functions of the bulk solvent as where ρ = ∑ s ρ s is the total number density of the bulk solvent mixture of molecular species s.To apply 3D-RISM-KH theory in calculating absolute solvation free energy, one has to be careful, as the theory overestimates solvation by a large margin, although it produces correct trends.To avoid overestimation, one can use the so-called universal correction (UC) scheme [56]: The chemical potential adjusted with the UC is obtained from the Gaussian Fluctuation (GF) chemical potential, corrected with partial molar volume (PMV).The correlation coefficients, a and b, in this expression are obtained from multiple linear regression analysis carried out against benchmarking results.Such applications were reported recently [57,58].
From a computational viewpoint, the most efficient way of solving the 3D-RISM-KH equation is via the modified direct inversion in the iterative subspace (MDIIS) accelerated numerical solver [59], a simple and less computer memory demanding protocol to solve the 3D site direct correlation function for solvent-solute site-site interactions, to zero the residuals of the equations, which are nonlocal functionals of 3D-DCF and arise via a difference between the 3D-site distribution function generated from Equation ( 1) and the closure relation.MDIIS is simple, robust, and stable, has relatively small memory usage, and provides substantial acceleration with quasiquadratic convergence in the whole range of root mean square residual.It reliably converges for complex charged systems with strong associative and steric effects, which is challenging for 3D integral equations.MDIIS is closely related to Pulay's DIIS solver to accelerate and stabilize the convergence of the Hartree-Fock self-consistent field equations [60].Other similar algorithms include the generalized minimal residual (GMRes) solver [61] which was also coupled with a Newton-Raphson-like approach on a coarse grid [62], and is limited to the solute repulsive core [63].
A direct improvement in the treatment of the required solvation shells for correct 3D-RISM-KH calculations was achieved by splitting the 3D-cubic box into a core region encapsulating the excluded volume of the solute macromolecule, followed by about three solvation shells around the solute and the remainder of the box region [64].The core region is built by analogy to solute cavity formation in the classical continuum models [65], via rolling a solvent site around each solute site.Another recent advance in closure relations is the Kobryn-Gusarov-Kovalenko closure, which is built on the KH closure as The performance of this newest closure relation has been under investigation for use in various types of systems [66,67].

Electrical Double Layer in Nanoporous Materials
Electric double layers (EDLs) in nanoporous carbon electrodes are unlike those in planar electrochemical capacitors because of the overlap of EDLs.At the inner surface of nanopores, the EDL is distorted compared to a planar electrode, and has the specific capacitance about 1-2 orders of magnitude less than that of a planar electrode.Furthermore, another EDL at the outer macroscopic surface of nanoporous material dramatically contributes to the specific capacitance of nanoporous electrodes.Because of the interplay of long-range electrostatic and short-range steric interactions and the chemical and mechanical balance between the sorbed electrolyte solution and that in the bulk, consistent molecular simulation of these systems is practically unfeasible.Conventional modeling either uses analytical description and MD simulation for an EDL of finite size ions in slit-like pores with no molecular solvent and chemical specificities [68], or runs an all-atom MD simulation of ions and the solvent in confinement with simplified geometry [69].
The replica RISM-KH-VM (modified Verlet) theory comes from generalizing RISM molecular solvation theory to solutions sorbed in disordered nanoporous material [32][33][34][35][36][37].It makes predictive modeling of sorption and supercapacitance of electrolyte solutions in functionalized nanoporous carbon electrodes feasible.In particular, it describes solventspecific wetting and water depletion in hydrophobic carbon nanopores, asymmetry in solvation and adsorption of cations and anions, specific adsorption in functionalized carbon nanopores, desalination of ions in hydrophobic nanopores, and the removal of ionic impurities from an aqueous waste stream in a nanoporous electrosorption cell.
The replica formalism of statistical-mechanical, integral equation theoriy of quenchedannealed systems treats "annealed" fluid with equilibrium temperature T 1 (species 1), sorbed in a porous matrix of "quenched" particles with a spatial distribution corresponding to an equilibrium ensemble with temperature T 0 (species 0).The mean free energy of the sorbed annealed fluid is obtained as a statistical average of the free energy with the canonical partition function Z 1 (q 0 ) of the fluid sorbed in the matrix with a particular spatial configuration of quenched particles q 0 over the ensemble of all realizations of matrix configurations q 0 , where A id 1 is the ideal gas free energy.The statistical average of a logarithm is not amenable to standard evaluation and is obtained by using the so-called replica identity, similar to the theory of spin glasses, that relates the logarithm to the analytical continuation of moments Z s as follows: ln Z 1 = lim s→0 dZ s /ds.The statistical average of the moments is given by the equilibrium canonical partition function of a fully annealed (s + 1)-component liquid mixture of matrix species 0 and s equivalent replicas of fluid species 1 not interacting with each other.The average free energy of the annealed fluid is obtained, assuming no replica symmetry breaking in the analytical continuation of the annealed replicated free energy A rep (s), The replica Ornstein-Zernike integral equations for a quenched-annealed atomic system [70][71][72] were extended to the replica DRISM integral equations for annealed associating molecular liquid sorbed in a quenched matrix [33][34][35][36][37], where ρ 1 γ and ρ 0 γ are the site number densities of liquid species and matrix nanoparticles, h ij αγ (r)' and c ij αγ (r) are the replica total and direct correlation functions between interaction sites α and γ of species i and j (1 for liquid molecules and 0 for matrix nanoparticles).The replica liquid-liquid total and direct correlation functions have the connected (c) and disconnected, or blocking (b) portions, The connected correlations marked with superscript follow from the correlations between the same replica of the liquid, and the blocking ones from those between different replicas of the liquid in the analytical continuation limit s → ∞ .The blocking correlations are a subset of liquid-liquid Mayer diagrams with all paths between two root vortices passing through at least one field vortex of matrix that are completely blocked by matrix vortices (the indirect, matrix-mediated part of the liquid-liquid correlations).The remaining part is the connected liquid-liquid correlations.The KH closure is applied to the matrix-matrix, liquid-matrix, and liquid-liquid correlations, whereas the modified Verlet closure is applied to the blocking part The VM bridge correction (25b) is expressed in terms of the nodal correlation function αγ (r), with the parameter a = 0.8 same as in the original Verlet correction.The blocking correlation closure (25) does not have an interaction potential, as different replicas of fluid do not interact with each other in the limit s → ∞ .This blocking approximation (25) adequately represents nonlinear blocking correlations in polar solvents as well as in electrolyte solutions in nanoporous materials, either neutral or charged.With the analytical treatment of the electrostatic asymptotics is similar to the bulk DRISM-KH Equations ( 3)-( 13), the replica DRISM-VM integral equations are converged using the MDIIS accelerated numerical solver [59], The excess chemical potential of liquid sorbed in disordered nanopores is decomposed into the host matrix contributions from to the liquid-matrix (ij = 10) correlations and the sorbed liquid part from the liquid-liquid (ij = 11) correlations, minus an additional term from the blocking (b) correlations [33,34], Like Equations ( 13)-( 18), the replica DRISM integral Equation ( 23) with the KH closure (24) yields the terms of the excess chemical potential in a closed analytical form of the liquid-matrix and liquid-liquid correlation functions for liquid (j = 1) and matrix (j = 0) species.The blocking correlations term in the chemical potential are obtained from the replica DRISM Equations (23e)-(23g) with the VM closure (25) using thermodynamic integration as ∆µ with the star function derived assuming unique functionality of the correlations [33] s The compressibility and other thermodynamic derivatives of the sorbed annealed solution have the exact analytical forms of the connected terms of the correlations.
An EDL forms in the electrolyte solution sorbed at the nanopores surface even in the absence of external electric charge because of the orientation of polar solvent molecules at the surface of the nanopores as well as the asymmetry of the cationic and anionic density distributions.For interaction site charges q 0 γ of chemical functional groups grafted on the surface of matrix nanoparticles, the statistical charge density distribution τ 0 c of interaction site charges q 1 γ of sorbed solution around a matrix nanoparticle is where matrix nanoparticles and chemical functional groups grafted on the nanoparticles surface are of sort c ∈ 0 and f ∈ 1, respectively.The statistical charge density distribution τ 0 c corresponds to the statistically averaged full local electrostatic potential ψ 0 c (r) around a matrix nanoparticle c, according to the Poisson equation: This local electrostatic potential also includes the term from an externally induced charge on the matrix nanoparticle, The external charge density on the electrode as well as the opposite external charge −q ex on the other electrode of the supercapacitor cause separation of electrolyte cations and anions that diffuse across the electrode separator to the electric double layers in their nanopores.The diffusion occurs until the ionic concentration bias in each electrode satisfies the condition of electroneutrality in the whole system, The diffusion exchange between the electrodes and bulk solution bath adjust the bias between the densities of cations and anions until they satisfy electroneutrality in each electrode.For connected carbon nanospheres of sorts α and sizes R 0 c , the external charge q ex is distributed among charges q 0 c on each sort of nanosphere, provided the electrostatic potential is the same inside all carbon nanospheres, On the other hand, the electrostatic potential far from each nanosphere of sort c gives the average electrode potential with both the carbon nanospheres and sorbed solution parts, The electrostatic potential of the carbon electrode is given by the potential of each carbon sphere φ c (q ex ).The potential change from the carbon nanoparticle to the solution outside the electrode consists of the following parts: (i) voltage across the electric field from φ c (q ex ) at the nanopores' intrinsic surface to the average level φ av (q ex ) of all other carbon nanoparticles in the nanoporous material, and (ii) via the electric field φ exEDL of the external EDL at the macroscopic surface of the electrode to the solution bulk outside the electrode.The "zero" potential level is therefore shifted from the potential in vacuum by the external EDL φ ex EDL .
For sorbed species s with density ρ s , the chemical potential comprises the ideal gas term of spheres with molecular weight m s , ∆µ id s = k B T ln(ρ s Λ s ), where Λ s = h 2 /(2πm s k B T) is the de Broglie thermal wave length, the excess chemical potential ∆µ s for liquid-liquid and liquid-matrix interactions in the nanopores, and the electrostatic potential of species charges q s in the average electrostatic field between the electrodes, These chemical potential terms arise, respectively, from the osmotic effect, the interaction in the nanopores, and the "zero" level of the nanoporous electrode relative to the vacuum potential following from the Nernst equation [73].Inside a charged electrode, the excess chemical potentials of ionic species strongly differ from the bulk solution, which causes diffusion of ions across the separator to have the electric field of ionic dipoles at the electrodes boundaries counterbalance the difference of the "interior" chemical potential µ id s + ∆µ s of electrodes I and II, and to equalize the chemical potential in the electrodes, µ I s = µ II s .The bias between the statistically averaged "zero" levels of the two electrodes φ II av (q ex ) and φ I av (−q ex ) with opposite external charges ±q ex is obtained as The same bias of the electrostatic potential levels φ II av − φ I av must satisfy Equation ( 37) for all solution species.Diffusion of the ionic species to the corresponding electrode with the lower chemical potential occurs due to a difference between the bias values necessary to counterbalance the 'interior" chemical potential µ id s + ∆µ s until the chemical equilibrium (37) is satisfied for both cations and anions.Additionally, neutral solvent molecules diffuse between the two electrodes until their osmotic term and excess chemical potentials satisfy the bias condition (37), which reduces to the equality of their "interior" parts µ id s + ∆µ s .The supercapacitor device voltage obtained by summing up the statistically averaged electrostatic potential changes across the device consists of the potential changes in the electrode I intrinsic EDL φ I av (−q ex ) − φ I c (−q ex ), from the "zero" level of electrode I to the solution bulk and then to the "zero" level of electrode II φ II av (q ex ) − φ I av (−q ex ), and in the electrode I intrinsic EDL φ II av (q ex ) − φ I av (q ex ).With the relation (37) for the average electrostatic potentials in terms of the chemical potentials and number densities of sorbed ions, the supercapacitor voltage assumes the form The chemical equilibrium conditions (37) have to be converged for fluid densities ρ s , with the replica DRISM-KH-VM integral Equations ( 23)-( 25) being converged at each outer loop of iterations for ρ s .This is followed by calculating the potential changes φ I,II av − φ I,II c at converged ρ s by solving the Poisson equation (30).
The purification efficiency relation for the electrosorption cell which holds sorbed electrolyte with ionic concentrations inside electrodes I and II at a much lower electrolyte concentration ρ blk s and excess chemical potential ∆µ blk s in the bulk solution efflux is obtained from the chemical equilibrium condition ∆µ blk s − ∆µ I s (−q ex ) − q s φ I av (−q ex ) for q s > 0 (cations) , (39a) Select applications of the developments in the 3D-RISM-KH molecular solvation theory are described briefly in the following sections.

(Macro)Molecular Simulations with the 3D-RISM-KH Theory
The system size and associated intricacies in the material and biological sciences requires a complex combination of scales for correctly ascribing structure and property combination.The use of quantum mechanics is simply impossible because of the large number of primitive functions needed to construct a wave function for the system.Molecular dynamics simulations, on the other hand, are limited by the timescale.The multiscale approach to handling chemical problems spanning across different size and time domains has been addressed in the literature and specific methods have been developed to address such issues.The applicability of the 3D-RISM-KH molecular solvation theory across multiple such domains have been demonstrated via effectively coupling this theory with electronic structure calculation algorithms and MD-engines (Figure 1) [74,75].
Select applications of the developments in the 3D-RISM-KH molecular solvation theory are described briefly in the following sections.

(Macro)Molecular Simulations with the 3D-RISM-KH Theory
The system size and associated intricacies in the material and biological sciences requires a complex combination of scales for correctly ascribing structure and property combination.The use of quantum mechanics is simply impossible because of the large number of primitive functions needed to construct a wave function for the system.Molecular dynamics simulations, on the other hand, are limited by the timescale.The multiscale approach to handling chemical problems spanning across different size and time domains has been addressed in the literature and specific methods have been developed to address such issues.The applicability of the 3D-RISM-KH molecular solvation theory across multiple such domains have been demonstrated via effectively coupling this theory with electronic structure calculation algorithms and MD-engines (Figure 1) [74,75].Electrolyte solutions are of paramount interest in the fields of energy as well as health sciences.The traditional NaCl-aqueous solution was analyzed using the 3D-RISM-KH molecular solvation theory to understand the solvation structure for a wide range of concentrations; as a benchmark test set for the validity of the theory at hand [20,21].Analysis of the 3D site distributions of water oxygen (O) and hydrogen (H) around a pair of Na + and Cl − ions at infinite dilution presents the solvation structure of the ions correctly with bridging water molecule(s) between the ions embedded in a weak second solvation shell.The features of the H-site of the water molecule show interaction with the chloride ion as well as with the outer solvation shell while positioned away from the sodium ion.This corresponds to water dipole-like water oriented around a cation and hydrogen-bonded to an anion.Water molecules form a bridge of strongly associated water molecules located in a ring between the ions, a situation known as contact ion pair arrangement.At a high salt concentration, many salt bridges form in addition to water hydrogen bonding, and like ion pairs stabilize in both contact ion pair (CIP) and solvent-separated ion pair (SSIP) Electrolyte solutions are of paramount interest in the fields of energy as well as health sciences.The traditional NaCl-aqueous solution was analyzed using the 3D-RISM-KH molecular solvation theory to understand the solvation structure for a wide range of concentrations; as a benchmark test set for the validity of the theory at hand [20,21].Analysis of the 3D site distributions of water oxygen (O) and hydrogen (H) around a pair of Na + and Cl − ions at infinite dilution presents the solvation structure of the ions correctly with bridging water molecule(s) between the ions embedded in a weak second solvation shell.The features of the H-site of the water molecule show interaction with the chloride ion as well as with the outer solvation shell while positioned away from the sodium ion.This corresponds to water dipole-like water oriented around a cation and hydrogen-bonded to an anion.Water molecules form a bridge of strongly associated water molecules located in a ring between the ions, a situation known as contact ion pair arrangement.At a high salt concentration, many salt bridges form in addition to water hydrogen bonding, and like ion pairs stabilize in both contact ion pair (CIP) and solvent-separated ion pair (SSIP) arrangements.In particular, the Cl − -Cl − ion pair at a high concentration in water is bridged by several Na + ions and waters, and forms a cluster.
The replica RISM-KH-VM calculations for KOH aqueous electrolyte solution in a nanoporous carbon supercapacitor that we carried out in this work showed that the solvation structure and thermodynamics of this device are dramatically different from planar electrode-electrolyte interfaces.Figure 2 presents the solvation structure of the aqueous solution of KOH electrolyte in the nanoporous carbon electrode with a positive, neutral, and negative charge.The external charge strongly changes the radial distribution functions (RDFs) between carbon nanospheres and ions, attracting counterions and repelling co-ions.On the other hand, the RDFs between carbon nanoparticles and water oxygen and hydrogen almost do not change with the nanoporous carbon charge.Figure 3 shows the electrostatic potential in the area from the surface of a carbon nanoparticle, through the intrinsic ELD at the nanoparticle surface, and to the bulk value φ C inside the nanoporous electrode averaged over the nanoporous carbon material.It is obtained from the Poisson equation (Equation ( 30)) with the charge density taken from    Figure 4 illustrates a diagram of the distributions of Li + and K + cations as well as OH - anions at the nanoporous carbon electrode surface.The dipole electric field formed by the Figure 3. Statistically averaged electrostatic potential φ 0 (r) around nanoparticles in carbon nanoporous electrode, with respect to the "zero" level φ C .The sorbed solution is in equilibrium with the bulk ambient aqueous solution of KOH electrolyte at concentration 120 ppm.Specific charge of the nanoporous electrode: q ex = 0 (black line); q ex = +16 [C/cm 3 ] (yellow line); q ex = +80 [C/cm 3 ] (red line); q ex = −16 [C/cm 3 ] (green line); q ex = −80 [C/cm 3 ] (blue line).Inset: statistical-mechanical average (red circle and distance vector) around a carbon nanoparticle (red ball) over carbon nanoparticles (black balls) and nanopores (white voids).
Figure 4 illustrates a diagram of the distributions of Li + and K + cations as well as OH − anions at the nanoporous carbon electrode surface.The dipole electric field formed by the solution at the electrode macroscopic boundary due to the chemical equilibrium conditions (37) shifts the bulk potential level q ex φ C in the nanopores.The difference between the "interior" part of the of ions in the nanoporous electrodes k B T ln(ρ s Λ s ) + ∆µ s is counterbalanced by an additional EDL appearing at the macroscopic boundaries of the two electrodes.This brings about a major portion of the supercapacitor voltage U(q ex ) caused by the chemical equilibrium between the solutions inside the nanoporous electrodes and the bulk solution outside them [75].For nanoporous carbon in an ambient aqueous solution of KOH electrolyte at a concentration of 1M, it was found that the electrochemical mechanism of the supercapacitor is driven mostly by the Nernst-Planck equation, determining the chemical balance of sorbed ions in the whole electrode rather than just the EDL potential change at a planar electrode [35][36][37].The purification efficiency of a nanoporous electrosorption cell (39) is determined with the same molecular forces [69].The sorption capacity and specific capacitance are an interplay of the EDL potential change in the Stern layer at the nanoparticles' surface and the Gouy-Chapman layer statistically averaged over the nanopores, the osmotic term coming from the difference of the ionic concentrations in the nanopores and in the bulk electrolyte solution, and the solvation chemical potentials of ions in the nanopores.Chemical specificities of ions, solvent, surface functional groups, and steric effects for sorbed ions strongly affect their solvation chemical potentials in nanopores.Note that the specific capacitance is strongly affected by enlarged effective sizes of sorbed ions, with strong implications on supercapacitor devices.A major factor affecting the specific capacitance is that two extra EDLs at the macroscopic boundaries of the nanoporous electrodes offset the differences of their chemical potentials from the solution bulk.This substantially changes the supercapacitor voltage.chemical potentials in nanopores.Note that the specific capacitance is strongly affected by enlarged effective sizes of sorbed ions, with strong implications on supercapacitor devices.A major factor affecting the specific capacitance is that two extra EDLs at the macroscopic boundaries of the nanoporous electrodes offset the differences of their chemical potentials from the solution bulk.This substantially changes the supercapacitor voltage.Another popular application of the 3D-RISM-KH molecular solvation theory is the calculation of solvation free energy and, hence, molecular partition coefficients.The successful application has been demonstrated in predicting solvation free energies in water, cyclohexane, and n-octanol solvents, using proper benchmark datasets.For all these solvents, a universal correction scheme is adopted for the final prediction of solvation free energy [57,58,76,77].
The formation of nanostructures at the interface of biology and chemistry is aggregationdriven by the so-called hydrophobic interaction.The best example of the 3D-RISM-KH theory addressing such issue is the prediction of such aggregates in the mixture of water with tert-butyl alcohol [78,79].The solvation structure of this chemical system changes with the concentration of tert-butyl alcohol in water, varying from infinite dilution to pure alcohol, as first a single alcohol molecule gets encaged in a tetrahedral hydrogen bonding network of water at infinite dilution, changing subsequently to micromicelles with a size of four to six alcohol units incorporated in the "head-to-head" arrangement in a hydrogen bonding cage of water with an increasing alcohol molar percentage, subtly disturbing the tetrahedral hydrogen bonding network of water, and finally giving away to a zig-zag network of alcohol aggregates with water embedded in it for a very high concentration of the alcohol.The RISM-KH molecular solvation theory predicted both the structure and thermodynamics of these alcohol-water mixtures in full agreement with experiment, providing thermodynamic properties of the systems like compressibility dependences on concentration and temperature as a function of temperature and pressure, the isosbestic point, and the minimum showing formation of alcohol micromicelles.
The macromolecular hierarchical self-assembly in the solution of the rosette nanostructures was analyzed via the 3D-RISM-KH molecular theory of solvation providing a molecular level description of wetting of the nanotube surface and channel by water, which imparts stability to the entire assembly.It predicts the molecular mechanism supramolecular chirality driven by a solvent in helical rosette nanotubes.For instance, structural solvent molecules localized in the pockets behave as molecular switches that result in the formation of a right-hand supramolecular helix in water as a thermodynamically controlled product in the nanotube formation; whereas a left-hand supramolecular helix is formed in methanol because of a kinetic barrier for the right-hand helix when initiating self-assembly in methanol, which is kinetic formation, in agreement with [80,81].
In biological systems, the application of the 3D-RISM-KH molecular solvation theory dates back to the inception of the KH closure.The applications range from protein-ligand binding, binding site mapping, molecular docking, and fast simulation techniques for protein folding.In this theory, the potential of mean force is defined as a 3D map of the effective potential between solvent species.The solute biomolecule determines the binding strength and most probable locations of structural solvent molecules, averaged over the statistical mechanical ensemble of the system [32].This mixture includes ions, cosolvent, and ligand(s), along with the solvent [82].Applying this protocol, one ends up with a novel method of mapping multiple probes on a receptor surfaces useful for protein-ligand binding and for fragment-based drug design.
It was realized in the ligand mapping based on 3D-RISM theory, or 3D-RISM-LM algorithm [39,83,84], and also in the ligand docking algorithm (3D-RISM-Dock) [85] implemented in the AutoDock package.It was validated by experiments with the modes and free energy of the binding of the antiprion compound, GN8, to mouse PrP protein and of thiamine to the extra-cytoplasmic-thiamine-binding lipoprotein, MG289 [86].
The Langevin approach includes artificial friction and random forces to stabilize the solutions systems [87][88][89][90][91].The target temperature is not assured, even though these forces satisfy the fluctuation-dissipation theorem, and a large viscosity coefficient can diverge the system from the real canonical state.Additionally, there is no conserved property in Langevin dynamics, which makes an analysis of the trajectories problematic.This drawback is overcome in the Nosé-Hoover (NH) chain method [92][93][94][95], which reproduces the canonical behavior by including extra variables in phase space that are associated with a thermostat.Its chain counterparts ensure that the ergodicity, and the temperature is controlled without using random numbers.
Quasidynamics is a multi-time-step molecular dynamics (MTS-MD), or optimized isokinetic Nosé-Hoover (OIN) thermostat, of biomolecules steered with 3D-RISM-KH mean solvation forces at long outer time steps (picoseconds) and extrapolated forward using generalized solvation force extrapolation (GSFE) at short inner time steps (femtoseconds) [38].It changed from the reference system propagator algorithm (RESPA) run with the microcanonical ensemble, to the isokinetic Nose-Hoover (NH) chain RESPA (INR) with heat baths coupled individually to each degree of freedom [96].Next, MD/3D-RISM was spread to the canonical ensemble by employing the Langevin (LN) thermostat [97].In the most advanced formulation, the canonical-isokinetic NH chain method consists of grouped degrees of freedom associated with a chain of individual thermostats with a given length and relaxation time.This leads to the OIN method with optimized and increased efficiency of the integration of motion.
Figures 5-7 show the tertiary structures of miniprotein 1L2Y and protein G, simulated by the OIN/3D-RISM-KH/GSFE quasidynamics [38].This quasidynamics keeps the secondary and tertiary structures at huge outer time steps of about 1 ps.In this excellent result with a number of numerical techniques used, the approximate character and the extrapolation of the force fields, and the 3D-RISM-KH integral equations are the main sources of possible uncertainties.Compared to the experiment, the tertiary structures of protein G show that these uncertainties do not affect the ability of the method to reproduce the conformational behavior.There is some difference between the theoretical results and experiment, as they correspond to liquid and crystal states, respectively.For comparison, the folded conformations obtained in explicit solvent conventional MD simulations are very similar to the quasidynamics results (the right parts of Figures 6 and 7).
secondary and tertiary structures at huge outer time steps of about 1 ps.In this e result with a number of numerical techniques used, the approximate character extrapolation of the force fields, and the 3D-RISM-KH integral equations are t sources of possible uncertainties.Compared to the experiment, the tertiary struc protein G show that these uncertainties do not affect the ability of the method to re the conformational behavior.There is some difference between the theoretical res experiment, as they correspond to liquid and crystal states, respectively.For com the folded conformations obtained in explicit solvent conventional MD simulat very similar to the quasidynamics results (the right parts of Figures 6 and 7).protein G show that these uncertainties do not affect the ability of the conformational behavior.There is some difference between th experiment, as they correspond to liquid and crystal states, respe the folded conformations obtained in explicit solvent conventio very similar to the quasidynamics results (the right parts of Figu    The up-to-100-fold time scale compression in the OIN/3D dynamics of protein G in water using 3D-RISM-KH mean solva accelerates the productivity of protein conformational samplin mental dynamics as well as conventional MD with an explicit so dynamics, folding miniprotein 1L2Y from a denatured state took time in explicit-solvent MD was similar to 4−9 μs in experiment.
Accurate calculation of equilibrium and conformational pr in solution can be performed using the stabilizing effect of an O dynamics OIN/3D-RISM-KH/GSFE simulations with gigantic ou The up-to-100-fold time scale compression in the OIN/3D-RISM-KH/GSFE quasidynamics of protein G in water using 3D-RISM-KH mean solvation forces significantly accelerates the productivity of protein conformational sampling compared to experimental dynamics as well as conventional MD with an explicit solvent [38].In our quasidynamics, folding miniprotein 1L2Y from a denatured state took 60 ns, while the folding time in explicit-solvent MD was similar to 4−9 µs in experiment.
Accurate calculation of equilibrium and conformational properties of biomolecules in solution can be performed using the stabilizing effect of an OIN thermostat in quasidynamics OIN/3D-RISM-KH/GSFE simulations with gigantic outer time steps of picoseconds.With this new protocol, one can sample the phase space for rare statistical events, including Thermo 2023, 3 391 the exchange and localization of solvent and ligand molecules in pockets and at the binding sites of biomolecules.This is distinct from explicit solvent MD, which involves huge computational time in such cases.

Conclusions
Nanoscale processes are very different from the macroscopic laws in continuous media, and predictive modeling has to adapt to the needs of coupling different techniques spanning over different time steps.It is therefore necessary to use multiscale methods that couple: (i) electronic structure theories for building blocks, (ii) molecular dynamics simulations for characteristic aggregates, (iii) statistical-mechanical theories for large assemblies of the aggregates and for mean properties over proper size and time scales and macroscopic properties, and finally, (iv) molecular solvation theory to address the issues related to solvation.A powerful platform to predict the solvation structure and thermodynamics of complex chemical and biomolecular systems in real-life solutions is offered by the integral equation theory of molecular liquids based on the first principles of statistical mechanics.An essential part of the multiscale methodology is 3D-RISM-KH molecular solvation theory, which operates, rather than with the trajectories of molecules, with their spatial distributions via analytical summation of the free energy diagrams.This yields the solvation structure and thermodynamics in the statistical-mechanical ensemble.3D-RISM-KH molecular theory of solvation has been coupled in a self-consistent field description with ab initio CASSCF, Kohn-Sham DFT, and orbital-free embedded (OFE) DFT quantum chemistry methods [98][99][100][101][102]. 3D-RISM-KH molecular solvation theory self-consistently coupled with Kohn-Sham DFT in the Amsterdam Density Functional (ADF) package [45,103,104] is also implemented in the Software for Chemistry & Materials suite [105].
For biomolecular simulations, this theory has been coupled with multiple time step molecular dynamics (MTS-MD) steered by 3D-RISM-KH-effective solvation forces, based on a novel algorithm of the optimized isokinetic Nosé-Hoover chain (OIN) thermostat and generalized solvation force extrapolation (GSFE) at short inner time steps.In the energy application, the high specific capacitance of the electric double-layer (EDL) supercapacitor with nanoporous materials has been explored with a new theoretical model, yielding theoretically significant results for aqueous potassium hydroxide electrolyte on a carbon nanoshpere matrix.This new model is based on chemical potentials of solvated ions in nanoporous confinement, bringing about two extra EDLs at the outer boundaries of the nanoporous electrodes to offset the chemical potential difference between the electrodes and solution bulk, which substantially contribute to the supercapacitor voltage.The present status of the 3D-RISM-KH molecular solvation theory has a very broad scope of application from a single molecule in solution to biomolecules in a cellular condition and to nanomaterials in a laboratory setup.Further progress is on the way to increase the scalability of the computational algorithms over various computing architectures.

Figure 2 .
Figure 2. Solvation structure of the aqueous solution of KOH electrolyte sorbed in the nanoporous carbon electrode.Specific charge of the nanoporous electrode: q ex = 0 (solid black lines); q ex = +80 [C/cm 3 ] (long-dashed red lines); q ex = −80 [C/cm 3 ] (short-dashed blue lines); RDFs of water O and H sites, and of K + and OH − ions around carbon nanoparticles.

Figure 2 .
Figure 2. Solvation structure of the aqueous solution of KOH electrolyte sorbed in the nanoporous carbon electrode.Specific charge of the nanoporous electrode: qex = 0 (solid black lines); qex = +80 [C/cm 3 ] (long-dashed red lines); qex = −80 [C/cm 3 ] (short-dashed blue lines); RDFs of water O and H sites, and of K + and OH − ions around carbon nanoparticles.

Figure 4 .
Figure 4. Illustration of Li + and K − cations and OH − anions in aqueous solution in a nanoporous carbon electrode.