The Poisson–Boltzmann model for implicit solvation of electrolyte solutions: Quantum chemical implementation and assessment via Sechenov coefficients

In contrast to the ab initio calculation of chemical processes in the gas phase, the investigation of such processes in solution adds an additional layer of tremendous complexity. The electrostatic interaction between solvent and solute directly alters the enthalpic landscape, whereas osmotic pressure and reorientation of the solvent molecules based on the electrostatic interactions has a profound influence on the free energy. Although these effects can in principle be calculated with molecular dynamics (MD) simulations, the amount of sampling required to converge the free energy combined with the large number of structures involved in chemical processes for which such a sampling would have to be performed renders this approach unfeasible at present.

Hence, implicit solvation models, rather than explicit ones, are usually the method of choice in computational quantum chemistry. Famous examples of the former are the polarizable continuum model (PCM)^1–4 or the conductor like screening model (COSMO)⁵ that exists in various flavors and specializations, are available in most multipurpose electronic structure packages, and have proven suitable in countless applications.^6,7

For the description of electrochemical processes, however, the effect of the electrolyte ions needs to be taken into account. In the realm of electrocatalysis, the phenomenological description of the solvent/electrolyte interaction with charged surfaces has a long history. Helmholtz was the first to describe the formation of a layer of counterions in a certain distance to a charged surface—the Helmholtz layer—which resembles a capacitor and stores electrostatic energy.⁸ A little later, Gouy⁹ and independently Chapman¹⁰ introduced the concept of a diffuse electrochemical double-layer based on the notion of an exponential decay of the electric potential away from the charged surface. A combination of both approaches, a rigid inner layer close to the charged surface, called the Stern layer, and the Gouy–Chapman type diffuse layer, was introduced in 1924 by Stern.¹¹ The Poisson–Boltzmann (PB) model^12–22 for implicit electrolyte solvation is a computational realization of the Gouy–Chapman model with the possibility to include a Stern-layer correction.²³ In this model, the electrostatic potential of the system is calculated from the Poisson equation appended by an additional term describing the charge density of the mobile electrolyte ions. The charge density of the mobile ions is modeled by a Boltzmann distribution mediated by the electrostatic potential, hence the term Poisson–Boltzmann (PB) model. Since the electrolyte charge density depends on the total electrostatic potential that is in turn determined by the PB equations, a self-consistent solution, often on a real-space grid, is required. The PB model is especially popular in MD simulations and biochemical applications^24–28 but has also found increasing interest in the modeling of chemical processes on surfaces and electrode-electrolyte interfaces,^29–31 mostly in software packages for extended solid systems employing periodic boundary conditions.^23,32,33 We would further like to mention that the “integral equation formalism” of PCM (IEF-PCM)^34–36 can account for electrolyte solutions but assumes a linear approximation and does not include a finite-ion size correction.³⁷ 

In this article, we present a three-dimensional real-space grid implementation in the general purpose electronic structure package Q-Chem.³⁸ This implementation paves the way to investigate chemical processes and reactions in electrolyte solutions for isolated molecules or small cluster models. Especially for the study of electrocatalytic reactions, such a solvent model is imperative in order to calculate qualitatively meaningful free energy profiles.

We will first describe the details of the theory and the hierarchy of approximations, before we analyze the choice of parameters of the model and optimize those with respect to experimental data. In Sec. V, we discuss the effect of an electrolyte concentration dependent Stern-layer thickness on the PB model and suggest future improvements on the model informed by either experiments or MD simulations.

The theory of the Poisson–Boltzmann implicit solvation model has been introduced and reviewed in numerous publications (see, e.g., Refs. 13, 19, 21–23, and 39). However, the specific parameterization of our implementation and the hierarchy of approximations we introduce here requires some theoretical background that we provide in Secs. II A–II D.

The total electrostatic potential $ϕtotr$ of a molecular system embedded in an electrolyte solution can be obtained as the solution to Poisson’s equation appended by a term describing the charge density of the electrolyte ions ρ^ions(r),

where ϵ(r) is the spatially dependent dielectric permittivity, ρ^sol(r) is the charge distribution of the solute, and the explicit expressions for f(ϕ^tot(r), r) define the different approximations we will discuss in the following. Following the implementation of the Poisson equation solver in the Q-Chem program package by one of us,⁴⁰ and equivalent to the soft-sphere model by Fisicaro et al.,⁴¹ we define the spatially dependent dielectric constant as a product of error functions for each atom,

with the dielectric permittivity in vacuum ϵ₀, the dielectric permittivity of the solvent ϵ_solv and

where d_α is a given atom-specific radius (vide infra) and Δ defines the interpolation length between the vacuum permittivity in the solvent region and the bulk solvent permittivity at larger distances. The error function interpolates smoothly between these two values over a region of about 4Δ.

If ρ^ions(r) is assumed to be defined by a Boltzmann distribution of the electrolyte charges, Eq. (1) becomes the nonlinear Poisson–Boltzmann (PB) equation with

where m is the number of electrolyte ion species, q_i is the charge of electrolyte i, $cib$ is its bulk concentration, k_B is the Boltzmann constant, and T is the absolute temperature. The ion exclusion function λ(r) ensures that the electrolyte ion concentration tends to zero inside the solute cavity and is defined in a similar way to ϵ(r) with an additional parameter a that corrects for the Stern layer by excluding an even larger region around the solute,

For a 1:1 electrolyte (q₁ = e, q₂ = −e, $c1b=c2b=cb$⁠, with the elementary charge e), the charge density of the electrolyte ions can be written as

The function f_PB(ϕ^tot(r), r) then reduces to

The Poisson–Boltzmann equation can be recast into an Euler–Lagrange equation¹³ for a Lagrangian L(ϕ^tot, r, ϕ_x, ϕ_y, ϕ_z) (with ϕ_α = ∂ϕ/∂α) that is still to be defined,

It follows from a theorem of the calculus of variations that the solution to this set of differential equations defines the minimum of the integral of L over the independent variables r = x, y, z. Since the minimum condition defines the thermodynamic equilibrium state, this integral describes the electrostatic free energy of the system,

In the case of the nonlinear Poisson–Boltzmann equation, integration gives the functional L for a 1:1 electrolyte as

The second term describes the excess osmotic pressure due to the presence of the electrolyte ions. The Euler-Lagrange equations are fulfilled for all L′ = C₁L + C₀, with the constants C₁ and C₀ that can be chosen to fulfill certain criteria for the units and boundary condition, respectively. We therefore fix C₁ = (4π)⁻¹ which is essentially a transformation to MKS units and C₀ = 0 such that in the absence of ρ^sol(r) [and hence when ϕ^tot(r) = 0] only the excess osmotic pressure term deviates from zero.

The electrostatic free energy can then be written as

where ΔΠ_PB is the excess osmotic pressure, E = ∇ϕ^tot(r) is the electric field, and D = ϵ(r)E is the electric displacement field.

Integrating Eq. (1) by parts, we obtain³⁹ 

and the electrostatic free energy can then be written in the more common form,²¹ 

The electrostatic free energy expression is drastically simplified in the linear Poisson–Boltzmann (LPB) equation where the exponential term in the electrolyte ion charge density distribution is approximated by the first two terms of the Taylor expansion,

Restricting ourselves again to the case of a 1:1 electrolyte, this reduces to

The functional L_LPB in the linear case can then be written as

With the same arguments as above, we can fix the constants C₀ = 0 and C₁ = 1/4π. The electrostatic free energy for the linear Poisson–Boltzmann equation then reads

which is the same expression as for a pure solvent without electrolyte and only the electrostatic potential is modified by the presence of the electrolyte ions.

The standard Poisson–Boltzmann equation assumes point-like ions, which leads to a charge accumulation of the electrolyte that is unphysically large. Finite ion size can however be included in the theory leading to the modified Poisson–Boltzmann (MPB) equation.^17,21,23,42 This can be achieved by a modification of the expression for the ionic concentrations,

The highest possible concentration for electrolyte species j assuming closest-packing $cjmax$ can be written as

with the closest-packing factor p = 0.74, the Avogadro number N_A, and the effective ion radius of the electrolyte ions R_j. Restricting ourselves again to a 1:1 electrolyte, we obtain for the charge density of the electrolyte ions

with

Integration with respect to ϕ^tot(r) gives the functional,

Proceeding in the same manner as above, we obtain the constants C₁ = 1/4π and C₀ = 0 leading to the following expression for the electrostatic energy:

or alternatively Eq. (13), where the excess osmotic pressure is now replaced by ΔΠ_MPB.

The finite ion size also alters the expression for the electrolyte charge density for the linearized version of the modified Poisson–Boltzmann (LMPB) equation,

Note that the last term in the denominator in the first line is usually very small and even vanishes for equal ionic radii. The second line of Eq. (24) offers a simple interpretation of the finite size effect. The charge density is scaled by the fraction of the maximum electrolyte density (c^b/c₁₊₂) when the ion exclusion function λ(r) < 1 meaning that far away from the solute, the charge density is unaffected by the finite size effects. As the electrolyte charge density approaches zero close to the solute anyway, the finite ion size effect is most important close to the Stern layer.

The LMPB expression for the electrostatic energy [Eq. (17)], however, is the same as for the LPB case.

As the total electrostatic free energy is rarely of interest, we will focus now on expressions for the electrostatic contribution to the free energy of solvation in an electrolyte solution. In general, this quantity can be expressed as

where the first term is the electrostatic free energy of the solute density optimized in the electrolyte solution (⁠$ρessol$⁠), the second term describes the electrostatic free energy of the same solute density in a vacuum environment and the last term corresponds to the pure electrolyte solution. The absence of the charge density of the solute [ρ^sol(r) = 0] in the last term results in ϕ^tot(r) = 0 and λ(r) = 1 such that the individual terms of Eq. (25) can be written as

Starting from Eq. (26), we can now write down the expressions for the electrostatic contributions to the free energies of solvation for the different approximations. In the case of the L(M)PB, the free energy of solvation is simply

because $∫ΔΠPBλ(r)=1,ϕtot(r)=0dr=0$ and $∫ΔΠPBdr=12∫ρions(r)ϕtot(r)dr$⁠. For the PB model without size modification of the ions, we have

and finally the expression for the MPB solvation free energy reads

We note here that this expression is identical to the expression for the electrostatic part of the free energy of solvation in Refs. 19, 23, and 43 when the packing factor for a simple cubic lattice (p ≈ 0.52) is assumed instead of the closest-packing factor. This corresponds to the lattice-gas model assumed in that work.

In order to calculate the full free energy of solvation, the electronic energy difference for the charge density distributions optimized in vacuum and in implicit solvent has to be added, as well as an additional term to account for nonelectrostatic effects,

Correction terms for the nonelectrostatic part of the free energy of solvation are usually calculated from the surface area and volume of the solute cavity.^7,44–47 Since the cavity is fixed in our implementation and we are only concerned in the following with differences in the free energy of solvation for varying electrolyte ion concentrations (ΔΔG_ion), the nonelectrostatic corrections cancel out and we will not consider them.

In order to solve Eq. (1), we follow the algorithm described by Andreussi et al. in Ref. 47 and solve iteratively for the electrostatic potential with a multigrid solver as described in Ref. 40, rather than applying a preconditioned conjugate gradient method as described in Ref. 39. We first isolate the effect of the dielectric response into a polarization charge density ρ^pol(r),

with ρ^tot(r) = ρ^sol(r) + ρ^ions(r). Starting from ρ^ions(r) = 0 and ϕ^tot(r) = 0, we converge the total electrostatic potential as shown in Algorithm 1. We observed that a much stronger damping (κ = 0.2) is necessary to converge ρ^ions(r) than is necessary for ρ^iter(r) (η = 0.6).

We will now discuss some of the properties of our parameterization on the example of 4-nitroaniline in a 15 Å cubic box (see Fig. 1). Figure 1(a) shows the error function for the dielectric permittivity ϵ(r) (blue line, left y-axis) and the ion exclusion function λ(r) (red line, right y-axis) for a hypothetical atom with d_α = 1.0 Å and a Stern-layer thickness a = 0.5 Å. In Fig. 1(b), these two functions are displayed for a cut through the main σ symmetry plane of 4-nitroanisole. Clearly, the solvent (blue)—characterized solely by its dielectric permittivity—is closer to the solute molecule than the electrolyte ions (red) because of the Stern-layer. The converged polarization charge density ρ^pol(r) and the ion charge density ρ^ions(r) are displayed in Fig. 1(c). As can be seen from Eq. (33), the polarization charge density resides in the domain where ϵ(r) varies, meaning in the onset of the error function. For the electrolyte ion charge density, a clear separation of the positively and negatively charged ions is observed, resulting in an accumulation of positive charge close to the negatively polarized part of the solute molecule and vice versa.

All calculations were carried out in a development version of Q-Chem 5.2.³⁸ Unless otherwise noted, we employed the ωB97X-V density functional⁴⁸ with a def2-TZVPP basis set⁴⁹ and standard integration grids. Our implementation builds up on a previous multigrid solver for the Poisson equation.^40,50 Unlike the Poisson solver that is reported in Ref. 40, we do not apply Gaussian blurring to the nuclear contribution to the electrostatic potential. In Sec. IV B, we place the molecule in a cubic box with an edge length of 20 Å. With 97 grid points in each dimension, the grid spacing amounts to 0.206 Å. After some initial tests, the convergence threshold for the conjugate gradient step on each multigrid level was set to 10⁻⁵ a.u., the convergence of the polarization charge density ρ^pol(r) was set to 10⁻³ a.u. and the convergence for the electrolyte ion charge density ρ_ion(r) was set to 10⁻⁴ a.u. for each grid point.

The initial implementation described here is not yet as efficiently parallelized as other highly optimized multigrid solver implementations described in Refs. 21, 22, and 42. Solving for the electrostatic potential hence takes considerably longer than the electronic structure calculation. In the case of caffeine, the Poisson–Boltzmann calculation takes roughly five times longer than the electronic structure calculation. The benefit of the real space grid implementation is certainly the independence of the calculations for individual grid points that enable excellent scaling over several threads. In future work, we will hence focus on a more efficient message passing interface (MPI) parallelization. The uniform Cartesian grid employed here might also not be an optimal choice since the electrostatic potential, as well as the charge density distributions vary more close to the solute cavity then they do further away from the solute. Adaptive grids²¹ or more generally adaptive mesh refinement methods^51,52 might be a more economical choice and can be combined with multigrid approaches^53–55 and help to focus the computational resources to the areas where they are needed most. For a moderate number of grid points, typical schemes for convergence acceleration of self-consistent approaches such as the direct inversion of the iterative subspace (DIIS) method^56,57 can further be applied to efficiently reduce the number of iterations. This approach can be understood as an advancement over the damping with the parameters η and κ in our current algorithm.

Having introduced the hierarchy of PB models with their respective approximations, we will now analyze the numerical accuracy, find an optimal choice for the parameters of the model, and discuss the implications upon inclusion of finite ion size.

In order to assess the numerical accuracy we can expect from our model, we first calculate the electrostatic free energies of solvation for a singly charged cation with a radius of a = 2 Å. If the linear Poisson–Boltzmann equation is to be employed, the result for the free energy of solvation can be compared to the analytical solution provided by the Debye–Hückel theory,^37,58,59

where a is the ion radius, ϵ_r is the dielectric permittivity of the solvent, and κ is the inverse Debye screening length given by

where z_i is the integer charge of the electrolyte species i, $cib$ is its bulk concentration, e is the elementary charge, k_B is the Boltzmann constant, and T is the absolute temperature. As in Ref. 37, we calculate the effect of the mobile ions on the free energy of solvation,

for three dielectric constants (ϵ_r = 4, 20, 80) and Debye screening lengths (λ = 25, 5, and 3 Å), with corresponding bulk concentrations given in Table I.

The spherical cavity has an additional parameter, the interpolation length l that defines a range of smooth interpolation between the vacuum and solvent permittivity by means of a hyperbolic tangent function⁵⁰ and is hence similar to the Δ parameter of the error function in Eq. (3) that we will use for all molecular solutes. To obtain converged solvation free energies, we solved the linear PB equation for cubic boxes of eight different sizes ranging from 15 Å to 50 Å, adjusting the grid points to a uniform grid spacing of l/3. We chose the grid spacing after a set of preliminary test calculations where this spacing proved to provide a proper interpolation in the transition region between vacuum and solvent permittivity. The solute ion was placed in the center of these boxes. Charge neutrality—cancellation of the integrated charge of the electrolyte ions and the solute ion—is not guaranteed for finite box sizes but serves as a criterion for extrapolation. We extrapolate the solvation free energy to that of a hypothetical box where charge neutrality is achieved. The explicit extrapolation scheme differs slightly for each Debye length and is detailed in the supplementary material.

The results are summarized in Table I. All in all, the mean absolute error (MAE) decreases with decreasing interpolation length l as expected, i.e., when the numerical model approaches the analytical model where no interpolation region is necessary. A charged molecule in a solvent of low dielectric permittivity combined with a low concentration of mobile ions is certainly the most challenging case for our implementation and extrapolation procedure. Consequently, relative errors of almost 25% are observed in these cases. We note here that a spherical, charged solute is the worst case for our Cartesian grid with necessarily finite simulation box size and we expect the results obtained in this section to be an upper limit to the inaccuracies we expect for our implementation.

Another model system that is more comparable to the uncharged solutes that we discuss in the remainder of this article is a point dipole in the spherical cavity that can analytically be described by the Kirkwood–Onsager model. In the original article by Onsager,⁶⁰ only a point dipole in a spherical cavity surrounded by a dielectric was considered but Kirkwood⁶¹ extended this model to multipoles in electrolyte solution and provided analytical expressions for the solvation free energy. The latter can be written as a sum of an electrolyte independent term ΔG₀ and an electrolyte dependent term ΔG(κ),

For a pure point dipole, these terms read

which is the Onsager result with the electric dipole moment μ, and

We model the point dipole by two charges of ±8e separated by d = 0.1 Å such that a ≫ d and μ = 1.512 D. Since the solute does not contain charges, charge neutrality is fulfilled even for small box sizes and we can calculate the free energy of solvation directly from the largest box size considered (30 Å) without extrapolation. If no convergence was achieved with the standard damping parameters introduced above, results were taken from calculations with smaller box size (20 Å or 15 Å).

Table II summarizes results for all combinations of the parameters ϵ = {20, 80}, c^b = {0.5, 1.0} mol/l, and a = {1, 2} Å. Again, the calculations converge to the analytic solution with decreasing interpolation length, which brings the numerical model closer to the analytical one where such a region is absent. Certainly, this comes at the price of a denser grid which can be optimized in future work by a more variable grid that is very fine in this interpolation region and more coarse everywhere else. The need for a dense grid is more pronounced for solvents with a large dielectric permittivity because the electrolyte is more concentrated close to the boundary region between solute and solvent. From the results for the Kirkwood–Onsager model, we conclude that in our implementation the numerical accuracy can be as accurate as 0.01 kcal/mol without the need for extrapolation.

We have examined the numerical accuracy we can expect from our implementation on two cases for which analytical solutions are known and will ignore the linear PB approximation from now on and exclusively present results that were obtained with the full nonlinear PB model.

The PB solvent models introduced above include a number of parameters that need to be carefully adjusted in order to obtain accurate results. These are the atom-specific radii d_α that define the solute cavity, the width parameter of the error function Δ, the finite ion radii R, and the Stern-layer thickness a. Motivated by the result of Sec. IV A that the interpolation width is not a crucial parameter as long as it is not too small, we keep Δ = 0.265 Å as in previous work.³⁹ We use experimental data on hydrated ions to set the finite ion radii R. Since we are only concerned with sodium chloride solutions in this study, we employ an ion size of R = 4.3 Å, which is the average of the hydrated sodium ion size (R = 4.7 Å)⁶² and that of the hydrated chloride ion (R = 3.9) Å.⁶² 

As is the case for other implicit solvation models such as the polarizable continuum model (PCM), the proper selection of atom-specific radii d_α to construct the solute cavity is crucial. This cavity determines the fraction of space where the solvent, characterized solely by its dielectric permittivity [Eq. (2)], resides, and, after inclusion of the additional Stern-layer parameter a, also the location of the mobile ions [Eq. (5)]. Two choices are commonly employed for these atom-specific radii d_α: van der Waals radii scaled by a factor of 1.2 (sVDW)^5,63,64 and the solvent-accessible-surface (SAS),^7,65 which is constructed by adding a solvent-dependent probe radius to the van der Waals radius,

where in the case of water $dαprobe=1.4$ Å, being half the distance to the maximum of the O-O radial distribution function of pure water.^66,67 The cavity created by the solvent-accessible surface is therefore usually larger than that of the scaled van der Waals radii which affects the total free energy of solvation but also its dependence on the electrolyte concentration ΔΔG_ion. The results in Fig. 2 exemplify the differences observed for the two sets of atom-specific radii on the example of the acetic acid molecule. As will be discussed in detail below, a linear dependence of the solvation free energy with the electrolyte concentration is observed experimentally, with the slope being described by the so-called Sechenov coefficient. This linear behavior is observed for the standard nonlinear PB equation, but finite ion size annihilates the linearity. Following previous authors,⁴³ we concentrate on the standard PB equation in the following when determining a suitable parameter for the Stern-layer thickness but will discuss this in more detail in Sec. IV D. For both the SAS and the sVDW cavity, hardly any difference is observed between the results obtained from the full PB expression compared to that from the osmotic pressure term only. In the following, we will employ the SAS in all further calculations. We have now determined all parameters of our model but one, the Stern-layer thickness a, which is optimized in Sec. IV C.

Since there are scarce experimental data on the exact thickness of the Stern-layer, especially around molecules rather than charged surfaces, we optimize this parameter in order to obtain the best possible agreement of our implicit solvent model with reliable experimental results.⁴³ The experimental data we can compare to are the Sechenov coefficients k_s introduced above^43,68,69 that are directly calculated from the slope of the free energy of solvation with the electrolyte concentration. The relation is simply

Experimental data frequently show a slight deviation from these linear data that are accounted for by an additional quadratic term.⁷⁰ However, the corresponding second-order Sechenov coefficient is usually two orders of magnitude smaller and is therefore of minor importance for electrolyte concentrations in the range of up to 2 mol/l. Additionally, the temperature dependence of the Sechenov coefficients is not accounted for in Eq. (41), as the experimental data are obtained at around room temperature and the temperature dependence of the Sechenov coefficients is negligible between 5 and 30 °C.⁷¹ We discussed already in Sec. IV B that this linear relation is only observed for the PB equation with point charges, whereas finite ion size leads to more pronounced salting-out effects for larger electrolyte concentrations.

In order to find an optimal Stern-layer thickness parameter a, we calculated the Sechenov coefficients for 39 molecules for which experimental data were available in Ref. 72. This set of molecule is our training set, whereas the 43 molecules from Ref. 73 serve as our validation set. The optimized molecular structures were taken from Ref. 74.

While it is a general shortcoming of implicit solvent models that they cannot capture explicit solute-solvent interaction, this is especially severe when electrolyte solutions are considered. Ion-pairing with electrolyte ions in zwitterionic molecules, structural changes upon cation-π interaction, and strong bonding to very polar groups are either not at all or inadequately described. Since we calculate the free energy of solvation for isolated molecules, any kind of solute-solute interaction is also neglected. We therefore expect to observe a number of outliers in our training set (and, less problematic, validation set) that need to be identified and excluded to not bias the optimization procedure. The random sample concensus (RANSAC) method⁷⁵ allows us to identify outliers and exclude them from the linear regression. We used the scikit-learn implementation⁷⁶ of RANSAC in python with 1000 trials in order to ensure that all combinations of two data points required for the linear regression model were included, effectively eliminating the randomness of the sampling for our rather small dataset. The residual threshold that defines the maximum residual for a data point to be counted as an inlier is another crucial parameter of the algorithm. We decided to scan several values for the threshold parameter (0.1, 0.075, 0.05, 0.025, 0.01, 0.005, and 0.001), where the loose thresholds usually identified few or no outliers, whereas for the tightest thresholds, almost all data points were identified as outliers. Figure 3(a) shows an example of the linear regression for the Stern-layer thickness a = 0.5 Å. The black line indicates perfect correlation between experimentally determined and calculated Sechenov coefficients, the light blue line is the result of an ordinary least squares fit (OLS), and the dark blue line results from the RANSAC algorithm with the threshold set to 0.075. Clearly, the identification of four outliers significantly increases the coefficient of determination (RANSAC: R² = 0.79 and OLS: R² = 0.48) and is also closer to ideal correlation.

After some initial trial calculations, we decided on two criteria for the final RANSAC threshold selection (and, hence, the outlier detection):

1. The maximum number of outliers should not exceed one third of the full training set to allow for a meaningful fit.

2. A new set of outliers identified with a tighter threshold is only accepted if accompanied by a substantial increase in the coefficient of determination R². This is to avoid a large number of outliers being identified with only a marginal increase in the fit quality. We required ΔR² = 0.02 per newly identified outlier.

These criteria were applied for calculated sets of Sechenov coefficients for eight Stern-layer thicknesses (a = 1.4, 1.2, 1.0, 0.8, 0.6, 0.5, 0.4, 0.3 Å). In all cases, the same four outliers were detected: bisphenol A, caffeine, 4-nitroanisole, and 4-nitroaniline. Bisphenol A was identified as an outlier also in a previous study,⁷⁴ where the deviation was attributed to an explicit interaction of the π-system with a sodium ion that induces structural changes. In the case of caffeine, a (partial) dimerization might cause the large deviation between the experimental and calculated Sechenov coefficient.⁷⁷ Why the aromatic nitrocompounds are outliers is difficult to speculate, especially since the agreement is much better for the two nitroalkanes: 1-nitropentane and 1-nitrohexane.

We identified the optimal Stern-layer thickness parameter a_opt, corresponding to the minimal root mean squared error (RMSE), by fitting a polynomial of fourth order to all data points obtained from the calculations on the training set with varying a [see Fig. 3(b)]. For the curve (red) without outlier detection—the OLS results—no minimum can be observed, whereas for the curve with outlier detection (blue)—the RANSAC results—a minimum at a_opt = 0.44 Å can be identified, corresponding to an RMSE of 0.027 l mol⁻¹. We hence recommend the thus optimized Stern-layer thickness in conjunction with the SAS for sodium chloride electrolyte solutions but note that the experimental data on the Sechenov coefficients are only of ml mol⁻¹ accuracy such that an RMSE that varies less than this contains no significant information. The range of 0.4 < a < 0.5 Å is therefore of indistinguishable accuracy and an equally suitable choice.

In Fig. 4, we show the results obtained with the PB model employing the optimized Stern-layer thickness for the training (green markers) and validation set (purple markers). The black line again corresponds to perfect correlation between theory and experiment, and the blue area indicates a deviation band of 0.1 l/mol. The RMSE for the validation set alone is 0.047 l mol⁻¹, and for the combined training (without outliers) and validation data it amounts to 0.039 l mol⁻¹.

In order to assess the quality of this result with other models, we rely on previous studies that certainly include different datasets but have substantial overlap with our data in all cases. The COSMO-RS model^78,79 shows a very high RMSE of 0.315 l mol⁻¹,⁷⁴ which is, however, due to a systematic overestimation because it can be reduced to 0.050 l mol⁻¹ when the empirical correction

is applied.⁷² The most meaningful comparison, however, is certainly the work of Ringe et al.,⁴³ since they use a similar PB model with a slightly different parameterization. Their calculations with sodium chloride as an electrolyte result in an RMSE of 0.068 l mol⁻¹, but it has to be noted that they do not employ any kind of outlier detection. Excluding outliers will certainly reduce the RMSE significantly and bring their accuracy in the realm of our results. To put these results into perspective, we compare to the polyparameter linear free energy relationships (pp-LFER) model.^80,81 These highly parameterized calculations result in an RMSE of 0.047 l mol⁻¹ and are therefore comparable in accuracy to our model and that of Ringe et al. but of no use for quantum-chemical calculations because the Sechenov coefficients are calculated directly without the electrostatic potential.

For the optimization of the Stern-layer thickness parameter in Sec. IV C, we compromised on the PB model with point-like ions instead of considering the ion size or the influence of the hydration shell. The linear dependence of the solvation free energy with the electrolyte dependence was only retained in this approximation and allowed us to fit to experimental data. This, however, is a disappointing result from a conceptual perspective: a physically more involved model is abandoned in favor of a much more coarse approximation because experimental relations are not reproduced by the more elaborate model. This nonlinear behavior cannot be explained by a quadratic term in Eq. (41) because experimentally determined quadratic Sechenov coefficients are usually orders of magnitude smaller than the linear one and might even have a negative sign.⁷⁰ In this section, we argue that this lack of agreement is likely due to a neglect of another effect: the dependence of the Stern-layer thickness on the electrolyte concentration.

In a recent study, the dependence of the Stern-layer thickness on the electrolyte concentration was probed by X-ray photoelectron spectroscopy on colloidal silica nanoparticles in sodium chloride solution.⁸² The authors observed a steep decrease in the Stern-layer thickness at low concentrations, followed by a domain where it decreases linearly. Taking into account both an electrolyte-concentration dependent Stern-layer thickness and finite size for the electrolyte ions in the MPB model can then yet again result in the linear Sechenov relation with the added benefit of less drastic approximations.

In Fig. 5(a), we show the results of MPB calculations on 2-octanone for seven electrolyte concentrations and eight Stern-layer thicknesses between 0.1 Å and 0.8 Å. The data points are fitted to a polynomial for proper interpolation (details on the polynomial fit are described in the supplementary material). The corresponding experimental curve (black line) obtained from the measured Sechenov coefficient is also shown as the points (red) where the calculated lines intersect with the experimental one. Assuming an otherwise adequate model—most importantly an accurate parameter for the finite ion size—these points mark the correct Stern-layer thickness for a given electrolyte concentration. We then plot the Stern-layer thicknesses of these points as a function of the electrolyte concentration [see Fig. 5(b)] and obtain a very similar pattern as the experimentally observed one. At very low concentrations, the concept of a Stern-layer is invalid, but in the physiological saline (∼0.15 mol l⁻¹) and seawater regime (∼0.6 mol l⁻¹), a linear dependence is observed. For even higher concentrations, the linear curve flattens out, which is to be expected because negative values for the Stern-layer thickness would be unphysical. In conclusion, our results obtained from the comparison of the MPB calculations with the experimental curve are in line with the experimentally observed electrolyte concentration dependence of the Stern-layer thickness.

The slopes (m) and intercepts (a_max) of these diagrams are compared for five molecules in Table III. While the slopes roughly correlate with the intercepts as might be expected, a correlation with the Sechenov constant or the size of the molecules is not evident.

The analysis of this section offers an explanation for the nonlinear behavior of the MPB model. In practical calculations however, the application of the MPB model is hampered by the additional parameters that are hard to be determined from first principles. These parameters are the electrolyte concentration dependence of the Stern-layer thickness and the finite ion size (pure ions or inclusion of some hydration). If one of these parameters can be deduced from either experimental data (as for the Stern-layer thickness of the silica nanoparticles) or from MD simulations, it will be possible to refit the remaining parameter with the aid of the experimental Sechenov coefficients to arrive at a more accurate implicit solvation model. For now, however, in the absence of either experimental or MD data, we recommend the standard PB model with the optimized parameters described above.

In this article, we introduced our multigrid implementation of an implicit solvation model that includes the effect of electrolytes. We describe the underlying Poisson–Boltzmann model with a hierarchy of approximations, namely, linear vs nonlinear Poisson–Boltzmann equations, inclusion of a Stern-layer correction and finite ion-size corrections. Our derivation of the free energy of solvation including the effect of the electrolyte is based on the notion that these equations can be recast into an Euler-Lagrange equation.¹³ We then challenged our implementation of the linear PB model by comparing to analytical solutions for the Born ion model and the Kirkwood–Onsager model and observed that the relative error for polar solvents and moderate to high salt concentrations is below 20%. We then discussed the parameters of the model and optimized the Stern-layer thickness a based on a fit to an experimentally observed linear relationship between the free energy of solvation and the electrolyte ion concentration, which is described by so-called Sechenov coefficients. For a cavity constructed as a solvent-accessible surface with a probe radius of 1.4 Å for water, we arrived at a_opt = 0.44 Å, by fitting to experimental data for 39 molecules and identifying outliers with the RANSAC algorithm. The RMSE of 0.039 l/mol that we obtained by including results for a validation set of 43 molecules is comparable to that of previous work by Ringe et al. and as good as that obtained with the highly parameterized pp-LFER model. In the final section, we attributed the failure of the size-modified PB model to a neglect of the electrolyte concentration dependence of the Stern-layer thickness. We then demonstrated that in a certain concentration range, knowledge about this linear dependence can restore the experimentally observed linearity of the solvation free energy with electrolyte concentration. We argued that this knowledge might be obtained either from an experiment or from MD simulations, where the latter can also be informative on the finite ion size, i.e., the size of the hydration shell around the electrolyte ions.

To our knowledge, this is the first implementation of an implicit solvation model for electrolyte solutions in a Gaussian orbital based multipurpose quantum chemistry program. A related implementation²³ based on numeric atom-centered orbitals can be found in the FHI-aims⁸³ package and for plane-wave basis sets, for example, in VASPsol (only linearized PB).³³ 

This enables us to investigate electrochemical processes for molecules or small cluster models with the inclusion of the effect of electrolyte solutions. Future work will therefore focus on the implementation of nuclear gradients for the grid-based PB model so that solvation induced structural changes can be analyzed in addition to the change in the solvation free energy. This will then also allow us to systematically study the effect of the linear vs nonlinear PB model and the influence of the finite ion size on electrocatalytic reactions.

See the supplementary material for a table demonstrating the convergence with grid density and examples of the extrapolation schemes for the Debye–Hückel model calculations.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing, and Office of Basic Energy Sciences, via the Scientific Discovery through Advanced Computing (SciDAC) program.

C.J.S. gratefully acknowledges financial support in the form of an Early Postdoc.Mobility fellowship from the Swiss National Science Foundation.

Work by J.M.H. on the Poisson solver was supported by the National Science Foundation (Grant No. CHE-1665322).

The authors declare the following competing financial interests: J.M.H. and M.H.-G. are part owners of Q-Chem, Inc.