Structure and thermodynamics of nondipolar molecular liquids and solutions from integral equation theory

ABSTRACT Solvent-induced solute polarisation of nondipolar solvents originates mainly from specific directional interactions and higher electrostatic multipole moments. Popular continuum solvation models such as the polarisable continuum models ignore such interactions and, therefore, cannot adequately model solvation effects on electronic structure in these environments. Important examples of nondipolar solvents that are indistinguishable by continuum methods are benzene and hexafluorobenzene. Both substances have very similar macroscopic properties, while solutes dissolved in either benzene or hexafluorobenzene behave differently due to their inverted electrostatic quadrupole moments and slightly different size. As a first step towards a proper and computationally feasible description of nondipolar molecular solvents, we present here integral equation theory results based on various forms of the reference interaction site model coupled to quantum-chemical calculations for benzene and hexafluorobenzene solutions of small molecules. We analyse solvation structures, also in comparison with molecular dynamics simulations, and show that predictions of transfer Gibbs energies, which define partition constants, benefit substantially from considering the exact, wave function-derived electrostatic field distribution beyond a simple point charge solute model in comparison with experimental data. Moreover, by constructing artificial uncharged and charge-inverted toy models of the solvents, it is possible to dissect the relative importance of dispersion and quadrupolar electrostatic effects on the partitioning equilibria. Such insight can help to design specifically optimised solvents to control solubility and selectivity for a wide range of applications.


Introduction
Since a large fraction of relevant solvents has dipolar character, thermodynamic solvation effects on solutes can frequently be described well by reducing solute-solvent interactions to the dipolar level, possibly augmented by phenomena attributed to the 'hydrophobic effect' or 'hydrophobic interactions' [1,2] and packing effects for CONTACT Stefan M. Kast stefan.kast@tu-dortmund.de Supplemental data for this article can be accessed at http://dx.doi.org/./... weakly interacting solvents [3]. In a broader context, molecular solvents can exhibit similar molecular packing, but lack dipolar interactions. A good example for such a solvent pair is benzene (C 6 H 6 )/hexafluorobenzene (C 6 F 6 ). Both molecules are highly symmetric (D 6h ) and not easily polarisable, as indicated by their low relative dielectric permittivities (ε r = 2.27 for C 6 H 6 and 2.02 for C 6 F 6 [4]). In the condensed phase and in the crystal structure, both exhibit a characteristic so-called 'herringbone structure' [5,6] that is preferred due to beneficial quadrupolar electrostatic interactions and which contains T-shape like aromatic ring dimers (see Figure 1) [7][8][9]. Furthermore, the molecular quadrupole moments, = (−33.3 ± 2.1) 10 -40 C m 2 for C 6 H 6 [10] and (31.7 ± 1.7) 10 −40 C m 2 for C 6 F 6 [10], have nearly the same magnitude, but opposite signs, favouring a π-π stacking between C 6 H 6 and C 6 F 6 molecules [5]. Yet, there are fundamental differences between both solvents regarding thermodynamic and kinetic properties, such as solvation free energies G solv [11] or rate constants of stereoselective reactions [12,13].
For a theoretical rationalisation of these phenomena, methods are required that reach beyond dielectric effects and incorporate subtle structural solvent structuring differences. Classical force field-based molecular dynamics (MD) simulations [14] treat solvation effects in a detailed atomistic manner, but lack a quantum-chemical description of the solute, which is required, for example, for the computational investigation of chemical reactions. This deficiency is tackled by ab initio MD [15] simulations, however, at much higher computational expense. Unlike ab initio MD, dielectric continuum solvation methods like the polarisable continuum model (PCM) [16] or the conductor-like screening model (COSMO) [17] are computationally cheap in comparison. Several extensions to PCM that account for quadrupolar solvation are available as described by Jeon and Kim [18] that, however, need many experimental parameters and are not easily extended to higher multipoles.
One way to fill the gap between a computationally feasible solvent model, which accounts for a structured solvent and the resulting electronic polarisation of a solute, is to use integral equation theory in combination with a quantum-chemical solute description. A computationally balanced and reasonably accurate way of implementing integral equation theory is the reference interaction site model (RISM) [19] approach. Starting from the onedimensional (1D) level of theory, the so-called 1D RISM method, a solvent model is defined, characterised by its solvent distribution functions and its site-site susceptibilities. Based on these, it is possible to calculate spatial solute-solvent site distribution functions for arbitrarily shaped solutes by solving the three-dimensional form (3D RISM) [20,21]. Our approach to couple solvent and electronic solute structure consists of mapping the solvent charge distribution, which results from the 3D RISM solvent distribution, onto a discrete point charge grid that polarises the electronic solute Hamiltonian. The electrostatic potential (ESP) from the new wave function repolarises and modulates the solvent structure in a subsequent 3D RISM calculation. In this way, by reiterating these steps, a self-consistent procedure is established, which we have termed 'embedded cluster reference interaction site model' (EC-RISM) [22]. Up to now, this approach was mostly applied to and successfully validated for polar solvents like water [22] or dimethyl sulfoxide (DMSO) [23].
Here, we apply the EC-RISM methodology to C 6 H 6 and C 6 F 6 solutions, which allows us to determine the impact of the quantum-mechanical level of theory, the integral equation approach and the solvent model on structure and thermodynamics. Various 1D RISM solvent susceptibility functions, that we will simply call 'solvent models, ' are generated and evaluated by comparing integral equation results to MD simulations. These solvent models are then used to evaluate free energy differences of transferring small molecules between hexafluorobenzene and benzene solutions, more specifically the transfer Gibbs energy in the Ben-Naim standard state [24], The transfer Gibbs energy is connected to the ratio of Henry law constants k H of dissolved species via (ρ is the solvent density and β is the inverse temperature) which is, therefore, an important quantity to estimate and understand the solubility of gases [25]. Note that, for computation of absolute Henry law constants from with standard pressure p 0 , one applies the usual standard state concentration and the corresponding Gibbs energy of solvation by for which G * →0 corresponds to the free enthalpy change from the 1 M infinitely diluted solution to the usual standard state [26] with p 0 being standard pressure and G * vac the vacuum Gibbs energy of the solute. Since we focus solely on transfer Gibbs energies here, the standard state superscripts will be dropped throughout this paper.
In the first part, we focus on pure solvent properties. In contrast to earlier integral equation-theoretical attempts to describe C 6 H 6 and C 6 F 6 by Lowden and Chandler [27] and by Steinhauser et al. [28], this work reaches beyond radial distribution functions by applying 3D RISM theory. We will show that significant structural detail is lost by radial averaging which has substantial influence on the thermodynamics. The connection between RISM solvation patterns and electrostatic quadrupole moments is rationalised by examination of the 'like' pairs (C 6 H 6 in C 6 H 6 and C 6 F 6 in C 6 F 6 ) and 'unlike' pairs (C 6 H 6 in C 6 F 6 and C 6 F 6 in C 6 H 6 ) with 1D, 3D and EC-RISM. In the subsequent part, the analysis of pure solvent properties is followed by predictions of transfer Gibbs energies of small gas molecules between the solvents using various levels of theory and approximations, augmented by an analysis of the reliability of simple point charge models for electrostatic interactions in contrast to explicitly accounting for 'exact' electrostatics originating from the wave function. An additional benefit for understanding predicted numbers is gained from constructing artificial uncharged and reversed-charge (i.e. charge-inverted) toy models of the solvents. In this way, it is possible to dissect the relative importance of dispersion and quadrupolar electrostatic effects on the transfer Gibbs energies and related partitioning equilibria.

RISM
Integral equation theory connects the interaction potential between fluid particles and their spatial distribution in the fluid phase analytically. For the simple case of atomic fluids the Ornstein-Zernike (OZ) equation couples the radial total correlation function h(r) that is related to the radial distribution function g(r) = h(r) + 1, to the radial direct correlation function c(r) with r (ij) being distances. Generalisation to molecular systems in a site-site picture yields the site-site Ornstein-Zernike (SSOZ) or 1D RISM equation [19] h = ω * c * ω + ω * c * ρh, with h = (h αγ ) and c = (c αγ ) now representing total and direct correlation matrices for atomic sites α and γ , and ρ = (ρ γ ) being the site density matrix; the star symbol denotes convolution products. ω = (ω αγ ) is the intramolecular correlation function matrix that consists of the elements for rigid molecular models. The 3D solvent structure around solutes is approximately accessible by 3D RISM theory [20,21]. Here, the 3D spatial solvent distribution function g γ (r) = h γ (r) + 1 of a solvent site γ around a solute molecule on spatial point r can be computed by solving the 3D RISM (solutesolvent) equations A similar form of the solute-solvent equation can be formulated within 1D RISM theory by including another sum over solute sites. Here, the pure solvent susceptibility matrix elements are χ γ γ ' (r) = ρ γ ω γ γ ' (r) + ρ γ h γ γ ' (r)ρ γ ' which can be pre-computed by solving the pure solvent 1D RISM equations. In the following, we will refer to the solvent susceptibilities simply as the 'solvent model. ' OZ-type integral equations need an additional relation, the so-called 'closure' relation that connects the solvent distribution to the solutesolvent interaction potential u. The bridge function B (which only has a formal meaning in site-based theories of molecular liquids) is not easily computationally accessible, which leads to the requirement to apply approximate closure relations such as the hypernetted-chain closure (HNC) [29] where B is simply discarded, or approximations thereof, like the 'partial series expansion' of order k (PSE-k) [30] The PSE-n closures gradually approach the HNC closure with increasing order, bypassing convergence difficulties with the HNC closure. Both closure approximations have the advantage that analytical expressions for the excess chemical potentials μ ex are available for 1D and 3D RISM approaches [30,31].

EC-RISM
In order to combine quantum-chemical calculations for the solute and 3D RISM theory, the solute-solvent interactions is split into the sum of apolar, i.e. for instance Lennard-Jones interactions, and electrostatic terms. For the latter, the ESP of the solute is calculated from the solute's wave function and employed directly, or alternatively a simple point charge representation of the solute's ESP can be determined to reduce computational costs; see [32,33] for procedural details. Background charges at position r i that polarise the solute Hamiltonian are derived from the solvent-charge density (with solvent site charges q γ ) ρ q (r) = γ q γ ρ γ g γ (r) (12) by integrating over grid point volume V i to yield These charges perturb the solute's wave function, change the electronic energy E sol and polarise the molecular ESP. The polarised ESP of the solute is then used to compute an update of the solvent distribution. This cycle is iterated towards self-consistency within the EC-RISM approach [22]. Finally, the corresponding per particle Gibbs energy of the solute in a given conformation is approximated by

... Solvent susceptibilities and D RISM calculations
Standard force fields were employed for the solvent sitesite interactions in the form of a Coulomb potential for electrostatic interactions and a Lennard-Jones 12-6-potential with Lorentz-Berthelot mixing rules for all dispersion/repulsion interactions. Force field parameters for C 6 H 6 and C 6 F 6 , namely the solvent site charge q, the Lennard-Jones potential depth ε and the corresponding contact distance σ were taken from a hybrid Amber/OPLS (Optimized Potentials for Liquid Simulations) force field [34,35]. For comparison, a second force field model for C 6 H 6 developed by Cornell et al. [36] was examined. Note that, consequently, transfer Gibbs energy calculations with the Cornell et al. benzene model have to be used in conjunction with the Amber/OPLS model for hexafluorobenzene. Geometric mixing rules that are common for the OPLS force field were not used to stay consistent with the otherwise applied treatment of mixing non-bonded parameters within 3D RISM calculations described below. Ideal D 6h -symmetric C 6 H 6 and C 6 F 6 structures were employed. For an overview of all solvent parameters and bond distances, see Table 1. Furthermore, two additional sets of solvent models were created, one where the partial charges of C 6 H 6 and C 6 F 6 were set to zero (referred to as 'q 0 -C 6 H 6 ' and 'q 0 -C 6 F 6 ') and another where the partial charges between C 6 H 6 and C 6 F 6 were interchanged, thus, essentially inverting the quadrupole moments of both molecules (referred to as 'q rev -C 6 H 6 ' and 'q rev -C 6 F 6 '). The Lennard-Jones parameters and the molecular structures of the Amber/OPLS models have been used for these two sets of artificial toy models.
The density was set to 873.4 kg/m 3 for benzene and to 1606.3 kg/m 3 for hexafluorobenzene (corresponding to particle densities of 6.785 and 5.218 nm −3 , respectively) [37] and the temperature was set to 298.15 K. The 1D RISM equations for pure solvents were solved with HNC or PSE-k closures, k ranging from 1 to 4 on a logarithmically spaced grid with 512 points analogously to earlier work [38]. The grid ranged from 0.0059Å to a maximum distance of 164.02Å. Iterations were converged to a maximum residuum norm of 1 × 10 −6 for the direct correlation functions. The dielectrically consistent RISM approach [39] was not applicable due to the lack of molecular dipole moments.
The force field models for all MD simulations were chosen corresponding to the 1D RISM calculations. In this context, cubic simulation cells were generated with a volume of 5.3 3 nm 3 containing 1010 benzene molecules and 5.8 3 nm 3 containing 1018 hexafluorobenzene molecules. The Packmol software (version 15.133) [40] was used to create initial simulation cells. Gromacs (version 4.6.3) [41,42] was used for MD equilibration and production runs with a time step of 1 fs. Cut-off radii of 1.2 nm for real space interactions were selected for all calculations. Long-range electrostatics were treated with the smooth particle mesh Ewald technique [43,44]. After equilibration, first within the NVT ensemble (T = 298.15 K) over 2 ns at the experimental densities, followed by NpT simulations (p = 1 bar, T = 298.15 K) over 2 ns production runs were performed in both ensembles over 10 ns. The isothermal-isobaric simulations were conducted to check the force field-specific solvent density while canonical MD was conducted in order to obtain pair distribution functions at the experimental density. Temperature coupling was realised through the Nosé-Hoover thermostat [45,46] and the desired pressure was adjusted with the Parrinello-Rahman coupling scheme [47,48], with relaxation times of 1 and 5 ps, respectively.

... Spatial solvent distributions and transfer Gibbs free energies
Excess chemical potentials and 3D solvent structures were determined from 3D RISM calculations on a cubic grid with 120 3 grid points and a grid spacing of 0.3Å. The 3D RISM equations were solved with the PSE-1 closure iteratively to a maximum residuum norm of 5×10 −6 for the direct correlation functions. EC-RISM calculations were performed with a convergence criterion of 10 −3 kcal mol −1 for the total Gibbs energy of the solute. The transfer Gibbs energies of nine small molecules, namely He, Ne, Ar, N 2 , O 2 , CO, CO 2 , CH 4 and CF 4 , for which experimental values have been reported by Wilhelm and Battino [11], were estimated by computing the difference of results from applying Equation (14) to EC-RISM-derived correlation functions for the small molecules in the respective solvents. A range of different quantum-chemical approaches was examined  Table 2). The investigated levels of theory were limited to B3LYP with 6-311G(d,p) and augcc-pVTZ basis sets for the q 0 models since solute polarisation is not possible for these artificial solvent systems. The solute-solvent interaction was either described by point charges, indicated by EC-RISM q , or by the exact ESP, abbreviated by EC-RISM ϕ . The quantum-chemical part of EC-RISM was performed by the Gaussian 03 program package (Rev. D.02) [64] also using its implementation of the ChelpG partial charge fitting procedure [65]. For all but the final EC-RISM step at the final desired level of theory, HF theory was applied during EC-RISM iterations. Exact ESP were computed from the respective self-consistent wave function before eventual perturbation corrections. Additionally, 1D RISM/PSE-1 solute-solvent calculations were performed exemplarily with the MP2/aug-cc-pVTZ structures and vacuum partial charges, using the same settings as for pure solvent 1D RISM evaluations. The EC-RISM results were compared to [16], to experimental values of Wilhelm and Battino [11] and, furthermore, to 1D and 3D RISM results utilising unpolarised HF vacuum point charges. 3D RISM solvent distribution functions for the 'like' and 'unlike' pairs of C 6 H 6 and C 6 F 6 were obtained with the same procedural parameters as for the small molecules. The solute structures, point charges and Lennard-Jones parameter were chosen in agreement with the 1D RISM solvent models for pure 3D RISM calculations. To obtain a qualitative picture of the quantumchemical polarisation influence on the 'like' and 'unlike' solute-solvent pairs, additional EC-RISM calculations were performed with the small 6-31G(d) basis and exact electrostatic solute-solvent electrostatics.

Pure solvent radial pair distribution functions
All solvent susceptibility functions of C 6 H 6 and C 6 F 6 were derived from pure solvent 1D RISM calculations in combination with all closure approximations applied in this work, PSE-(1-4) and HNC, with the exception of the combination of Amber/OPLS with reversed charges for C 6 F 6 due to convergence issues with HNC. The solvent site-site correlation functions for the whole set of possible force field and closure combinations are depicted in Figures S1-S9 in the Supporting Information. In contrast to earlier results for water [30,66], the pure solvent pair distribution functions of C 6 H 6 and C 6 F 6 are practically closure-insensitive.
In comparison with MD results, the 1D RISM pair distribution functions predominantly exhibit shorter sitesite contact distances as illustrated in Figure 2. This is consistent with 1D RISM calculations on other fluids which stimulated the development of the so-called 'repulsive bridge correction' [67] that broadens short distance peaks at least in the case of water. Beyond the contact distances, the MD-derived structural characteristics of C 6 H 6 are well reproduced by 1D RISM theory, sometimes slightly shifted to shorter distances also for the second and third solvation shells; examples are the two maxima of the carbon-carbon pair distribution function. A comparison of the 1D RISM/AMBER/OPLS results of benzene with the Cornell force field reveals that both models yield nearly identical results with slightly different heights of the first two maxima for g HH and g CC within 1D and MD.
Focusing on g CC and g CF from MD, the hexafluorobenzene structure has similar features compared to the corresponding benzene distributions although they are not as well reproduced by 1D RISM as in the benzene case. The major difference between the radial solvation structures is exhibited by g FF and g HH . While g FF shows a maximum at around 3Å, there is none for g HH .
The simulated solvent densities within the NpT ensemble show that both force fields underestimate the density by 6.9%-12.4% as reported in Table 3. Here, Amber/OPLS performs slightly better than Cornell. It is conceivable that the pair distribution functions from NpT simulations are slightly higher with a force field model that reflects the densities more appropriately, as is the case for the NVT MD results.
The underlying cause behind the difference of g HH and g FF seems not to be the difference in quadrupole moments but to be the difference in the dispersion interactions. As shown in Figure 3, where the results without electrostatic moments (q 0 ) and reversed quadrupole moments (q rev ) are illustrated, the drastic change in the electrostatics between both solvents leads to nearly the same 1D RISM results for radial solvation patterns. Note that especially the exchange of partial charges has only minor impact on pair distribution functions, while notably the results with the uncharged model are slightly flatter.
These distribution functions can be compared to the radially averaged results from 3D RISM calculations on the solute-solvent pairs C 6 H 6 /C 6 H 6 and C 6 F 6 /C 6 F 6 by integrating over the angular coordinates using Lebedev-Laikov grids [68]. The comparison is illustrated in Figure 4. These radial functions are in much better agreement with the simulated results for short distances than the 1D RISM pair distributions. Furthermore, the shifts of the maxima observed for 1D RISM are corrected by 3D RISM, matching the MD-derived distance distributions better. However, the heights of the radially averaged 3D RISM distribution functions are much lower, probably due to the use of PSE-1 as 3D RISM closure. In conclusion, the radial solvation structure is appropriately represented by 3D RISM, but the peak height is underestimated. Nonetheless, particularly the better behaviour for small distances not captured by a priori radially averaged 1D RISM theory has great impact on the solvation free energies [67], see also the following.
Regarding the 3D solvent distributions of the like and unlike pairs illustrated in Figures 5 and 6, a fundamental difference in the solvation pattern is obtained that is not reproduced by 1D RISM theory. The 3D RISM results with the AMBER/OPLS benzene and hexafluorobenzene models are reasonable in the context of matching quadrupole interactions between both ring systems and indicate π-stacking for the unlike pairs C 6 H 6 /C 6 F 6 and C 6 F 6 /C 6 H 6 , and T-shaped solvation structures for the like pairs C 6 H 6 /C 6 H 6 and C 6 F 6 /C 6 F 6 . Our results exhibit similarities with the crystal structure packing of C 6 H 6 and C 6 F 6 [5] and are in very good agreement with the 'reconstructed 3D solvent structures estimated by the empirical potential structural refinement' analysis technique [69,70] carried out by Headen et al. [71]. An inversion of the solvent molecule's quadrupole moment yields an inverted 3D solvation pattern, as shown in the first and third columns of Figure 6. Furthermore, turning off the partial charges of the solvent species totally disrupts the T-shape like solvation patterns. After reduction of the spatial distribution functions of the results from q 0 and q rev models to radial solvent distribution functions, a levelling effect on the solvation distribution is observed, as illustrated in Figure 4. While not as strong as the pure 1D RISM results but still noticeable, this means that the 3D RISM radial functions of q rev and the original solvent models are very similar to each other, despite showing nearly opposing features in three dimensions. Therefore, the C 6 H 6 /C 6 H 6 solvent distribution of the regular solvent model matches the distributions of the q rev results for the pair C 6 F 6 /C 6 H 6 and so forth.
As illustrated in Figure 5, there exist just minor differences in the solvation patterns between the pure 3D RISM and the EC-RISM ϕ results, although the solute's exact ESP is applied and the empirical Amber/OPLS force field partial charges are used for the pure 3D RISM evaluations. The largest difference between the generic 3D RISM and EC-RISM data is exhibited in the C 6 F 6 /C 6 F 6 case. Nonetheless, the overall solvation patterns are very similar here. The consequences of the radial averaging and the solvent model on Gibbs energies of small molecules are presented and discussed in the next section.

Gibbs energies and solvation structures of small molecules
We now report and discuss transfer Gibbs energies of small molecules for all C 6 H 6 /C 6 F 6 solvent models determined from EC-RISM calculations with exact and point charge-approximated electrostatics, EC-RISM ϕ and EC-RISM q , in comparison with 1D/3D RISM results obtained with vacuum charges, i.e. without accounting for electronic polarisation. The results are summarised in Tables S1-S19 of the Supporting Information. Figure 7 distinguishes between different solvent models and solute-solvent interaction approaches averaged over all quantum-chemical levels of theories. The influence of the quantum-chemical part will be discussed below.  All integral equation methodologies clearly differentiate between C 6 H 6 and C 6 F 6 solvation for the process X(C 6 H 6 ) → X(C 6 F 6 ), where X represents a gas species, in the sense that the RISM transfer Gibbs energies show the same trend as the experimental data and exhibit negative values. Notably, this is not the case for PCM (see Table 4) which ignores the directionality of the interactions accompanying the transfer process, leading to small positive transfer Gibbs energies for every level of theory and species. Therefore, we skip the PCM model results in the figures to focus more on integral equation approaches. The EC-RISM results with the original solvent models exhibit a difference to the experiments below 1.2 kJ mol −1 , except for CO 2 . Especially the results for argon, oxygen and carbon monoxide are very close to experimental results. In the cases of the uncharged (q 0 ) and the reversed-charge (q rev ) models, nearly all values are by far not negative enough and, therefore, underestimate the transfer Gibbs energies. The sole exception is CO 2 , for which G trans with the uncharged solvent model has the smallest deviation from experiments. All realistic solvent models predict G trans values lower than −3.9 kJ mol −1 for CO 2 , which is off by 1.9 kJ mol −1 or more. 1D RISM calculations perform quite well for the three noble gases as well as for nitrogen, but are quite bad for the other compounds compared to the 3D and EC-RISM results. Figure 8 (see also Figure S10 in the Supporting Information) illustrates the root mean squared deviation (RMSD) between experiment and results for EC-RISM with a certain quantum-chemical level of theory as well as for 3D RISM. It is evident that the basis set and the level of theory have no systematic influence on the pure 3D RISM results. In contrast, the EC-RISM q and even more so the EC-RISM ϕ results are improved considerably with a larger basis. The results for the q rev solvent models also improve, but not as stringent as the realistic models. The smallest deviation is achieved by MP2/augcc-pVTZ in combination with EC-RISM ϕ , which is just slightly better than MP2/aug-cc-pVQZ/EC-RISM ϕ . Furthermore, the influence of the basis set is more pronounced than the influence of the ab initio or density functional theory (DFT) level. Figure 8 illustrates that the results are more sensitive to a change of the solvation model than to the quantumchemical level of theory. Comparing the overall mean signed errors (MSEs) and RMSDs from experimental values as shown in Table 5 and Figure 9, it is obvious that the combination of the Amber/OPLS force field with the EC-RISM methodology and exact ESP is the best choice for transfer Gibbs energy calculations. The Cornell model has the smallest MSE, but a much larger RMSD, indicating larger scatter and, therefore, less reliability. Overall, the results are improved by EC-RISM in comparison with generic 3D RISM and are generally superior with exact electrostatics in comparison with ChelpG partial charges. Surprisingly, the second best RMSD result is obtained with the uncharged (q 0 ) EC-RISM q Amber/OPLS model calculations. However, here the G trans values are always slightly larger than experimental values, resulting in a larger MSE which, therefore, indicates an unbalanced, systematically shifted methodology. The 1D RISM results deviate nearly as much as the artificial q rev model, but with a much better MSE. The q rev model results exhibit the largest difference to experimental values and the largest MSE between all RISM approaches.
The uncharged q 0 solvent model results represent the thermodynamic effect of solutes immersed in a purely dispersive/repulsive solvent. Therefore, the results in Figure 9 and Table 5 demonstrate that dispersive effects are important for the transfer Gibbs energy of small Table . EC-RISM-derived transfer Gibbs energies for the process X(C  H  ) → X(C  F  ) for small molecules resulting from averaging over HF, PBE, BLYP and MP calculations combined with different Pople and Dunning basis sets on respective PCM structures in comparison with experimental data, and averaged PCM and D (MP/aug-cc-pVTZ geometry/vacuum charges)/D RISM (all geometries/vacuum charges) results.
molecules, but also show that a realistic representation of the electrostatics improves these results that are in general not negative enough in their absence. In contrast, in the reversed-charge q rev case, the inclusion of electrostatic solute-solvent interactions worsens the results by further increasing G trans values. To quantify the influence of dispersive solute-solvent interactions, we define the fraction of the purely dispersive interaction of the transfer Gibbs energy as with M representing any solvation method other than EC-RISM q with the q 0 solvation model. Due to small   field, reveals that approximately 87%-89% of the transfer free energy is due to dispersion (see Table 6). The fraction x disp in this table is evaluated and averaged over the same levels of theory that are used in combination with the q 0 solvent model in order to make a fair comparison. The transfer Gibbs energies of CO 2 and CF 4 show the largest difference to those from the q 0 model and, thereby, with 63%-64% and 79%-82%, the smallest amount of dispersive influence. There is nearly no difference between the 3D RISM and EC-RISM q results with B3LYP/6-311G(d,p) and B3LYP/aug-cc-pVTZ. In addition, the EC-RISM results for the dispersive fraction with exact electrostatics and with ChelpG point charges are very similar, except for N 2 , O 2 and CO, where x disp is lower for the EC-RISM ϕ results.
To explain the influence of the electrostatic representation on the transfer Gibbs energies and in particular on the differences between quadrupolar and dispersive solvation effects, it is instructive to focus on the N 2 and CO results. Figure 10 shows that the EC-RISM q solvation structure around N 2 is dominated by two solvent density 'rings. ' A comparison of these densities with the solvent structure around N 2 from the q 0 solvation model reveals that these rings are already obtained with purely dispersive solvation. Figure 10(m) and 10(n) illustrate the ESP originating from the solute that is experienced by the solvent in the corresponding EC-RISM calculations. Within the partial solute site charge approach EC-RISM q , the solvent does not experience any electrostatic solutesolvent interaction at all for N 2 due to symmetry, while, in contrast, the solute interacts through higher electrostatic moments with the solvent within EC-RISM ϕ , i.e. with exact electrostatics. Consequently, the solvent structure with exact electrostatic solute-solvent interactions differs drastically from the EC-RISM q and EC-RISM (q 0 ) results. Thus, N 2 and O 2 can be polarised by EC-RISM ϕ , although both molecules have atom-centred partial charges of zero. Only higher-order multipolar interactions can cause polarisation for N 2 and O 2 , which is not possible with regular 3D RISM and EC-RISM q models. These findings are supported by calculations of Mohan et al. [72] for the C 6 F 6 /N 2 dimer. Another example for which the x disp values are dissimilar between EC-RISM q and EC-RISM ϕ is CO, although it is a dipolar molecule and multipolar influences appear not to be substantial. But again, as for N 2 , the EC-RISM ϕ solvent site densities around CO are clearly different from those computed by EC-RISM q , which is in fact again explicable by the complex molecular ESP of CO, see Figure 10(o) and 10(p).
To further differentiate between energetic contributions and to investigate polarisation effects, we split the and are the unpolarised contributions of the electronic energy (Equation (17)) and the excess chemical potential (Equation (18)) of the solute that corresponds to the 3D  RISM results with vacuum charges. The other two contributions that describe the influence of self-consistent polarisation are described by E pol,trans = E sol (C 6 F 6 ) − E sol (C 6 H 6 ) − E 0,trans μ ex pol,trans = μ ex (C 6 F 6 ) − μ ex (C 6 H 6 ) − μ 0,trans .
The estimated contributions, averaged over all quantum-mechanical levels of theory, are illustrated in Figure 11 for the EC-RISM q results with the AMBER/OPLS solvent model. The calculated transfer Gibbs energy is dominated by the unpolarised excess chemical potential. The changes of E or μ ex that are caused by polarisation are timid, and play a minor role for the G trans predictions.

Concluding remarks
Benzene and hexafluorobenzene solvent models were developed and assessed in the context of liquid structure analyses and of transfer free energy predictions by an integral equation approach in order to demonstrate the relevance of higher multipolar electrostatic effects and the methodical implementation to account for them adequately. In terms of structure, both solvent models yield reasonable radial pair distribution functions and spatial solvent distributions, as is evident by comparing MD and 3D RISM data. Conversely, radial averaging tends to level important details for distinguishing between these solvents. In terms of thermodynamics, transfer Gibbs energies can be predicted for small molecules with a total RMSD of about 0.9-1.2 kJ mol −1 in good agreement with experiments and better than with other theoretical estimates like, e.g. scaled particle theory [11,73]. This quality is only possible since both electronic polarisation and electric multipoles beyond the dipolar approximation are properly accounted for.
More physical insight could be gained from investigating uncharged and charge-inverted solvent models by separating dispersion and electrostatic contributions and by quantifying their importance for the calculated transfer Gibbs energies. In this context, it was observed that dispersion interactions already capture a large amount of the Gibbs energy differences (and therefore of Henry law constants) of C 6 H 6 and C 6 F 6 solutions. In particular, the q 0 models account for roughly 90% of the solvation discrimination measured by transfer Gibbs energies and already show the correct partitioning tendency. Adding electrostatic interactions only amplifies the effect, while it is diminished for reversed-charge models. Thus, electrostatic interactions turn out to be responsible for a differential solvation effect in the sense that they act as an additional factor over the baseline solvent discrimination resulting from predominantly dispersive interactions which are of course coupled to packing effects within the Lennard-Jones model potential.
Consequently, the use of exact solute-solvent electrostatics, which yields a more realistic solvent description and slightly but systematically better thermodynamic predictions than point charge models, is emphasised and supported by more reasonable solvation patterns that are not properly represented by the latter. The combination of an exact electrostatic representation for the solutesolvent interactions allows to incorporate higher multipole moments in a simple way and to treat phenomena like sigma hole effects [74] between solute and solvent molecules without specifically parameterised force fields. As a perspective, the methodology presented in this work has the potential to allow for a rational design of solvent properties to modulate and control chemistry in the solution state.