#  Dielectric continuum methods for quantum chemistry

John M<EMAIL_ADDRESS>
Department of Chemistry and Biochemistry, The Ohio State University, Columbus,
Ohio, USA

###### Abstract

This review describes the theory and implementation of implicit solvation
models based on continuum electrostatics. Within quantum chemistry this
formalism is sometimes synonymous with the polarizable continuum model, a
particular boundary-element approach to the problem defined by the Poisson or
Poisson-Boltzmann equation, but that moniker belies the diversity of available
methods. This work reviews the current state-of-the art, with emphasis on
theory and methods rather than applications. The basics of continuum
electrostatics are described, including the nonequilibrium polarization
response upon excitation or ionization of the solute. Nonelectrostatic
interactions, which must be included in the model in order to obtain accurate
solvation energies, are described as well. Numerical techniques for
implementing the equations are discussed, including linear-scaling algorithms
that can be used in classical or mixed quantum/classical biomolecular
electrostatics calculations. Anisotropic models that can describe interfacial
solvation are briefly described.

###### Contents

  1. 1 Overview
  2. 2 Continuum Electrostatics
    1. 2.1 Poisson’s equation
    2. 2.2 Solute cavity
    3. 2.3 Polarizable continuum models
    4. 2.4 Comparison of ASC-PCMs
  3. 3 Implementation
    1. 3.1 Matrix equations for PCMs
    2. 3.2 Discretization
    3. 3.3 Isodensity and self-consistent cavity surfaces
    4. 3.4 Linear-scaling algorithms
  4. 4 Solvation Energies
    1. 4.1 SM$x$ and SASA-based models
    2. 4.2 Physics-based models
  5. 5 Nonequilibrium Solvation
    1. 5.1 Conceptual overview
    2. 5.2 State-specific approach
    3. 5.3 Linear response approach
  6. 6 Anisotropic Solvation
  7. 7 Closing remarks
  8. Conflict of Interest

##  1 Overview

The use of dielectric continuum models in quantum chemistry dates to the
mid-1970s,1, 2, 3, 4, 5, 6 to when the field itself was still in its infancy.
In their simplest form, these models describe the solvent in terms of a single
parameter $\varepsilon_{\text{s}}$, the (static) dielectric constant, a
dimensionless quantity equal to the electric permittivity relative to vacuum
and ranging (for simple liquids) from $\varepsilon_{\text{s}}\approx 2$ for
non-polar solvents such as benzene and hexane, up to
$\varepsilon_{\text{s}}=78$ for water and $\varepsilon_{\text{s}}=110$ for
formamide. This constant describes the solvent’s ability to screen a charge,
and the Coulomb interaction between charges $Q_{1}$ and $Q_{2}$ separated by a
distance $r$ is modified from $V(r)=Q_{1}Q_{2}/4\pi\varepsilon_{0}r$ in the
gas phase, to an attenuated version
$V(r)=Q_{1}Q_{2}/4\pi\varepsilon_{0}\varepsilon_{\text{s}}r$ within the
dielectric medium. The continuum description represents the ultimate in
coarse-graining, reducing the solvent to a single parameter, which has obvious
utility in a quantum chemistry calculation whose cost rises steeply with
system size. Within a continuum description, there is no need for sampling
over solvent degrees of freedom (e.g., solvent reorganization in response to
an electron transfer event that modifies the solute’s charge distribution),
because this averaging is implicitly encoded into the value of
$\varepsilon_{\text{s}}$. While advantageous from the standpoint of cost,
limitations of the continuum description are equally apparent: “specific”
solvation effects such as hydrogen bonding are not captured, and dielectric
continuum theory alone does not describe nonelectrostatic interactions
including dispersion and Pauli repulsion. Without the latter, there is nothing
to imbue the molecules with finite size, necessitating ad hoc introduction of
a “solute cavity” to define the interface between the atomistic solute and the
continuum solvent, as depicted in Fig. 1.

Fig. 1:  (a) Zwitterionic tautomer of glycine (-O2CCH2NH${}_{3}^{+}$) in a
molecular cavity constructed from atom-centered spheres. Coloring reflects the
sign and magnitude of the molecular electrostatic potential evaluated at the
cavity surface, $\varphi^{\rho}(\mathbf{s})$ for $\mathbf{s}\in\bm{\Gamma}$.
(b) Schematic illustration of the same molecular cavity (in green) embedded in
a dielectric medium (in blue), illustrating how the continuum polarizes in
response to the solute’s electrostatic potential. The orange probe sphere
illustrates how the atomic radii that define the van der Waals (vdW) cavity
surface might be augmented to afford a “solvent-accessible surface” (SAS). The
region interior to the solute cavity is designated as $\bm{\Omega}$, and for a
sharp dielectric interface one sets
$\varepsilon(\mathbf{r})\equiv\varepsilon_{\text{in}}$ for
$\mathbf{r}\in\bm{\Omega}$. If the solute is described using quantum chemistry
then the natural choice is $\varepsilon_{\text{in}}=1$. Outside of the cavity,
the permittivity function $\varepsilon(\mathbf{r})$ takes the value
$\varepsilon_{\text{out}}$, which is usually the static dielectric constant of
the solvent, $\varepsilon_{\text{s}}$. Panel (b) is adapted from Ref. 7;
copyright 2008 John Wiley & Sons.

Some of the aforementioned limitations can be overcome, in principle, by
admission of a small number of explicit solvent molecules into the atomistic
part of the calculation, in what is often called a “semicontinuum” or a
“cluster + continuum” approach.8 As such, the continuum description serves as
a flexible and highly useful starting point for the description of solvation
effects in quantum chemistry. Introduction of continuum solvation models into
quantum chemistry was led by the group of Tomasi in Pisa,9, 10, 11, 12, 13,
14, 15, 16 who refer to this approach as the polarizable continuum models
(PCM). That terminology will be introduced in Section 2.3, to refer to a
particular class of continuum solvation models that replace the three-
dimensional partial differential equations of continuum electrostatics with a
two-dimensional boundary element problem, defined on the surface of a cavity
representing the interface between atomistic solute and continuum solvent.
Although PCMs are ubiquitous in quantum chemistry, they are not the only
continuum solvation models that are used in the field and not the only ones
discussed herein. In any case, the combination of a quantum-mechanical (QM)
description of the solute sets up a self-consistent reaction-field (SCRF)
problem in which the solute’s charge distribution both polarizes, and is
polarized by, its environment. The two effects must be iterated to self-
consistency.

The remainder of this review is organized as follows. Section 2 provides the
elementary specification of the continuum electrostatics problem, starting
from the Poisson and Poisson-Boltzmann equations. The mechanics of turning
that formalism into a computationally-tractable model are discussed in Section
3, with an emphasis on the smooth discretization approach developed by this
author’s group.17, 18, 19, 20, 21 The focus here is on continuum solvation
models in quantum chemistry but the formalism in Sections 2 and 3 is perfectly
applicable to biomolecular implicit solvent calculations, in which a
macromolecular solute is described using a classical force field.19 Section 3
also introduces the various flavors of PCM that can be found in the literature
and discusses how they can be understood in one another’s context. It should
be noted that the solution of Poisson’s equation, or equivalently the apparent
surface charge PCM problem, specifies only the electrostatic contribution to
the solvation energy. Other contributions including cavitation, dispersion,
Pauli repulsion, and hydrogen bonding must be included in order to predict
free energies of solvation ($\Delta_{\text{solv}}\mathcal{G}$) in reasonable
agreement with experiment. Section 4 introduces models for these
nonelectrostatic contributions and provides an overview of the accuracy that
can be expected for $\Delta_{\text{solv}}\mathcal{G}$. Section 5 introduces
several “nonequilibrium” formulations of the continuum electrostatics problem,
which are used to describe the continuum’s response to a sudden change in the
solute’s charge density, as in photoexcitation or photoionization. This
provides the machinery to compute solvent effects on vertical excitation,
ionization, or fluorescence energies. Finally, Section 6 discusses
modifications to the isotropic continuum model that are necessary in order to
describe anisotropic solvation environments, such as the air/water interface
or the solid-state/aqueous interface.

This review is focused on the theoretical framework and computational
mechanics of continuum solvation models, not on applications. While some
limited data to describe the performance of the models are provided, for a
complete overview of continuum solvation methods in action, the reader is
directed to several general reviews,13, 14, 15, 22, 23, 24, 8 as well as
specialized ones concerning the application of PCMs to specific types of
spectroscopy.25, 26, 27, 28, 29, 30, 31 In lieu of a great deal of data, the
present work provides copious references to the primary literature.

##  2 Continuum Electrostatics

This section reviews the basic electrostatic formalism that underlies
continuum solvation theory. The physical model is defined by Poisson’s
equation in three-dimensional space (Section 2.1) but is not fully specified
without a surface to demarcate the boundary of the atomistic region (Section
2.2). The polarizable continuum models that were introduced into quantum
chemistry by Tomasi and co-workers9, 6, 13 are introduced in Section 2.3, and
stem from a reformulation of the Poisson problem as a boundary-element of
“apparent surface charge” problem. There are several variants, and these are
compared side-by-side in Section 2.4.

###  2.1 Poisson’s equation

The basic tenet of dielectric continuum theory is an assumption that the
electric response of matter can be coarse-grained in the form of a dipole
density $\mathbf{P}$ that defines the polarization of the medium. In the
presence of this dielectric medium, the role of the electric field
$\mathbf{E}(\mathbf{r})$ in vacuum is supplanted by the electric displacement
field (or electric induction),

$\mathbf{D}(\mathbf{r})=\mathbf{E}(\mathbf{r})+4\pi\mathbf{P}(\mathbf{r})=\varepsilon(\mathbf{r})\,\mathbf{E}(\mathbf{r})\;.$
(2.1)

The electric permittivity $\varepsilon(\mathbf{r})$ is defined by the manner
in which polarization $\mathbf{P}(\mathbf{r})$ is induced by the external
field $\mathbf{E}(\mathbf{r})$, and Eq. (2.1) amounts to the definition of a
linear dielectric material whose polarization is proportional to field
strength. Nonlinear susceptibilities are harder to describe within a continuum
formalism and have received less attention.13 Whereas a fully general
discussion of (linear) dielectric materials would allow for a permittivity
that is a function of frequency also (or even a nonlocal function of space
and/or time, in some formulations),32, 33, 34 the ground-state SCRF problem
does not require this generality, and unless otherwise specified,
$\varepsilon$ will mean the static (zero-frequency) dielectric constant,
$\varepsilon_{\text{s}}$. (The continuum electrostatics community has
stubbornly defied the suggestion35 that “dielectric constant” is obsolete and
should be replaced by “relative electric permittivity”.) For an anisotropic
medium, $\varepsilon$ would take the form of a $3\times 3$ tensor rather than
scalar, which could be used to model a liquid crystal in which the electric
susceptibility depends on the orientation of the applied field.36, 37 This
review will not consider such cases, although a different form of anisotropic
solvation is considered in Section 6. Herein, $\varepsilon$ is a scalar.

That said, Eq. (2.1) does express the permittivity as a scalar-valued
function, $\varepsilon(\mathbf{r})$, rather than simply a dielectric constant.
This allows for a situation such as that depicted in Fig. 1(b), wherein a
“solute cavity” (two-dimensional surface $\bm{\Gamma}$) defines an interface
between the continuum solvent and an atomistic region, the latter of which
shall be described using quantum chemistry. Within the cavity, inter-particle
Coulomb interactions are included explicitly in the Hamiltonian and thus
$\varepsilon=1$ in this region. Outside of the cavity, the permittivity
function takes on a value equal to the (static) dielectric constant of the
solvent, e.g., $\varepsilon=78$ for water. Given a charge density
$\rho(\mathbf{r})$ for the solute, including both nuclei and electrons,
Maxwell’s equation for the displacement field $\mathbf{D}(\mathbf{r})$ is

$\hat{\bm{\nabla}}\bm{\cdot}\mathbf{D}(\mathbf{r})=4\pi\rho(\mathbf{r})\;.$
(2.2)

Recognizing that the electric field is
$\mathbf{E}(\mathbf{r})=-\hat{\bm{\nabla}}\varphi(\mathbf{r})$, where
$\varphi(\mathbf{r})$ is the electrostatic potential, Eq. (2.2) can be written
in the more familiar form of Poisson’s equation,

$\hat{\bm{\nabla}}\bm{\cdot}\left[\varepsilon(\mathbf{r})\,\hat{\bm{\nabla}}\varphi(\mathbf{r})\right]=-4\pi\rho(\mathbf{r})\;.$
(2.3)

All of these equations are expressed in Gaussian electrostatic units, where
$4\pi\varepsilon_{0}=1$.38

Poisson’s equation is the mathematical starting point for continuum
electrostatics. Given the charge density $\rho(\mathbf{r})$ from an electronic
structure calculation, Eq. (2.3) is solved for $\varphi(\mathbf{r})$
throughout space, including both the atomistic region and the surrounding
dielectric medium. This potential can be separated into two parts,

$\varphi(\mathbf{r})=\varphi^{\rho}(\mathbf{r})+\varphi_{\text{rxn}}(\mathbf{r})\;,$
(2.4)

where the first term is the electrostatic potential generated by the solute’s
charge density,

$\varphi^{\rho}(\mathbf{r})=\int\frac{\rho(\mathbf{r}^{\prime})}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}\;d\mathbf{r}^{\prime}\;,$
(2.5)

whereas the “reaction field” contribution $\varphi_{\text{rxn}}(\mathbf{r})$
arises from polarization of the continuum, resulting in an additional charge
density $\rho_{\text{pol}}(\mathbf{r})$. Having obtained $\rho(\mathbf{r})$
from Schrödinger’s equation and then $\varphi(\mathbf{r})$ by solving Eq.
(2.3), the electrostatic solvation energy can be expressed variously as18, 39

$\mathcal{G}_{\text{elst}}=\frac{1}{2}\int\varphi_{\text{rxn}}(\mathbf{r})\;\rho(\mathbf{r})\;d\mathbf{r}=\frac{1}{2}\int\varphi^{\rho}(\mathbf{r})\;\rho_{\text{pol}}(\mathbf{r})\;d\mathbf{r}\;.$
(2.6)

This is sometimes alternatively called the polarization energy
($\mathcal{G}_{\text{pol}}$),18 but we will use the term electrostatic energy
($\mathcal{G}_{\text{elst}}$) as that is the nomenclature typically
encountered in the literature on dielectric materials.38 Unlike the theory of
intermolecular interactions,40 it makes little sense in the present context to
consider electrostatics separate from polarization but the reader may, if
desired, substitute the phrase “electrostatics + polarization” whenever
“electrostatics” is used herein. The quantity $\mathcal{G}_{\text{elst}}$ is a
free energy insofar as the dielectric formalism implicitly accounts for the
averaging over solvent degrees of freedom, and the factor of $1/2$ in Eq.
(2.6) reflects the fact that the interaction energy is reduced, by precisely
half its value, on account of the work required to polarize the
environment.41, 9, 42, 18

From the point of view of electronic structure theory,
$\mathcal{G}_{\text{elst}}[\varepsilon,\rho]$ is a functional of both the
permittivity function $\varepsilon(\mathbf{r})$ and the solute’s charge
density $\rho(\mathbf{r})$. The total (free) energy is

$\mathcal{G}_{0}[\Psi]=\big{\langle}\Psi\big{|}\hat{\mathcal{H}}_{\text{vac}}\big{|}\Psi\big{\rangle}+\mathcal{G}_{\text{elst}}[\varepsilon,\rho]\;,$
(2.7)

where the first term represents the electronic energy in vacuo. (It is written
here in a form that suggests wave function quantum mechanics but could be
replaced by a density functional.) This can equivalently be expressed as a
total energy functional

$\mathcal{G}_{0}[\Psi]=\big{\langle}\Psi\big{|}\hat{\mathcal{H}}_{\text{vac}}+\tfrac{1}{2}\hat{\mathcal{R}}_{0}\big{|}\Psi\big{\rangle}$
(2.8)

where $\hat{\mathcal{R}}_{0}$ is a reaction-field operator that generates the
integral in Eq. (2.6). Minimization of $\mathcal{G}_{0}[\Psi]$, or
$\mathcal{G}_{0}[\rho]$ in density functional theory (DFT), in conjunction
with Poisson’s equation to obtain the electrostatic potential that defines
$\hat{\mathcal{R}}_{0}$, defines the SCRF problem. Provided that the
electronic structure model satisfies a variational principle, as it does for
self-consistent field (SCF) models, then the total energy defined by Eq. (2.7)
satisfies a variational principle as well,43, 18

Equation (2.3) is sometimes called the “generalized” form of Poisson’s
equation, with the “ordinary” form being

$\varepsilon\,\hat{\nabla}^{2}\varphi(\mathbf{r})=-4\pi\rho(\mathbf{r})\;.$
(2.9)

The distinction is that the permittivity function $\varepsilon(\mathbf{r})$ in
Eq. (2.3) is replaced in Eq. (2.9) by a scalar, $\varepsilon$. It is the
latter equation that is often taken to define the continuum electrostatics
problem, but this requires additional specification because $\varepsilon=1$
for the atomistic interactions. Some sort of molecular surface is needed to
delineate the boundary with the continuum, as shown in Fig. 1 where the cavity
is defined by a union of atom-centered spheres. Given a cavity surface, Eq.
(2.9) is shorthand for Eq. (2.3) with the permittivity function

$\varepsilon(\mathbf{r})=\begin{cases}\varepsilon_{\text{in}},&\mathbf{r}\in\bm{\Omega}\\\
\varepsilon_{\text{out}},&\mathbf{r}\notin\bm{\Omega}\end{cases}\;.$ (2.10)

Note that $\varphi(\mathbf{r})$ is continuous across the cavity surface but
its derivative is not.9 In general it makes sense to take
$\varepsilon_{\text{in}}=1$ although larger values (commonly
$\varepsilon_{\text{in}}=2$–$4$,44, 7, 45 but in some cases
$\varepsilon_{\text{in}}=10$–$20$46, 47, 48, 45) have sometimes been used in
biomolecular electrostatics calculations, in an effort to approximate a
“dielectric constant of protein”. (However, the very concept that such a
“constant” exists has been vociferously criticized.49, 44, 50, 51, 52) If the
solute is described using quantum mechanics, however, then any choice other
than $\varepsilon_{\text{in}}=1$ represents an inconsistent treatment of the
Coulomb interactions unless the Coulomb operators that define
$\hat{\mathcal{H}}_{\text{vac}}$ are modified, which is seldom done.

Although we have introduced it as a quantum-mechanical SCRF problem that must
be solved in tandem with the electronic structure problem, Eq. (2.3) is also
used in classical electrostatics calculations for biomolecules,53, 54, 55, 56
where $\rho(\mathbf{r})$ is comprised of classical point charges (or higher-
order multipoles57, 58) that come from a force field, e.g.,

$\rho(\mathbf{r})=\sum_{A}^{\text{atoms}}Q_{A}\,\delta(\mathbf{r}-\mathbf{R}_{A})\;.$
(2.11)

For biomolecular applications, one is often interested in an aqueous solvent
containing some concentration of dissolved ions. The continuum analogue of
that situation is described by the Poisson-Boltzmann equation,59, 60, 61, 62,
7

$\hat{\bm{\nabla}}\bm{\cdot}\left[\varepsilon(\mathbf{r})\,\hat{\bm{\nabla}}\varphi(\mathbf{r})\right]=-4\pi\big{[}\rho(\mathbf{r})+\rho_{\text{ions}}(\mathbf{r})\big{]}\;,$
(2.12)

in which the right side of Eq. (2.3) has been augmented with a term that
accounts for a thermal distribution of dissolved ions.60, 7 Whereas the
solute’s charge density $\rho(\mathbf{r})$ reflects atomistic modeling, using
either a classical force field or else an electronic structure calculation,
the density $\rho_{\text{ions}}(\mathbf{r})$ of “mobile” ions is treated
statistically.63, 60 For an electrolyte with dissolved ion concentrations
$\\{c_{i}\\}$, whose ionic charges are denoted $\\{q_{i}\\}$, the statistical
charge density for the ions is7

$\rho_{\text{ions}}(\mathbf{r})=\sum_{i}^{\text{ions}}q_{i}c_{i}\lambda_{i}(\mathbf{r})\exp\left(\frac{-q_{i}\varphi(\mathbf{r})}{k_{\text{B}}T}\right)\;.$
(2.13)

The ion accessibility function $\lambda_{i}(\mathbf{r})$ represents some type
of step function that serves to exclude the ions from the atomistic region.
With this form for the density of mobile ions, Eq. (2.12) is sometimes known
as the size-modified version of the (nonlinear) Poisson-Boltzmann equation.7,
64 In the case of a 1:1 electrolyte with monovalent ions ($q_{1}=e=-q_{2}$),
Eq. (2.13) reduces to

$\rho_{\text{ions}}(\mathbf{r})=-2c\lambda(\mathbf{r})\,\text{sinh}\left(\frac{e\,\varphi(\mathbf{r})}{k_{\text{B}}T}\right)\;.$
(2.14)

At the physiological ionic strengths, the hyperbolic sine function can be
linearized without significant error,65, 66 resulting in a linearized Poisson-
Boltzmann equation,62, 7

$\hat{\bm{\nabla}}\bm{\cdot}\left[\varepsilon(\mathbf{r})\,\hat{\bm{\nabla}}\varphi(\mathbf{r})\right]=-4\pi\rho(\mathbf{r})+\kappa^{2}\lambda(\mathbf{r})\,\varphi(\mathbf{r})\;,$
(2.15)

where

$\kappa=\left(\frac{8\pi e^{2}c}{k_{\text{B}}T}\right)^{1/2}\;.$ (2.16)

The dissolved ions screen electrostatic interactions over a length scale
$\sim\kappa^{-1}$ known as the Debye screening length, such that the potential
that appears in Debye-Hückel theory is the attenuated Coulomb potential
$e^{-\kappa r}/(\varepsilon_{\text{s}}r)$.67, 60, 7, 68

Within the biomolecular electrostatics community there has been significant
discussion regarding the accuracy of the linearization approximation, with
studies noting that the nonlinear form [Eq. (2.12), with Eq. (2.13) for
$\rho_{\text{ions}}(\mathbf{r})$] affords better agreement with explicit
solvent simulations when the ionic strength is high.66, 69 Deficiencies in the
Poisson-Boltzmann model itself—even even in its nonlinear form and especially
for polyvalent ions—have also been pointed out.70 It is therefore worth noting
that for the small solutes that characterize most quantum chemistry
applications, the effect of the mobile ions on $\mathcal{G}_{\text{elst}}$ is
quite modest,68, 64 although there are effects on activity coefficients.64, 71
These effects are presumably magnified for a solute the size of a protein, but
the intermediate size regime has hardly been explored.

Methods for solving the partial differential equations introduced in this
section will be introduced below. First, however, we need to discuss one more
aspect of the model problem itself, namely, the definition a surface that
establishes the boundary between atomistic solute and continuum solvent.

###  2.2 Solute cavity

For the case of a sharp dielectric boundary, Eq. (2.3) has an analytic
solution if the cavity surface is spherical and contains the entire charge
density $\rho(\mathbf{r})$. For a solute consisting of a single point charge,
$Q$, centered in a spherical cavity of radius $\bar{R}$ in a medium with
dielectric constant $\varepsilon$, this solution affords the well-known “Born
model” for ion solvation,72, 73

$\Delta\mathcal{G}_{Q}=-\frac{Q^{2}}{2\bar{R}}\left(\frac{\varepsilon-1}{\varepsilon}\right)\;.$
(2.17)

Here, $\Delta\mathcal{G}$ indicates the change in $\mathcal{G}_{\text{elst}}$
from its gas-phase value of zero to the solution-phase value obtained from Eq.
(2.6). Replacing the point charge $Q$ by a point dipole $\mu$, the solvation
energy is

$\Delta\mathcal{G}_{\mu}=-\frac{(\varepsilon-1)\mu^{2}}{(2\varepsilon+1)\bar{R}^{3}}\;.$
(2.18)

(This result is often attributed to Onsager,74 although it was derived a few
years earlier by Bell.75) The dipole solvation energy can alternatively be
written
$\Delta\mathcal{G}_{\mu}=-\tfrac{1}{2}(\bm{\mu}\bm{\cdot}\mathbf{E}_{\text{rxn}})$,
where

$\mathbf{E}_{\text{rxn}}=\underbrace{\frac{1}{\bar{R}^{3}}\left(\frac{2(\varepsilon-1)}{2\varepsilon+1}\right)}_{g_{1}(\varepsilon,\bar{R})}\bm{\mu}=g_{1}(\varepsilon,\bar{R})\,\bm{\mu}$
(2.19)

is the “reaction field” induced by $\bm{\mu}$.41 If $\bm{\mu}$ is polarizable
($\bm{\mu}=\bm{\mu}_{0}+\bm{\alpha}\bm{\cdot}\mathbf{E}$), then Eq. (2.19)
provides the earliest example of a SCRF model.41 Historically, this model was
used to formulate a microscopic theory for the dielectric constant of a polar
liquid,74, 76, 77 going beyond the mean-field description contained in the
Clausius-Mossotti equation.78, 79 These attempts were not particularly
quantitative,77 and for modern purposes a dipolar description of the solute
constitutes a needless approximation.

The aforementioned results for $\Delta\mathcal{G}$ were quickly generalized by
Kirkwood,80, 81, 82 and later by others,83, 84, 85, 86, 87, 57 to an arbitrary
multipole in a spherical cavity, with the general result

$\mathcal{G}_{\text{elst}}=-\frac{1}{2}\sum_{\ell\geq
0}\underbrace{\left[\frac{(\ell+1)(\varepsilon-1)}{\ell+(\ell+1)\varepsilon}\right]\frac{1}{\bar{R}^{2\ell+1}}}_{g_{\ell}(\varepsilon,\bar{R})}\sum_{m=-\ell}^{\ell}\big{|}\big{\langle}\Psi\big{|}\hat{T}_{\ell
m}\big{|}\Psi\big{\rangle}\big{|}^{2}$ (2.20)

where the $\hat{T}_{\ell m}$ are spherical multipole operators.88, 89 Since
then, analytic results for ellipsoidal cavities,82, 84 off-center charges,81
and higher-order off-center multipoles87 have also been derived, along with
formulas for the interaction of multipoles contained in disjoint spheres with
a dielectric medium in between.90 Insofar as an arbitrary charge distribution
$\rho(\mathbf{r})$ has a single-center multipole expansion, these results
provide an analytic solution for a charge density of arbitrary complexity in a
spherical or ellipsoidal cavity, assuming there is no penetration of
$\rho(\mathbf{r})$ into the continuum region. (The latter effect, known as
volume polarization,91, 92, 93, 94, 95, 96, 97 is discussed below.) The
multipole expansion method has been called a “generalized Kirkwood” solvation
model,57 and using multipolar formulas for $\hat{\mathcal{R}}_{0}$ in Eq.
(2.8), it can be turned into a generalized Kirkwood SCRF for quantum
chemistry. Multipolar methods are reviewed elsewhere,9, 98 but have
effectively been rendered obsolete by the polarizable continuum methods
described in Section 2.3. In the absence of volume polarization, these
approaches also afford an exact (albeit numerical) solution to the continuum
electrostatics problem, yet can used in conjunction with molecule-shaped
cavities.

A spherical cavity may make sense if a large number of explicit solvent
molecules are included in the atomistic solute region, and such approaches are
known as solvent boundary potential methods.99, 100, 101 These have been
developed as replacements for periodic boundary conditions in both QM/MM
simulations102, 103, 104, 105, 106, 107, 108 and in fully quantum-mechanical
calculations.109, 110, 111, 112 Spherical cavities make little sense for
single molecules, however, and it is clear from Eq. (2.17) that
$\Delta\mathcal{G}$ will be quite sensitive to the cavity radius. This is
especially problematic if that cavity does not correspond closely to the size
of the molecule itself, and results in Section 3.3 demonstrate that even for
realistic “molecule-shaped” cavities, solvation energies are quite sensitive
to cavity construction.

Fig. 2:  Illustration of cavity surface definitions using a set of atomic
spheres (in gray) whose envelope defines the vdW surface (in black). The
solvent-accessible surface (in green) is defined either by augmenting the
atomic radii by a probe radius ($R_{\text{probe}}$, taken here to be smaller
than the vdW radii), or equivalently as the center point of the probe sphere
as it rolls over the vdW surface. The solvent-excluded surface (SES, in
magenta, sometimes called the “molecular surface”) is traced out by arcs of
the probe that connect points of contact between the probe and the vdW
surface. The SES eliminates the cusps in the vdW surface that otherwise appear
along seams of intersection between atom-centered spheres.

Examples of some molecule-shaped cavity constructions are provided in Fig. 2.
A simple union of atom-centered spheres is often called a van der Waals (vdW)
cavity surface. The atomic radii might simply be empirical parameters of the
model,113 or alternatively covalent radii extracted from crystal
structures,114, 115 or based on calculations of atomic size.116 In the latter
two cases these radii represent close-contact distances, meaning that an
implicit or continuum solvent should not be allowed to approach all the way to
the vdW radii of the solute atoms. This exclusion effect has been treated in
several ways, most commonly by scaling the atomic vdW radii by a factor
$\alpha>1$. A scaling factor $\alpha=1.2$ was chosen early on,117 and factors
$\alpha\approx 1.1$–$1.2$ have since become standard choices, albeit with
little theoretical justification to choose one value over another within a
modest range. Alternatively, and with somewhat better justification, the
atomic vdW radii can be augmented by a “probe radius” representing the assumed
size of a solute molecule, which can be extracted from the liquid structure of
the solvent. For example, $R_{\text{probe}}=1.4$ Å is often used for water,
representing half the distance to the first peak in the oxygen–oxygen radial
distribution function.118 Values for water as small as $R_{\text{probe}}=0.2$
Å have also been used, however, in an effort to match solvation energy
benchmarks from simulations with explicit solvent.119 The cavity surface
generated using atomic radii $R_{\text{vdW}}+R_{\text{probe}}$ is known as the
solvent-accessible surface (SAS), which was first introduced as a means to
measure accessible surface area in proteins and is often used to define the
ion accessibility function $\lambda(\mathbf{r})$ in Eq. (2.15). Note that it
does not make sense to augment a scaled vdW radius with a probe radius, as
this double counts the size of the exclusion layer.

Both the vdW surface and the SAS exhibit cusps where atomic spheres intersect.
These cusps are eliminated by solvent-excluded surface (SES) that is generated
by the probe sphere as it rolls over the vdW surface; see Fig. 2. (In
principle, this procedure could be applied to eliminate cusps in the SAS as
well, but that is less often done because cusps in the SAS are less
problematic, numerically speaking.21) The center of the probe sphere traces
out the SES while its points of contact with the vdW surface, combined with
the concave arcs of the probe sphere that smooth over the cusps, constitute
the SES. (The SES is also known as the Connolly surface,120 or sometimes
simply the “molecular surface”,121 as it is probably the closest approximation
to the true shape of the molecule. However, these names have sometimes been
used interchangeably or ambiguously in the literature, and in particular the
term “Lee-Richards surface”122 has been used to mean either the SAS and the
SES.123, 119, 21 (The notation used here is standard in the quantum chemistry
literature.124) The SES has an analytic construction,120, 21 although it has
most often been constructed numerically, for visualization purposes.

To a greater or lesser degree each of these cavity definitions seems
physically reasonable, but beyond that there is little theoretical
justification for any of them, and no guarantee that small changes in one or
more atomic radii will not have a significant impact on computed
observables.125, 97 It has been suggested that the “optimal” atomic radius for
a given atom likely ought to vary as a function of its partial atomic
charge,126, 127 but this effect is generally handled empirically, if at all.
“United-atom” models, in which hydrogen atoms are not given explicit spheres
and the radii of the other atoms are increased to compensate,126 are also in
widespread use. A more justifiable definition of the solute cavity surface
might involve an isocontour of the molecular electron density
$\rho(\mathbf{r})$,128, 94, 129, 130, 96 as it is the density that is
ultimately responsible for molecular size and shape.131 Isodensity cavity
definitions are considered in Section 3.3.

With the introduction of a molecule-shaped cavity, one must forego analytic
solution of Poisson’s equation. A variety of numerical algorithms have been
introduced, either for classical biomolecular applications,132, 133, 134, 135,
136, 137, 138, 139, 140, 141, 142 or else for quantum electronic structure
calculations.143, 144, 145, 39, 64 Numerically, Eq. (2.3) and its Poisson-
Boltzmann analogue are partial differential equations (for $\varphi$) that
require discretization of three-dimensional space, which must extend well into
the continuum region due to the $(\varepsilon_{\text{s}}r)^{-1}$ asymptotic
decay of $\varphi(\mathbf{r})$. This necessitates the use of multiresolution
techniques for efficiency.132, 134, 62, 55, 142, 39 Worth mentioning also is
the Langevin dipoles model,146, 147, 148 which discretizes the continuum
solvent using a three-dimensional grid of point dipoles, and is therefore a
direct realization of the conceptual notion that $\mathbf{P}$ is a dipole
density.

As compared to these approaches, each of which requires discretization of
three-dimensional space, a much more efficient two-dimensional formulation of
the continuum electrostatics problem is often feasible. The transformation of
Eq. (2.3) from a partial differential equation in three dimensions to an
integral equation in two dimensions (over the cavity surface $\bm{\Gamma}$)
defines the class of methods known as polarizable continuum models (PCMs),13,
15, 19 which we next describe.

###  2.3 Polarizable continuum models

In quantum chemistry, PCMs are so widely used as to be nearly synonymous with
continuum solvation itself. This popularity stems from efficiency, which in
turn derives from a transformation of the volumetric polarization theory of
the Poisson or Poisson-Boltzmann equation (Section 2.1) into a surface charge
problem that can be solved with numerical methods that are described in
Section 3.2. That transformation, and the basic working equations of PCMs, are
introduced in this section. Matrix analogues of these integral equations,
which define the various PCMs in practice, are introduced in Section 3.1.

Physically speaking, polarization of the continuum extends beyond the solute
cavity surface, as indicated pictorially in Fig. 1(b). Transformation of the
three-dimensional polarization problem into a two-dimensional problem defined
on the cavity surface $\bm{\Gamma}$ relies on a characteristic feature of a
sharp dielectric interface, namely, a discontinuity in the electric field
resulting from a buildup of charge at the interface. This occurs at any
dielectric interface, including the one that defines the boundary between
atomistic solute and continuum solvent. Let $\mathbf{s}\in\bm{\Gamma}$ denote
a point on the interface (cavity surface), and let $\mathbf{n}_{\mathbf{s}}$
be the outward-pointing unit vector normal to the cavity surface at the point
$\mathbf{s}$. At a sharp interface in the permittivity function
$\varepsilon(\mathbf{r})$, as in Eq. (2.10), the normal electric field
satisfies a “jump” boundary condition,9, 92

$\varepsilon_{\text{out}}\;\bigl{(}\mathbf{n}_{\mathbf{s}}\bm{\cdot}\hat{\bm{\nabla}}\bigr{)}\,\varphi(\mathbf{s})\Big{|}_{\mathbf{s}=\mathbf{s}^{+}}=\varepsilon_{\text{in}}\;\bigl{(}\mathbf{n}_{\mathbf{s}}\bm{\cdot}\hat{\bm{\nabla}}\bigr{)}\,\varphi(\mathbf{s})\Big{|}_{\mathbf{s}=\mathbf{s}^{-}}\;.$
(2.21)

This ensures that the electric displacement
$\mathbf{D}(\mathbf{r})=\varepsilon(\mathbf{r})\,\mathbf{E}(\mathbf{r})$ is
continuous across the interface.41, 38 The notation $\mathbf{s}^{\pm}$
indicates that these are one-sided derivatives to be evaluated immediately
inside ($\mathbf{s}^{-}\in\bm{\Omega}$) or outside
($\mathbf{s}^{+}\notin\bm{\Omega}$) of the solute cavity, as $\varphi$ is only
semi-differentiable at the interface.

Polarization of the continuum manifests as a surface charge that accumulates
at the boundary between the atomistic and the continuum region in order to
satisfy Eq. (2.21), and whose magnitude is proportional to the outward-
pointing normal electric field at the dielectric interface.38 Let us call that
charge $\sigma(\mathbf{s})$, where $\mathbf{s}\in\bm{\Gamma}$, so as to
distinguish it from a volume charge such as $\rho(\mathbf{r})$, where
$\mathbf{r}\in\mathbb{R}^{3}$. Introducing the notation
$\hat{\partial}_{\,\mathbf{s}}=\mathbf{n}_{\mathbf{s}}\bm{\cdot}\hat{\bm{\nabla}}$
to indicate the normal derivative, the surface charge at the dielectric
boundary can be expressed in several equivalent ways,92, 14 two of which are

$\sigma(\mathbf{s})=\frac{1}{4\pi}\left(\frac{\varepsilon_{\text{out}}-\varepsilon_{\text{in}}}{\varepsilon_{\text{in}}}\right)\hat{\partial}_{\,\mathbf{s}}\varphi(\mathbf{s})\Big{|}_{\mathbf{s}=\mathbf{s}^{+}}=\frac{1}{4\pi}\left(\frac{\varepsilon_{\text{out}}-\varepsilon_{\text{in}}}{\varepsilon_{\text{out}}}\right)\hat{\partial}_{\,\mathbf{s}}\varphi(\mathbf{s})\Big{|}_{\mathbf{s}=\mathbf{s}^{-}}\;.$
(2.22)

These differ depending on whether the normal electric field is evaluated
immediately inside or outside of the cavity. In the usual quantum chemistry
case where $\varepsilon_{\text{in}}=1$ and
$\varepsilon_{\text{out}}=\varepsilon_{\text{s}}$, the first form of this
equation (in which the field is evaluated within the continuum region) is
merely the definition of the polarization,

$\mathbf{P}=\left(\frac{\varepsilon_{\text{s}}-1}{4\pi}\right)\mathbf{E}\;,$
(2.23)

as obtained from Eq. (2.1). If the entirety of the solute charge
$\rho(\mathbf{r})$ is confined within the cavity surface, then Eq. (2.22) is
simply a reflection of Gauss’ Law. For cavities that are realistic
approximations to the size and shape of a molecular solute, however, the tails
of a quantum-mechanical charge distribution will penetrate into the continuum
region.91, 92, 93, 94, 95, 96, 97 This outlying or “escaped” charge is
discussed in more detail below but will be ignored for now.

It has been argued that the second form of Eq. (2.22), in which the derivative
is evaluated inside of the cavity, should be used in order to avoid “self-
polarization” of the medium.9 Then taking $\varepsilon_{\text{in}}=1$ and
$\varepsilon_{\text{out}}=\varepsilon_{\text{s}}$ and recognizing that
$\mathbf{E}=-\hat{\bm{\nabla}}\varphi$, one obtains a model

$\sigma(\mathbf{s})=-\frac{1}{4\pi}\left(\frac{\varepsilon_{\text{s}}-1}{\varepsilon_{\text{s}}}\right)\underbrace{\big{[}\mathbf{E}^{\rho}(\mathbf{s})+\mathbf{E}^{\sigma}(\mathbf{s})\big{]}\bm{\cdot}\mathbf{n}_{\mathbf{s}}}_{E_{\perp}(\mathbf{s})}=\frac{1}{4\pi}\left(\frac{\varepsilon_{\text{s}}-1}{\varepsilon_{\text{s}}}\right)\left(\frac{\partial\varphi}{\partial\mathbf{n}_{\mathbf{s}}}\right)_{\mathbf{s}=\mathbf{s}^{-}}$
(2.24)

that corresponds to the original PCM introduced by Tomasi and co-workers.3, 4,
5, 9 Here, the induced surface charge density $\sigma(\mathbf{s})$ is
proportional to the outward-directed normal electric field,
$E_{\perp}(\mathbf{s})=\mathbf{E}(\mathbf{s})\bm{\cdot}\mathbf{n}_{\mathbf{s}}$,
which is separated into two contributions in Eq. (2.24), analogous to how
$\varphi(\mathbf{r})$ is partitioned in Eq. (2.4). The contribution
$\mathbf{E}^{\rho}=-\hat{\bm{\nabla}}\varphi^{\rho}$ comes from the solute,
whose electrostatic potential $\varphi^{\rho}(\mathbf{r})$ is defined in Eq.
(2.4), whereas $\mathbf{E}^{\sigma}=-\hat{\bm{\nabla}}\varphi^{\sigma}$ can be
identified as the reaction-field contribution and originates in the
electrostatic potential

$\varphi^{\sigma}(\mathbf{r})=\int_{\mathbf{s}\in\bm{\Gamma}}\frac{\sigma(\mathbf{s})}{\|\mathbf{s}-\mathbf{r}\|}\;d\mathbf{s}$
(2.25)

that is generated by the surface charge $\sigma$. To compute
$\varphi^{\sigma}(\mathbf{r})$, it is only necessary to discretize the cavity
surface $\bm{\Gamma}$ rather than the whole of three-dimensional space.
Historically, Eq. (2.24) was the first example of what has been called an
apparent surface charge (ASC) formulation of the continuum electrostatics
problem. Nowadays the term “ASC model” is essentially synonymous with PCM, as
multipolar expansions and other simplified treatments see little use in
contemporary quantum chemistry because the model defined by Eq. (2.24) makes
it easy to use the exact solute density $\rho(\mathbf{r})$ along with a cavity
or arbitrary shape. (Discretization of the surface is considered in Section
3.) For a spherical cavity, Eq. (2.24) is equivalent to the use of the
Kirkwood multipolar expansion formulas if the latter are carried to arbitrary
order, as can be proven explicitly.149

In early literature, the model defined by Eq. (2.24) is often called “the”
ASC-PCM,9 whereas in contemporary literature it is usually called D-PCM,13, 15
or sometimes the Miertuš-Scrocco-Tomasi model.150, 98 The somewhat arbitrary
decision to use the right side of Eq. (2.22) can be avoided, and a somewhat
different version of this model obtained, by noting that the discontinuity in
$E_{\perp}$ at the cavity surface can be expressed as92

$\displaystyle\hat{\partial}_{\,\mathbf{s}}\,\varphi^{\sigma}(\mathbf{s})\Big{|}_{\mathbf{s}=\mathbf{s}^{-}}$
$\displaystyle=2\pi\sigma(\mathbf{s})+\hat{\partial}_{\,\mathbf{s}}\,\varphi^{\sigma}(\mathbf{s})$
(2.26a)
$\displaystyle\hat{\partial}_{\,\mathbf{s}}\,\varphi^{\sigma}(\mathbf{s})\Big{|}_{\mathbf{s}=\mathbf{s}^{+}}$
$\displaystyle=-2\pi\sigma(\mathbf{s})+\hat{\partial}_{\,\mathbf{s}}\,\varphi^{\sigma}(\mathbf{s})\;.$
(2.26b)

Adding these two equations and combining them with Eq. (2.22) affords a
different expression for the surface charge,151, 149, 92

$\sigma(\mathbf{s})=\left(\frac{f_{\varepsilon}}{2\pi}\right)\hat{\partial}_{\,\mathbf{s}}\big{[}\varphi^{\rho}(\mathbf{s})+\varphi^{\sigma}(\mathbf{s})\big{]}\;,$
(2.27)

in which the normal derivative is evaluated at (rather than near) the point
$\mathbf{s}\in\bm{\Gamma}$. The permittivity-dependent prefactor in this
expression is

$f_{\varepsilon}=\frac{\varepsilon_{\text{out}}-\varepsilon_{\text{in}}}{\varepsilon_{\text{out}}+\varepsilon_{\text{in}}}\;.$
(2.28)

(For a QM solute the only sensible choice is $\varepsilon_{\text{in}}=1$, but
the more general notation is retained for now, in order to accommodate
different values of $\varepsilon_{\text{in}}$ that are sometimes used in
biomolecular electrostatics calculations.46, 47, 44, 48, 7, 45)

The model of Eq. (2.27) can be recast in a convenient form by defining an
operator $\hat{D}^{\dagger}$ that acts on surface functions
($\mathbf{s}\in\bm{\Gamma}$) according to94, 95, 96

$\hat{D}^{\dagger}\sigma(\mathbf{s})=\int_{\bm{\Gamma}}d\mathbf{s}^{\prime}\sigma(\mathbf{s}^{\prime})\frac{\partial}{\partial\mathbf{n}_{\mathbf{s}}}\left(\frac{1}{\|\mathbf{s}-\mathbf{s}^{\prime}\|}\right)=-E_{\perp}^{\sigma}(\mathbf{s})\;.$
(2.29)

The first equality defines $\hat{D}^{\dagger}$ and the second follows upon
realization that
$\partial/\partial\mathbf{n}_{\mathbf{s}}=\mathbf{n}_{\mathbf{s}}\bm{\cdot}\hat{\bm{\nabla}}$
can be pulled outside of the integral, leaving
$(\mathbf{n}_{\mathbf{s}}\bm{\cdot}\hat{\bm{\nabla}})\varphi^{\sigma}(\mathbf{s})=-E_{\perp}^{\sigma}(\mathbf{s})$.
The operator $\hat{D}^{\dagger}$ is often called $\hat{D}^{\ast}$ in the
literature,94, 95, 96, 152 but is the adjoint of the double-layer operator
$\hat{D}$ that is introduced below, and the notation used here is selected to
reflect that fact. Using $\hat{D}^{\dagger}$ to rewrite Eq. (2.27) affords an
alternative to Eq. (2.24), namely

$\left[\left(\frac{2\pi}{f_{\varepsilon}}\right)\hat{1}-\hat{D}^{\dagger}\right]\sigma(\mathbf{s})=-E_{\perp}^{\rho}(\mathbf{s})\;.$
(2.30)

This equation makes clear that the sole ingredient needed to determine the
induced surface charge is

$E_{\perp}^{\rho}(\mathbf{s})=-(\mathbf{n}_{\mathbf{s}}\bm{\cdot}\hat{\bm{\nabla}})\varphi^{\rho}(\mathbf{s})=-\left(\frac{\partial\varphi^{\rho}}{\partial\mathbf{n}_{\mathbf{s}}}\right)\;,$
(2.31)

which is the normal electric field due to the solute, evaluated at the cavity
surface. Chipman refers to Eq. (2.30) as the surface polarization for
electrostatics (SPE) method,95 but others have called it D-PCM,152 or simply
PCM.153

The D-PCM approach, which requires explicit calculation of the electric field
at the cavity surface, has largely been superseded by alternative ASC-PCMs
that determines $\sigma(\mathbf{s})$ using only electrostatic potentials and
not their derivatives, which might be more sensitive to discretization error.
The modern approach is known as the integral equation formalism (IEF-)
PCM,154, 152 and is based on a reformulation of the continuum electrostatics
problem as a boundary-value problem.155, 156, 157 This reformulation is exact
provided that the escaped charge is zero (e.g., for a classical solute), and
we continue to defer a discussion of the escaped charge problem. IEF-PCM is
formulated in terms of integral operators $\hat{S}$ and $\hat{D}$ that act on
surface functions (defined on $\bm{\Gamma}$) to generate the single- and
double-layer potentials, respectively. These operators are defined by

$\hat{S}\sigma(\mathbf{s})=\int_{\bm{\Gamma}}d\mathbf{s}^{\prime}\frac{\sigma(\mathbf{s}^{\prime})}{\|\mathbf{s}^{\prime}-\mathbf{s}\|}=\varphi^{\sigma}(\mathbf{s})\;,$
(2.32)

which generates the electrostatic potential associated with the surface charge
distribution $\sigma(\mathbf{s})$, and

$\hat{D}\sigma(\mathbf{s})=\int_{\bm{\Gamma}}d\mathbf{s}^{\prime}\sigma(\mathbf{s}^{\prime})\underbrace{\frac{\partial}{\partial\mathbf{n}_{\mathbf{s}^{\prime}}}\left(\frac{1}{\|\mathbf{s}^{\prime}-\mathbf{s}\|}\right)}_{D(\mathbf{s},\mathbf{s}^{\prime})}\;,$
(2.33)

which generates the double-layer potential, in the language of integral
equations.152 [The operator $\hat{D}$ is the adjoint of $\hat{D}^{\dagger}$ in
Eq. (2.29),13 as becomes clear upon reversing the indices of the kernel
$D(\mathbf{s},\mathbf{s}^{\prime})$.] Using $\hat{S}$ and $\hat{D}$, the
continuum electrostatics problem can be recast as an integral equation on the
surface of the cavity:13

$\left[\left(\frac{2\pi}{f_{\varepsilon}}\right)\hat{1}-\hat{D}\right]\hat{S}\sigma(\mathbf{s})=\big{(}-2\pi\hat{1}+\hat{D}\big{)}\,\varphi^{\rho}(\mathbf{s})\;.$
(2.34)

Equation (2.34) is the basic working equation of IEF-PCM. In early papers the
working equation was formulated somewhat differently, requiring
$E_{\perp}^{\rho}$ in addition to $\varphi^{\rho}$.155, 156, 157 That form is
sometimes called simply “IEF”,95, 152 to distinguish it from the IEF-PCM that
requires $\varphi^{\rho}$ but not its derivative, and should therefore be more
stable numerically. Equivalence of the two forms is demonstrated in Ref. 158,
and in fact Eq. (2.34) can be cast in a variety of equivalent forms,94, 95,
152 by taking advantage of the fact that $\hat{S}=\hat{S}^{\dagger}$ and

$\hat{D}\hat{S}=\hat{S}\hat{D}^{\dagger}\;.$ (2.35)

However, except for spherical cavities (for which
$\hat{D}=\hat{D}^{\dagger}$),94 the operator identity in Eq. (2.35) is not
generally preserved upon discretization,153, 20 with the practical result that
various forms of Eq. (2.34) are inequivalent as finite-dimensional matrix
equations.20, 19 These equations are discussed in Section 3.1.

For now, we simply note that Eq. (2.34) is an exact reformulation of the
classical continuum electrostatics problem, meaning the problem that is
defined by Poisson’s equation with a sharp dielectric boundary and where the
solute’s charge density $\rho(\mathbf{r})$ is contained entirely within the
cavity. That caveat is satisfied if the solute consists of atomic point
charges (or higher-order multipoles, including polarizable ones) from a force
field. In that case, solution of Eq. (2.34) for $\sigma(\mathbf{s})$
constitutes an exact solution to the electrostatics problem, and the
electrostatic solvation energy is18

$\mathcal{G}_{\text{elst}}=\frac{1}{2}\int_{\mathbb{R}^{3}}\varphi^{\sigma}(\mathbf{r})\;\rho(\mathbf{r})\;d\mathbf{r}=\frac{1}{2}\int_{\bm{\Gamma}}\varphi^{\rho}(\mathbf{s})\;\sigma(\mathbf{s})\;d\mathbf{s}\;.$
(2.36)

These are analogous to the two forms of $\mathcal{G}_{\text{elst}}$ that are
given in Eq. (2.6), however the second form in Eq. (2.36) requires only
surface integration. Using this ASC formulation, in lieu of discretizing
three-dimensional space, has significant advantages over the traditional
approach to biomolecular electrostatics, which require discretization far
enough into the continuum such that $\varphi^{\rho}(\mathbf{r})$ has decayed
to zero. Most contemporary biomolecular electrostatics calculations are also
based on finite-difference evaluation of the Laplacian
$\hat{\nabla}^{2}\varphi(\mathbf{r})$,159, 160 which leads to problems in
obtaining smooth forces for molecular dynamics.161, 162, 163 Discretization of
the action of $\hat{D}$ and $\hat{S}$ on surface functions can be accomplished
in a manner that affords inherently smooth forces;17, 18, 19 see Section 3.2.
Especially for biomolecular applications it is worth noting that the IEF-PCM
has been adapted to provide a solution to the linearized Poisson-Boltzmann
problem,155, 156, 164, 68 including the “size-modified” version that accounts
for finite size of the mobile ions.68 Large biomolecular solutes have been
tackled in this way,165, 166, 167, 19 but this requires iterative solvers for
the matrix equations that define the PCM. Linear-scaling implementations that
can handle biomolecular solutes are discussed in Section 3.4.

For QM solutes, however, there is always escaped charge for realistic cavity
sizes, and therefore IEF-PCM is not a fully equivalent substitute for
Poisson’s equation. The extent to which this is a problem is unclear from the
original derivation of IEF-PCM provided by Cancès et al.,155, 156, 157, 152
which does not provide much physical insight, nor does it emphasize the
assumption (inherent in the derivation) that there is no outlying charge. That
issue was addressed directly by Chipman,92, 93, 94, 95 who assumes from the
start that there is outlying charge and that as a consequence the correct
reaction-field potential is not $\varphi^{\sigma}(\mathbf{r})$ but rather

$\varphi_{\text{rxn}}(\mathbf{r})=\varphi^{\sigma}(\mathbf{r})+\varphi^{\beta}(\mathbf{r})\;,$
(2.37)

where $\varphi^{\sigma}(\mathbf{r})$ arises from the accumulation of charge
$\sigma(\mathbf{s})$ at the dielectric interface but is accompanied by an
additional potential $\varphi^{\beta}(\mathbf{r})$ due to volume polarization.
The latter originates with the tail of the solute’s charge density that
extends beyond the cavity surface. Including an additional term for
$\varphi^{\beta}$ on the right side of Eq. (2.27), and recognizing that
$\hat{\partial}_{\,\mathbf{s}}\varphi^{\sigma}(\mathbf{s})=\hat{D}^{\dagger}\sigma(\mathbf{s})$,
an exact equation for the surface charge that includes the effects of volume
polarization is94

$\displaystyle\begin{aligned}
\left[\hat{1}-\left(\frac{f_{\varepsilon}}{2\pi}\right)\hat{D}^{\dagger}\right]\sigma(\mathbf{s})&=\frac{f_{\varepsilon}}{2\pi}\big{[}\hat{\partial}_{\,\mathbf{s}}\,\varphi^{\rho}(\mathbf{s})+\hat{\partial}_{\,\mathbf{s}}\,\varphi^{\beta}(\mathbf{s})\big{]}\\\
&=-\frac{f_{\varepsilon}}{2\pi}\big{[}E_{\perp}^{\rho}(\mathbf{s})+E_{\perp}^{\beta}(\mathbf{s})\big{]}\;.\end{aligned}$
(2.38)

The potential $\varphi^{\beta}(\mathbf{r})$ can be modeled as the
electrostatic potential generated by a charge density

$\beta(\mathbf{r})=\begin{cases}0&\text{for $\mathbf{r}\in\bm{\Omega}$}\\\
\left(\varepsilon_{\text{out}}^{-1}-\varepsilon_{\text{in}}^{-1}\right)\rho(\mathbf{r})&\text{for
$\mathbf{r}\notin\bm{\Omega}$}\end{cases}$ (2.39)

that satisfies a vacuum-like Poisson equation,

$\hat{\nabla}^{2}\varphi^{\beta}(\mathbf{r})=-4\pi\beta(\mathbf{r})\;.$ (2.40)

Equation (2.38) can be solved numerically,91 to afford an exact solution to
the continuum electrostatics problem, in a method that Chipman calls surface
and volume polarization for electrostatics (SVPE).95 This approach is
challenging in practice because the volume charge density $\beta(\mathbf{r})$
is discontinuous at the cavity surface, however Eq. (2.38) can be recast into
a form that requires only the normal electric field
$E_{\perp}^{\rho}(\mathbf{s})$ at the cavity surface, along with the solution
of the ASC-PCM that is introduced next.96 This provides a practical means to
access exact electrostatics, even in the presence of outlying charge, while
staying within a two-dimensional surface integral formalism.

A simplified (if approximate) treatment is possible, which eliminates the
normal electric field in Eq. (2.38) and ultimately connects back to IEF-PCM.
This model is obtained by introducing an additional surface charge
$\alpha(\mathbf{s})$, distinct from $\sigma(\mathbf{s})$, that is defined by
the condition

$\hat{S}\alpha(\mathbf{s})=\varphi^{\beta}(\mathbf{s})\;,$ (2.41)

meaning that its electrostatic potential $\varphi^{\alpha}=\hat{S}\alpha$ must
reproduce $\varphi^{\beta}$ on the cavity surface.91, 92, 93, 94 This also
ensures that $\varphi^{\alpha}(\mathbf{r})=\varphi^{\beta}(\mathbf{r})$ for
all interior points $\mathbf{r}\in\bm{\Omega}$, and while the two potentials
may differ outside of the cavity, those contributions are scaled by
$\varepsilon_{\text{s}}^{-1}$ and therefore less important. (This is confirmed
in numerical tests.91, 97) Assuming that the true surface charge, augmented to
reflect volume polarization, is
$\tilde{\sigma}(\mathbf{s})=\alpha(\mathbf{s})+\sigma(\mathbf{s})$, the term
$\hat{\partial}_{\,\mathbf{s}}\,\varphi^{\beta}(\mathbf{s})$ in Eq. (2.38) can
be manipulated into a form that is consistent with the ASC formalism.94 The
result is a model that Chipman has called surface and simulation of volume
polarization for electrostatics [SS(V)PE],94, 95, 96

$\underbrace{\hat{S}\left(\hat{1}-\frac{f_{\varepsilon}}{2\pi}\hat{D}^{\dagger}\right)}_{\hat{K}_{\varepsilon}}\tilde{\sigma}(\mathbf{s})=\underbrace{f_{\varepsilon}\left(\frac{1}{2\pi}\hat{D}-\hat{1}\right)}_{\hat{Y}_{\varepsilon}}\varphi^{\rho}(\mathbf{s})\;.$
(2.42)

Using the identity in Eq. (2.35), this equation is easily rearranged to afford
Eq. (2.34) that defines IEF-PCM, and therefore the two models are equivalent
at the level of integral operators.94, 95, 158 (They differ in practice, as
described in Section 3.1.) Importantly, what the derivation of SS(V)PE makes
clear is that the surface charge that is determined by solving Eq. (2.34)
implicitly contains the (approximate) effects of volume polarization, which
was not evident from the derivation by Cancès and co-workers.155, 156, 157,
152 Chipman’s derivation clarifies that both approaches constitute an exact
treatment of electrostatic interactions in the limiting case that there is no
escaped charge ($\alpha\equiv 0$).

###  2.4 Comparison of ASC-PCMs

Table 1 presents electrostatic solvation energies
($\mathcal{G}_{\text{elst}}$) for several small molecules and ions, in both a
low-dielectric solvent (toluene, $\varepsilon_{\text{s}}=2.4$) and a high-
dielectric solvent (water, $\varepsilon_{\text{s}}=78$).95 The SVPE method
[Eq. (2.38)] affords the exact result but SS(V)PE solvation energies are
within 0.1 kcal/mol. In comparison, the SPE method of Eq. (2.30), which is
equivalent to D-PCM, exhibits noticeable differences, especially for ions. The
amount of outlying charge ($Q_{\text{out}}$) is also quantified in Table 1.
For future reference we note the obvious definitions

$\displaystyle Q_{\text{in}}$
$\displaystyle=\int_{\mathbf{r}\in\bm{\Omega}}\rho(\mathbf{r})\;d\mathbf{r}$
(2.43a) $\displaystyle Q_{\text{out}}$
$\displaystyle=\int_{\mathbb{R}^{3}}\rho(\mathbf{r})\;d\mathbf{r}-Q_{\text{in}}\;.$
(2.43b)

The charge density $\rho(\mathbf{r})$ includes both nuclei and electrons, so
$Q_{\text{in}}+Q_{\text{out}}=0$ for a neutral molecule. The escaped charge is
generally $|Q_{\text{in}}|\approx 0.1$–$0.2e$ for small solutes.91, 97

Table 1:  Electrostatic solvation energies in toluene
($\varepsilon_{\text{s}}=2.4$) and in water ($\varepsilon_{\text{s}}=78.3$),
computed with various approaches.a The SVPE method [Eq. (2.38)] affords the
exact result and SPE is the method in Eq. (2.30).

Solute | $\varepsilon_{\text{s}}$ |  $\mathcal{G}_{\text{elst}}$ (kcal/mol) | $Q_{\text{out}}$
---|---|---|---
|  | SVPE | SS(V)PE | SPE | COSMOb | (a.u.)c
|  |  |  |  | $\zeta=0$ | $\zeta=1/2$ |
H2O | $2.4$ | $-3.9$ | $-3.9$ | $-4.0$ | $-4.8$ | $-3.9$ |
CH3CONH2 | $2.4$ | $-5.3$ | $-5.0$ | $-5.2$ | $-5.9$ | $-4.8$ |
NO+ | $2.4$ | $-52.2$ | $-52.2$ | $-55.3$ | $-52.5$ | $-43.4$ |
CN- | $2.4$ | $-39.4$ | $-39.4$ | $-35.0$ | $-39.4$ | $-32.5$ |
H2O | $78.3$ | $-8.6$ | $-8.6$ | $-8.7$ | $-8.6$ | $-8.6$ | $-0.06$
CH3CONH2 | $78.3$ | $-10.9$ | $-10.8$ | $-11.1$ | $-10.9$ | $-10.8$ | $-0.15$
NO+ | $78.3$ | $-89.5$ | $-89.5$ | $-94.7$ | $-89.5$ | $-88.9$ | $-0.07$
CN- | $78.3$ | $-67.4$ | $-67.3$ | $-56.8$ | $-67.3$ | $-66.9$ | $-0.17$

  * •

aFrom Ref. 95; all calculations used an isodensity cavity with
$\rho_{0}=0.001$ a.u.. bUsing a renormalization factor
$\tilde{f}_{\varepsilon}(\zeta)$, Eq. (2.45). cFrom Ref. 97.

By arbitrarily dropping the $\hat{D}$\- and $\hat{D}^{\dagger}$-dependent
terms in Eq. (2.42), one obtains a model
$\hat{S}\sigma=-f_{\varepsilon}\varphi^{\rho}$ that we rewrite as

$\hat{S}\sigma(\mathbf{s})=-\tilde{f}_{\varepsilon}(\zeta)\;\varphi^{\rho}(\mathbf{s})\;.$
(2.44)

This introduces a parameter $\zeta$ into the permittivity factor that was
defined in Eq. (2.28), rewriting it as

$\tilde{f}_{\varepsilon}(\zeta)=\frac{\varepsilon_{\text{s}}-1}{\varepsilon_{\text{s}}+\zeta}$
(2.45)

for the usual QM case where $\varepsilon_{\text{in}}=1$ and
$\varepsilon_{\text{out}}=\varepsilon_{\text{s}}$. The model in Eq. (2.44) has
a long history and a variety of names, one of which is the conductor-like
screening model (COSMO).168 Since the neglected double-layer operator embodies
the electric field discontinuity at the cavity surface, this condition is not
satisfied when these terms are neglected to obtain Eq. (2.44),152, 68 and
rescaling the surface charge by $\tilde{f}_{\varepsilon}(\zeta)$ can be seen
as an attempt to mimic the effect of the jump boundary condition.

The name COSMO hints at the original derivation of the model in Eq. (2.44):
for a conductor ($\varepsilon_{\text{s}}=\infty$), the total electrostatic
potential vanishes at the cavity surface and the ASC formulation of the
Poisson problem is simply $\hat{S}\sigma=-\varphi^{\rho}$.93 A scaling factor
$\tilde{f}_{\varepsilon}(\zeta)$ is then introduced to account for the fact
that $\varepsilon_{\text{s}}$ is finite, and that there is charge penetration
into the medium. The scaling factor [Eq. (2.45)] can be justified based on the
normalization condition for the total surface charge, which is92, 156, 169

$\int_{\bm{\Gamma}}\sigma(\mathbf{s})\;d\mathbf{s}=-\left(\frac{1}{\varepsilon_{\text{in}}}-\frac{1}{\varepsilon_{\text{out}}}\right)Q_{\text{in}}\;.$
(2.46)

For $\varepsilon_{\text{in}}=1$ and
$\varepsilon_{\text{out}}=\varepsilon_{\text{s}}$, the total surface charge
$Q_{\text{surf}}=-[(\varepsilon_{\text{s}}-1)/\varepsilon_{\text{s}}]Q_{\text{in}}$,
and it was originally suggested to take $\zeta=1/2$ in order to renormalize
for outlying charge,168 although the value $\zeta=0$ was later recommended for
ions,170 and values $0\leq\zeta\leq 2$ have also been used.10 The choice
$\zeta=0$ has variously been called the “generalized COSMO” (GCOSMO)
model,171, 172, 173, 174 or the conductor-like PCM (C-PCM).175, 176 In fact,
the normalization condition in Eq. (2.46) forms the basis of various ad hoc
attempts to rescale the surface charge,3, 5, 177, 178, 179 but tests against
exact results from Poisson’s equation suggest that none of these schemes is
very satisfactory.91 Charge rescaling also complicates analytic energy
gradients.180

At a practical level, the choice $\zeta=1/2$ for neutral solutes and $\zeta=0$
for ions works remarkably well in comparison to the IEF-PCM and SS(V)PE
methods, even in low-dielectric solvents. This is suggested by the smattering
of data in Table 1 and confirmed by calculations on a much larger data set,
for which the statistical difference between these methods is $\lesssim 0.1$
kcal/mol for neutral solutes and $\approx 0.5$ kcal/mol for ions, even at
$\varepsilon_{\text{s}}=2$.170 Some justification for these parameter choices
is found in the fact that $\tilde{f}_{\varepsilon}(0)$ looks like the
$\varepsilon$-dependent factor in the Born ion model [Eq. (2.17)], whereas
$\tilde{f}_{\varepsilon}(1/2)$ affords the prefactor in the Bell-Onsager model
of dipole solvation [Eq. (2.18)]. With appropriate choice of $\zeta$, COSMO is
therefore little different from SS(V)PE or IEF-PCM in practice, and
considerably simpler. It can be extended in a straightforward way to solvents
with nonzero ionic strength, with or without size exclusion.68 For large
biomolecular applications in aqueous solvent, there seems little reason not to
use this approach, in lieu of exact but more complicated models, and as an
alternative to finite-difference solution of Poisson’s equation.

##  3 Implementation

Poisson’s equation with a sharp dielectric interface makes for a cute problem
in applied mathematics and there has sometimes been a tendency in that
community to analyze its properties, or those of the various cavity surfaces
discussed in Section 2.2, at fixed molecular geometries only. This is not
useful for chemistry, which requires a potential energy surface. Exploration
of that surface requires gradients of the total energy functional
$\mathcal{G}[\Psi]$, which in turn requires discretization of the integral
equations introduced above. In addition, many spectroscopic observables can be
formulated as analytic energy derivatives,181, 182, 183 so a solution-phase
theory that encompasses molecular properties requires a differentiable model.
Numerical implementation of the models introduced in Section 2 is the topic of
this section. We deal specifically with PCMs in Sections 3.1 and 3.2, but
admit more general Poisson equation models in Section 3.3. Linear-scaling
implementations, which are useful for hybrid QM/MM simulations with PCM
boundary conditions, are discussed in Section 3.4.

###  3.1 Matrix equations for PCMs

In practice, the integral equation that defines any PCM must be discretized in
order to obtain a finite-dimensional matrix equation. As such, it is
convenient to rewrite the equation that defines SS(V)PE and IEF-PCM as

$\hat{K}_{\varepsilon}\sigma(\mathbf{s})=\hat{Y}_{\varepsilon}\,\varphi^{\rho}(\mathbf{s})\;,$
(3.1)

which is simply a restatement of Eq. (2.42). Using various definitions for
$\hat{K}_{\varepsilon}$ and $\hat{Y}_{\varepsilon}$, this equation can
encompass a whole family of ASC-PCMs including the simplified C-PCM and COSMO
methods (Table 2).95, 129, 20, 19

Table 2:  Matrices that define various PCMs according to
$\mathbf{K}_{\varepsilon}\mathbf{q}=\mathbf{Y}_{\\!\varepsilon}\mathbf{v}^{\rho}$
[Eq. (3.2)].

Method | Matrix $\mathbf{K}_{\varepsilon}$ | Matrix $\mathbf{Y}_{\\!\varepsilon}$
---|---|---
C-PCMa,b | $\mathbf{S}$ | $-\tilde{f}_{\varepsilon}(0)\,\bm{1}$
COSMOb | $\mathbf{S}$ | $-\tilde{f}_{\varepsilon}(\zeta)\,\bm{1}$
DESMOc | $\mathbf{S}$ | $-\bm{1}+(1/\varepsilon)\mathbf{M}$
SS(V)PEd | $\mathbf{S}-(f_{\varepsilon}/4\pi)(\mathbf{DAS}+\mathbf{SAD}^{\dagger})$ | $-f_{\varepsilon}\big{[}\bm{1}-(1/2\pi)\mathbf{DA}\big{]}$
IEF-PCMd | $\mathbf{S}-(f_{\varepsilon}/2\pi)\mathbf{DAS}$ | $-f_{\varepsilon}\big{[}\bm{1}-(1/2\pi)\mathbf{DA}\big{]}$

  * •

aAlso known as GCOSMO.
b$\tilde{f}_{\varepsilon}(\zeta)=(\varepsilon-1)/(\varepsilon+\zeta)$.
c$M_{ij}=\delta_{ij}\,\varphi^{\rho}_{\kappa}(\mathbf{s}_{i})/\varphi^{\rho}_{0}(\mathbf{s}_{i})$.
d$f_{\varepsilon}=(\varepsilon-1)/(\varepsilon+1)$.

To discretize Eq. (3.1), one must first generate a surface grid of points
$\mathbf{s}_{i}\in\bm{\Gamma}$. The surface charge distribution
$\sigma(\mathbf{s})$ is thereby replaced by a set of point charges
$\\{q_{i}\\}$ at the discretization points $\\{\mathbf{s}_{i}\\}$. The details
of this procedure are discussed in the next section but for now it suffices to
introduce a matrix notation

$\mathbf{K}_{\varepsilon}\mathbf{q}=\mathbf{Y}_{\\!\varepsilon}\mathbf{v}^{\rho}$
(3.2)

for the discretized form ofEq. (3.1), where $\mathbf{q}$ is a vector of
surface charges and $\mathbf{v}^{\rho}$ denotes the molecular electrostatic
potential evaluated at the cavity surface,
$\upsilon_{i}^{\rho}=\varphi^{\rho}(\mathbf{s}_{i})$. In this form, surface
integrals are replaced by scalar products, e.g.,
$\mathcal{G}_{\text{elst}}=\tfrac{1}{2}\mathbf{q}\bm{\cdot}\mathbf{v}^{\rho}$
replaces Eq. (2.36). An alternative form of Eq. (3.2) that is sometimes
encountered is $\mathbf{q}=\mathbf{Q}_{\varepsilon}\mathbf{v}^{\rho}$ where
$\mathbf{Q}_{\varepsilon}=\mathbf{K}_{\varepsilon}^{-1}\mathbf{Y}_{\\!\varepsilon}$,
but whereas the corresponding operator
$\hat{Q}_{\varepsilon}=\hat{K}_{\varepsilon}^{-1}\hat{Y}_{\\!\varepsilon}$ is
self-adjoint, this property is generally not preserved by discretization
except in the special case of C-PCM and COSMO.18 This means that the mapping
from Eq. (3.1) to Eq. (3.2) is not unique, because discretization fails to
preserve the condition $\hat{D}\hat{S}=\hat{S}\hat{D}^{\dagger}$.20 In matrix
form, this implies that $\mathbf{DAS}\neq\mathbf{SAD}^{\dagger}$ (except for
spherical cavities),153, 20 where $\mathbf{A}$ is a diagonal matrix consisting
of the surface area $a_{i}$ associated with each discretization point
$\mathbf{s}_{i}$. This leads to an ambiguity in the matrix representation of
the operator
$\hat{K}_{\varepsilon}=\hat{S}-(f_{\varepsilon}/2\pi)\hat{S}\hat{D}^{\dagger}$
[Eq. (2.42)], since it can be argued that any matrix of the form

$\mathbf{K}_{\varepsilon}=\mathbf{S}-\left(\frac{f_{\varepsilon}}{4\pi}\right)\underbrace{\big{(}c_{1}\mathbf{DAS}+c_{2}\mathbf{SAD}^{\dagger}\big{)}}_{\mathbf{X}}\;.$
(3.3)

is an equally valid representation, provided that $c_{1}+c_{2}=1$.20 The
various choices for $\mathbf{X}$ are inequivalent in matrix form.
Historically, IEF-PCM has been implemented using $\mathbf{X}=\mathbf{DAS}$
(that is, $c_{1}=1$ and $c_{2}=0$), whereas SS(V)PE was implemented using the
symmetrized matrix ($\mathbf{DAS}+\mathbf{SAD}^{\dagger})/2$.95, 19 This is
indicated in Table 2, which provides matrix definitions for these and other
commonly-encountered PCMs. Precise definitions of the matrices $\mathbf{S}$
and $\mathbf{D}$ that represent the operators $\hat{S}$ and $\hat{D}$ can be
found elsewhere.129, 124, 18, 20, 19 These depend somewhat upon the
discretization algorithm that is selected, but generally $S_{ij}$ represents
the Coulomb interaction between $q_{i}$ and $q_{j}$ (which is straightforward
to discretize except when $i=j$), whereas $D_{ij}$ incorporates the effects of
the outward-pointing electric field.

Fig. 3:  (a) Comparison of aqueous solvation energies,
$|\mathcal{G}_{\text{elst}}(\mbox{PCM})-\mathcal{G}_{\text{elst}}(\mbox{APBS})|$,
computed using either an ASC-PCM or else by numerical solution of Poisson’s
equation using the apbs software. The data set of solutes consists of amino
acids described using atomic partial charges from a force field, so that there
is no outlying charge. The traditional implementation of IEF-PCM corresponds
to $\mathbf{X}=\mathbf{DAS}$ (see Table 2) but results are also shown for the
transpose matrix, $\mathbf{SAD}^{\dagger}$, and SS(V)PE is essentially the
average of these two choices. (b) Convergence of $\mathcal{G}_{\text{elst}}$
for classical histidine as a function of the solvent’s dielectric constant,
using the SwiG discretization scheme described in Section 3.2. Dielectric
constants for a few common solvents are indicated along the top axis. Adapted
from Ref. 20; copyright 2011 Elsevier.

In the absence of outlying charge, the IEF-PCM and SS(V)PE models are exact up
to discretization errors, which can be driven to zero in a controlled way.
Although this fact follows from the derivation in Section 2.3, it is worth
emphasizing via numerical calculations. For classical solutes with no outlying
charge, which are described using atomic partial charges from a force field,
Fig. 3(a) presents a comparison of PCM solvation energies to those obtained by
numerical solution of Poisson’s equation,20, 19 using a standard
multiresolution algorithm as implemented in the apbs software.184 To examine
the ambiguity regarding the choice of $\mathbf{X}$, these calculations test
both $\mathbf{X}=\mathbf{DAS}$ (i.e., IEF-PCM) and also
$\mathbf{X}=\mathbf{SAD}^{\dagger}$; the SS(V)PE model should be essentially
the average of these two choices. Calculations in Fig. 3 use dense but finite
discretization grids, and a small systematic discrepancy is evident between
the two choices of $\mathbf{X}$, indicating a systematic discrepancy between
SS(V)PE and IEF-PCM. Numerical values of $\mathcal{G}_{\text{elst}}$ obtained
using the IEF-PCM choice of $\mathbf{X}$ agree to within $\lesssim 0.1$
kcal/mol of apbs results, demonstrating the operational equivalence of the ASC
and the volumetric implementation of continuum electrostatics. This renders
the Kirkwood (multipolar expansion) formulas obsolete, since the requisite
formulas are only valid in the absence of escaped charge anyway, and in that
case the ASC-PCM formalism furnishes an exact numerical solution to the
continuum electrostatics problem, for cavities of arbitrary shape, while also
providing a correction for outlying charge.94, 169, 96

The conductor-like model (C-PCM) represents the high-dielectric limit of IEF-
PCM and the two models do indeed become equivalent as
$\varepsilon\rightarrow\infty$,93, 68, 20 as demonstrated numerically in Fig.
3(b). In practice, there is little difference between the two models already
for moderately polar solvents. For $\varepsilon>10$, there seems little
justification for the increased complexity of IEF-PCM relative to C-PCM. That
said, for non-spherical cavities only the $\mathbf{X}=\mathbf{DAS}$ form of
the $\mathbf{K}_{\varepsilon}$ matrix achieves the correct conductor limit for
finite discretization grids, as demonstrated by both formal and numerical
arguments.20 In general, $\hat{D}^{\dagger}$ proves more challenging to
implement as compared to $\hat{D}$,20, 185 so the absence of
$\hat{D}^{\dagger}$ is one reason to prefer the IEF-PCM form with
$\mathbf{X}=\mathbf{DAS}$.

It is illustrative to note that C-PCM can be derived in an alternative way,68
which also generalizes the method to the Poisson-Boltzmann case in which the
solvent contains a dissolved electrolyte. This can be done by introducing the
following ansatz for the electrostatic potential:68, 19

$\varphi(\mathbf{r})=\begin{cases}\varphi_{0}^{\rho}(\mathbf{r})+\varphi_{0}^{\sigma}(\mathbf{r})&\text{for
$\mathbf{r}\in\bm{\Omega}$}\\\
\varphi_{\kappa}^{\rho}(\mathbf{r})/\varepsilon&\text{for
$\mathbf{r}\notin\bm{\Omega}$}\\\ \end{cases}\;.$ (3.4)

Within the solute cavity, this looks like the solute’s electrostatic potential
($\varphi^{\rho}_{0}\equiv\varphi^{\rho}$) plus that arising from the apparent
surface charge, whereas outside of the cavity the potential is modified to
include screening by dissolved ions:

$\varphi_{\kappa}^{\rho}(\mathbf{r})=\int\rho(\mathbf{r}^{\prime})\;\frac{\exp(-\kappa\|\mathbf{r}-\mathbf{r}^{\prime}\|)}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}\;d\mathbf{r}^{\prime}\;.$
(3.5)

The screened Coulomb potential $e^{-\kappa r}/(\varepsilon_{\text{s}}r)$ is
the form encountered in Debye-Hückel theory,67, 60, 7, 68 with screening
length $\kappa^{-1}$ that is defined in Eq. (2.16). C-PCM is immediately
recovered by requiring $\varphi(\mathbf{r})$ in Eq. (3.4) to remain continuous
across the solute cavity surface.68 This simple ansatz cannot, however, be
made to satisfy the jump boundary condition in Eq. (2.21) and incurs an error
of ${\cal O}(1/\varepsilon)$.68, 19 (This is consistent with the fact that
C-PCM can be “derived” from IEF-PCM simply by dropping the $\hat{D}$\- and
$\hat{D}^{\dagger}$-dependent terms from IEF-PCM and rescaling the surface
charge to compensate,152 as described in Section 2.4.) The simple ansatz in
Eq. (3.4) is easily modifiable to incorporate the effects of an ion exclusion
layer (Stern layer) around the solute cavity, as in the size-modified Poisson-
Boltzmann equation of Eq. (2.12), where the ion accessibility function
$\lambda(\mathbf{r})$ serves the same purpose. In homage to GCOSMO, and in
recognition of the fact that this approach generalizes Debye-Hückel theory to
cavities of arbitrary shape, this approach has been named the Debye-Hückel-
like screening model (DESMO).68, 19

###  3.2 Discretization

Having introduced various PCMs in matrix form, we now turn to the details of
discretizing the cavity surface. Historically, this has been accomplished
using various “tessellation” schemes,124 in which small, flat surface elements
approximate the curved surface of the cavity. The “gepol” algorithm186, 187,
188, 189 is a popular version of this finite-element approach, and an example
of a molecular cavity discretized in this way is presented in Fig. 4(a). This
approach has several limitations, however, including the fact that the number
of tesserae per atomic sphere cannot be increased arbitrarily and thus the
discretization error cannot be systematically driven to zero.20 Furthermore,
the solid geometry of the tessellation procedure is complicated, leading to
very complex formulas for surface areas190 and analytic energy gradients.191
(In fact, second derivatives of the tesserae areas $a_{i}$ were considered
sufficiently complicated that they were not originally formulated, and the PCM
Hessian was implemented in a semi-analytic way, via finite-difference
evaluation of $\partial^{2}a_{i}/\partial x\partial y$.153) More recently,
these complexities have been overcome by discretizing the surface using atom-
centered Lebedev grids,192, 193, 194, 17, 18, 20, 19, 21 which are already
widely used in density functional theory,195, 196, 197, 198 and therefore
widely available in quantum chemistry codes. An example is depicted in Fig.
4(b), where the outline of the vdW surface is evident even though only the
discretization grid points (and the underlying nuclear framework) are shown.
Relative to gepol and other tessellation schemes, Lebedev grids have the
advantage of being systematically improvable so that results can be converged
to the infinite-grid limit.18, 20 Fully analytic Hessians have been formulated
and implemented.199

Fig. 4:  (a) Example of a molecular surface (top) and its tessellation using
the gepol algorithm (bottom). (b) Surface discretization grid for a double-
stranded segment of DNA. Panel (a) is reprinted from Ref. 15; copyright 2012
John Wiley & Sons.

An important issue faced by both tessellation and quadrature is ensuring that
the discretization algorithm produces a smooth potential energy surface as the
atoms are displaced. The appearance of discontinuities in molecular surface
area algorithms used in biomolecular electrostatics calculations was pointed
out long ago,200 and discontinuities are likely the cause of anecodotal
complaints about slow convergence of geometry optimizations using PCMs. These
discontinuities arise because grid points may disappear into (or emerge from
within) the interior of the solute cavity, as displacement of the nuclei
modifies the extent of interpenetration amongst the atomic spheres that define
the cavity. An example is shown in Fig. 5(a), which plots convergence of the
energy during geometry optimization of (adenine)(H${}_{\text{2}}$O)52 in
implicit solvent. Two discretization algorithms, the variable tesserae number
(VTN) method201 and the fixed points with variable areas (FixPVA) approach,202
are shown to exhibit repeated spikes in the energy.17 The VTN algorithm uses a
fixed surface grid that unceremoniously discards surface elements that are
swallowed by the cavity, and it is unsurprising that the corresponding
potential surface exhibits discontinuities, although their magnitude ($>20$
kcal/mol in one case) is disturbing. The FixPVA algorithm, however,
specifically introduces a switching function to attenuate the surface area of
each atomic tesserae within the cavity’s interior. In fact, the sharp changes
in energy along the FixPVA optimization in Fig. 5(a) are not discontinuities
per se but rather near-singularities induced by the switching function, which
allows surface discretization charges to approach one another much more
closely as compared to the VTN scheme.17, 20

Fig. 5:  Demonstration of the SwiG-PCM discretization approach. (a) Geometry
optimization trajectories of (adenine)(H${}_{\text{2}}$O)52 in a C-PCM
representation of bulk water, using several different algorithms to discretize
the vdW cavity surface. (The surface itself is not shown, but the atomistic
region appears in the blue inset.) Optimizations are performed in Cartesian
coordinates so the total number of steps is large. (b) Schematic of the SwiG
discretization algorithm, in which the surface charges $\\{q_{i}\\}$ are
subject to Gaussian blurring and also to a switching function that attenuates
the quadratures weights near the cavity surface. (c) Nonelectrostatic
solvation energy ($\mathcal{G}_{\text{nonelst}}$) as spheres $A$ and $B$ are
pulled part. The value of $\mathcal{G}_{\text{nonelst}}$ is related to the
solvent-exposed surface area and inherits any discontinuities in the surface
area function. Panels (a) and (c) are adapted from Ref. 17; copyright 2010
American Chemical Society.

To avoid this, one can use a combination of a switching function with Gaussian
blurring of the surface charges, as illustrated in Fig. 5(b), in what has been
called the switching/Gaussian (SwiG) discretization procedure. This ensures
that Coulomb interactions between surface discretization elements remain
finite even as the distance between them approaches zero. Such a procedure was
originally developed by York and Karplus to obtain a “smooth COSMO”,192 then
later extended by Lange and Herbert to the complete family of PCMs.17, 18, 20,
21 This scheme uses Lebedev points to discretize the surface rather than
tesserae of finite area. Nevertheless, the solvent-accessible surface area
(SASA) for atom $B$ (whose radius is $R_{B}$) is easily defined:

$\text{SASA}(B)=R_{B}^{2}\sum_{i\in B}^{\text{grid}}w_{i}F_{i}\;.$ (3.6)

Here, $w_{i}$ is the quadrature weight for Lebedev discretization point
$\mathbf{s}_{i}$, and $F_{i}$ is the switching function associated with that
point ($0\leq F_{i}\leq 1$), so $a_{i}=w_{i}F_{i}$ is the surface area that is
assigned to the point $\mathbf{s}_{i}$.17 Equation (3.6) is considerably
simpler than geometric algorithms for determining the exposed surface
areas.190 As discussed in Section 4.1, models of the nonelectrostatic
contributions to the solvation energy usually include terms proportional to
the solvent-exposed surface area so continuity of the potential energy surface
also demands that the surface area be a continuous function of the nuclear
coordinates. For SwiG discretization, it is evident from Fig. 5(c) that this
is indeed the case.

Analytic gradients of SwiG-PCMs are greatly simplified relative to those of
the corresponding gepol-discretized models.18 SwiG-PCM potential energy
surfaces are provably continuous and differentiable,192, 18 and are free of
the unwanted oscillations that plague the FixPVA approach (see Fig. 5).18, 20
SwiG discretization is well-behaved enough to be used for ab initio molecular
dynamics simulations involving bond-breaking, as shown in Fig. 6 for
intramolecular proton transfer in glycine. SwiG-PCM provides forces that are
stable enough to afford good energy conservation despite the significant
deformation of the solute cavity as it transforms between two tautomeric forms
of glycine. This is demonstrated by the energy profile in Fig. 6(b), which
provides a closeup view of $\approx 250$ fs of dynamics during which a bond-
breaking event occurs. On the other hand, it appears that SwiG discretization
exacerbates differences between the $\mathbf{X}=\mathbf{DAS}$ and
$\mathbf{X}=\mathbf{SAD}^{\dagger}$ forms of $\mathbf{K}_{\varepsilon}$ [Eq.
(3.3)], at least in comparison to gepol results where these differences are
small.153 This is not a major problem insofar as the $\mathbf{X}=\mathbf{DAS}$
form provides correct formal properties and good numerical agreement with
Poisson’s equation, as demonstrated in Fig. 3.

Fig. 6:  (a) Electrostatic solvation energy along an ab initio molecular
dynamics trajectory of glycine (PBE0/6-31+G* level) in implicit water
(SwiG/C-PCM). The simulation starts at $t=0$ from the amino acid tautomer
(energy data in orange), which is the most stable form of gas-phase glycine,
but this species spontaneously transfers a proton to form the zwitterionic
tautomer (energy data in blue), during the time window highlighted in the gray
box and shown in detail in (b). The time step is 0.97 fs. (b) Close-up view of
$\mathcal{G}_{\text{elst}}$ in the region where the proton transfer occurs.
Energy fluctuations are smooth despite the bond-breaking even that occurs
around 2.3 ps. Data are from Ref. 18.

The SwiG implementation of PCMs first appeared in the Q-Chem program203 and
related discretization schemes have since been adopted in other software.204,
205 Other discretization approaches have been described206, 207, 208, 209 but
it is unclear that they have been formulated with gradients in mind. Given the
simplicity and success of the SwiG approach, it is also unclear what is to be
gained from these formulations.

###  3.3 Isodensity and self-consistent cavity surfaces

Smooth SwiG-PCM discretization schemes been developed for both vdW cavities18
(including the SAS) and for the “molecular surface” or SES,21 thus providing a
solution to the practical problem of discontinuities in potential energy
surfaces. Unfortunately, these methods do not alter an arguably more
fundamental problem, namely, that construction of the solute cavity remains a
nontrivial source of arbitrariness in both PCMs and in biomolecular Poisson-
Boltzmann electrostatics calculations. Solvation energies119 and other
properties210 can be quite sensitive to the particular atomic radii that are
used, and often the atomic radii that work well for small-molecule solvation
energies do not work well for proteins.211

Both the vdW and SAS cavity constructions consist of atom-centered spheres
with radii

$R_{A}=\alpha R_{\text{vdW},A}+R_{\text{probe}}\;.$ (3.7)

The atomic vdW radii $\\{R_{\text{vdW},A}\\}$ might be taken from crystal
structure data (e.g., Bondi’s radii set of radii and its subsequent
extensions),114, 115, 116 or might simply be parameters of the model.113 For
vdW cavities, $R_{\text{probe}}=0$ and a typical scaling factor is
$\alpha=1.2$,9 whereas one generally does not scale the vdW radii for SAS
cavities. (As discussed in Ref. 9, the choice $\alpha$ = 1.2 does not result
from any kind of elaborate fitting procedure and is intended only as a rough
guide; common choices range from $\alpha=1.1$–$1.4$.) As an example of just
how sensitive $\mathcal{G}_{\text{elst}}$ is to cavity construction, Table 3
reports calculations for the two tautomers of glycine from Fig. 6, using
various cavity definitions. A change from $\alpha=1.2$ to either $\alpha=1.1$
or $\alpha=1.3$ results in changes of anywhere from 3–9 kcal/mol in
$\mathcal{G}_{\text{elst}}$.

Table 3:  Aqueous solvation energies computed with the SS(V)PE model for two
tautomers of glycine using atomic radii in Eq. (3.7). Electronic structure
calculations were performed at the B3LYP/6-31+G* level.

Cavity | $\alpha$ | $R_{\text{probe}}$ |  $\mathcal{G}_{\text{elst}}$ (kcal/mol)
---|---|---|---
|  | (Å) | amino acid | zwitterion
vdWa | 1.0 | 0.0 | $-26.1$ | $-68.2$
vdWa | 1.1 | 0.0 | $-20.2$ | $-56.0$
vdWa | 1.2 | 0.0 | $-16.1$ | $-46.8$
vdWa | 1.3 | 0.0 | $-12.9$ | $-39.0$
vdWa | 1.4 | 0.0 | $-10.6$ | $-32.6$
SASa | 1.0 | 0.2 | $-18.0$ | $-51.4$
SASa | 1.0 | 1.4 | $-4.3$ | $-13.9$
isodensityb |  |  | $-16.4$ | $-48.1$

  * •

aUsing $R_{\text{vdW}}$ = 1.10 Å (H), 1.70 Å (C), 1.55 Å (N), and 1.52 Å (O),
discretized using SwiG with 302 points per atom. bUsing an isocontour
$\rho_{0}=0.001$ a.u. and 1,202 grid points.

When the solute is described using electronic structure theory, a more
satisfying cavity definition uses the solute’s own electron density.128, 129,
130, 96 It is possible to settle on a numerical isocontour value (generally
$\rho_{0}\sim 0.001$ a.u.) that appears to have some universal validity,125,
210 and results using $\rho_{0}=0.001$ a.u. are shown in Table 3. It is
noteworthy that the standard vdW scaling factor $\alpha$ = 1.2 affords the
best agreement with the isodensity result for both tautomers of glycine,
despite their very different solvation energies. It is also noteworthy that
the two SAS cavities in Table 3, using values of $R_{\text{probe}}$ that are
common in biomolecular Poisson-Boltzmann calculations,119 afford very
different solvation energies. The probe radius $R_{\text{probe}}$ = 1.4 Å is a
realistic measure of the size of a water molecule (based on radial
distribution functions for liquid water),118 but affords rather small values
of $\mathcal{G}_{\text{elst}}$. In contrast, $R_{\text{probe}}$ = 0.2 Å is
much too small to represent an actual water molecule but affords solvation
energies much closer to values obtained using the isodensity and canonical vdW
surface definitions. Values $R_{\text{probe}}=0.2$–$0.3$ Å are consistent with
typical solvent probe radii used since the early days of continuum solvation
models in quantum chemistry.3, 4, 5

Although the isodensity cavity construction is an appealing choice on physical
grounds, existing algorithms to compute this surface are subject to occasional
failure for certain molecular geometries.129, 130 In principle, these
difficulties could likely be overcome using an implementation based on the
“marching cubes” algorithm that is well known in computer graphics.212, 213 A
more fundamental problem is that for an isodensity cavity the surface normal
vector $\mathbf{n}_{\mathbf{s}}$, which is needed to define the ASC-PCM
double-layer operator $\hat{D}$, depends on the density:128, 129, 96

$\mathbf{n}_{\mathbf{s}}=-\frac{\hat{\bm{\nabla}}\rho(\mathbf{s})}{\|\hat{\bm{\nabla}}\rho(\mathbf{s})\|}\;.$
(3.8)

The surface area associated with the discretization point $\mathbf{s}_{i}$
inherits a density dependence as well, and this significantly complicates the
formulation of analytic energy gradients. To date, these have never been
published for the combination of PCM (or related models) with an isodensity
cavity construction. This limitation could be overcome using an analytically-
differentiable pseudo-density, as has sometimes been used in biomolecular
electrostatics calculations,214, 19 perhaps using a superposition of frozen
atomic densities. Although this might remove some arbitrariness from the
selection of atomic radii, however, it would not represent a self-consistent
determination of the cavity surface that could deform to reflect changes in
the molecular electronic structure. Conversely, a pseudo-density that is
determined in order to reproduce the molecular electrostatic potential
$\varphi^{\rho}(\mathbf{r})$ may afford better solvation energies as compared
to vdW radii that do not respond to the electronic structure,215 but
reintroduces the problem of how to compute the analytic gradient. As such, it
is unclear whether such constructions offer advantages over the simplicity of
the vdW cavity. As a simpler workaround, united-atom radii (in which hydrogen
atoms are not given atomic spheres) have been parameterized in an effort to
reproduce results obtained with an isodensity cavity.126

Within the context of solving Poisson’s equation in three-dimensional space,
this idea has more merit, as one can define “soft” atomic spheres that
interpolate the permittivity $\varepsilon(\mathbf{r})$ between limiting values
inside ($\varepsilon=1$) and outside ($\varepsilon=\varepsilon_{\text{s}}$) of
the cavity.216 Related to this, and offering a more compelling scheme for
self-consistent determination of the solute/continuum interface, is a scheme
that takes $\varepsilon(\mathbf{r})$ to be a functional of the solute’s charge
density, $\varepsilon[\rho](\mathbf{r})$, with limiting values $\varepsilon=1$
near the nuclei and $\varepsilon=\varepsilon_{\text{s}}$ farther away. In
practice, “near” and “far” are determined not by distance but by comparison of
$\rho(\mathbf{r})$ to a pair of parameters $\rho_{\text{max}}$ and
$\rho_{\text{min}}$, the latter of which establishes what constitutes the
“tail” of the density. This idea was originally developed by Fattebert and
Gygi217, 218, 219 then later refined by others;220, 221, 222, 223 see Ref. 24
for a review. In its modern incarnation,221 it uses a permittivity functional

$\varepsilon[\rho](\mathbf{r})=\begin{cases}1&\rho(\mathbf{r})>\rho_{\text{max}}\\\
\exp\big{[}t(\ln\rho(\mathbf{r}))\big{]}&\rho_{\text{min}}<\rho(\mathbf{r})<\rho_{\text{max}}\\\
\varepsilon_{\text{s}}&\rho(\mathbf{r})<\rho_{\text{min}}\\\ \end{cases}$
(3.9)

in which $t(x)$ is a switching function that interpolates smoothly between
values $t(\ln\rho_{\text{min}})=\ln\varepsilon_{\text{s}}$ and
$t(\ln\rho_{\text{max}})=0$, so that $\varepsilon(\mathbf{r})$ achieves the
limits indicated in Eq. (3.9). Inserting this ansatz into Poisson’s equation
[Eq. (2.3)] affords a model in which the dielectric interface is smooth,
rather than sharp as it is in PCMs, yet one where the definition of the
interface is updated self-consistently as the density $\rho(\mathbf{r})$ is
iterated to convergence. For that reason, this approach has been called simply
the self-consistent continuum solvation (SCCS) model.24 The dependence of
$\varepsilon_{\text{s}}(\mathbf{r})$ on the density does mean that the Fock
operator $\delta\mathcal{G}/\delta\rho$ acquires an extra term relative to
what was discussed in Section 2.1, namely221

$\upsilon_{\varepsilon}[\rho](\mathbf{r})=-\frac{1}{8\pi}\big{\|}\hat{\bm{\nabla}}\varphi(\mathbf{r})\big{\|}^{2}\left(\frac{\delta\varepsilon[\rho]}{\delta\rho(\mathbf{r})}\right)\;.$
(3.10)

The SCCS model is increasingly being used in ab initio simulations of
materials, e.g., to model the aqueous electrolyte/solid-state interfaces
relevant in electrochemistry.223, 224, 225, 24, 226, 227, 228, 71, 229, 230
Some of that work points to limitations of the linear dielectric model itself
(i.e., the assumption that $\mathbf{P}\propto\mathbf{E}$), because the
rotational response of the water molecules saturates at the electrode
interface and consequently the susceptibility is smaller than in bulk
water.227, 228 Limitations in the linearized Poisson-Boltzmann description of
electrolyte effects have also been demonstrated.226, 227, 228

###  3.4 Linear-scaling algorithms

As noted above, the electrostatic solvation energy obtained from IEF-PCM
should be exactly equivalent to that obtained by solving Poisson’s equation,
in the case of a classical solute for which there is no outlying charge. (This
equivalence holds only up to discretization errors, but those are controllable
and can be driven to zero if systematically-improvable grids are employed.) It
is therefore surprising that the PCM formulation of the classical continuum
problem has seen very little use in biomolecular electrostatics calculations.
Such calculations are almost always performed in water, meaning that the
simpler C-PCM should be essentially exact, and modifications have been
proposed to treat the linearized Poisson-Boltzmann problem,155, 156, 164, 68
including modifications to simulate an ion exclusion layer.68 The Poisson-
Boltzmann equation in three dimensions [Eq. (2.12)] is typically solved using
a finite-difference scheme for the Laplacian operator,231, 39 and this poses
problems for molecular dynamics simulations with implicit solvent, because the
forces need not be continuous when discretized in this way. Although
significant effort has gone into obtaining high-quality forces,232, 233, 234,
235, 236, 161, 162, 163 from a certain point of view this looks like
engineering an elaborate means of escape from the very deep hole created by a
numerical framework that admits discontinuities. A better strategy is not to
get trapped in that hole in the first place. The SwiG discretization procedure
for PCMs (Section 3.2) was designed as a starting point that is free of
discontinuities.

All of this suggests that SwiG-PCMs could make very attractive replacements
for the finite-difference Poisson-Boltzmann solvers that are commonplace in
biomolecular electrostatics calculations. Efforts to do so will quickly run up
against the size of the PCM matrix equation, which is equal to the number of
discretization grid points, itself proportional to the number of solvent-
exposed atoms. This means that a straightforward solution of Eq. (3.2), based
on constructing $\mathbf{K}_{\varepsilon}^{-1}\mathbf{Y}_{\varepsilon}$ or its
equivalent, incurs a CPU cost of ${\cal O}(N^{3}_{\text{atoms}})$ and a memory
footprint of ${\cal O}(N^{2}_{\text{atoms}})$, albeit with significant
prefactors in both cases that reflect the number of discretization points per
atom. In QM/MM calculations with continuum boundary conditions, it is the size
of the MM region that dictates the matrix dimension and for small QM regions
with large MM environments, one can easily encounter scenarios wherein the
cost of classical electrostatics (PCM part of the calculation) exceeds the
cost of doing the quantum mechanics!19

A straightforward solution to that problem is to introduce iterative solvers
that do not require storage or inversion of the matrix
$\mathbf{K}_{\varepsilon}$,237, 238, 19 such as conjugate gradient (CG) or
biconjugate gradient algorithms.19 The bottleneck in these methods is
computing the Coulomb interactions between surface discretization points, and
this can be accelerated using the fast multipole method (FMM),238 in either
its original formulation or using a simpler tree-code approach.141, 19
Parallelization strategies have also been discussed.19 Data for polyalanine
helices (Fig. 7) show that wall times for iterative solution of the PCM
equation (to obtain surface charges) can be reduced to a few seconds, even for
(Ala)4000 with $5.6\times 10^{6}$ surface discretization points. Proof-of-
concept MM/PCM molecular dynamics simulations have been reported in which the
atomistic MM region consists of a segment of DNA bound to a histone (21,734
classical atoms from a force field), with $\approx$124,000 point charges used
to discretize the surface.19

Fig. 7:  Timings and parallel scalability data for linear-scaling PCM solvers
applied to polyalanine helices, using a classical force field description of
the solute. (a) Strong-scaling data for a CG-FMM algorithm applied to
(Ala)250, running on 12 shared-memory cores per node. (b) Weak-scaling data
for (Ala)n helices of increasing length, versus the number of Lebedev grid
points used to discretize the cavity surface. (c) Comparison of timing data
for CG-FMM versus ddCOSMO for (Ala)n. The CG-FMM data in (c) are the same as
those in (b), but all ddCOSMO calculations were run on a single 12-core node.
Data in (a) and (b) are from Ref. 19 and ddCOSMO data in (c) are from Ref.
165.

An alternative algorithm for fast iterative solution of the PCM equation is
based on a domain decomposition procedure.165, 166, 167, 239, 240, 241 Here,
the solute cavity is divided into overlapping domains, each consisting of a
single atomic sphere, for which the PCM equation can be solved analytically.
These single-sphere solutions form the basis for an iterative solution of the
PCM equation for the full domain, which can be formulated as sparse matrix
equations that can be solved in ${\cal O}(N_{\text{atoms}})$ time. This method
was originally developed for C-PCM and is thus named “ddCOSMO”,165, 166, 167,
239, 240 although it has now been extended to arbitrary PCMs (“ddPCM”).241
Timing data on just a single compute node [Fig. 7(c)] show that the method is
competitive with a parallel CG-FMM solver. The ddCOSMO algorithm has been
implemented for MM/PCM molecular dynamics in the Tinker-HP code,242 using a
switching function is used to guarantee smooth forces.243 Such an approach
seems like a practical replacement for finite-difference Poisson-Boltzmann
solvers that is inherently free of problems with discontinuous forces.

##  4 Solvation Energies

Perhaps the single most important property afforded by a solvation model is
the free energy of solvation, $\Delta_{\text{solv}}\mathcal{G}$. Using a
solvation model to compute $\Delta_{\text{solv}}\mathcal{G}^{\circ}$ for the
reactant and product species in a chemical reaction, combined with a gas-phase
calculation of the reaction energy $\Delta_{\text{rxn}}\mathcal{U}$ and free
energy

$\Delta_{\text{rxn}}\mathcal{G}^{\circ}\mbox{[gas]}=\Delta_{\text{rxn}}\mathcal{U}-RT\ln\left(\frac{Z_{\text{products}}}{Z_{\text{reactants}}}\right)\;,$
(4.1)

where $Z$ represents the partition function, one may obtain a value for the
solution-phase reaction energy,
$\Delta_{\text{rxn}}\mathcal{G}^{\circ}\mbox{[solv]}$. This procedure is
illustrated by the thermodynamic cycle shown in Fig. 8. However, the continuum
electrostatics problem discussed heretofore defines only the electrostatic
contribution to $\Delta_{\text{solv}}\mathcal{G}^{\circ}$, and this is
generally insufficient to predict accurate solvation energies. This point is
demonstrated vividly by statistical errors in IEF-PCM solvation energies,244
summarized for aqueous solutes in Table 4. Errors with respect to experimental
solvation energies average about 6 kcal/mol for charge-neutral solutes, 8
kcal/mol for anions, and 13 kcal/mol for cations, despite much smaller errors
with respect to exact Poisson electrostatics (see Table 1). In nonaqueous
solvents, IEF-PCM errors with respect to experimental values of
$\Delta_{\text{solv}}\mathcal{G}^{\circ}$ are 4.9 kcal/mol for neutral solutes
and 12.4 kcal/mol for ions.244 Similarly poor results confirmed using other
electrostatics-only PCM approaches including C- and D-PCM.245, 220

Fig. 8:  Thermodynamic cycle connecting gas- and solution-phase reaction
energies for $\rm A+B\rightarrow C$. The gas-phase value of
$\Delta_{\text{rxn}}\mathcal{G}^{\circ}$ contains the vibrational entropy
change in addition to the reaction enthalpy, and a continuum solvation model
is used to compute $\Delta_{\text{solv}}\mathcal{G}^{\circ}$ for each reactant
and product species. This completes the thermodynamic cycle to compute the
solution-phase value of $\Delta_{\text{rxn}}\mathcal{G}^{\circ}$. Changes in
shape signify that geometries of A, B, and C may be different in solution than
they are in the gas phase, in which case the solvation energies
$\Delta_{\text{solv}}\mathcal{G}^{\circ}$ should include a term representing
the gas-phase deformation energy.  Table 4:  Mean unsigned errors (MUEs) for
aqueous hydration energies ($\Delta_{\text{hyd}}\mathcal{G}^{\circ}$), using
continuum solvation models, for small molecules in the Minnesota Solvation
Database.246, 247, 248, 249 Estimated uncertainties in the experimental data
are $\pm 0.2$ kcal/mol for neutral solutes246, 247 and $\pm 3$ kcal/mol for
ions.247

Data | $N_{\text{data}}$ | MUE (kcal/mol)
---|---|---
Set |  |  |  | PCM-based |  | Poisson-based
|  | IEF-PCMa | SM12b | SMDc | CMIRSd |  | SCCSe | soft-spheref
neutrals | 274 | $5.7$ | $0.7$ | $0.8$ | $0.8$ |  | $1.1$ | $1.1$
cations | 52 | $12.7$ | $3.5$ | $3.2$ | $1.8$ |  | $2.3$ | $2.1$
anions | 60 | $8.0$ | $3.8$ | $6.2$ | $2.8$ |  | $5.5$ | $3.0$
all ions | 112 | $9.7$ | $3.7$ | $4.7$ | $2.4$ |  | $4.0$ | $2.6$

  * •

aB3LYP/6-31G* level using united-atom radii,126 from Ref. 244. bB3LYP/6-31G*
level using ChElPG charges, from Ref. 248, cmPW1PW/6-31G* level, from Ref.
250, dB3LYP/6-31G* level using isodensity cavity with $\rho_{0}=0.001$ a.u.,
from Ref. 251. ePBE in a plane-wave basis, from Ref. 252. fPBE in a wavelet
basis, from Ref. 216.

To do better, these electrostatic solvation models must be augmented to
include nonelectrostatic interactions ($\mathcal{G}_{\text{nonelst}}$),253 a
catch-all term that includes

  * •

cavitation, meaning the energy required to carve out a space in the continuum
solvent;

  * •

Pauli repulsion, i.e., short-range repulsive interactions with the solvent
molecules;

  * •

dispersion, which is non-specific but attractive, and finally

  * •

hydrogen-bonding between solute and solvent.

Harder-to-define entropic changes to the solvent make a contribution to
$\Delta_{\text{solv}}\mathcal{G}^{\circ}$ as well,244 although this effect
could be classified as part of the cavitation (free) energy. Indeed, the
phenomena on the list above are not wholly independent, e.g., short-range
repulsion is also related to cavitation, and hydrogen bonding has
contributions from electrostatics, repulsion, and dispersion. In protic
solvents, and for solutes where hydrogen bonding is important, one or more
explicit solvent molecules could be included as part of the atomistic solute,
thereby removing the largest part of the last item on the list. Semicontinuum
or “cluster + continuum” approaches of this kind are common, e.g., for
calculating hydration energies of ions.8, 254 An important special case is
p$K_{\text{a}}$ calculations,255, 256, 257, 8 corresponding to the ionization
reaction

$\mbox{HA(aq)}\rightarrow\mbox{H}^{+}\mbox{(aq)}+\mbox{A}^{-}\mbox{(aq)}\;.$
(4.2)

Even with one or more explicit solvent molecules to capture hydrogen bonding,
however, the remaining contributions to $\mathcal{G}_{\text{nonelst}}$ must be
added to the electrostatic model in order to obtain accurate reaction
energetics. It has occasionally been noted that semicontinuum predictions of
hydration energies $\Delta_{\text{hyd}}\mathcal{G}^{\circ}$ may fail to
converge as the number of explicit water molecules is increased.258, 259, 260,
261 In truth, there is no reason why such a calculation should converge to
$\Delta_{\text{hyd}}\mathcal{G}^{\circ}$, at least in the absence of sampling
over the explicit water molecules, because such an approach treats entropic
contributions in an unbalanced way.

Nonelectrostatic models appropriate for use with PCMs have been developed over
the years on an ad hoc basis.9, 13 These have often appeared under the name of
“Miertuš-Scrocco-Tomasi” solvation model,98, 262, 263, 253 a term that in
practice means a PCM used in conjunction with a SASA-type parameterization of
$\mathcal{G}_{\text{nonelst}}$ along the lines of what is discussed in Section
4.1. These models have retained something of a do-it-yourself aesthetic and
have not been packaged into “canned” or “black-box” solvation models, a fact
that has occasionally led to hard feelings when errors in IEF-PCM solvation
energies (such as those in Table 4) are pointed out.253, 264 Section 4.1
examines one particular class of black-box solvation model, namely, the SM$x$
models developed by Cramer, Truhlar, and co-workers.247, 265, 266, 244, 250,
248, 249 An even blacker-box approach is “COSMO-RS”,267, 268 which is an
extension of the COSMO method168 designed for “real solvents”. The formulation
and parameters of COSMO-RS are proprietary, however, so the black box cannot
be opened and this method’s performance cannot be independently implemented,
verified, or tested. As such, COSMO-RS is not discussed any further.

###  4.1 SM$x$ and SASA-based models

For computing free energies of solvation, the most successful and popular
models within the quantum chemistry community are the “SM$x$” models.244 These
have version numbers through $x=12$,248 and later $x=\text{D}$ (for
“density”).249 Only the most recent versions (SM12 and SMD) are included in
Table 4.

To understand how these models work, we first take a step back to discuss
implicit solvation in classical biomolecular simulations. Owing in no small
part to difficulties in obtaining stable numerical forces for Poisson-
Boltzmann electrostatics,236, 161, 162, 163 as discussed in Section 3.4, the
most popular implicit solvation models in that context are not based directly
on solution of the Poisson-Boltzmann equation but instead are methods known as
generalized Born (GB) models,269, 270 whose name is a nod to the Born ion
formula of Eq. (2.17). GB models use the Born model as a pairwise ansatz to
compute the total electrostatic solvation energy according to

$\mathcal{G}^{\text{GB}}_{\text{elst}}=-\frac{1}{2}\left(\frac{1}{\varepsilon_{\text{in}}}-\frac{1}{\varepsilon_{\text{out}}}\right)\sum_{A,B>A}^{\text{atoms}}\frac{Q_{A}Q_{B}}{f_{AB}}\;.$
(4.3)

This expression allows for the possibility of a dielectric constant
$\varepsilon_{\text{in}}$ that is different from unity inside of the solute
cavity. The quantity $f_{AB}^{-1}$ in Eq. (4.3) is a parameterized, effective
Coulomb potential, the canonical example of which is271, 269, 270

$f_{AB}=\left[R_{AB}^{2}+\bar{R}_{A}\bar{R}_{B}\exp\big{(}-R_{AB}^{2}/4\bar{R}_{A}\bar{R}_{B}\big{)}\right]^{1/2}\;.$
(4.4)

Here, $R_{AB}=\big{\|}\mathbf{R}_{A}-\mathbf{R}_{B}\big{\|}$ is the
interatomic distance and the quantity $Q_{A}^{2}/f_{AA}$ is the Coulomb self-
energy in the Born ion model, with an effective cavity radius
$f_{AA}\equiv\bar{R}_{A}$. As such, the quantities $\\{\bar{R}_{A}\\}$ are a
set of effective radii that measure the “electrostatic size” of each atom in
the molecule + continuum environment.272, 273, 274 For a classical solute, in
which $\rho(\mathbf{r})$ is a collection of atom-centered point charges
$\\{Q_{A}\\}$, these radii can be computed exactly by solving either Poisson’s
equation,273 or else its PCM analogue,275 once per atomic charge $Q_{A}$, in a
cavity representing the entire molecule. This procedure effectively forces the
pairwise GB ansatz in Eq. (4.3) to reproduce the exact electrostatic energy
defined by the continuum electrostatics problem, and these values of
$\bar{R}_{A}$ have been called the “perfect” Born radii.273, 270 The procedure
just described is not a practical way to determine $\bar{R}_{A}$ but can be
useful to generate data sets that may suggest new analytical forms for GB
models.276, 277, 275, 278, 19

In practice, the radii $\\{\bar{R}_{K}\\}$ are typically determined using
various surface integration procedures.274 The most popular of these is272,
269, 274, 270

$\bar{R}_{K}=\left(\frac{1}{R_{\text{vdW},K}}+\frac{1}{4\pi}\int_{\mathbf{r}\in\bm{\Omega}}\frac{\theta\big{(}\|\mathbf{r}-\mathbf{R}_{K}\|-R_{\text{vdW},K}\big{)}}{\|\mathbf{r}-\mathbf{R}_{K}\|^{4}}\;d\mathbf{r}\right)^{-1}\;.$
(4.5)

Here, $\theta$ is a step function with $\theta(x)=1$ for $x>0$ and
$\theta(x)=0$ for $x\leq 0$; this limits the integration domain to the region
of the solute cavity ($\mathbf{r}\in\bm{\Omega}$) that lies outside of the vdW
sphere for atom $K$. The expression for $\bar{R}_{K}$ in Eq. (4.5) comes from
the so-called Coulomb-field approximation,272, 269, 270 a charge-in-a-sphere
model for the reaction-field potential. This is thought to overestimate
electrostatic size, especially for charges close to the cavity surface, and
alternatives have been developed,274, 270 but Eq. (4.5) is the form used in
the SM$x$ models that are described below. The requisite integral can be
evaluated numerically using concentric atomic radical shells,279 but its
gradient with respect to nuclear displacements, which is needed to obtain
$d\mathcal{G}^{\text{GB}}_{\text{elst}}/dx$, is complicated when defined in
terms of the solid geometry of interpenetrating vdW spheres.190 Both the
integral in Eq. (4.5) and its gradient would be straightforward to evaluate
using SwiG discretization; see Fig. 5(b).

The GB ansatz in Eq. (4.3) specifies only the electrostatic component of the
solvation energy. The model is completed by adding a nonelectrostatic term
$\mathcal{G}_{\text{nonelst}}$, for which two forms are common. The simpler of
the two consists of a term proportional to the volume of the solute cavity, to
model the cavitation energy, and another term proportional to the SASA, which
serves to model dispersion:

$\mathcal{G}_{\text{nonelst}}=\beta
V_{\text{cavity}}+\gamma\times(\text{SASA})\;.$ (4.6)

The quantities $\beta$ and $\gamma$ are empirical parameters. An alternative
form that is also widely used is

$\mathcal{G}_{\text{nonelst}}=\sum_{K}^{\text{atoms}}\gamma_{K}\times\text{SASA}(K)\;,$
(4.7)

where SASA($K$) is the solvent-exposed surface area of atom $K$ and the
parameter $\\{\gamma_{K}\\}$ have units of surface tension. Note that there is
no volumetric term in Eq. (4.7). It has been argued, based on the scaled-
particle theory of hard-sphere fluids,280 that for small molecules the
cavitation energy ought to be parametrizable in terms of the solvent-exposed
surface area,149, 281, 98, 282 although the form in Eq. (4.7) has been used
for macromolecules as well.283 Models that combine GB electrostatics [Eq.
(4.3)] with force-field charges to represent the solute, using either Eq.
(4.6) or Eq. (4.7) for $\mathcal{G}_{\text{nonelst}}$, are known as “MM/GBSA”
methods.284, 285, 286, 287, 288, 289, 290, 291, 292 There is an analogous set
of “MM/PBSA” methods that substitute Poisson-Boltzmann electrostatics in place
of the GB model.284, 285, 286, 287, 288, 289, 290, 291, 292, 293 The MM/GBSA
method is the most widely-used implicit solvation model in biomolecular
simulation,270 but both approaches are popular in drug-discovery applications,
e.g., for calculating ligand–protein binding affinities.294, 293, 295, 296 In
that context, these methods are often used as low-resolution screening tools
representing a level of sophistication that is intermediate between
“knowledge-based” (but largely physics-free) docking or scoring-function
procedures, and much more expensive atomistic free energy simulations in
explicit solvent.

In quantum chemistry contexts, it is not uncommon to use different definitions
of the cavity surface for the electrostatic and nonelectrostatic interactions.
Models for $\mathcal{G}_{\text{nonelst}}$ introduced above suggest a SAS
cavity but numerical results for $\mathcal{G}_{\text{elst}}$ in Table 3
suggest that atomic radii $R=1.2R_{\text{vdW}}$ correspond more closely to
isodensity results as compared to SAS radii
($R=R_{\text{vdW}}+R_{\text{probe}}$), at least when $R_{\text{probe}}$
corresponds to a realistic estimate of the size of a solvent molecule. The
SM$x$ models,244 which are the most popular methods for predicting
$\Delta_{\text{solv}}\mathcal{G}$ in quantum chemistry calculations, exemplify
this distinction. The SM8265 and SM12248 variants, for example, employ GB
electrostatics based on radii obtained from Eq. (4.5), along with atom-
centered charges derived from the quantum-mechanical charge density by means
of certain charge models, CM$x$.297, 298, 247, 299, 300 The latter are
empirically-parameterized improvements upon standard atomic partial charge
prescriptions (Mulliken, Hirshfeld, etc.), with additional parameters designed
to improve the dipole moments obtained from the partial charges, since
Mulliken and Hirshfeld charges known to do a poor job of reproducing the
dipole moment corresponding to the SCF charge density from which they were
obtained.301 CM$x$ charges inherit the basis-set sensitivity of the underlying
wave function-derived charges, therefore the SM$x$ models through SM12 are
each parameterized for use with particular (small) basis sets, e.g. 6-31G*.
The SMD model249 was introduced to overcome this limitation, by substituting
IEF-PCM electrostatics in place of a GB model, with the result that SMD can be
used with arbitrary basis sets whereas other SM$x$ models should only be used
in conjunction with the basis sets for which they were parameterized.

The SM$x$ models are completed by adding a nonelectrostatic term of the form
in Eq. (4.7). Atomic “surface-tension” parameters $\\{\gamma_{K}\\}$ are
themselves empirically-fitted functions of certain “solvent descriptors”,
which include the surface tension of the solvent itself, its refractive index
(which is related to polarizability), and certain Lewis acidity parameters.302
The result is a “universal” model,244 in the sense that once the fitting
procedure is completed, parameters $\\{\gamma_{K}\\}$ are available for any
solvent for which the various solvent descriptors are known. The SM8 model
additionally parameterizes $\mathcal{G}_{\text{nonelst}}$ as a function of
temperature.303, 304

Error statistics for small-molecule hydration free energies
($\Delta_{\text{hyd}}\mathcal{G}^{\circ}$), obtained using SM12 and SMD, are
listed in Table 4. (The SM8 model has slightly larger errors for ions; see
Ref. 244.) The comparison between IEF-PCM and SMD is especially revealing
because these two approaches use the same approach for
$\mathcal{G}_{\text{elst}}$: IEF-PCM is just SMD without
$\mathcal{G}_{\text{nonelst}}$. Evidently, the nonelectrostatic contribution
dramatically reduces errors in $\Delta_{\text{hyd}}\mathcal{G}^{\circ}$ for
both cations and for charge-neutral solutes. Errors for anions are reduced as
well although the impact of $\mathcal{G}_{\text{nonelst}}$ is more muted in
that case, which may simply reflect the difficulty of the anion solvation
problem. In any case, it is clear that nonelectrostatic interactions must be
included to obtain accurate solvation energies. For properties that depend on
only a single geometry (e.g., for some spectroscopic applications), the
nonelectrostatic contributions may be less important but to compute a reaction
profile—reactant, product, and transition state energies in solution—one might
worry that changes in geometry along the reaction coordinate could affect the
nonelectrostatic contribution differently for the various chemical species
involved.

Error statistics for solvents other than water are on par with results for
aqueous solutes, e.g., mean unsigned errors (MUEs) of 0.6 kcal/mol for neutral
solutes and 4.5 kcal/mol for ions, in the case of SM12,248 versus 0.7 kcal/mol
for neutrals and 4.2 kcal/mol for ions using SMD.249 Prediction of
$\Delta_{\text{solv}}\mathcal{G}^{\circ}$ in nonaqueous solvents is necessary
in order to predict partition coefficients between different solvents. In
particular, the octanol/water partition coefficient (equilibrium constant
$K_{\text{ow}}$) is a common measure of lipophilicity (or conversely,
hydrophobicity),305, 306, 307 and is related to solvation energies according
to

$\Delta_{\text{solv}}\mathcal{G}^{\circ}[\text{octanol}]-\Delta_{\text{solv}}\mathcal{G}^{\circ}[\text{water}]=-RT\ln
K_{\text{ow}}\;.$ (4.8)

The value of $K_{\text{ow}}$ is widely used in drug-discovery
applications,308, 309, 307 and atomic decomposition of terms contributing to
$\Delta\Delta\mathcal{G}^{\circ}$ in Eq. (4.8) has been used to determine
similarity indices for predicting quantitative structure–activity
relationships.307 For environmental toxicology purposes, $K_{\text{ow}}$ is an
important physical parameter to determine for any new compound.310 The
octanol/water partition coefficient was the subject of a recent blind
challenge for theoretical methods,311 with the SM$x$ models emerging as
amongst the best performers with a root-mean-square error of 0.44 in units of
$\log_{10}K_{\text{ow}}$.312

Note that by using Eq. (4.7) for $\mathcal{G}_{\text{nonelst}}$, the SM$x$
models do not contain any term proportional to the cavity volume, as is often
present in MM/GBSA and MM/PBSA models, and cavitation effects are therefore
included by means of area-dependent parameters. Although formally justified by
appeal to scaled-particle theory,98 that argument is probably most convincing
for small-molecule solutes. As such, and the success of SM$x$ may partly
reflect the fact that it was parameterized using experimental solvation
energies of mostly small molecules. (The largest molecules in the training set
are $n$-hexadecane and ethyloctadecanoate, at 51 and 63 atoms, respectively,
but most of the solutes are much smaller.250, 249) Small solutes have limited
conformational flexibility and there may not be too much difference in the
volumes of different conformers, therefore a term that explicitly accounts for
cavitation may be largely redundant and not required to obtain good solvation
energies. In contrast, MM/GBSA and MM/PBSA methods are usually parameterized
for (or at least tested on) proteins and other macromolecular solutes, where
folded and unfolded structures likely have rather different volumes. The
benchmarks for MM/GBSA and MM/PBSA are usually solvation energies obtained
from molecular dynamics simulations in explicit solvent. Since the SM$x$
models have not been scaled up to macromolecules, it remains unclear how they
would perform in that context.

SM$x$ approaches have historically been the go-to models for computing
solvation energies in quantum chemistry but two newer approaches based
directly on Poisson’s equation seem very promising.252, 216 One of these is
the SCCS approach that was described in Section 3.3, which uses a density-
dependent permittivity functional $\varepsilon[\rho](\mathbf{r})$ to determine
the solute cavity self-consistently alongside the charge density
$\rho(\mathbf{r})$, with Poisson’s equation solved in three-dimensional space
to obtain $\mathcal{G}_{\text{elst}}$. Combined with the simple two-parameter
model for $\mathcal{G}_{\text{nonelst}}$ in Eq. (4.6), the SCCS approach
achieves a statistical accuracy for aqueous solvation energies that is
comparable to SM12 or SMD (see Table 4).252 This is despite using just two
empirical parameters in $\mathcal{G}_{\text{nonelst}}$, along with two
parameters [and an ansatz, Eq. (3.9)] for the cavity determination. The SM$x$
models have a considerably larger number of empirical parameters, although
they are designed to be “universal” models for all solvents, whereas the SCCS
model is parameterized only for water. Improvement upon SCCS results for
anions is obtained using a “soft-sphere” model,216 which constructs the cavity
from vdW spheres but then interpolates the dielectric function (from
$\varepsilon=1$ to $\varepsilon=78$) over a narrow switching region centered
on the vdW radius.216 When combined with the two-parameter model for
$\mathcal{G}_{\text{nonelst}}$, the soft-sphere solvation model affords a MUE
of 3.0 kcal/mol in $\Delta_{\text{hyd}}\mathcal{G}^{\circ}$ for anions (Table
4), which is on par with the estimated uncertainty (2–3 kcal/mol) in the
experimental data themselves.313, 247

###  4.2 Physics-based models

Occasionally there have been attempts to put the nonelectrostatic
contributions to the solvation energy on a more rigorous footing,314, 315,
316, 317, 318, 319, 320, 321, 322, 323, 324 or at least to develop
parameterized models that are more closely connected to the physics of
intermolecular interactions as compared to the wholly empirical SASA-type
approaches.325, 326, 327, 328, 329 As a simple example of the latter variety,
one might borrow from QM/MM methodology and assume classical interaction
potentials, e.g., of Lennard-Jones type that are centered on the solute atoms,
with an $r^{-6}$ distance dependence for dispersion. This affords a model of
the form326, 329

$\mathcal{G}_{\text{disp}}=\sum_{A}^{\text{solvent}}\bar{\rho}_{A}\sum_{B}^{\text{solute}}\gamma_{\\!AB}\sum^{\text{grid}}_{i\in
B}a_{i}\left(\frac{\mathbf{r}_{iB}\bm{\cdot}\mathbf{n}_{\mathbf{s}_{i}}}{3r_{iB}^{6}}\right)\;,$
(4.9)

where $\bar{\rho}_{A}$ is the average number density of solvent atom $A$, and
the parameters $\gamma_{\\!AB}$ come from a force field. A similar model,
depending for example on $r_{iB}^{-12}$, can be developed to describe the
repulsive contribution, $\mathcal{G}_{\text{rep}}$.327, 329 The quantity
$\mathcal{G}_{\text{disp}}+\mathcal{G}_{\text{rep}}$ is usually then combined
with a cavitation energy of the form

$\mathcal{G}_{\text{cav}}=\sum_{B}^{\text{solute}}\left(\frac{\mbox{SASA}(B)}{4\pi
R_{B}^{2}}\right)\Delta G_{\text{HS}}(R_{B})\;,$ (4.10)

in which $\Delta G_{\text{HS}}(R_{B})$ is the solvation energy for a hard
sphere of radius $R_{B}$, obtained from scaled-particle theory.280 Typically
the atomic radii used in these nonelectrostatic terms are SAS radii, i.e.,
they include a probe radius representing the assumed size of a solvent atom;
see Eq. (3.7). This is the case even if $\mathcal{G}_{\text{elst}}$ is
evaluated using a vdW cavity surface.

This is a classical model for the nonelectrostatic interactions but there have
also been attempts to derive fully quantum-mechanical dispersion–repulsion
models within a continuum framework. One such model, due to Amovilli,314, 315
starts from the generalized Casimir-Polder expression for the dispersion
energy of a supramolecular complex A$\cdots$B,330, 331, 332, 333, 334, 335

$\mathcal{U}_{\text{disp}}^{\text{\scriptsize A}\cdots\text{\scriptsize
B}}=-\frac{\hbar}{2\pi}\left(\frac{e^{2}}{4\pi\varepsilon_{0}}\right)^{\\!2}\int_{0}^{\infty}d\omega\int_{\mathbb{R}^{3}}d\mathbf{r}_{1}d\mathbf{r}_{2}d\mathbf{r}_{1}^{\prime}d\mathbf{r}_{2}^{\prime}\;\left(\frac{\chi^{\text{\scriptsize
A}}(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime}|\dot{\iota}\omega)\;\chi^{\text{\scriptsize
B}}(\mathbf{r}_{2},\mathbf{r}_{2}^{\prime}|\dot{\iota}\omega)}{\|\mathbf{r}_{1}-\mathbf{r}_{1}^{\prime}\|~{}\|\mathbf{r}_{2}-\mathbf{r}_{2}^{\prime}\|}\right)\;,$
(4.11)

in which $\chi(\mathbf{r},\mathbf{r}^{\prime}|\omega)$ is the frequency-
dependent density susceptibility (also known as the polarization
propagator),334, 335 evaluated in Eq. (4.11) at imaginary frequencies. When
the monomer separation $R_{\text{AB}}$ is large, second-order perturbation
theory affords the “uncoupled” approximation,336, 337, 338, 339, 331, 340
which was first derived by London:341, 342

$\mathcal{U}_{\text{disp}}^{\text{\scriptsize A}\cdots\text{\scriptsize
B}}(R_{\text{AB}})=-\left[\frac{3\hbar}{\pi}\int_{0}^{\infty}\bar{\alpha}^{\text{\scriptsize
A}}(\dot{\iota}\omega)\;\bar{\alpha}^{\text{\scriptsize
B}}(\dot{\iota}\omega)\;d\omega\right]\frac{1}{R_{\text{AB}}^{6}}\;.$ (4.12)

The quantity $\bar{\alpha}(\omega)$ is the frequency-dependent isotropic
polarizability. Components of the polarizability tensor, in the spectral
representation, are

$\alpha_{ab}(\omega)=\frac{1}{\hbar}\sum_{n>0}\left[\frac{\langle
0|\hat{\mu}_{a}|n\rangle\langle
n|\hat{\mu}_{b}|0\rangle}{\omega-\omega_{n0}}+\frac{\langle
0|\hat{\mu}_{b}|n\rangle\langle
n|\hat{\mu}_{a}|0\rangle}{\omega+\omega_{n0}}\right]\;,$ (4.13)

where $\omega_{n0}=(E_{n}-E_{0})/\hbar$.343 The term in brackets in Eq. (4.12)
provides a microscopic formula for the Lennard-Jones $C_{6}$ coefficient.

Returning to the more general expression in Eq. (4.11), a sum-over-states
expression (Lehmann representation) can be introduced for the requisite
density susceptibilities, writing them in terms of transition densities
$\rho_{0n}(\mathbf{r})$:333, 335

$\chi(\mathbf{r},\mathbf{r}^{\prime}|\dot{\iota}\omega)=-\frac{2}{\hbar}\sum_{n>0}\frac{\omega_{0n}\;\rho_{0n}(\mathbf{r})\;\rho_{0n}(\mathbf{r}^{\prime})}{\omega_{n0}^{2}+\omega^{2}}\;,$
(4.14)

Inserting this formula for $\chi^{\text{\scriptsize
A}}(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime}|\dot{\iota}\omega)$ in Eq. (4.11)
affords an expression for the dispersion energy of “A-in-B”, i.e., solute A
interacting with an implicit representation of B:314

$\mathcal{G}_{\text{disp}}^{\text{A-in-B}}=\frac{1}{\pi}\int_{0}^{\infty}d\omega\sum_{n>0}\frac{\omega_{n0}^{\text{\scriptsize
A}}}{(\omega_{n0}^{\text{\scriptsize
A}})^{2}+\omega^{2}}\int_{\mathbb{R}^{3}}d\mathbf{r}\int_{\bm{\Gamma}}d\mathbf{s}\left(\frac{\rho^{\text{\scriptsize
A}}_{0n}(\mathbf{r})\;\sigma_{\text{\scriptsize
B}}[\varepsilon^{\text{\scriptsize
B}}(\dot{\iota}\omega),\rho^{\text{\scriptsize
A}}_{0n}](\mathbf{s})}{\|\mathbf{r}-\mathbf{s}\|}\right)\;.$ (4.15)

In writing this equation, the susceptibility $\chi^{\text{\scriptsize B}}$ in
Eq. (4.11) has been eliminated by first introducing a polarization density

$\rho^{\text{\scriptsize
B}}_{\text{pol}}(\mathbf{r}_{2})=\int_{\mathbb{R}^{3}}d\mathbf{r}_{1}^{\prime}d\mathbf{r}_{2}^{\prime}\left(\frac{\chi^{\text{\scriptsize
B}}(\mathbf{r}_{2},\mathbf{r}_{2}^{\prime}|\dot{\iota}\omega)\;\rho^{\text{\scriptsize
A}}(\mathbf{r}_{1}^{\prime})}{\|\mathbf{r}_{1}^{\prime}-\mathbf{r}_{2}^{\prime}\|}\right)\;,$
(4.16)

induced on B by the presence of A. This quantity is then replaced with a
surface charge $\sigma_{\text{\scriptsize B}}(\mathbf{s})$, in the spirit of
the PCM formulation of the continuum electrostatics problem, and which depends
on the frequency-dependent dielectric function
$\varepsilon(\dot{\iota}\omega)$ evaluated at imaginary frequencies. The
function $\varepsilon(\dot{\iota}\omega)$ is the central quantity in the
Lifshitz’s general theory of the Casimir force and dispersion energies.344,
330, 345, 346 In a separate approach, Ninham and coworkers323, 324 use the
frequency-dependent dipole- and higher-order (hyper)polarizabilities for the
solute, $\bar{\alpha}(\dot{\iota}\omega)$ etc., in conjunction with models for
$\varepsilon(\omega)$, to model solute–solvent dispersion.

Models for $\mathcal{G}_{\text{disp}}$ can now be constructed based on Eq.
(4.15) by modeling the dielectric function $\varepsilon(\dot{\iota}\omega)$
and the surface charge density $\sigma_{\text{\scriptsize B}}(\mathbf{s})$
that is induced by various excited states of solute A, via transitions
densities $\rho_{0n}^{\text{\scriptsize A}}(\mathbf{r})$. These transition
densities, along with the corresponding excitation frequencies
$\omega^{\text{\scriptsize A}}_{n0}$ that appear in Eq. (4.15), might be
computed explicitly (since the solute A is described by quantum chemistry),314
or perhaps modeled using SCF eigenvalues.315, 316 A model suggested for the
surface charge is315

$\sigma_{\text{\scriptsize B}}[\rho^{\text{\scriptsize
A}}_{0n}](\mathbf{s})=-\frac{1}{4\pi}\left(\frac{\Omega^{2}}{\Omega^{2}+\omega^{2}}\right)\left(\frac{\varepsilon_{\infty}-1}{\varepsilon_{\infty}}\right)E_{\perp}[\rho^{\text{\scriptsize
A}}_{0n}](\mathbf{s})\;,$ (4.17)

where $E_{\perp}[\rho^{\text{\scriptsize A}}_{0n}](\mathbf{s})$ is the normal
electric field generated by the density $\rho^{\text{\scriptsize
A}}_{0n}(\mathbf{r})$. Modulo a factor of $\Omega/(\Omega+\omega^{2})$, this
is precisely the apparent surface charge of the D-PCM approach [Eq. (2.24)],
albeit generated by the transition density $\rho^{\text{\scriptsize A}}_{0n}$
rather than the ground-state density and with the “optical” dielectric
constant $\varepsilon_{\infty}$ replacing the static dielectric constant
$\varepsilon_{\text{s}}$. (As discussed in Section 5, the value
$\varepsilon_{\infty}$ is the appropriate one for re-polarization upon sudden
or vertical excitation, without orientational contributions from the solvent.)
Regarding the factor of $\Omega^{2}/(\Omega^{2}+\omega^{2})$, the quantity
$\hbar\Omega$ is the characteristic ionization energy of the solvent. [This
comes from the approximation of setting every $\omega_{n0}$ equal to
$\Omega$;342 In practice,
$\hbar\Omega=\varepsilon_{\infty}\times(\text{IE})_{\text{solvent}}$ has been
used.315, 316] The integral over $\omega$ in Eq. (4.15) remains to be
evaluated and this factor interpolates the apparent surface charge between