Effect of the Solute Cavity on the Solvation Energy and its Derivatives within the Framework of the Gaussian Charge Scheme

The treatment of the solvation charges using Gaussian functions in the polarizable continuum model results in a smooth potential energy surface. These charges are placed on top of the surface of the solute cavity. In this article, we study the effect of the solute cavity (van der Waals‐type or solvent‐excluded surface‐type) using the Gaussian charge scheme within the framework of the conductor‐like polarizable continuum model on (a) the accuracy and computational cost of the self‐consistent field (SCF) energy and its gradient and on (b) the calculation of free energies of solvation. For that purpose, we have considered a large set of systems ranging from few atoms to more than 200 atoms in different solvents. Our results at the DFT level using the B3LYP functional and the def2‐TZVP basis set show that the choice of the solute cavity does neither affect the accuracy nor the cost of calculations for small systems (< 100 atoms). For larger systems, the use of a vdW‐type cavity is recommended, as it prevents small oscillations in the gradient (present when using a SES‐type cavity), which affect the convergence of the SCF energy gradient. Regarding the free energies of solvation, we consider a solvent‐dependent probe sphere to construct the solvent‐accessible surface area required to calculate the nonelectrostatic contribution to the free energy of solvation. For this part, our results for a large set of organic molecules in different solvents agree with available experimental data with an accuracy lower than 1 kcal/mol for both polar and nonpolar solvents.

over the surface of the solute cavity. There exist two main types of solute cavities used to treat the contribution of solvation to the electrostatics of a system: (a) the vdW-type cavity, (b) the solventexcluded surface (SES)-type cavity. [10] The first kind of cavity is the one contained inside the external surface of the vdW spheres centered on each atom (see Figure 1b). The latter, on the other hand, is obtained by following the inward part of a solvent sphere rolling over the solute (see Figure 1c), which makes its outer surface smoother than the vdW surface. The SES-type cavity is, in most of the cases, generated through the GEPOL (GEnerating POLyhedra) algorithm. [11][12][13] Other strategies to generate the SES are based on an implicit function derived from equivalence statements and Voronoitype diagrams. [14] The use of solvation point charges leads, however, to discontinuities in the Coulomb potential that affect the selfconsistent field (SCF) energy and its derivatives with respect to nuclear displacements. This is firstly due to the fact that the interaction between surface elements depends on r −1 ij , r ij being the distance between charges i and j. Then, when r ij tends to zero (solvation charges overlap), this directly affects the SCF energy, as it is computed from the solvation charges. A very efficient way to overcome this problem is to represent the effect of the solvent in our system by a charge density made of spherical Gaussians centered at the point charge coordinates. This strategy was first proposed by York and Karplus [15] and further refined by Scalmani and Frisch [16] as well as Lange and Herbert. [17] In all cases, together with the Gaussian treatment of the solvation charges, a switching function is used to attenuate the contribution of the charges to the solvation properties. This last point is crucial in order to avoid discontinuities in the solvation potential due to the charges that may appear or disappear when optimizing the geometry of a molecule. This whole strategy, that from now we call the Gaussian charge scheme (GCS), results in a rigorously smooth potential energy surface possessing continuous derivatives.
The use of Gaussian charges together with the use of a switching function permits to model the contribution of solvation to the electrostatic energy of the system. However, a rigorous treatment of solvation should also include nonelectrostatic contributions to the energy. There exist different ways to consider the nonelectrostatic solvation energy, ΔG nel . A simple way to do so is to assume this contribution to linearly depend on the solvent-accessible surface area (SASA) and on the contribution of each atom to the SASA. [18,19] The SASA is obtained by following the center of a solvent probe sphere rolling over the vdW surface (see Figure 1d). [10,20] A well-known continuum solvation model that involves an accurate description of ΔG nel is the Solvation Model based on Density (SMD). [21] In this case, ΔG nel is modeled as a sum of contributions that are proportional to the SASA of the atoms forming the solute via a set of geometrydependent proportionality constants called atomic surface tensions.
In this model, the SASA is generated by rolling a solvent sphere with a fixed radius R S of 0.4 Å over the vdW surface of the solute. In principle, the fact of not considering R S to be solvent-dependent is somehow hidden in the careful parametrization of the atomic surface tensions. Regarding this last point, there is a lack of studies in the literature that consider a different value for R S depending on the solvent of interest. Although the main goal of our study is to report the effect of the cavity used for electrostatics on the properties of a system using the GCS, we also want to prove that the use of a solventdependent R S to construct the SAS leads to the accurate prediction of solvation energies. The present work is divided into two parts. In the first part, we study the effect of the solute cavity used to compute the solvation contribution to electrostatics within the C-PCM using a GCS for the polarization charges. For that purpose, the GCS has been implemented into the ORCA quantum chemistry suite. [22,23] We focus on the results for the first derivatives of the SCF energy for small and big systems (up to 300 atoms). In order to verify the correctness of our implementation we compare our results using the GCS with those for systems where the point charge scheme results in singularities in the SCF energy and its derivatives. In the second part of this study, we have parametrized the nonelectrostatic contribution to the free energy of solvation in order to predict experimental data for this last quantity for organic solutes in different solvents. Our aim in this last part is to show that a strategy that considers a solventdependent R S for nonelectrostatic solvation effects together with a Gaussian treatment of the solvation charges leads to accurate results for solvation data.
This article is organized in the following way. In the theory section, we present the basics of the C-PCM and the GCS, we discuss the different types of solute cavities and we introduce the equations to compute the free energy of solvation. The computational details of our calculations, together with the procedure used to parametrize the free energy of solvation are provided afterward. In the results section, we show the outcome regarding the efficiency and computational cost associated to the GCS, the performance of this scheme for large systems, and the results for the parametrization of the solvation free energy. Finally, the last section is devoted to the main conclusions that can be drawn from this study.
F I G U R E 1 The different type of cavity surfaces for a given solute (a), van der Waals surface (b), solvent-excluded surface (c), and solvent-accessible surface (d). The solvent molecule is represented as a dashed sphere labeled as "S" 2 | THEORY 2.1 | C-PCM equations For a solvated system within the PCM, the Hamiltonian of the solute in vacuum, H 0 , is perturbed by a solvation term that depends on the N q charges, q, spread out over the surface of the solute cavity, with r i the position of the ith point charge.
Within the C-PCM, the charge vector q is calculated by solving the following system of equations: where V stands for the electrostatic potential at the position of the charges, A is the LHS matrix, and f(ε) is a scaling function that depends on the dielectric permittivity, ε, of the solvent The parameter x is equal to 0 in the case of the C-PCM. The elements of the matrix A have a different functional form depending on how the charges on the surface of the cavity are treated, as we discuss later on.
The SCF energy for the solvated system reads as where P μνσ = P i c * μiσ c νiσ is the solute density matrix for spin σ, η = α, β, with c μiσ being the molecular orbital (MO) coefficients. We adopt the Mulliken notation for the two-electron integrals over the atom-centered basis functions {μ}. The term V NN accounts for the nuclear repulsion energy while the fourth term describes the interaction between the solute and the solvent. Finally, E XC ρ σ r ð Þ ½ is the exchange-correlation energy with c x equal to the fraction of Hartree-Fock exchange and ρ σ r ð Þ = P μν P μν μ r ð Þν r ð Þ being the electron density for spin σ at position r.
The first derivative of E with respect to an infinitesimal nuclear displacement X, is equal to The matrices S and W are the overlap and energy-weighted density matrices, respectively. If the superscript X is enclosed in parenthesis, then the derivative does not involve the calculation of ∂P/∂X. Notice that, by definition, W = PFP with F being the Fock matrix. When we adopt the C-PCM approach, a solvation contribution is also added to the Fock matrix elements over basis functions, F μνσ Þdr is the exchange-correlation potential and V μνσ is the potential due to shells χ μ and χ ν at the position of the charges.

| Discontinuities in the SCF energy
The left-hand side (LHS) of Equation (2) (5)) that diverge faster (A X depends on r −3 ij ). Then, even if we have a small system, when optimizing it, it is very likely that some of the charges get very close to each other causing numerical instabilities in the SCF energy that will lead to erroneous optimized geometries. If one calculates vibrational frequencies, this systematic error is accumulated again, due to two main reasons. The first one is associated to the fact that we are using an optimized geometry that can be biased. Second, the Hessian, if calculated analytically, does not only depend on r −1 ij and r −3 ij but also on r −5 ij . If we compute it numerically through finite differences, we have to generate the cavity as many times as we compute the gradient. However, although using Gaussian charges produces finite solvation contributions to the SCF energy, it does not completely remove the discontinuity in the derivatives of the solvation potential with respect to nuclear displacements. This issue arises from those charges changing discontinuously in the intersection between cavity spheres from one molecular configuration to the other during geometry optimizations. This problem, which is discussed in detail in Reference [15], is solved by the use of a switching function that continuously turns on and off the contributions of charges to the energy, its gradient and further derivatives.

| Gaussian Charge Scheme
Within the GCS, the effect of the solvent is represented by a charge density made of spherical Gaussian functions centered at positions r i on the surface of the solute cavity [15] Here, q i is the magnitude of the charge located at r i and ζ i is the width of the Gaussian. If we consider that the surface of each of the spheres that compose the solute cavity is discretized following a Lebedev quadrature scheme, then each charge i with width ζ i is represented by a position r i , its Lebedev weight w i and the radius R I of the sphere on which the charge is placed. In this case, ζ i is computed as where the exponent ζ is a dimensionless parameter adjusted to obtain the exact Born solvation energy for a conductor and a uniform surface charge distribution. This parameter depends on the Lebedev grid chosen. [15] The elements of the LHS matrix A in Equation (2) have the following functional form under the GCS.
. Regarding the diagonal matrix elements, A ii , they are calculated from the so-called switching function F i . This function weights the contribution of each charge to the solvation properties and can be calculated as a product of elementary switching where N sph is the number of spheres that define the solute cavity and R J is the position vector of sphere J. We adopt the improved switching-Gaussian (ISWIG) function proposed in Reference [17].
with R J being the radius of sphere J, and r iJ the distance between charge i and the Jth sphere center.
It is especially interesting to analyze the functional form of the nondiagonal elements of the LHS A matrix in Equation (10). Differently than for the point charge scheme, where A ij = r − 1 ij , we have this quantity weighted by erf(ζ ij r ij ). Then, That is, A ij does not diverge when two charges overlap, as it is the case for the point charge scheme. This fact prevents discontinuities in the SCF energy. Indeed, with the choice of A ii , A ij and F i we have a strictly continuous potential energy surface.
The contribution to the surface area of the solute cavity by the charge i, a i , is calculated as Then, the total surface area of the chosen cavity is just The total volume of the cavity, V T , is In order to perform geometry optimizations within the C-PCM using the GCS we need the derivatives of A ii and A ij with respect to nuclear displacements. This involves calculating the gradient of F i . If we consider the nuclear coordinates of sphere A, On the other hand, The terms that depend on ∂ζ i ∂RA and ∂ζ ij ∂RA are equal to 0 for a vdWtype cavity because they are proportional to the derivatives of the radius of the cavity spheres with respect to nuclear displacements. In the case of the SES, the equations used in the present study to compute these terms are provided in Data S1. The quantities ∂ri ∂RA in Equations (17) and (18) are calculated slightly different depending on the solute cavity chosen, as we discuss in the following subsection.

| Type of cavities used for the solvation contribution to electrostatics
We consider two types of solute cavities to be used together with the GCS to compute the solvation contribution to electrostatics, (a) a vdW-type cavity and (b) a SES-type cavity.

| VdW-type cavity
This type of cavity is generated by scaling the vdW radius of each atom I, R vdW,I , by a factor f, that is, R I = fR vdW,I and taking the external surface of these spheres. We use Bondi's vdW radii for those elements where this quantity is available, [24] except for hydrogen where we adopt a value of 1.1 Å, as suggested in Reference [25]. For 16 of the main-group elements in the periodic table for which Bondi's radii are not defined, we adopt the radii proposed by Mantina et al. [26] This is the case of the following elements: Be, B, Al, Ca, Ge, Rb, Sr, Sb, Cs, Ba, Bi, Po, At, Rn, Fr, and Ra. For those elements that are covered neither by Bondi nor by Reference [26], we consider a value of 2 Å.
Although previous studies have shown that the use of a solventdependent f improves the electrostatic interaction between solute and solvent, [27][28][29][30] in our approach, f = 1.2, as done in other studies. [17] The Gaussian charges are finally spread over the scaled vdW surface using a Lebedev grid as pointed out in the previous subsection.

| SES-type cavity
The second type of solute cavity is the one that has the SES as its external surface. In this case, we use the GEPOL algorithm to generate it but not to discretize its surface. The GEPOL algorithm iteratively considers all possible pairs of existing spheres and checks if there is space in between them that should be inaccessible to the solvent. If it is the case, a new sphere (N) is added in between the pair of spheres.
There exist three different situations depending on the distance between the pair of spheres, their radii and the radius of the solvent probe sphere. These different cases, A, B, and C, are shown in Figure 2. The equations used to compute the position vector of the new spheres, r N , are provided in Data S1. Regarding the radii of the spheres for the SES-type cavity, we adopt the same strategy that for the vdW-type cavity, which is scaling the vdW radius of each atom by a scaling factor f, with f being equal to 1.2. The radii of the new GEPOL spheres (those that are not associated to atoms) are calculated from the radii of the pair of spheres considered to create it, as well as from the distance between them and from the radius of the solvent molecule R S . In order to determine R S for each solvent, we start from the Stokes-Einstein formula, [31] that gives the radius, R SE , of a sphere that diffuses at the same rate as the solvent under consideration, where D is the self-diffusion coefficient of the solvent, η its viscosity, k B the Boltzmann's constant, and T is the temperature. In the case of water at ambient conditions (T = 298.15 K), D = 2.26 × 10 −9 m 2 /s, η = 8.9 cP, which results in R w SE = 1:09 Å . This value is slightly lower than the experimentally determined radius for a water molecule, which is R w S = 1.38 Å. In the present study, we adopt the later value of R S for water. If we consider a different solvent than water, we calculate the radius of the solvent molecules with the following equation, where R w S =R w SE ≈ 1:27. This prefactor slightly corrects for the fact that R SE underestimates the value of R S . The values of R S for the solvents F I G U R E 2 The different pairs of spheres considered within the GEPOL algorithm to generate the SES: case A (a), case B (b), and case C (c). The solvent molecule is represented as a dashed sphere labeled as "S" considered in this study are shown in Table 1. Although an estimation of R S based on R SE can seem a crude approximation, this strategy permits us to take into account the size of the solvent molecule, something that is ignored in most of the implicit solvation models.
Once the GEPOL spheres are generated (position and radius) we spread the Gaussian charges over the SES using a Lebedev grid. A quantity that requires special attention are the derivatives of the charge positions, r i , with respect to nuclear displacements. Let us consider, for instance, the first derivative of r i with respect to the infinitesimal displacement of atom A along the x direction, that is ∂r i /∂x A . If the charge i belongs to atom A, then ∂r i /∂x A = (1, 0, 0). If it is not the case, and i belongs to a GEPOL sphere N (not associated to any atom) then we should check if the atom A is one of the two spheres that was used to create N. If it is the case, due to the fact that r i depends on the coordinates of atom A, R A , then ∂r i /∂x A 6 ¼ 0. The equations used to compute the gradient for cases A, B, and C in Figure 2 are provided in Data S1.

| Free energy of solvation
The solvation free energy, ΔG solv , is the free energy of transfer of a solute from the gas phase to the condensed phase. If we denote by R v and R l the equilibrium coordinates of a solute M in the gas phase and in solution, respectively, the fixed-concentration free energy of solvation of M can be written as, [32] The first term, E solv (R v , R l ), accounts for the difference between the liquid-phase expectation value of the gas-phase Hamiltonian, E(R l ), and the gas-phase potential energy surface E(R v ) The second contribution to Equation (21), ΔG c (R v ), represents the changes in the free energy due to the coupling between the solute and the solvent. This term can be approximated by with ΔG el and ΔG nel being the electronic-polarization and nonelectrostatic contributions to ΔG solv , respectively. Regarding ΔG nel , it can be split into two terms: (a) a cavitation energy, ΔG cav , defined as the reversible work to create a cavity inside the solvent to accommodate the solute, and (b) a term, ΔG disp , that accounts for the changes in the dispersion energy when solvating the solute.
With respect to the term G solv,lib − tr in Equation (21), it accounts for the difference between the liberational energy of the system in solution, G s,lib , and the translational free energy of the solute in the gas phase, G v,tr , If we choose standard states in which the concentration of the solute in the gas-phase is the same as that in solution, G s,lib = G v,tr and then G solv,lib − tr = 0. Finally, the term G solv,rv is defined as the difference between the rotational-vibrational free energy of the solvated system, G s,rv , and that for the system in vacuum, G v,rv .
In this case, it is also customary to consider G solv,rv = 0. This approximation is based on the fact that the contribution of G solv,rv to ΔG solv is in general small with respect to the other terms in Equation (21) and expected to be smaller than the mean intrinsic error of the continuum solvation model.
If we then consider Equations (25) and (26) to be equal to zero, and combine Equations (22), (23), and (24) with Equation (21), then ΔG solv reads as, The term ΔG el is equal to the fourth term in Equation (4), and is calculated from the charges spread over the surface of the vdW-type cavity or the SES. Regarding the last two terms in Equation (27), they are calculated using the SASA (see Figure 1d). The radius of the solvent is calculated using the approach described in the previous subsection. For the cavitation energy ΔG cav , we adopt the equation proposed by Pierotti, [33] based on scaled particle theory for a hardsphere solute I in a hard-sphere fluid of solvent molecules. In this case, ΔG cav,I is expanded in powers of the ratio R between the radius of the solute sphere R I and that of the solvent spheres R S .
with y = 8πρ S R 3 S =6, being ρ S the density of the solvent, and R = R I =R S . The last term in Equation (28) depends on the thermodynamic conditions, that is, temperature and pressure. Equation (28)  spheres. [34] In this case, the contribution to the total ΔG cav by sphere I is weighted by a prefactor that depends on the exposed SASA I of that sphere.
Here, the summation over the total number of spheres can be replaced by a summation over the total number of surface elements (N el ) in which the SASA is divided into. This number is, in general, larger than the number of charges placed on top of either the vdW and the SES surfaces. Finally, I stands for the index of the sphere to which the surface element i belongs to.
Regarding the last term in Equation (27), ΔG disp , we assume it to depend linearly on the contribution of each atom to the SASA [18,19] ΔG disp = X Nsph Here, the factor σ I is the atomic surface tension of the sphere I.
The equations for the gradient of ΔG cav and ΔG disp are provided in Data S1.

| Computational details
All the calculations in this study have been performed using a development version of the ORCA 4.2 electronic structure package at the DFT level of theory. The GCS has been implemented in ORCA, as described in the theory section. We denote by GVDW and GSES the cases where we use the scaled vdW-type cavity and the SES-type cavity for the contribution of solvation to the electrostatics of the system, respectively. When the nonelectrostatic contributions to the free energy of solvation are included in the SCF energy and its derivatives, we label the solvation approach as GVDW_nel and GSES_nel. In all calculations using the GCS, we place 110 charges per sphere. Regarding the SASA, which is used in the calculation of ΔG cav and ΔG disp , we consider 434 Lebedev points, due to the fact that the SASA is larger than the vdW surface or the SES. We have also performed a set of calculations using point charges in order to prove the validity of the GCS to treat the contribution to the electrostatics. In this case, we use a SES-type cavity generated through the GEPOL algorithm [11][12][13] with 60 triangles per sphere, discarding surface points that are closer than 0.1 a 0 from each other. The other GEPOL settings are NDIV = 5 and RMIN = 0.5 Å. We refer to this last method as PSES. The different approaches used in this study, with the type of charges and cavities employed are shown in Table 2.
In all calculations, we have employed the B3LYP functional [35][36][37] with the def2-TZVP basis set [38] using the resolution of the identity and the chain of spheres approximation (RIJCOSX). [39] With regard to the DFT integration grid, we consider the Grid7 settings that correspond to N ang = 770 and ξ = 5.67, being N ang the number of points of the largest angular grid, and ξ the number of radial shells for a given atom. In order to speed-up the calculation of molecular integrals, ORCA uses the Libint library, version 2.0.2-stable. [40] The energy convergence tolerance for the SCF cycle is equal to 10 −8 E h (VeryTightSCF). Geometry optimizations were done using the Ver-yTightOpt settings, which involve an energy tolerance equal to calculations. [41][42][43][44][45] To parametrize ΔG solv we use part of the SMD training set, [21] that involves experimental free energies of solvation, ΔG exp solv , for 318 solutes in 91 different solvents. The total number of considered solutes per solvent in our study is listed in Table 3, while the information on the type of solutes considered is provided in the Table S1.
Our training set involves organic solutes that just contain H, C, N, and O atoms. In order to derive the values for the atomic surface tensions in Equation (30), we proceed in the following way: 1. We perform a geometry optimization of the system in the gasphase. We get a final energy E 0 (it is equal to the second term in Equation (27)).
2. We perform a geometry optimization of the same system in solution. In this case, we include the contribution of ΔG cav to the SCF energy (Equations (28) and (29)) and its gradient but the effect of ΔG disp is not included neither in the SCF nor its gradient. We get a final energy E no− disp

We compute ΔG
5. We fit ΔG exp disp vs SASA I in order to get the different atomic surface tensions σ I via multivariate regression (Equation (30)).
6. Once we get the different σ I , we perform Step 2 but now including the effect of ΔG disp in the SCF energy and its gradient and calculate ΔG solv as shown in Equation (27). If it is needed, we refine slightly σ I . In addition to the systems listed in Table S1, we have tested our approach for ΔG solv for a set of cations and anions in water (see Table S2). We have also carried geometry optimizations of two large systems, the oxygen evolving complex of photosystem II, and vancomycine in protein and water, respectively. The corresponding computational details are provided in the results section.
where k runs over N leb , E i,k is the SCF energy of molecule i for the sub-  Figure 5, we show the ratio between t SCF , t grad , and t cycle grad , and the same quantity for N leb = 770 for both the GVDW and GSES models. All three quantities increase with N leb . If we focus on t/t 770 for the gradient, we observe a linear increase with respect to N leb . In this case, t grad,110 /t grad,770 is equal to 0.718 and 0.705 for the GVDW and the GSES models (t cycle grad,110 =t cycle grad,770 is equal to 0.695 and 0.692, respectively). When N leb increases from 110 to 194, there is an increase of 5% of t grad and t cycle grad with respect to the case of N leb = 770. For the SCF the same increase is of 2 and 4% for the GVDW and the GSES model. If we increase N leb a bit more, from 110 to 302, then the

| Efficiency
The efficiency of the GVDW and the GSES models has been stud- However, the next level of tessellation for the PSES model would involve 60 × 4 point charges per sphere, which is quite large according to other studies that employ this same methodology. The ratio between t SCF , t grad , and t cycle grad for one of the methods with respect to the same quantity for the other is provided in Table 4 for the same systems studied in the previous subsection in water. If we focus first on t cycle grad we observe that this quantity is slightly larger for the GVDW and the GSES models with respect to the PSES model (1.035 and 1.049, respectively). This very small difference is due to the reason explained above, that is, N q is larger when the GCS is employed and then there are more operations to compute per each gradient cycle.
Due to the fact that N GSES q > N GVDW q , t cycle grad is slightly larger for the GSES model as compared to the GVDW model. If we focus on t SCF and t grad , these quantities are slightly smaller when adopting the GVDW and GSES models with respect to the PSES model. In this case, although the cost per SCF iteration and gradient cycle should increase when adopting the GCS scheme, the total amount of cycles is slightly smaller in total. This is also the case when we compare the GVDW and the GSES models, being t SCF and t grad slightly smaller for the first model. From the results shown in Table 4 we do not observe a significant difference in terms of computational cost between the three schemes. This is especially interesting, as it shows that adopting the GCS does not make calculations more expensive as compared to the point charge scheme.
For completeness, we have also computed the ratio between t cycle grad and t cycle SCF (the average time spent in the SCF iterations before each optimization cycle) for the GVDW and GSES with respect to the vacuum case. The corresponding results, presented in Table 5, show that the GCS only adds an overhead of 7-8% relative to the gas-phase calculations. Regarding t cycle grad , we have studied how this quantity increases with the size of the system. In Figure 6, we plot t cycle grad as a function of the product N atoms × N q for the GVDW and the GSES. In both cases,

| Performance
In the present subsection, we study the effect of the GCS on the convergence of the SCF energy and its gradient. For that purpose, we focus on the oxygen evolving complex (OEC) of photosystem II (PSII).
The cluster model considered has 238 atoms (see Figure 7)  functional with the Zeroth-order regular approximation (ZORA) [48][49][50] and ZORA-contracted basis sets [51] (TZVP for Mn, Ca, O, and N atoms and the SARC-SVP for C and H atoms). The atom-pairwise dispersion correction with Becke-Johnson damping is used. [41][42][43] The system is implicitly solvated (ε = 8.0). For the SCF convergence we choose the tight criterium, while for the gradient optimization we adopt the opt criterium, according to the nomenclature used in ORCA. The resolution of the identity (RI) is used to accelerate the SCF convergence.
The results for E − E conv for the GVDW, GSES, and PSES models are shown in Figure 8. If we focus first on the results for the GVDW, we observe a very smooth convergence of the energy, without oscillations. This situation is slightly different when using a SES-type cavity, where E − E conv has some jumps, much smaller for the case of the GSES model than for the PSES model. If we focus on the PSES model,  This situation is worse in calculations where the initial structure is far away from its final counterpart (data not shown for simplicity) or for solvents where the geometry of the solute is prone to change substantially during its optimization. For these cases, the SCF energy and the gradient present several oscillations of large amplitude for the PSES. This is something that occurs more frequently for large systems but is also observed for small systems.
Regarding, the GSES model, Grad GSES MAX ≈ 3 × Grad GVDW MAX . In our implementation of the GSES, all contributions to the gradient are calculated, with no approximations for any of the terms in Equations (17) and (18). Therefore, the results shown in Figures 8 and   9 indicate that when performing a geometry optimization of a big systems, the use of the GSES is not as adequate as using the GVDW model. To prove that we have also optimized the structure of  Figure 10. This structure is used as a starting point for the calculations using the GVDW and the GSES models.
The convergence of the SCF energy with respect to the geometry optimization cycle is shown in Figure 11. Due to the fact that the calculation using GSES converges much slower than that using the GVDW model, for the GSES model we plot E − E, being  From the results exposed in Figures 8, 9, 11

| Free energy of solvation
Following the steps described in the computational details section, we have parametrized the dispersion contribution to ΔG solv for organic molecules containing H, C, N, and O atoms. In order to make the analysis of the results easier, we divide our results into two subsets: (a) solvents with low ε (ε < 10) and (b) solvents with mid and large values of ε (ε > 10). For completeness, the performance of our approach to compute ΔG solv has also been tested for a set of cations and anions in water.

| Solvents with low ε
This subset considers four nonpolar solvents: benzene, chloroform, cyclohexane, and toluene. The results for the fittings of the atomic surface tensions using the GVDW_nel and the GSES_nel models are shown in Table 6. In all cases the R Square coefficients for the multivariate regressions are higher than 0.95, being in most of the cases equal or larger than 0.99. There is almost no difference between the results for the GDW_nel and the GSES_nel models. Indeed, if we take as a reference the value for the GSES_nel model, then the average relative error between both sets of data is equal to 1.7%. In Figure 13, we show the correlation plot between ΔG calc solv using the atomic surface tensions in Table 6 and ΔG exp solv . The level of agreement between the experimental and the calculated free energies of solvation is very good, being the mean absolute error (MAE) lower than 1.0 kcal/mol (see Table 7). In general, MAE SES < MAE vdW (except for toluene where the number of solutes is slightly smaller than for the other solvents); however, the difference is too small to extract any trend. It is interesting to see the improvement of the calculated energies when adopting the GDVW_nel and GSES_nel with respect to the case of the GVDW and GSES models (see data in parenthesis in Table 7). For this last two models, the MAE is around three times larger than when considering nonelectrostatic contributions to the free energy of solvation showing that ΔG nel represents an important percentage of ΔG solv for nonpolar solvents.
Although the agreement between calculated and experimental ΔG solv is good, this quantity should not be the only guideline for the parametrization of a solvation model. Indeed, large polarization effects can occur within a continuum model but be counterbalanced by the term ΔG nel in Equation (27), giving a good prediction of ΔG solv . Such a large electrostatic contribution to ΔG solv is due in general, to the choice of a too small cavity for ΔG el in Equation (27). A variable in our model that controls the polarization of the solute is the parameter f that scales the vdW radii of the solute atoms. For simplicity, we have considered a solventindependent value for f equal to 1. dipole. [28][29][30] For instance, values of f = 1.6 and 1.8 are suggested for neutral solutes in carbon tetrachloride (ε = 2.24), and chloroform (ε = 4.90), respectively. [28,29] In order to evaluate our assumption of a solvent-independent f, we compare the value of the solute dipole moment in solution, μ S , to its counterpart in vacuum, μ vac . We restrict the analysis to the GVDW_nel model, as the obtained results are almost the same than those for the GSES_nel model. The corresponding plots μ S versus μ vac are shown in Figure 14, together with the linear fits μ S = a + b × μ vac . For solutes in cyclohexane, we observe an increase of around 16% in the dipole moment. Similar results are observed for benzene and toluene (increase of 19%). These results agree with the increase observed in Reference [29]. A larger enhancement is observed, however, for chloroform where b = 1.281 (increase of 28%). This last value does not agree with the value provided in Reference [29] for a series of solutes in chloroform (increase of 12% in μ).

| Solvents with ε > 10
For this subset, we consider two solvents: octanol and water. In the case of octanol, the value of ε is still low (ε = 10.30). Indeed, octanol (C 8 H 18 O) can be considered as nonpolar due to the fact that the polar nature of the molecule only comes from the O─H bond (being the rest of bonds nonpolar). Then, although ε is slightly larger than for the case of pure nonpolar solvents, there should not be significant effects from the nonconsideration of a first solvation explicit shell on the computed ΔG solv . The situation for water is very different. Water is a highly polar molecule with a large value of ε, being able to form hydrogen bonds with other molecules, forming quite rigid solvation shells around solutes and consequently screening out the interaction between them. Then, the derivation of a set of σ I for water that are independent of the protic or aprotic character of the solute is not so clear. The multivariate fittings of σ I for octanol and water for the GVDW_nel and the GSES_nel models are shown in Table 8, while the correlation plots ΔG calc solv vs ΔG exp solv are provided in Figure 15. The agreement between calculated and experimental solvation energies for octanol is slightly worse than that for the case of nonpolar solvents, as seen in Figure 15a. Nevertheless, most of the points lie in the 1 kcal/mol error region, being the MAE equal to 1.09 kcal/mol and 1.10 kcal/mol for the GVDW_nel and the GSES_nel models, respectively (see Table 9). Again, the effect of the type of solute cavity is negligible on ΔG solv . The electrostatic contribution to ΔG solv plays a larger role for octanol than for nonpolar solvents. This can be seen from the value of the MAE for the GVDW and GSES models, which equals 1.58 kcal/mol and 1.54 kcal/mol, respectively, that is, 1.5 times T A B L E 6 Atomic surface tension σ vdW and σ SES for H, C, N, and O for solvents with low ε larger than the same quantity for the GVDW_nel and the GSES_nel (while it was three times larger for nonpolar solvents).
For water, the number of calculated points that deviate more than 1 kcal/mol from their experimental counterpart is larger than for octanol (see Figure 15b), being the MAE equal to 1.27 kcal/mol for both type of solute cavities (see Table 9). This value is not so far from the MAE for the GVDW and GSES models (1.68 and 1.64 kcal/mol), showing that for water, the nonelectrostatic contributions play a smaller effect than for solvents with low ε. Indeed, if we plot μ S versus μ vac for water using the GVDW model (see Figure 16), we observe an increase of around 36% in μ. This enhancement agrees with the value of 35% reported in Reference [30]. If we assume ΔG nel to be linearly proportional to the SASA, as it is suggested from experimental transfer data for hydrocarbons into water, [52] then we obtain a MAE equal to 1.67 kcal/mol. This confirms that the strategy described in Equations (29) and (30) is more suitable than simply taking ΔG nel to depend on the total SASA. However, as mentioned before, part of the MAE shown in Table 9 for water may come from the fact that our approach considers the same values for σ H , σ C , σ N , and σ O , whatever is the content of the solute. In order to prove that, we split out our training set for water into four different subsets: solutes that just con-  Table 10. More information about the solutes used per subset is provided in Table S1. After having defined the subsets, we fit the σ I for each of them and compare ΔG calc solv with ΔG exp solv . The corresponding results for the σ I are shown in Table 11, the correlation plots for the solvation energy are provided in Figure 17, and finally the MAE for each subsets and solute cavity are shown in Table 12.
The improvement in the results is rather important, as seen from Figure 17 and Table 12, with all subsets having now a MAE lower than 1 kcal/mol, regardless of the solute cavity used. In order to understand the reasons behind the improvement, we analyze the values of σ I per each subset (Table 11). First, we observe that σ H for solutes that contain N atoms is twice the same quantity for solutes that do not contain N atoms. Indeed, molecules with N atoms have N─H bonds, F I G U R E 1 4 Correlation plot between the dipole moment in solution, μ S , and the same quantity in vacuum, μ vac , for solvents with low ε for the GVDW_nel model. The slope b of the linear fit (μ S = a + b × μ vac ) and the corresponding correlation coefficient r are shown for completeness With the results shown above, we believe that both points are sufficiently addressed. for ions than for neutrals. [21,27] [21] A further reduction of f = 1.10 gives very accurate results for ΔG solv with a MAE lower than 3 kcal/mol and a percentage error lower than 5%. Although decreasing f also improves the agreement between calculated and experimental data for the case of anions, the magnitude of the MAE and the percent error is larger than for cations. For f = 1.10, both the MAE and the percent error for anions using the GVDW_nel and the GSES_nel models are around 2.5 times larger than the same quantities for cations. A further decrease of f (f = 1.00) would improve the results for anions but it is not recommended as it can lead to unrealistic estimates of the electrostatic part of ΔG solv . [27] Refining the atomic surface tensions in Equation (30) for anions should not change the accuracy substantially as electrostatic effects dominate over nonelectrostatic effects in ions. A better approach would be to refine slightly the vdW radii for the solute atoms used to construct the vdW-type or the SES-type cavity for the treatment of solvation electrostatic effects. However, such a change would involve a careful T A B L E 1 1 Atomic surface tension σ vdW and σ SES for solutes with different atomic content in water  When the SES-type cavity is adopted within the GCS, the SCF energy converges much slower. In this case, the maximum component of the gradient vector oscillates almost indefinitely with a small amplitude.

| Ions in water
Therefore, when SES-type cavities are used, it is recommended to soften slightly the thresholds for the gradient in order to reach convergence. In any case, the use of the GCS prevents sudden jumps in the SCF energy observed when adopting the point charge scheme, and which are on the order of 100 kcal/mol.