Spicing up continuum solvation models with SaLSA: The spherically averaged liquid susceptibility ansatz

A common ingredient in continuum solvation models is the formation of a cavity that excludes the solvent from a region of space occupied by the solute. Atom-based parametrizations typically exclude the centers of solvent molecules from a union of spheres centered on each solute atom with radius equal to the sum of the atomic and solvent van-der-Waals (vdW) radii, and then construct a dielectric cavity that is smaller by an empirical solvent-dependent radius.⁷ In contrast, the density-based approaches adopt a smoothly varying dielectric constant which is a function of the local solute electron density that switches from the vacuum to bulk solvent value at a solvent-dependent critical electron density n_c.^11,12

Here, since we explicitly account for the nonlocal response, we require only the distribution of the solvent molecule centers. Because the distance of these centers from the solute atoms corresponds to the distance of nearest approach of two non-bonded systems, vdW radii, which are defined in terms of non-bonded approach distances,¹⁹ provide a reasonable guess. However, directly using the vdW radii does not account for changes in the electronic configuration between the isolated atom and molecules or solids. A description based on the electronic density would instead naturally capture this dependence.

The interaction of non-bonded systems is dominated by Pauli repulsion at short distances, which depends on the overlap of the electron densities of the two systems. Indeed, Fig. 1 shows that the atom separation, R₁₂, at which the electron density overlap $n ̄ ( R 12 ) =∫d r → n 1 ( r → ) n 2 ( r → )$ crosses a threshold value of $n ̄ c =1.42×1 0 − 3 a 0 − 3$ correlates well with the sum of vdW radii¹⁸ R₁ + R₂ of the two atoms. Here, $n 1 ( r → )$ and $n 2 ( r → )$ are spherical electron densities of the two atoms (calculated using OPIUM²⁰), centered R₁₂ apart, and we obtain $n ̄ c$ by minimizing $∑ i j ( R i + R j − R i j ) 2$⁠. This result is insensitive to the choice of exchange-correlation approximation used to calculate the densities: at the optimized $n ̄ c$⁠, the RMS relative error in R₁₂ compared to R₁ + R₂ is 8.1% for the local-density approximation,²¹ 8.3% for the Perdew-Burke-Ernzerhof (PBE) generalized-gradient approximation,²² and 9.9% for Hartree-Fock theory.

The above analysis provides a universal threshold on the density product that can estimate the approach distance of non-bonded systems. We utilize this capability to determine the distribution of solvent molecule centers around a solute with electron density $n ( r → )$⁠. Our ansatz requires a spatially varying but isotropic initial distribution of molecules. Therefore, we compute overlaps with the spherically averaged electron density $n lq 0 ( r ) =∫ d n ˆ 4 π n lq ( r n ˆ )$⁠, where $n lq ( r → )$ is the electron density of a single solvent molecule and $n ˆ$ is a unit vector. Following Ref. 14, we describe the spatial variation of the solvent distribution by the smooth “cavity shape” functional

which smoothly transitions from vacuum (s = 0) to bulk fluid (s = 1) as the overlap of the solute and solvent electron densities (readily calculated as a convolution) crosses the universal overlap threshold $n ̄ c$⁠.

We begin with the in-principle exact joint density-functional description¹⁶ of the free energy of a solvated electronic system

where A_HK is the Hohenberg-Kohn functional¹ of the solute electronic density $n ( r → )$⁠, Φ_lq is a free energy functional for the solvent in terms of nuclear site densities ${ N α ( r → ) }$⁠, and ΔA captures the free energy of interaction between the solute and solvent. (See Ref. 16 for details about the theoretical framework.)

We adopt the Kohn-Sham prescription² with an approximate exchange-correlation functional for A_HK[n] and the polarizable “scalar-EOS” free energy functional approximation²³ for Φ_lq[{N_α}]. We separate ΔA in (2) as the mean-field electrostatic interaction and a remainder that is dominated by electronic repulsion and dispersion interactions. We then assume that the remainder is responsible for determining the initial isotropic distribution $N 0 ( r → ) =N bulk s ( r → )$ of solvent molecules, where N_bulk is the bulk number density of solvent molecules and $s ( r → )$ is given by (1). Substituting the free energy functional from Ref. 23, we can then write (2) as

Above, analogously to the Kohn-Sham approach, the liquid free energy functional employs the state of the corresponding non-interacting system specified by two sets of independent variables. The first, $p ω ( r → )$⁠, is the orientation probability density of finding a solvent molecule centered at $r →$ with orientation ω ∈ SO(3). The solvent site densities $N α ( r → )$ are dependent variables that are calculated from $p ω ( r → )$⁠. The second variable is the polarization density, $P α ( r → )$⁠, for each solvent site.

The second term in (3), Φ₀[N₀], collects all the free energy contributions due to the initial isotropic distribution, so that all the subsequent terms are zero when $p ω ( r → ) = N 0 ( r → )$ and $P α ( r → ) =0$⁠. The third term is the non-interacting rotational entropy at temperature T, the fourth term is a weighted-density correlation functional for dipole rotations, and the fifth term is the potential energy for molecular polarization with site susceptibilities χ_α. C_rot and C_pol parametrize correlations in the rotations and polarization, respectively, and are constrained by the bulk static and optical dielectric constants. The final term is the mean-field interaction between the solute charge density $ρ el ( r → )$ and the induced charge density in the liquid $ρ lq ( r → ) −ρ lq 0 ( r → )$ (where $ρ lq 0$ is the charge in the initial isotropic distribution), and the self-energy of that induced charge. Here,

where ρ_α(r) and w_α(r), respectively, specify decomposition of the solvent molecule’s charge density and nonlocal susceptibility into spherical contributions at the solvent sites. Bulk experimental properties of the liquid and ab initio calculations of a single solvent molecule constrain all involved parameters. See Ref. 23 for details; the terms above are identical, except for the inclusion of solute-solvent interactions in the final electrostatic term and for the separation of the initial isotropic contributions into Φ₀[N₀].

Next, we treat the orientation-dependent pieces perturbatively by expanding $p ω ( r → ) = N 0 ( r → ) ( 1 + ∑ l m m ′ x m m ′ l ( r → ) D m m ′ l ( ω ) )$⁠, where $D m m ′ l ( ω )$ are the Wigner D-matrices (irreducible representations of SO(3)).²⁴ We then expand free energy (3) to quadratic order in the independent variables $x m m ′ l ( r → )$ (rotation) and $P → α ( r → )$ (polarization), formally solve the corresponding linear Euler-Lagrange equations, and substitute those solutions back into the quadratic form. After some tedious but straightforward algebra involving orthogonality of D-matrices, addition of spherical harmonics and their transformation under the D-matrices, as well as Fourier transforms to simplify convolutions, we can show that the resulting free energy to second order is exactly

Here, $K ˆ$ is the Coulomb operator and $χ ˆ$ is the nonlocal “spherically averaged liquid susceptibility” (SaLSA), expressed conveniently in reciprocal space as

where $f ̃ ( G → )$ is the Fourier transform of $f ( r → )$ for any f. The first term of (6) captures the polarization response, where $N α 0 ( r → ) = N 0 ( r → ) *δ ( r − R α ) / ( 4 π R α 2 )$ is the site density at the initial configuration $p ω = N 0 ( r → )$ of a solvent site at a distance R_α from the solvent molecule center. The second term of (6) captures the rotational response of solvent molecules with charge distribution $ρ mol ( r → )$⁠, decomposed in Fourier space as $ρ ̃ mol ( G → ) = ∑ l m ρ ̃ mol l m ( G ) Y l m ( G ˆ )$⁠. The prefactor $C rot l =C rot$⁠, the dipole rotation correlation factor²³ for l = 1, and it equals unity for all other l.

For practical calculations, we rewrite the last term of (5) as ∫(ϕ − ϕ₀) ρ_el/2, where $ϕ 0 = K ˆ ρ el$ is the electrostatic potential in vacuum and ϕ is the total (mean-field) electrostatic potential which solves the modified Poisson-like equation $( ∇ 2 + 4 π χ ˆ ) ϕ=−4πρ el$⁠. Note that this term is analogous to the solvent polarization energy in the polarizable continuum model (PCM) approach⁷ and other local solvation models,^11,12 except that the nonlocality of $χ ˆ$ (given by (6)) would produce bound charges within a finite range of the cavity surface, instead of only at the surface. The l = 1 rotational and polarization terms in $χ ˆ ϕ$ have the structure $∇⋅w*N ( r → ) w*∇ϕ$ which resembles the Poisson equation for an inhomogeneous dielectric, except for the convolutions that introduce the nonlocality. For neutral solvent molecules, the l = 0 term captures the interaction of the solute with a spherical charge distribution of zero net charge, and is zero except for small contributions from non-zero but negligible overlap of the solute and solvent charges. However, note that the SaLSA response easily generalizes to mixtures, and for ionic species in the solution, the l = 0 terms convert the Poisson-like equation to a Helmholtz-like equation that naturally captures the Debye-screening effects of electrolytes.¹² The l > 1 terms resemble (nonlocal versions of) higher-order differential operators and capture interactions with higher-order multipoles of the solvent molecule, which decrease in magnitude with increasing l. We find that including terms up to l = 3 is sufficient to converge the solvation energies to 0.1 kcal/mol.

The nonlocality of the SaLSA response allows the fluid bound charge to appear at a distance from the edge of the cavity. For example, for a water molecule in liquid water, the SaLSA cavity transitions at about the first peak of the radial distribution function g_OO(r) of water (Fig. 2), whereas the bound charge is dominant at smaller distances. In contrast, local solvation models require a smaller unphysical cavity to produce bound charge at the appropriate distance to capture the experimental solvation energy. This key difference from the local models allows a non-empirical description of the electric response in SaLSA.

At this stage, solvated free energy (5) is fully specified except for the free energy of the initial configuration Φ₀[N₀], dominated by electronic repulsion, dispersion, and the free energy of forming a cavity in the liquid. We set Φ₀[N₀] = G_cav[s] + E_disp[N₀], where G_cav is a non-local weighted density approximation to the cavity formation free energy and E_disp empirically accounts for dispersion and the remaining contributions, exactly as in Ref. 15. (See Ref. 15 for a full specification.) Briefly, G_cav is completely constrained by bulk properties of the solvent including density, surface tension, and vapor pressure, and reproduces classical density functional and molecular dynamics predictions for the cavity formation free energy²³ with no fit parameters.

E_disp employs a semi-empirical pair potential dispersion correction²⁶ which includes a scale parameter s₆, which we fit to solvation energies below. In principle, it might be possible to eliminate this parameter as well by treating dispersion explicitly using the ab initio polarizability of the solute²⁷ or using nonlocal correlation functionals that include dispersion interactions.²⁸ However, these methods require significantly higher computational effort and are much more complicated to implement in a self-consistent solvation code which requires analytical derivatives for electronic and geometry optimization. For simplicity and computational expediency, we therefore retain the single-parameter empirical dispersion functional from Ref. 15 in this work. Note, however, that unlike previous continuum solvation models, the dominant electrostatic response includes no parameters that are fit to solvation energies.

We implement the SaLSA solvation model in the open source density-functional software JDFTx,²⁹ and perform calculations using norm-conserving pseudopotentials²⁰ at 30 E_h plane-wave cutoff and the revTPSS meta-generalized-gradient exchange-correlation functional.³⁰ For three solvents, water, chloroform, and carbon tetrachloride, we fit the sole parameter s₆ to minimize the RMS error in the solvation energies of a small but representative set of neutral organic molecules with a variety of functional groups and chain lengths (same set for each solvent as in Ref. 15). Table I summarizes the optimum values of s₆ and the corresponding RMS error in solvation energies. The RMS errors of SaLSA are only slightly higher than those of the local solvation model from Ref. 15 that includes additional fit parameters for the electric response.

Figure 3 compares the aqueous solvation energy predictions of SaLSA with those of the linear and nonlinear local-response models from Ref. 14. (Briefly, the linear PCM in Ref. 14 is an iso-density solvation model with an empirical surface tension term to describe cavity formation and dispersion, and the nonlinear PCM additionally accounts for dielectric saturation in the liquid response. See Ref. 14 for details.) All three models perform comparably for the neutral molecule set (Fig. 3(a)) used for the parameter fit, but SaLSA is substantially more accurate for inorganic ions (Fig. 3(b)). In particular, the local models severely over-predict the solvation energies of small cations, and the nonlocal SaLSA model reduces the error by a factor of three for Li+ and Na+. However, SaLSA does not correct the systematic over-solvation of cations compared to anions, a known deficiency of electron-density based solvation models.¹³ 

In order to directly compare the accuracy of SaLSA to existing models, Table II shows the mean absolute errors in the solvation energies of an identical set¹³ of 240 neutral molecules, 51 cations, and 55 anions in water, predicted by SaLSA and various parametrizations of the self-consistent continuum solvation (SCCS) model¹¹ and the integral equation formalism polarizable continuum model (IEF-PCM)^34,35 in GAUSSIAN.³⁶ SaLSA exhibits comparable accuracy to the empirical solvation models for the neutral solutes. The accuracy for charged solutes is only slightly worse, which is remarkable considering that the solvation energies for these solutes is dominated by the electrostatic term which includes no fit parameters in SaLSA.

Finally, we turn to solvation in diffusion quantum Monte-Carlo (DMC) calculations. Conventional density-based solvation models^11–14 are sensitive to the electron density in the low density regions of space (⁠$n ( r → ) ∼1 0 − 4 −1 0 − 3 a 0 − 3$⁠). This sensitivity imposes extremely stringent restrictions on the statistical noise in the DMC electron density, rendering self-consistent solvated DMC calculations impractical. A scheme correct to first order that combines a solvated density-functional calculation with a DMC calculation in an external potential provides reasonable accuracy without calculating DMC electron densities.¹⁷ However, full self-consistency would be particularly important for systems where density-functional approximations fail drastically.

The nonlocality of the SaLSA model, particularly the dependence of cavity shape (1) on a convolution of the density rather than the local density, significantly reduces its sensitivity to noise in the electron density, and enables self-consistent DMC solvation in a straightforward manner. We start with an initial guess for the density (from solvated DFT), compute the potential on the electrons from the SaLSA fluid model, and perform a DMC calculation (using CASINO³⁷) in that external potential while collecting electron density (using the mixed estimator in CASINO). We then mix the DMC electron density with the previous guess, update the fluid potential, and repeat till the density becomes self-consistent. The estimation of the solvated energy at each cycle proceeds exactly as in Ref. 17 with slight differences in the details of the DMC calculation: we use Trail-Needs pseudopotentials^38,39 with a DFT plane-wave cutoff of 70 E_h and a DMC time step of 0.004 $E h − 1$⁠.

Figure 4 shows that the solvated DMC energy converges to within 0.1 kcal/mol in just three self-consistency cycles. The solvation energies of methanol and carbon monoxide are essentially unchanged from the density-functional results, whereas the solvation energy of the fluoride anion, which suffers from strong self-interaction errors in DFT, is corrected by 3 kcal/mol. This change, although in the right direction, is small compared to the 20 kcal/mol error in the predicted SaLSA solvation energy using DFT (Fig. 3(b)). The error in the solvation of anions is therefore not predominantly caused by the inaccuracy of electronic density-functional approximations for anions. In fact, we expect small anions to form strong hydrogen bonds with water, an effect that is missed by the electrostatic response of all solvation models including SaLSA. We therefore expect that empirical corrections for hydrogen bonding,⁴⁰ or equivalently the charge asymmetry of solvation,⁴¹ will be necessary to remedy the accuracy for anions.

The linear-response limit of joint density-functional theory, combined with an electron-density overlap based estimate of the initial solvent molecule distribution, provides a nonlocal continuum solvation model with no empirical parameters for the electric response. Consequently, this “SaLSA” model extrapolates reliably from neutral organic molecules to ions and is an excellent candidate for describing highly polar and charged surfaces in solution; a comparison of the performance of various solvation models for such systems would be highly desirable and the subject of future work. Further, the nonlocality of this model enables self-consistent solvation in diffusion Monte Carlo calculations. This opens up the possibility of studying systems in solution for which standard density-functional approximations fail, such as the adsorption of carbon monoxide on transition metal surfaces. Additionally, SaLSA should facilitate the development of more accurate and perhaps more empirical models that, for example, account for the charge asymmetry in the solvation of cations and anions.

This work was supported as a part of the Energy Materials Center at Cornell (EMC²), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0001086.