limits of unity for $\omega=0$, for which the solute sees the full induced
polarization response, and zero as $\omega\rightarrow\infty$ because when
$\omega\gg\Omega$ the excitation frequency is so large that the response
averages to zero. These models have interesting possibilities for the
description of solute–environment dispersion in excited states, which are only
starting to be explored.316, 317

A rather different formulation of the solute–continuum dispersion interaction
has been put forward by Pomogaeva and Chipman,320 who borrow from the nonlocal
correlation energy functional developed by Vydrov and Van Voorhis.347, 348,
349, 350 This functional, usually known as VV10,351 represents an attempt to
incorporate dispersion interactions into density functional theory in a
rigorous way and is itself a simplified form of the nonlocal functional
developed by Langreth and Lundqvist.352, 353, 354, 355 (In a mildly annoying
bit of physicist reductionism, the Langreth-Lundqvist functional is often
known as “the” van der Waals functional, as if such a designation could
possibly be unique.) These nonlocal correlation functionals are already based
on a pairwise ansatz,351 and to use them in conjunction with a continuum
representation of the solvent one simply replaces the density of one
interacting partner with the bulk solvent density, $\bar{\rho}$. The
functional form is then320, 251

$\mathcal{G}_{\text{disp}}=A\int_{\mathbf{r}\in\bm{\Omega}}\left(\frac{\rho(\mathbf{r})\;I(\mathbf{r})}{w[\rho](\mathbf{r})\big{(}w[\rho](\mathbf{r})+\bar{\rho}^{1/2}_{\text{solvent}}\big{)}}\right)\;d\mathbf{r}$
(4.18)

where

$I(\mathbf{r})=\int_{\mathbf{r}^{\prime}\notin\bm{\Omega}}\left(\frac{d\mathbf{r}^{\prime}}{\|\mathbf{r}-\mathbf{r}^{\prime}\|^{6}+\delta^{6}}\right)$
(4.19)

and

$w[\rho](\mathbf{r})=\left(\rho(\mathbf{r})+\frac{3C}{4\pi}\frac{\big{\|}\hat{\bm{\nabla}}\rho(\mathbf{r})\big{\|}^{4}}{\rho(\mathbf{r})^{4}}\right)^{1/2}\;.$
(4.20)

The parameter $C$ in Eq. (4.20) is taken without modification from VV10,348,
349 but a parameter $\delta$ is introduced in Eq. (4.19) to prevent divergence
of the integral. Aside from that, the only additional parameter that has been
introduced (beyond those already present in VV10) is an overall scaling factor
$A$ in Eq. (4.18). The density $\rho(\mathbf{r})$ is the solute’s electron
density and the integral in Eq. (4.18) is evaluated over the solute cavity
($\mathbf{r}\in\bm{\Omega}$). However, the integration domain in Eq. (4.19) is
the region $\mathbf{r}^{\prime}\notin\bm{\Omega}$ that is outside of the
solute cavity, which requires integration of three-dimensional space. (In
practice the discretization need not extend very far beyond the cavity, since
the integrand decays rapidly.) It is worth noting that analytic gradients have
not been implemented for the quantum-mechanical dispersion models described in
this section, although they have been implemented for the QM/MM-style approach
of Eq. (4.9).329 However, the functional form for $\mathcal{G}_{\text{disp}}$
in Eq. (4.18) may better lend itself to analytic gradients as compared to Eq.
(4.5), insofar as the analytic gradient of VV10 has already been reported.349

In contrast to dispersion, a functional form for Pauli repulsion is rather
straightforward. This effect arises due to interpenetration of the tails of
two non-bonded densities, and one functional form that has been suggested is
simply

$\mathcal{G}_{\text{exch}}=B\int_{\mathbf{r}\notin\bm{\Omega}}\rho(\mathbf{r})\;d\mathbf{r}$
(4.21)

where $B$ is an empirical scaling factor.315 (Note also that the integration
domain is over the solvent, $\mathbf{r}\notin\bm{\Omega}$.) An alternative
is320

$\mathcal{G}_{\text{exch}}=B\int_{\mathbf{r}\notin\bm{\Omega}}\big{\|}\hat{\bm{\nabla}}\rho(\mathbf{r})\big{\|}\;d\mathbf{r}\;.$
(4.22)

This latter form is suggested by an exact result for the exchange-repulsion
energy of two hydrogen atoms,356 but both models have been used in practice.
For example, the model in Eq. (4.21) has been used to develop an “extreme
pressure” (XP-) PCM,357, 358 which is based on the thermodynamic relation
$p=-(\partial\mathcal{G}/\partial V)_{T}$ and calculation of analytic
derivatives of $\mathcal{G}_{\text{exch}}$ and $\mathcal{G}_{\text{elst}}$
with respect to the cavity volume. XP-PCM has been used to study how pressures
$p>1$ GPa affect both molecular geometries and the equilibrium positions of
chemical reactions.359, 360 The exchange-repulsion model in Eq. (4.22) has
been used as part of a black-box solvation model that is described next.

In an attempt to develop a first-principles implicit solvation model that can
compete with SM$x$, Pomogaeva and Chipman318, 319, 320, 321, 322 have combined
SS(V)PE electrostatics (based on an isodensity cavity construction) with a
“minimally parameterized” nonelectrostatic model of the form

$\mathcal{G}_{\text{nonelst}}=\mathcal{G}_{\text{disp}}+\mathcal{G}_{\text{exch}}+\mathcal{G}_{\text{FESR}}\;.$
(4.23)

Here, $\mathcal{G}_{\text{disp}}$ and $\mathcal{G}_{\text{exch}}$ are the
models in Eq. (4.18) and Eq. (4.22), respectively, and
$\mathcal{G}_{\text{FESR}}$ is “field-effect short-range” (FESR) term for
hydrogen bonding. The functional form for $\mathcal{G}_{\text{FESR}}$ takes as
inputs the maximum and minimum values of the normal electric field at the
cavity surface, and has the form319

$\mathcal{G}_{\text{FESR}}=C\big{(}\min|E_{\perp}|\big{)}^{\xi}+D\big{(}\max|E_{\perp}|\big{)}^{\xi}\;.$
(4.24)

(The correlation between hydrogen bond strength and local electric fields in
water has been made before.361, 362, 363) This ansatz is likely only
appropriate for small molecules since it would seem to accommodate only a
single hydrogen-bond donor site and a single acceptor site.

The combination of SS(V)PE electrostatics with $\mathcal{G}_{\text{nonelst}}$
in Eq. (4.23) has been called the composite method for implicit representation
of solvent (CMIRS).320, 321, 322, 251 This model contains five empirical
parameters for a given solvent: the linear coefficients $A$, $B$, $C$, and
$D$, along with the exponent $\xi$. (The FESR term is omitted for non-polar
solvents.) Pomogaeva and Chipman determined these parameters, by fitting to
experimental solvation energies, for benzene,320 cyclohexane,320 water,321
dimethyl sulfoxide,322 and acetonitrile,322 although all of those parameters
were later adjusted to fix an error in the original implementation.251
Parameters for methanol have also been reported based on the original
implementation.364 For any given solvent, this is is considerably fewer
parameters than are contained in any of the SM$x$ models, although the latter
are designed as “universal” models in which the nonelectrostatic parameters
$\\{\gamma_{K}\\}$ in Eq. (4.7) are determined only once, and then the model
is available for any solvent whose macroscopic descriptors are available.244
For example, the training set for SM12 contains 92 solvents,248 as compared to
the six for which CMIRS parameters are currently available.

Fig. 9:  Comparison of experimental hydration energies in the Minnesota
solvation database246, 247, 248, 249 to predictions from the CMIRS model: (a)
all solvation energies, including both neutral molecules as well as ions (with
the number of data points indicated in each case), versus (b) results for
charge-neutral solutes only. Reprinted from Ref. 251; copyright 2016 American
Chemical Society.

Of course, water is the solvent to have if you’re only having one, and error
statistics versus experimental hydration energies (Table 4) demonstrate that
CMIRS is somewhat more accurate (on average) than SM12 or SMD values,
especially for ions, despite fewer parameters in the model. The CMIRS error of
2.4 kcal/mol for ions is comparable to the estimated uncertainty (2–3
kcal/mol) in the experimental solvation energies for aqueous ions,313, 247 and
thus represents the practical lower limit for any solvation model trained on
these data. Correlation between CMIRS and experimental hydration energies is
excellent; see Fig. 9. That said, despite all of the physical considerations
that went into the CMIRS model for $\mathcal{G}_{\text{nonelst}}$, and despite
the limited number of parameters used per solvent, an error in the original
implementation $\mathcal{G}_{\text{disp}}$ went unnoticed despite the fact
that it modifies dispersion energies in the training set by up to 8
kcal/mol.251 This was able to escape notice because the $B$ parameter in Eq.
(4.22) absorbs the discrepancy, making $\mathcal{G}_{\text{exch}}$ more
repulsive in order to offset the too-attractive dispersion energy.251 A
similar cancellation between the repulsive cavitation term and the attractive
dispersion term has been noted elsewhere,324 and in fact
$\mathcal{G}_{\text{disp}}+\mathcal{G}_{\text{exch}}$ is often parameterized
together as a single entity in empirical models, including SM$x$ and others
that use the atomic surface tension ansatz of Eq. (4.7). This may hide certain
subtleties, such as the fact that cavitation effects are more important than
dispersion to explain binding affinities of rare-gas guest atoms to
cucurbituril host molecules,365 or that the unfavorable hydration energies of
small nonpolar polymers are well approximated by the cavitation energy
($\Delta_{\text{hyd}}\mathcal{G}\approx\mathcal{G}_{\text{cav}}$),366
suggesting near-cancellation of other effects.

##  5 Nonequilibrium Solvation

###  5.1 Conceptual overview

How does a continuum solvent respond to a sudden change in the solute’s charge
distribution? This question must be considered for electronic spectroscopies,
including absorption to (or emission from) an excited electronic state, or
photoelectron spectroscopy that removes an electron from the solute. The
general theory of time-dependent processes in dielectric materials introduces
a frequency-dependent dielectric function $\varepsilon(\omega)$, such that the
induction field $\mathbf{D}$ responds to a frequency-dependent electric field
$\mathbf{E}$ according to
$\mathbf{D}(\mathbf{r},\omega)=\varepsilon(\omega)\,\mathbf{E}(\mathbf{r},\omega)$.367,
368, 369, 370 In the presence of a time-dependent field, the permittivity
$\varepsilon(\mathbf{r},\omega)$ is complex-valued, in order to describe the
phase lag between $\mathbf{E}$ and $\mathbf{D}$.367, 370 Often, the complex-
valued permittivity is denoted as $\hat{\varepsilon}$ or $\varepsilon^{\ast}$
but we will not do so here. Where needed, we will simply indicate the
frequency dependence explicitly.

The permittivity is real-valued in two important limits, namely,
$\varepsilon_{\text{s}}\equiv\varepsilon(0)$, which is the static (zero-
frequency) limit, and also in the high-frequency limit. The limiting value
$\varepsilon_{\infty}=\lim_{\omega\rightarrow\infty}\varepsilon(\omega)$ is
known as the optical dielectric constant, for reasons that are described
below, and describes the fact that there is always some part of the
polarization that is able to remain in phase with the applied field. Switching
to the time domain and labeling that part of the medium’s response as “fast”
polarization,

$\mathbf{P}_{\text{fast}}(t)=\left(\frac{\varepsilon_{\infty}-1}{4\pi}\right)\mathbf{E}(t)\;,$
(5.1)

the remaining (“slow”) component is then defined by the time-dependent
analogue of Eq. (2.1):367

$\mathbf{D}(t)=\varepsilon_{\infty}\mathbf{E}(t)+4\pi\mathbf{P}_{\text{slow}}(t)\;,$
(5.2)

The frequency components of $\mathbf{P}_{\text{slow}}(t)$ depend on
$\varepsilon(\omega)$, but generally speaking the slow polarization response
consists of vibrational contributions with timescales of $10^{-12}$–$10^{-14}$
s, along with orientational components having timescales slower than
$10^{-12}$ s. For example, the primary relaxation timescale for liquid water
is 8–10 ps under ambient conditions.77, 371, 372, 373

It is possible to model the frequency dependence of $\varepsilon(\omega)$
directly.374, 375, 376, 377, 378 A phenomenological model is

$\varepsilon(\omega)=\varepsilon_{\infty}+(\varepsilon_{\text{s}}-\varepsilon_{\infty})\sum_{k}\frac{c_{k}}{1+\dot{\iota}\omega\tau_{k}}\;,$
(5.3)

where the $\tau_{k}$ are a set of characteristic time constants, representing
microscopic relaxation processes, and $\sum_{k}c_{k}=1$.367 (The version with
only a single timescale was originally introduced by Debye.379, 77, 370) When
such a model is used in the context of continuum solvation theory, the
polarization charge becomes explicitly time-dependent and the permittivity for
“fast” polarization is the real part of $\varepsilon(\omega)$ for frequencies
larger than the perturbation of interest. Such models have been used to
simulate the time-dependent Stokes shift in the fluorescence energy,374, 375,
376 and to simulate the combined response of the molecule and the medium to an
external field that is resonant with an excited state of the solute.378 The
latter application makes the most sense when combined with electronic
structure methods that simulate time-dependent electron dynamics,380, 381 but
these explicitly time-dependent approaches are not considered in this work.

Instead, this work is focused on the nonequilibrium response to a sudden
(vertical or Franck-Condon) change in the solute density. For a vertical
absorption, emission, or ionization process, the nuclei are fixed and the
continuum solvent “molecules” cannot vibrate or reorient. This limits the
continuum solvation response to the fast component of the polarization, which
is electronic in nature and remains in equilibrium with the sudden change in
$\rho(\mathbf{r})$. The slow component is dictated by the solute’s initial
electronic state and cannot respond on the timescale of a vertical excitation
or ionization process, and is therefore out of equilibrium with the final
electronic state. Phenomenologically, this is accomplished by reducing
$\varepsilon(\omega)$ to its limiting values, $\varepsilon_{\text{s}}$ and
$\varepsilon_{\infty}$.

The value of $\varepsilon_{\infty}$ can be related to the polarizability of
the solvent molecules (Lorenz-Lorentz equation),367 but in the present context
is simple a measured property of the solvent, determined from the solvent’s
refractivity $n(\omega)=\sqrt{\varepsilon(\omega)}$.367, 382 Formally this
should be done in the limit $\omega\rightarrow\infty$ but also needs to be
measured away from any resonances and therefore optical wavelengths are often
used, hence the “optical” dielectric constant. Values are often measured at
the sodium D-line, $\lambda=589$ nm,383 and for water the value obtained is
$\varepsilon_{\infty}(\lambda)=1.78$ at $\lambda=589$ nm, as compared to
$\varepsilon_{\infty}(\lambda)=1.95$ at $\lambda=200$ nm.384 In older
literature (and sometimes repeated more recently),385, 386 values of
$\varepsilon_{\infty}=4.0$–$5.5$ are occasionally reported for water,387, 77
the larger value originally being thought to agree better with predictions
from Onsager’s reaction-field theory.387 In fact, the larger value of
$\varepsilon_{\infty}$ is based on a false extrapolation using Debye’s model
of a single relaxation timescale ($\tau_{1}=8$–$10$ ps), whereas permittivity
data that extend to terahertz frequencies reveal at least two distinct
relaxation timescales,388, 389, 390, 373, 391 including a faster process
$\tau_{2}<1$ ps.388, 389, 373 (The microscopic explanation for these
timescales remains a topic of current debate.392, 393) The data therefore
appear to approach a limiting value $\varepsilon_{\infty}\approx 5$ in the
microwave regime,77, 394 but further decay to $\varepsilon_{\infty}\approx 2$
at terahertz frequencies,389, 391 the latter value being consistent with
$n^{2}(\lambda)$ measured at optical wavelengths.384

Table 5:  Static dielectric constantsa ($\varepsilon_{\text{s}}$) and optical
dielectric constantsb ($\varepsilon_{\infty}$) for some common solvents.
(These are given as dimensionless values relative to vacuum permittivity.)

Solvent | $\varepsilon_{\text{s}}$ | $\varepsilon_{\infty}$ | Solvent | $\varepsilon_{\text{s}}$ | $\varepsilon_{\infty}$
---|---|---|---|---|---
$n$-hexane | $1.9$ | $1.89$ | ethanol | $24.3$ | $1.85$
cyclohexane | $2.0$ | $2.03$ | ethylene glycol | $30.9$ | $2.05$
benzene | $2.3$ | $2.25$ | methanol | $33.0$ | $1.77$
toluene | $2.4$ | $2.24$ | nitrobenzene | $34.7$ | $2.41$
diethyl ether | $4.2$ | $1.83^{c}$ | acetonitrile | $36.0$ | $1.81$
chloroform | $4.7$ | $2.08$ | dimethyl acetamide | $39.6$ | $2.07$
dichloromethane | $8.9$ | $2.03$ | dimethyl sulfoxide | $46.6$ | $2.18$
2-propanol | $18.2$ | $1.92$ | water | $78.4$ | $1.78^{d}$
acetone | $20.8$ | $1.85$ | formamide | $109.6$ | $2.10$

  * •

aAt 25∘C, from Ref. 395. bAt 20∘C except where noted, based on refractive
indices $n(\lambda)$ measured at $\lambda=589$ nm, from Ref. 383. cAt 16.5∘C.
dFrom Ref. 396, valid for 20–25∘C.

Looking at a modern tabulation of the data for common solvents (Table 5),383
one finds little variation in refractive indices at visible wavelengths, which
generally range from $n\approx 1.3$–$1.5$. Despite having static dielectric
constants ranging from from $\varepsilon_{\text{s}}=2$–$78$, there is
considerable uniformity in $\varepsilon_{\infty}$, which ranges from
$\varepsilon_{\infty}=1.7$–$2.3$. This reflects the the fact that typical
solvents are closed-shell, small-molecule insulators with band gaps in the
vacuum ultraviolet, and therefore roughly similar molecular polarizabilities.
Inorganic solids may have considerably larger indices of refraction,397 e.g.,
$n=2.43$ at 589 nm (and therefore $\varepsilon_{\infty}=5.90$) for
BaTiO3(s),398 a material used in nonlinear optical applications. These larger
values are attributable to low-lying excited states that facilitate more
substantial electronic polarization and therefore significant dispersion of
light, but this behavior is simply not found in common solvents.

Before introducing a modern computational formalism for nonequilibrium
polarization, we first consider two historical examples. The first is a well-
known result in electronic spectroscopy that relates the Stokes shift
($\Delta\nu=\nu_{\text{abs}}-\nu_{\text{fluor}}$) to the change in the
solute’s dipole moment upon excitation ($\Delta\mu$), and which is known as
the Lippert-Mataga equation,399, 400, 41, 401

$hc\Delta\nu=\text{constant}+\frac{2(\Delta\mu)^{2}}{\bar{R}^{3}}\underbrace{\left[\frac{\varepsilon_{\text{s}}-1}{\varepsilon_{\text{s}}+1}-\frac{\varepsilon_{\infty}-1}{2\varepsilon_{\infty}+1}\right]}_{F(\varepsilon_{\text{s}},\varepsilon_{\infty})}\;.$
(5.4)

In practice this equation is used to determine excited-state dipole moments
(assuming that the ground-state dipole moment is known) by measuring the
Stokes shift in solvents of differing polarity. In fact, a variety of
alternative formulas for this purpose have been suggested,402, 403, 404, 405,
406, 407, 408, 409, 410 along the lines of Eq. (5.4) but differing somewhat in
their treatment of excited-state polarization, which leads to differences in
the form of the “solvent polarity function”
$F(\varepsilon_{\text{s}},\varepsilon_{\infty})$. These models (including the
Lippert-Mataga one) are so crude that often experimental data can be fit
equally well to any one of them.411, 412, 413, 414, 415, 416, 417 More
important is the basic molecular physics that underlies this approach.
Comparison to the model of a dipole in a spherical cavity [Eq. (2.18)] shows
that the physical content of of Eq. (5.4) is to take the difference dipole
moment $\Delta\bm{\mu}$ and solvate it using permittivity
$\varepsilon_{\infty}$ rather than $\varepsilon_{\text{s}}$.

The solvent parameter $\varepsilon_{\infty}$ also makes an appearance in
Marcus’ theory of electron transfer,418, 419, 420, 421, 422, 423 in which the
“outer-sphere” reorganization energy is given by

$\lambda_{\text{outer}}=(\Delta
Q)^{2}\left(\frac{1}{\varepsilon_{\infty}}-\frac{1}{\varepsilon_{\text{s}}}\right)\left(\frac{1}{2R_{\text{D}}}+\frac{1}{2R_{\text{A}}}-\frac{1}{\|\mathbf{r}_{\text{D}}-\mathbf{r}_{\text{A}}\|}\right)\;.$
(5.5)

This formula is derived from what is essentially a nonequilibrium formulation
of the Born ion model [Eq. (2.17)], combined with a Coulomb interaction
between charges centered in a donor sphere (radius $R_{\text{D}}$ centered at
$\mathbf{r}_{\text{D}}$) and an acceptor sphere (radius $R_{\text{A}}$
centered at $\mathbf{r}_{\text{A}}$). The electron transfer is assumed to
occur instantaneously—before the orientational degrees of freedom of the
solvent molecules can respond—hence the change in $\mathcal{G}_{\text{elst}}$
involves $\varepsilon_{\infty}$ in addition to $\varepsilon_{\text{s}}$. The
prefactor of $(\varepsilon_{\infty}^{-1}-\varepsilon_{\text{s}}^{-1})$ in Eq.
(5.5) is sometimes called the Pekar factor,424, 386 and replaces the prefactor
of $(1-\varepsilon_{\text{s}}^{-1}$) in the equilibrium version of the Born
model. For water, $\varepsilon_{\text{s}}$ is large due to the sizable H2O
dipole moment, leading to a large orientational component to the dielectric
screening effect, but for non-polar solvents most or all of the solvent
response is electronic since
$\varepsilon_{\infty}\approx\varepsilon_{\text{s}}$ (see Table 5). The
electronic contribution comes from the intrinsic polarizability of the solvent
molecules, and for this reason $\varepsilon_{\infty}$ has sometimes been
called the “dielectric constant for induced polarization”.367

###  5.2 State-specific approach

The phenomenology introduced above can be generalized to a rigorous
description of electrostatics beyond the Born or models,425, 426 affording a
continuum theory of nonequilibrium solvation.425, 426, 424, 385, 427, 422, 428
Several variations and implementations have been formulated for PCMs,429, 430,
431, 432, 185, 433 as well as for continuum solvation based on Poisson’s
equation.434, 39 Operationally, a charge density $\rho_{0}(\mathbf{r})$
corresponding to the initial electronic state of the solute is first
equilibrated with a continuum whose dielectric constant is
$\varepsilon_{\text{s}}$. The, following excitation (or ionization) from state
$|0\rangle$ to state $|k\rangle$, the difference density
$\Delta\rho_{k}(\mathbf{r})=\rho_{k}(\mathbf{r})-\rho_{0}(\mathbf{r})$ is
allowed to polarize a continuum whose dielectric constant is
$\varepsilon_{\infty}$. We now consider this in detail.

Consistent with the appearance of the zero- and infinite-frequency dielectric
constants in Eqs. (5.4) and (5.5), a general nonequilibrium theory of
continuum electrostatics is based upon a partition
$\mathbf{P}=\mathbf{P}_{\text{fast}}+\mathbf{P}_{\text{slow}}$ in which the
fast component of the polarization remains in equilibrium with the solute even
when the latter experiences a sudden change in its charge density. The slow
polarization is frozen on this timescale, fixed at the value determined by
equilibrium solvation of the initial state. This is accomplished via a
partition of the (linear) electric susceptibility,
$\chi(\omega)=[\varepsilon(\omega)-1]/4\pi$. Separating the susceptibility
$\chi=\chi_{\text{f}}+\chi_{\text{s}}$ into fast and slow contributions, the
form of $\chi_{\text{f}}$ is suggested by Eq. (5.1), namely410, 426, 435, 436,
185

$\displaystyle\chi_{\text{f}}$ $\displaystyle=(\varepsilon_{\infty}-1)/4\pi$
(5.6a) $\displaystyle\chi_{\text{s}}$
$\displaystyle=(\varepsilon_{\text{s}}-\varepsilon_{\infty})/4\pi\;,$ (5.6b)

This has been called the “Marcus partition” of the of the polarization
response,185 although it was actually formalized by Brady and Carr.410 It
embodies the phenomenological concepts introduced above, e.g., that
$\mathbf{P}_{\text{fast}}$ originates with molecular polarizability while
$\mathbf{P}_{\text{slow}}$ is orientational.

An alternative to Eq. (5.6) is the so-called “Pekar partition”,436, 185 which
was originally introduced to describe self-trapped polarons,437, 438 then
subsequently adapted for optical spectroscopy,439, 404, 400 and still later
adopted for use in ASC-type continuum solvation models.5, 440, 441, 442, 13
Its use is prevalent in older literature so the distinction with respect to
Eq. (5.6) is worth pointing out. Within the Pekar approach, the induced
surface charge is partitioned into “inertial” and “dynamic” components
($\mathbf{P}=\mathbf{P}_{\text{in}}+\mathbf{P}_{\text{dyn}}$). The total
surface charge $\sigma(\mathbf{s})$ is given by Eq. (2.24) and the dynamical
part by the analogous expression5

$\sigma_{\text{dyn}}(\mathbf{s})=\frac{1}{4\pi}\left(\frac{\varepsilon_{\infty}-1}{\varepsilon_{\infty}}\right)\left(\frac{\partial\varphi}{\partial\mathbf{n}_{\mathbf{s}}}\right)_{\mathbf{s}=\mathbf{s}^{-}}$
(5.7)

that is obtained with $\varepsilon_{\infty}$ replacing
$\varepsilon_{\text{s}}$. The inertial charge represents the difference,
$\sigma_{\text{in}}(\mathbf{s})=\sigma(\mathbf{s})-\sigma_{\text{dyn}}(\mathbf{s})$.

The difference between these two partitions can readily be understood using
the reaction-field for a ground-state dipole $\bm{\mu}_{0}$ in a spherical
cavity. The slow contribution to the reaction field is436

$\mathbf{E}_{\text{rxn}}^{\text{s}}=\begin{cases}\displaystyle
g_{1}(\varepsilon_{\text{s}},\bar{R})\left(\frac{\varepsilon_{\text{s}}-\varepsilon_{\infty}}{\varepsilon_{\text{s}}-1}\right)\bm{\mu}_{0}&\text{(Marcus-
Brady-Carr)}\\\
\displaystyle\big{[}g_{1}(\varepsilon_{\text{s}},\bar{R})-g_{1}(\varepsilon_{\infty},\bar{R})\big{]}\bm{\mu}_{0}&\text{(Pekar)}\end{cases}\;,$
(5.8)

where $g_{1}(\varepsilon,\bar{R})$ is the “Onsager factor” defined in Eq.
(2.19). The Marcus-Brady-Carr result follows from the fact that the slow
polarization constitutes a fraction
$\chi_{\text{s}}/\chi=(\varepsilon_{\text{s}}-\varepsilon_{\infty})/(\varepsilon_{\text{s}}-1)$
of the total polarization, according to Eq. (5.6), whereas the Pekar result is
set by fiat. Both partitions afford the same total reaction field,436
$\mathbf{E}_{\text{rxn}}^{\text{s}}+\mathbf{E}_{\text{rxn}}^{\text{f}}$, and
therefore the same nonequilibrium free energy,436, 443, 185 up to some minor
issues involving discretization along the lines of what was discussed in
Section 3.1.185 However, the partition into fast and slow components is
different. Noting that
$g_{1}(\varepsilon_{\text{s}},\bar{R})\approx\bar{R}^{-3}$ for high-dielectric
solvents, Brady and Carr410 noted that the Pekar result for
$\mathbf{E}_{\text{rxn}}^{\text{s}}$ seems oddly small (and also decoupled
from the value of $\varepsilon_{\text{s}}$) in this limit. Indeed, for water
one obtains $\mathbf{E}_{\text{rxn}}^{\text{s}}=0.97\bm{\mu}_{0}/\bar{R}^{3}$
for the Marcus-Brady-Carr case and
$\mathbf{E}_{\text{rxn}}^{\text{s}}=0.63\bm{\mu}_{0}/\bar{R}^{3}$ for the
Pekar partition. For that reason, the Marcus-Brady-Carr partition is the more
common choice in modern literature although the Pekar partition has not
vanished.443 (The Marcus partition has also been criticized recently, and some
alternatives have been suggested.386) Side-by-side expressions for the free
energy $\mathcal{G}_{\text{elst}}$ in the Marcus-Brady-Carr and the Pekar
partitions are provided by Tomasi et al..13 As those authors note, there is
some confusion in the literature regarding the names, e.g., Eq. (5.6) is
called Pekar partition by Chipman.431 As such, Tomasi et al. designate these
schemes as “partition I”, meaning Eq. (5.6), and “partition II”, meaning Eq.
(5.7). Although this notation has been adopted in a few places,427, 443 the
names “Marcus” (for partition I) and “Pekar” (for partition II) remain common.

Having settled on the partition given in Eq. (5.6), the basic idea of
nonequilibrium polarization is that the susceptibility $\chi_{\text{s}}$
should be used to induce polarization for the initial state (“0”), whose
solute charge density is $\rho_{0}(\mathbf{r})$, and then $\chi_{\text{f}}$
should be used in conjunction with the difference density
$\Delta\rho(\mathbf{r})$ in order to adjust the polarization to reflect the
final state. To realize this in practice, one first computes the surface
charge $\sigma_{0}(\mathbf{s})$ that is induced by $\rho_{0}(\mathbf{r})$ in a
medium whose dielectric constant is $\varepsilon_{\text{s}}$, according to a
normal equilibrium solvation calculation. Next, $\sigma_{0}(\mathbf{s})$ is
partitioned into fast and slow contributions,426, 385, 432

$\displaystyle\sigma^{\text{f}}_{0}(\mathbf{s})$
$\displaystyle=\left(\frac{\varepsilon_{\infty}-1}{\varepsilon_{\text{s}}-1}\right)\sigma_{0}(\mathbf{s})$
(5.9a) $\displaystyle\sigma^{\text{s}}_{0}(\mathbf{s})$
$\displaystyle=\left(\frac{\varepsilon_{\text{s}}-\varepsilon_{\infty}}{\varepsilon_{\text{s}}-1}\right)\sigma_{0}(\mathbf{s})\;.$
(5.9b)

The quantity $\sigma^{\text{s}}_{0}(\mathbf{s})$ is retained, whereas
$\sigma^{\text{f}}_{0}(\mathbf{s})$ is replaced by a surface charge induced by
the excited-state charge distribution, in a medium whose dielectric constant
is $\varepsilon_{\infty}$.

In order to derive rigorous formulas for the electrostatic free energy of an
excited state, introduce a Schrödinger equation of the form

$\underbrace{\left(\hat{\mathcal{H}}_{\text{vac}}+\hat{\mathcal{R}}_{0}^{\text{s}}+\hat{\mathcal{R}}_{k}^{\text{f}}\right)}_{\hat{\mathcal{H}}_{k}^{\text{SS}}}\big{|}\Psi_{k}\big{\rangle}=\mathcal{E}^{\text{SS}}_{k}\,\big{|}\Psi_{k}\big{\rangle}$
(5.10)

with $k=0$ for the ground state. The quantity $\hat{\mathcal{H}}_{\text{vac}}$
is the vacuum Hamiltonian and the reaction-field contribution
$\hat{\mathcal{R}}_{k}=\hat{\mathcal{R}}_{0}^{\text{s}}+\hat{\mathcal{R}}_{k}^{\text{f}}$
consists of a slow component $\hat{\mathcal{R}}_{0}^{\text{s}}$ that
originates with the ground-state density $\rho_{0}$ and susceptibility
$\chi_{\text{s}}$, along with a fast component
$\hat{\mathcal{R}}_{k}^{\text{f}}$ based on the final-state density
$\rho_{k}(\mathbf{r})$ and susceptibility $\chi_{\text{f}}$. Because
$\hat{\mathcal{R}}_{k}^{\text{f}}$ depends on the wave function
$|\Psi_{k}\rangle$ that is used to compute the excited state’s electrostatic
potential, the Hamiltonian
$\hat{\mathcal{H}}_{k}^{\text{SS}}=\hat{\mathcal{H}}_{\text{vac}}+\hat{\mathcal{R}}_{k}$
that is introduced in Eq. (5.10) is “state-specific” (SS), and straightforward
attempts to solve this equation encounter significant complications including
convergence difficulties and ambiguous formulas for transition moments.444
These problems can be circumvented by treating
$\hat{\mathcal{R}}_{k}^{\text{f}}$ as a perturbation on top of zeroth-order
states that are eigenfunctions of $\hat{\mathcal{H}}_{0}^{\text{SS}}$, as
discussed below.

First, let us consider an expression for the free energy in an excited state.
Note that Eq. (2.8) for the ground-state free energy $\mathcal{G}_{0}$ is
implicitly based on a Hamiltonian

$\hat{\mathcal{H}}_{0}\equiv\hat{\mathcal{H}}_{0}^{\text{SS}}=\hat{\mathcal{H}}_{\text{vac}}+\hat{\mathcal{R}}_{0}^{\text{s}+\text{f}}\;,$
(5.11)

however $\mathcal{G}_{0}$ differs from
$\langle\Psi_{0}|\hat{\mathcal{H}}_{0}|\Psi_{0}\rangle$ by an amount equal to
the work
$\mathcal{W}_{0}=\tfrac{1}{2}\langle\Psi_{0}|\hat{\mathcal{R}}_{0}|\Psi_{0}\rangle$
that is required to polarize the continuum. This expression for the work can
be generalized to an arbitrary state $|\Psi_{k}\rangle$:

$\mathcal{W}_{k}=\frac{1}{2}\big{\langle}\Psi_{k}\big{|}\hat{\mathcal{R}}_{k}\big{|}\Psi_{k}\big{\rangle}=\frac{1}{2}\int_{\bm{\Gamma}}\sigma_{k}(\mathbf{s})\;\varphi^{\rho_{k}}\\!(\mathbf{s})\;.$
(5.12)

Superscripts “s” or “f” will be added to $\sigma_{k}(\mathbf{s})$, and thus to
$\hat{\mathcal{R}}_{k}$ and $\mathcal{W}_{k}$, signifying the partition into
slow or fast charge according to Eq. (5.9). With this notation, the excited-
state generalization of $\mathcal{G}_{0}$ is185

$\mathcal{G}_{k}^{\text{SS}}=\mathcal{E}_{k}^{\text{SS}}-\mathcal{W}_{0}^{\text{s}}-\mathcal{W}_{k}^{\text{f}}+\mathcal{W}_{0,k}\;.$
(5.13)

where

$\mathcal{E}_{k}^{\text{SS}}=\big{\langle}\Psi_{k}\big{|}\hat{\mathcal{H}}_{k}^{\text{SS}}\big{|}\Psi_{k}\big{\rangle}=\big{\langle}\Psi_{k}\big{|}\hat{\mathcal{H}}_{\text{vac}}+\hat{\mathcal{R}}_{0}^{\text{s}}+\hat{\mathcal{R}}_{k}^{\text{f}}\big{|}\Psi_{k}\big{\rangle}$
(5.14)

and

$\mathcal{W}_{0,k}=\frac{1}{2}\int_{\bm{\Gamma}}\varphi^{\sigma_{0}^{\text{s}}}(\mathbf{s})\left[\sigma_{k}^{\text{f}}(\mathbf{s})-\sigma_{0}^{\text{f}}(\mathbf{s})\right]d\mathbf{s}\;.$
(5.15)

Equation (5.13) has a straightforward interpretation. To obtain the free
energy $\mathcal{G}_{k}^{\text{SS}}$ for state $k$, which includes the effects
of averaging over implicit solvent degrees of freedom, the energy
$\mathcal{E}_{k}^{\text{SS}}$ that is obtained from the Schrödinger equation
must be reduced by the work
$\mathcal{W}_{0}^{\text{s}}+\mathcal{W}_{k}^{\text{f}}$ that is required for
the ground- and excited-state polarization processes. The final term,
$\mathcal{W}_{0,k}$, accounts for the Coulomb interaction between initial- and
final-state surface charge. This term has sometimes been omitted from similar
treatments that start from a nonequilibrium free energy expression that is
otherwise analogous to Eq. (5.13),376, 445 however its presence is necessary
when the Marcus partition of the polarization is used.436, 426, 435, 13, 443,
185 The nonequilibrium expression for excitation energies is
$\hbar\omega_{k}=\mathcal{G}_{k}^{\text{SS}}-\mathcal{G}_{0}$, which is185

$\mathcal{G}_{k}^{\text{SS}}-\mathcal{G}_{0}=\Delta\mathcal{E}_{k}^{\text{SS}}-\mathcal{W}_{k}^{\text{f}}+\mathcal{W}_{0}^{\text{f}}+\mathcal{W}_{0,k}\;,$
(5.16)

where
$\Delta\mathcal{E}_{k}^{\text{SS}}=\mathcal{E}_{k}^{\text{SS}}-\mathcal{E}_{0}$.
Equation (5.16) also has a straightforward interpretation. The quantity
$\Delta\mathcal{E}_{k}^{\text{SS}}$ is the difference between ground- and
excited-state eigenvalues of the SS Hamiltonian [Eq. (5.10)], but must be
corrected by the difference in the work required to polarize either state,
which is restricted to the fast part of the response
($\mathcal{W}_{k}^{\text{f}}-\mathcal{W}_{0}^{\text{f}}$) since that is all
that changes upon vertical excitation.

As a simple example of the nonequilibrium formalism, we consider calculations
of vertical ionization energies (VIEs) in aqueous solution. These can be
measured using liquid microjet photoelectron spectroscopy,446, 447, 448 and
may be quite different from gas-phase values,446 especially for singly-charged
anions X-(aq) where the initial state is solvated much more strongly than the
final state, and while equilibrium solvation models might be appropriate for
computing adiabatic ionization energies, in which solvent has the opportunity
to re-equilibrate following ionization, these models do a poor job of
predicting VIEs.449, 446 From a computational perspective, the change in
charge upon photoionization means that long-range polarization effects are
significant, requiring hundreds of explicit solvent molecules (with
concomitant sampling) to obtain a converged result.450, 451, 452, 453, 454
Convergence is significantly accelerated by continuum boundary conditions,
using an atomistic solute X-(H${}_{\text{2}}$O)n that contains one or two
solvation shells of explicit water molecules.434, 39 The limited size of the
atomistic region not only reduces the cost of the quantum chemistry
calculation for any one structure, but also reduces the sampling that is
required in order to obtain converged averages. (A similar cluster + continuum
strategy is often used in p$K_{\text{a}}$ calculations, for essentially the
same reasons.255, 256, 257, 8) PCM boundary conditions have also been shown to
accelerate convergence of absorption spectra with respect to inclusion of
explicit water molecules,455 although the effects (in absolute energy shifts)
are not as dramatic as they are for ionization.

Table 6:  Vertical ionization energies (VIEs) for aqueous ions, comparing
experimental results to calculations using nonequilibrium MP2 + PCM
calculations.456 Each calculation contains $\approx 30$ explicit water
molecules and each calculated VIE represents an average over $\approx 100$
snapshots from a simulation.

Solute | Exptl. VIE | Computed VIE (eV)
---|---|---
| (eV) | Noneq. PCM |  | Equil. PCM | No PCM
|  | sphericala | SASb |  | sphericala | SASb |
Li+ | $60.4^{c}$ | 61.8 | 61.6 |  | 61.3 | 61.0 | 64.2
Na+ | $35.4^{c}$ | 36.5 | 36.3 |  | 36.0 | 35.8 | 38.9
H2O | $11.7^{d}$ | 11.6 | 11.6 |  | 11.1 | 10.9 | 13.8
$e^{-}$ | $3.7^{e}$ | 3.2 | 3.2 |  | 2.6 | 2.6 | 1.8
F- | $11.6^{f}$ | 11.4 | 11.5 |  | 10.8 | 10.9 | 10.0
Cl- | $9.6^{c}$ | 9.4 | 9.4 |  | 8.8 | 8.8 | 7.9

  * •

aSingle spherical cavity for the entire atomistic region, $R=7.525$ Å. bEq.
(3.7) with $\alpha=1.0$, $R_{\text{probe}}=1.4$ Å, and atomic radii
$\\{R_{\text{vdW},A}\\}$ from Ref. 115. cRef. 449. dRef. 457. eRef. 458. fRef.
448.

Table 6 shows aqueous VIEs computed for several small solutes using a cluster
+ continuum approach with an atomistic region extending to a radius of 5.5 Å
around the solute, containing $\approx 30$ explicit water molecules. Shown
side-by-side in Table 6 are VIEs computed using “gas phase” (vacuum) boundary
conditions that include the explicit water molecules but lack any continuum
model, versus results using equilibrium and nonequilibrium PCMs. In the latter
case, results are shown both using a SAS cavity but also using a cavity
consisting of a single sphere around atomistic region, which affords VIEs that
are essentially identical to the SAS values. The nonequilibrium PCM results
are $\approx 1$ eV too large for Li+(aq) and Na+(aq) but significantly more
accurate for the aqueous halide ions. For the halides, these calculations are
also significantly more accurate than previous QM/MM or equilibrium PCM
calculations.449 For neat liquid water, these calculations represent the most
accurate VIE to date, in line with new experiments,457 and are also one of the
most accurate VIE calculations to date for the challenging $e^{-}$(aq)
system.459

More germane to the present discussion is the comparison of various boundary
conditions. The “gas phase” results, meaning X-(H${}_{\text{2}}$O)30 and
M+(H${}_{\text{2}}$O)30, are in serious error relative to experiment, with
VIEs that are too small for the anions and too large for the cations. This is
consistent in both cases with understabilization of the state with larger
charge. Addition of equilibrium PCM boundary conditions modifies VIEs for
these systems by up to 3 eV for the cations, which is perhaps unsurprising
given that the divalent ion M2+ incurs very long-range polarization effects in
the final state, but more surprisingly even neutral H${}_{\text{2}}$O as the
solute sees a shift of 3 eV when continuum boundary conditions are activated.
While the application of an equilibrium PCM pushes the VIE substantially in
the right direction with respect to experiment, results remain far from
quantitative until the nonequilibrium correction is added. The latter ranges
in magnitude up to $\approx 0.6$ eV. Other calculations for aqueous
nucleobases find nonequilibrium corrections to VIEs of $\approx 1$ eV,460 and
the latter results suggest that for those particular systems at least, the
solvent response upon ionization is the most important difference between
vertical and adiabatic ionization energies, more so than geometric relaxation
of the ionized solute.460

The nonequilibrium formalism is relatively straightforward for ionization,
assuming that the final state is the ground state of the ionized species, but
is more complicated for excited states. For one, the state-specific nature of
the Hamiltonian in Eq. (5.10) may cause convergence problems in the presence
of (near-)degeneracies.444, 461 Even when the states are well-separated,
properties such as oscillator strengths are ill-defined because the various
final-state wave functions are not eigenfunctions of a common Hamiltonian, and
therefore are not mutually orthogonal.444 These problems are not unique to
continuum solvation methods and arise for any kind of polarizable model of the
environment, e.g., for QM/MM methods that employ polarizable force fields.462,
444 A solution to this conundrum is to treat
$\hat{\mathcal{R}}_{k}^{\text{f}}$ in Eq. (5.10) as a perturbation.463, 376,
445, 185, 432, 433 To do so, first introduce a set of zeroth-order states

$\hat{\mathcal{H}}_{0}\left|\Psi_{k}^{(0)}\right\rangle=\mathcal{E}_{k}^{(0)}\left|\Psi_{k}^{(0)}\right\rangle\;,$
(5.17)

such that the eigenvalue $\mathcal{E}_{k}^{(0)}$ includes the effects of the
ground-state reaction field, $\hat{\mathcal{R}}_{0}$. Eigenfunctions in Eq.
(5.17) are orthonormal, thus a convenient approximation for
$\mathcal{G}_{k}^{\text{SS}}$ in Eq. (5.13) uses $\mathcal{E}_{k}^{(0)}$ in
place of $\mathcal{E}_{k}^{\text{SS}}$ and furthermore uses
$\big{|}\Psi_{k}^{(0)}\big{\rangle}$ to evaluate the electrostatic potential
for state $k$. This avoids the complexities of the SS approach, and is
equivalent to first-order perturbation theory with respect to a perturbation
$\hat{W}=\hat{\mathcal{R}}_{k}^{\text{f}}-\hat{\mathcal{R}}_{0}^{\text{f}}$,
obtained from a partition

$\hat{\mathcal{H}}_{k}^{\text{SS}}=\underbrace{\hat{\mathcal{H}}_{\text{vac}}+\hat{\mathcal{R}}_{0}^{\text{s}+\text{f}}}_{\hat{\mathcal{H}}_{0}}+\underbrace{\hat{\mathcal{R}}_{k}^{\text{f}}-\hat{\mathcal{R}}_{0}^{\text{f}}}_{\hat{W}}$
(5.18)

of the SS Hamiltonian. This has been called the perturbation theory state-
specific (ptSS) approach to nonequilibrium solvation.185, 432, 433 (In
principle, it could be extended to higher-order perturbation theory but it is
not clear that is warranted.) Note that a widely-used “corrected linear
response” (cLR) procedure,376 introduced for excited-state PCM calculations at
the level of time-dependent (TD-)DFT is really a ptSS approach; the “linear
response” refers to TDDFT, not to the linear-response PCM formalism that is
discussed in Section 5.3. In the context of the ptSS or cLR approach, it is
best to view TDDFT as a form of configuration interaction with single
substitutions (CIS), which provides an eigenvalue equation of the form in Eq.
(5.17). This is consistent with the idea that the SS version of TDDFT, as
implemented by Improta et al.,464, 465 is a fully iterative version of Eq.
(5.10) within a CIS-style ansatz. Both the ptSS and the fully SS approaches to
TDDFT do require construction of the “relaxed” density466, 467, 468 for the
TDDFT excited state in question, in order to compute its electrostatic
potential.

To obtain practical formulas, e.g., for ASC-PCMs, let us introduce a vector
notation for surface integrals. As an example, we rewrite Eq. (2.36) for
$\mathcal{G}_{\text{elst}}$ in the form

$\mathcal{G}_{\text{elst}}=\frac{1}{2}\int_{\bm{\Gamma}}\sigma(\mathbf{s})\;\varphi^{\rho}(\mathbf{s})\;d\mathbf{s}=\tfrac{1}{2}\mathbf{q}\bm{\cdot}\mathbf{v}^{\rho}\;.$
(5.19)

The quantities $\mathbf{q}$ and $\mathbf{v}^{\rho}$ were introduced in Eq.
(3.2) and the dot-product notation represents how ASC-PCM surface integrals
are evaluated in practice, upon discretization of the cavity surface. Using
this notation, an expression for the nonequilibrium free energy for excited
state $k$ can be manipulated into the form376

$\mathcal{G}_{k}^{\text{neq}}=\mathcal{E}_{k}^{(0)}+\tfrac{1}{2}\mathbf{v}_{0}\bm{\cdot}\mathbf{q}_{0}+\tfrac{1}{2}(\mathbf{v}_{k}-\mathbf{v}_{0})\bm{\cdot}(\Delta\mathbf{q}_{k}^{\text{f}})+\mathcal{W}_{0,k}\;,$
(5.20)

where $\mathbf{v}_{0}$ and $\mathbf{v}_{k}$ represent the electrostatic
potentials for states $\big{|}\Psi_{0}^{(0)}\big{\rangle}$ and
$\big{|}\Psi_{k}^{(0)}\big{\rangle}$, and
$\Delta\mathbf{q}^{\text{f}}_{k}=\mathbf{q}_{k}^{\text{f}}-\mathbf{q}_{0}^{\text{f}}$
is the difference in the fast polarization charges in the two states. The
latter is computed according to

$\Delta\mathbf{q}^{\text{f}}_{k}=\mathbf{Q}_{\varepsilon_{\infty}}(\mathbf{v}_{k}-\mathbf{v}_{0})=\mathbf{Q}_{\varepsilon_{\infty}}\mathbf{v}^{\Delta\rho}$
(5.21)

for a reaction-field operator
($\mathbf{Q}_{\varepsilon}=\mathbf{K}_{\varepsilon}^{-1}\mathbf{Y}_{\varepsilon}$)
involving the optical dielectric constant. The second equality in Eq. (5.21)
recognizes that $\mathbf{v}^{\Delta\rho}=\mathbf{v}_{k}-\mathbf{v}_{0}$ is the
electrostatic potential corresponding to the difference density
$\Delta\rho_{k}(\mathbf{r})=\rho_{k}(\mathbf{r})-\rho_{0}(\mathbf{r})$.
Somewhat similar expressions to Eq. (5.20) can be found elsewhere,430, 464,
465 but the connection to the actual free energy of the excited state is most
explicit in the work of Caricato et al..376 For cases where the solvent
polarization has had time to fully equilibrate with respect to the excited-
state density, an analogous expression for the equilibrium free energy
$\mathcal{G}_{k}^{\text{eq}}$ is obtained from Eq. (5.20) by replacing
$\Delta\mathbf{q}_{k}^{\text{f}}$ with $\mathbf{q}_{k}-\mathbf{q}_{0}$, where
both ground- and excited-state charges are equilibrium ones, so that in the
equilibrium case
$\mathbf{q}_{k}=\mathbf{Q}_{\varepsilon_{\text{s}}}\mathbf{v}_{k}$.376

The solvent-corrected excitation energy is simply the difference between
ground- and excited-state free energies,376, 469

$\hbar\omega_{0k}^{\text{neq}}=\mathcal{G}_{k}^{\text{neq}}-\mathcal{G}_{0}=\Delta\mathcal{E}_{k}^{(0)}+\tfrac{1}{2}(\mathbf{v}_{k}-\mathbf{v}_{0})\bm{\cdot}(\Delta\mathbf{q}^{\text{f}}_{k})+\mathcal{W}_{0,k}\;.$
(5.22)

The quantity
$\Delta\mathcal{E}_{k}^{(0)}=\mathcal{E}_{k}^{(0)}-\mathcal{E}_{0}$ is the
excitation energy computed in the presence of the ground-state reaction field,
which is then corrected in Eq. (5.22) for the change in the fast polarization.
These equations remain valid for the Pekar partition if $\mathcal{W}_{0,k}$ is
omitted from Eqs. (5.20) and (5.22). As noted by Cammi et al.,463 Eq. (5.22)
is the detailed analogue of the heuristic theories of excited-state solvation
developed much earlier, e.g., by McRae,439, 404, 470 Lippert,400 and
Mataga.399, 401 This becomes clear upon considering the special case of a
polarizable dipole in a spherical cavity.463, 464

Although presented here in the notation of ASC-PCMs, an analogous theory of
nonequilibrium solvation can be developed based directly on Poisson’s
equation.434, 39 In that context, the total charge density
$\rho_{\text{tot}}(\mathbf{r})=\rho(\mathbf{r})+\rho_{\text{pol}}(\mathbf{r})$
includes an induced polarization $\rho_{\text{pol}}(\mathbf{r})$ [as in Eq.
(2.6)], in addition to the charge density $\rho(\mathbf{r})$ due to the
atomistic solute. The reaction-field potential is the electrostatic potential
arising from $\rho_{\text{pol}}(\mathbf{r})$ and surface integrals such as the
ones in Eqs. (5.12) and (5.15) are replaced by volumetric integration.
Conveniently, the dot product notation introduced in Eq. (5.19) is ambivalent
to this distinction and a formula analogous to Eq. (5.20) can be derived,39
where the dot product signifies volumetric integration.

In the interest of brevity, the notation introduced above omits a subscript
“elst” on both $\mathcal{G}_{k}^{\text{SS}}$ [Eq. (5.13)] and
$\mathcal{G}_{k}^{\text{neq}}$ [Eq. (5.20)], and for that matter on
$\mathcal{G}_{0}$ in Eq. (2.8) as well. In each case, these quantities
represent only the electrostatic contribution to the free energy. Some of the
earliest theoretical work on solvatochromatic shifts was concerned not just
with changes in the chromophore’s dipole moment, as in the Onsager-style
treatment leading to the Lippert-Mataga equation, but also with the role of
dispersion effects.471, 472, 473, 470 In modern formulations of continuum
theory, however, there have been only preliminary efforts to incorporate
nonelectrostatic interactions (as described in Section 4) into excited-state
calculations.316, 474, 475, 317 This is an interesting problem insofar as
solute–solvent dispersion is likely more attractive in an excited state, as
the excited-state wave function is probably more polarizable than the ground
state, but at the same time Pauli repulsion probably increases in the excited
state due to its larger size. To a significant extent, the favorable results
for solvatochromic shifts that are obtained with electrostatic-only models432,
443, 476, 461 likely rely on some error cancellation along these lines. Models
that introduce state-specific nonelectrostatic interactions also do well for
solvatochromic shifts.475

There is one remaining source of complexity when the nonequilibrium theory is
applied to excited states, in that the Schrödinger equation in Eq. (5.10)
leaves open the question of what level of self-consistency should be sought in
obtaining the excited-state reaction-field operator,
$\hat{\mathcal{R}}_{k}^{\text{f}}$, which depends on $|\Psi_{k}\rangle$. The
multiple-choice answer to this question leads to several categories of methods
that are mapped out in Fig. 10, and which are called “perturbation to energy”
(PTE), “perturbation to density” (PTD), and “perturbation with self-consistent
energy and density” (PTED). This nomenclature derives from efforts to use
perturbation theory to include correlation in the ground-state calculation
(e.g., MP2 + SCRF),477, 478, 479, 480 and the question of whether (and how)
electron correlation should be included in the density that is used to
polarize the continuum. The same notation has been adopted for PCM
calculations using non-perturbative models such as coupled-cluster theory,481,
482 and is used here in a discussion that is formulated specifically with
excited-state calculations in mind. As illustrated in Fig. 10(a), the PTE
scheme involves self-consistent solution of the SCF + SCRF problem followed by
a single-shot post-SCF calculation using solvent-polarized molecular orbitals
(MOs). This represents a kind of “zeroth-order” inclusion of solvation effects
in the correlated calculation,432 and has obvious advantages in terms of cost:
assuming that the post-SCF step dominates the cost of the gas-phase
calculation, then then addition of SCRF boundary conditions adds little to the
overall cost of the calculation. Alternatively, in the PTD scheme the
correlated calculation is performed in the gas phase and then the correlated
density (rather than the SCF density) is used to polarize the solvent. This
introduces solvation effects beyond zeroth order in perturbation theory,479 at
marginally increased cost: it is still a single-shot correlation calculation,
however the relaxed density is required, which entails computational effort
along the lines of a gradient calculation at the correlated level of theory. A
slightly better-performing variant of the traditional PTD approach is the PTE-
PTD scheme,433 in which the SCF + SCRF calculation is solved self-consistently
and those MOs are used in the post-SCF calculation, but then the correlated
density is used in a final, single-shot PCM calculation to compute the
solvation energy $\mathcal{G}_{\text{elst}}$. None of the aforementioned
approaches constitutes a fully self-consistent treatment of post-SCF
correlation effects, but this can be accomplished using the PTED scheme that
is mapped out in Fig. 10(a). Here, the correlated density is used to obtain
the PCM surface charge and this procedure is iterated to self-consistency.
This approach is significantly more expensive because the correlated
calculation is performed at each SCRF iteration.

Fig. 10:  Flowcharts of various state-specific pathways for combining a PCM or
other SCRF procedure with a quantum chemistry method requiring a post-SCF
step. (a) Illustration of the PTE and PTD schemes and two combinations
thereof. Forward-backward arrows ($\rightleftarrows$) indicate where solute
density ($\rho$) and polarization charge ($\sigma$) are iterated to self-
consistency, and downward arrows indicate the point at which various
contributions to the energy are extracted. (b) Schematic representation of the
PTE and PTED procedures for an excited-state (ES) calculation, along with the
PTES procedure designed as a lower-cost approximation to PTED. Panel (a) is
adapted from Ref. 433; copyright 2017 The PCCP Owner Societies. Panel (b) is
adapted from Ref. 428; copyright 2019 John Wiley & Sons.

These ideas have since been extended beyond their perturbation theory origins
and represent the various possible levels of self-consistency in any
calculation that combines an SCRF approach with a quantum chemistry method
requiring a post-SCF calculation,483, 484, 485, 432, 433, 428, 461 including
TDDFT. An alternative pictorial representation of the simplest method (PTE)
and most complete scheme (PTED) is provided in Fig. 10(b), which provides a
flowchart for an excited-state calculation, indicating which densities are
used to construct the various reaction-field operators $\hat{\mathcal{R}}$.
Because of the expense associated with the fully self-consistent PTED
approach, approximations have been developed in which both the ground- and
excited-state calculations are iterative, but those two iterative sequences
are uncoupled to one another.480, 483, 484, 485, 486 This scheme, which
Caricato calls “PTES” and has implemented at the coupled-cluster level of
theory,483, 484, 485 is analogous to a “vertical excitation model” introduced
for TDDFT.474 At the DFT level, the PTED scheme in Fig. 10 is essentially
equivalent to the SS-TDDFT + PCM method introduced by Improta et al..464, 465

Fig. 11:  Solvatochromic shifts in lowest ${}^{1}\pi\pi$ state for derivatives
of nitrobenzene (PhNO2) in different solvents, comparing experimental values
to ADC(2) + C-PCM calculations. Solvent effects are described using PTE and
PTD variants of the ptSS approach. Also shown are results for an empirically-
scaled version of the nonequilibrium PTD correction. Adapted from Ref. 432;
copyright 2015 American Chemical Society.

Figure 11 presents solvatochromic shifts for a set of nitrobenzene
derivatives,432 with excitation energies computed at the level of second-order
algebraic diagrammatic construction [ADC(2)],487 which is something of an
excited-state analogue of MP2. Solvent contributions in Fig. 11 are
incorporated using either the PTE or PTD variant of the ptSS approach.
Differences between the two variants are negligible, and both approaches show
good agreement with experimental shifts, without the need to invoke the more
expensive PTED scheme. For many of these molecules, the first-order ptSS
contribution to the solvatochromic shift (representing fast polarization) is
0.10–0.15 eV, in shifts ranging up to 0.6 eV. The remaining part comes from
the zeroth-order contribution of simply inserting solvent-polarized MOs into
the correlated part of the calculation.432 Tests on a more diverse set of
systems do reveal a small systematic error in the PTE approach,461 but the
mean error with respect to experiment remains $<0.1$ eV and the systematic
error can be eliminated by intermediate approaches that do not require the
full self-consistency of the PTED approach.461 In particular, the PTE-PTD
scheme [see Fig. 10(a)] works well in this regard.433 This approach requires
the correlated density but is not iterative at the correlated level of theory.
Other benchmark studies comparing continuum approaches to large QM
calculations with explicit solvent have suggested that QM/PCM excitation
energies may agree better with full-QM result as compared to QM/MM
calculations, but explicit water molecules in the QM/PCM calculation are
required to obtain good agreement for oscillator strengths.488 Simulation of
band shapes requires thermal sampling, which cannot be accomplished without at
least some explicit solvent molecules.

###  5.3 Linear response approach

Despite its computational complexities, the SS approach to excited-state
solvation is conceptually straightforward. An alternative to the state-by-
state approach, which has fewer moving parts at the computational level, is
based on linear-response (LR) quantum chemistry methods in which excitation
energies are computed from the poles of the frequency-dependent response to a
perturbation, rather than from the Schrödinger equation. Given that the PCM
electrostatic contribution to the energy is

$\mathcal{G}_{\text{elst}}=\frac{1}{2}\int_{\bm{\Gamma}}\int_{\bm{\Gamma}}\varphi(\mathbf{s})\;Q_{\varepsilon}(\mathbf{s},\mathbf{s}^{\prime})\;\varphi(\mathbf{s}^{\prime})\;d\mathbf{s}\;d\mathbf{s}^{\prime}\;,$
(5.23)

where $Q_{\varepsilon}(\mathbf{s},\mathbf{s}^{\prime})$ is the kernel of the
solvent response operator
$\hat{Q}_{\varepsilon}=\hat{K}_{\varepsilon}^{-1}\hat{Y}_{\varepsilon}$, the
solvent model contributes only a one-electron potential,
$\upsilon^{\text{PCM}}(\mathbf{r})=\delta\mathcal{G}_{\text{elst}}/\delta\rho(\mathbf{r})$.
The matrix elements of $\upsilon^{\text{PCM}}(\mathbf{r})$ are489

$\upsilon^{\text{PCM}}_{\mu\nu}=\int_{\bm{\Gamma}}\int_{\bm{\Gamma}}\varphi(\mathbf{s})\;Q_{\varepsilon}(\mathbf{s},\mathbf{s}^{\prime})\;\varphi_{\mu\nu}(\mathbf{s}^{\prime})\;d\mathbf{s}\;d\mathbf{s}^{\prime}=(\mathbf{v}^{\rho})^{\dagger}\mathbf{Q}_{\varepsilon}\mathbf{v}^{\mu\nu}\;,$
(5.24)

where $\varphi_{\mu\nu}(\mathbf{s}^{\prime})$ is the electrostatic potential
generated by the function pair $\mu\nu$ at the point $\mathbf{s}^{\prime}$.
The second equality in Eq. (5.24) demonstrates how the Fock matrix
contribution from $\upsilon^{\text{PCM}}(\mathbf{r})$ is evaluated in
practice, and analogous expressions exist for three-dimensional Poisson
approaches.39 The quantities $\mathbf{v}^{\rho}$ and $\mathbf{v}^{\mu\nu}$
involve only one-electron integrals, so incorporating the PCM contribution
into a LR calculation incurs negligible overhead with respect to the cost of
the gas-phase calculation, and meanwhile this approach is free of the
iterative complexities of the SS method. A general LR-PCM formulation has been
given by Cammi et al.,490, 463 and specific formulations for different excited
state methods are available as well, e.g., for TDDFT and other single-
excitation theories,490, 491, 489, 492, 493 following on earlier
implementations of the coupled-perturbed SCF + PCM procedure for response
properties;180, 494 for multiconfigurational SCF wave functions;88 for
ADC;445, 486 for the $GW$/Bethe-Salpeter equation formalism.495 Finally, it
has been implemented for coupled-cluster theory,496, 497, 498, 499, 500 based
on the coupled-cluster response formalism.482

For isolated-molecule quantum chemistry, the LR formalism for excitation
energies is generally equivalent to solving the corresponding Schrödinger
equation,501 but the LR- and SS-PCM formalisms are not equivalent.463, 502 The
general form of the LR-PCM result is463

$\hbar\omega_{0k}^{\text{neq,LR}}=\hbar\omega_{k}^{(0)}+\underbrace{\big{\langle}\Psi_{k}\big{|}\hat{\mathcal{V}}\big{|}\Psi_{0}\big{\rangle}\big{\langle}\Psi_{0}\big{|}\hat{\mathcal{Q}}^{\text{f}}\big{|}\Psi_{k}\big{\rangle}}_{\Delta\mathcal{R}^{\text{f}}(\bm{\mu}_{0k})}\;,$
(5.25)

where
$\big{\langle}\Psi_{k}\big{|}\hat{\mathcal{V}}\big{|}\Psi_{0}\big{\rangle}$ is
the electrostatic potential generated by the transition density,
$\rho_{k0}(\mathbf{r})$, and
$\big{\langle}\Psi_{0}\big{|}\hat{\mathcal{Q}}^{\text{f}}\big{|}\Psi_{k}\big{\rangle}$
is the apparent surface charge induced by $\rho_{k0}(\mathbf{r})$. For
comparison, the SS-PCM result in Eq. (5.22) can be rewritten in similar
notation:

$\hbar\omega_{0k}^{\text{neq,SS}}=\hbar\omega_{k}^{(0)}+\underbrace{\tfrac{1}{2}\Big{[}\big{\langle}\Psi_{k}\big{|}\hat{\mathcal{V}}\big{|}\Psi_{k}\big{\rangle}-\big{\langle}\Psi_{0}\big{|}\hat{\mathcal{V}}\big{|}\Psi_{0}\big{\rangle}\Big{]}\bm{\cdot}\Big{[}\big{\langle}\Psi_{k}\big{|}\hat{\mathcal{Q}}^{\text{f}}\big{|}\Psi_{k}\big{\rangle}-\big{\langle}\Psi_{0}\big{|}\hat{\mathcal{Q}}^{\text{f}}\big{|}\Psi_{0}\big{\rangle}\Big{]}}_{\Delta\mathcal{R}^{\text{f}}(\Delta\rho_{k})}\;.$
(5.26)

In both cases, the quantity
$\hbar\omega_{0k}^{(0)}\equiv\Delta\mathcal{E}_{k}^{(0)}$ is the zeroth-order
approximation to the solution phase excitation energy, calculated in the
static reaction field of the ground state. The quantity
$\Delta\mathcal{R}^{\text{f}}{}$ represents the change in the dynamical part
of the reaction field, which is a function of the transition dipole moment
$\bm{\mu}_{0k}=\langle\Psi_{0}|\hat{\bm{\mu}}|\Psi_{k}\rangle$ in the LR case
but a function of the difference density $\Delta\rho_{k}(\mathbf{r})$ in the
SS case. A detailed analysis of the two formalisms suggests that their
differences arise from the nonlinear nature of the SS Hamiltonian combined
with the lack of entanglement between the atomistic wave function and its
continuum environment.502

Whatever the origin of the discrepancy, the form of the LR-PCM correction in
Eq. (5.25) is problematic because the correction vanishes for optically-
forbidden transitions, as is readily seen from a model of a dipole in a
spherical cavity, for which
$\Delta\mathcal{R}^{\text{f}}=-g_{1}(\varepsilon_{\infty},\bar{R})\bm{\mu}_{0k}$.464
For the same reason, the LR-PCM correction
$\Delta\mathcal{R}^{\text{f}}(\bm{\mu}_{0k})$ will be rather small for any
excitation involving significant displacement of charge, whereas intuitively
(and in the SS formalism) one expects a significant solvent effect for a
charge-transfer excitation in a polar solvent. Indeed, SS-PCM results are
consistently superior to LR-PCM calculations for excited states with charge-
transfer character.503, 504, 505, 506, 507, 508 (As discussed in Section 5.2,
the cLR formalism encountered in some of these studies is really a ptSS-PCM
approach.) Even for states not dominated by charge transfer, the SS-PCM
approach generally affords smaller errors for solvatochromatic shifts in the
absorption spectrum as compared to LR-PCM calculations,443, 476, 432 although
it is worth bearing in mind that the experimental $\lambda_{\text{max}}$ need
not correspond to the origin (0–0) transition, due to vibrational
structure.509, 510, 511 The ptSS-PCM approach also affords more accurate
results for emission energies,512 although LR-PCM calculations can be used to
optimize the excited-state geometries, which simplifies the procedure. It has
also been argued that the LR correction
$\Delta\mathcal{R}^{\text{f}}(\bm{\mu}_{0k})$ constitutes a solute–continuum
dispersion interaction,502, 513, 428 insofar as it has the form of the solute
charge distribution oscillating at the Bohr frequency $\omega_{0k}^{(0)}$ and
coupling to the dynamical response of the environment. As such, some studies
have opted to include both the LR- and ptSS-PCM corrections to
$\omega_{0k}^{(0)}$.432

##  6 Anisotropic Solvation

Up to this point we have assumed that the continuum environment is isotropic,
which is usually the case for a bulk liquid although there are certain
exceptions (notably, liquid crystals) where polarization of the medium depends
upon the orientation of the electric field vector. This can be described by
allowing $\varepsilon$ to take the form of a $3\times 3$ matrix, with
orientation-dependent permittivities $\varepsilon_{xx}$, $\varepsilon_{yy}$,
and $\varepsilon_{zz}$. The ASC-PCM formalism, and in particular IEF-PCM, has
been formulated to handle a permittivity tensor.155, 156, 157, 514

Fig. 12:  Illustrations of anisotropic permittivity functions
$\varepsilon(\mathbf{r})$ for use in Poisson’s equation. (a) Cluster +
continuum description of chlorate ion at the air/water interface, in which the
atomistic solute is ClO${}_{3}^{-}$(H${}_{\text{2}}$O)30. The background color
shows the function $\varepsilon(\mathbf{r})$, interpolating between
$\varepsilon_{\text{out}}=1$ above the Gibbs dividing surface (GDS) and
$\varepsilon_{\text{out}}=78$ below it, with $\varepsilon_{\text{in}}=1$
inside of the solute cavity. The horizontal line indicates the position of the
dividing surface, $z_{\text{GDS}}=3.5$ Å. (b) Periodic water slab bounded on
either side by continuum water ($\varepsilon=78$, shown in purple), with
regions characterized by $\varepsilon>15$ shown in blue. From left to right,
the interpolating function is modified (using a “filling threshold” parameter)
in order to exclude pockets of high permittivity that encroach into the
interstices between the atomistic water molecules. Panel (b) is adapted from
Ref. 225; copyright 2019 American Chemical Society.

A more general class of anisotropic solvation problems are interfacial
phenomena. An example is shown in Fig. 12(a), in which an atomistic solute
consisting for a ClO${}_{3}^{-}$ ion with approximately two solvation shells
of explicit water molecules is situated at a dielectric representation of the
air/water interface. Here, the atomistic region is subject to a dielectric
environment characterized by $\varepsilon=78$ on one side but $\varepsilon=1$
on the other. The basic PCM formalism is not equipped to handle such a
situation, as it is predicted on a sharp dielectric interface between
$\varepsilon_{\text{in}}=1$ inside the cavity and a bulk solvent value
($\varepsilon_{\text{out}}$) outside, although it can be accomplished
piecewise if the medium is divided into separate domains, each with its own
value of $\varepsilon$.515 Alternatively, to describe the liquid/vapor
interface the PCM equations have been to interpolate matrix elements between
different values of $\varepsilon$.516, 517, 518, 519, 520, 521, 522, 523
Treatment of nonelectrostatic effects proves to be crucial. It has been
suggested that continuum models cannot describe specific-ion (“Hofmeister”)
effects,524, 525 in part because the Born model [Eq. (2.17)] cannot
distinguish between cations and anions, although that particular shortcoming
is better ascribed to an overly-simplified cavity shape.526 Another commonly
held view is that continuum models are incapable of describing the interfacial
affinity exhibited by soft ions,524, 525 because (so the logic goes) the ion
in a continuum solvent ought to be repelled from the interface by its own
image charge.525, 527 Despite this conventional wisdom, continuum models that
include nonelectrostatic interactions have been shown to predict interfacial
free energy minima for dragging a soft ion through the air/water
interface.516, 518, 528, 529

Fundamentally, however, the interfacial solvation problem seems to cry out for
a permittivity function $\varepsilon(\mathbf{r})$ that can take different
values in different regions of space, i.e., a method that solves the
generalized Poisson equation with an anisotropic permittivity function
$\varepsilon(\mathbf{r})$. Such a strategy has been pursued to describe the
interface between a solid-state electrode and an aqueous electrolyte,223, 224,
225, 24, 226, 227, 228, 71, 229, 230 as well as host/guest systems where the
guest experiences a low-dielectric environment despite the fact that the host
is dissolved in water.530 Finally, anisotropic models have been used to
compute VIEs of solutes at the air/water interface,434, 39 in order to connect
with liquid microjet photoelectron spectroscopy.

The setup for an interfacial calculation of this type is illustrated in Fig.
12(a), which depicts an atomistic model of ClO${}_{3}^{-}$(aq) at the
air/water interface and shows how the function $\varepsilon(\mathbf{r})$ is
defined. In this particular example, two solvation shells of explicit water
molecules are included in the atomistic region in order to account for
hydrogen-bonding effects; the continuum model takes care of long-range
polarization upon ionization of ClO${}_{3}^{-}$. VIE calculations based on a
nonequilibrium formulation of Poisson’s equation suggest that the VIEs of
common inorganic anions are very nearly the same at the air/water interface as
they are in bulk water.39 Even for an exotic anion like $e^{-}$(aq),531, 459
it appears that the interfacial VIE that is no more than 0.2–0.4 eV different
from its bulk value.434, 39

In continuum solvation calculations where the atomistic region contains
explicit solvent molecules, one must be careful to parameterize the function
$\varepsilon(\mathbf{r})$ to avoid artificial penetration of high-dielectric
regions into the interstices between molecules. This includes not just methods
that are directly based on Poisson’s equation but also PCMs, where explicit
solvent molecules might be also be added, for the purpose of describing
solute–solvent hydrogen bonds, as in the aforementioned p$K_{\text{a}}$
calculations. Oddly, the dielectric penetration problem in semicontinuum
calculations has received scant attention in the quantum chemistry
literature,39, 225 although there is an analogous problem in classical
Poisson-Boltzmann electrostatics calculations for biomolecules that is widely
discussed. In that context, the problem is that standard cavity construction
algorithms (based on intersecting atomic spheres) may leave pockets of high-
dielectric “solvent” in the hydrophobic interior of a protein.532, 211, 533,
534, 535, 52, 119 This can be mitigated in Poisson-based continuum methods by
appropriate adjustment of the interpolating function that defines the
dielectric boundary. A spatially-varying permittivity function has been
suggested as a solution to the problem that there is no single optimal value
for the dielectric “constant” inside of a protein.52

Figure 12(b) presents an example in which the SCCS model of Section 3.3, which
is based on a functional $\varepsilon[\rho](\mathbf{r})$, is used to remove
artificial high-dielectric regions. In these simulations, the goal is to
perform ab initio molecular dynamics simulations of explicit water but to use
continuum boundary conditions in order to limit the size of the atomistic QM
simulation cell that is required. The parameters that define
$\varepsilon[\rho](\mathbf{r})$ must be chosen carefully, lest high-dielectric
regions penetrate between the explicit water molecules, as they do on the left
side of Fig. 12(b). This situation is physically incorrect because the QM
calculation is based on Coulomb operators that assume vacuum permittivity.
This undesirable situation has been rectified by introducing “solvent
awareness” into the definition of the permittivity function, so that
$\varepsilon(\mathbf{r})$ depends on the coordinates in the atomistic region
directly, and not just implicitly through the functional
$\varepsilon[\rho]$.225 Moving from left to right in Fig. 12(b), this solvent
awareness is activated and removes the spurious high-dielectric regions.

##  7 Closing remarks

With the contents of this review serving as a lengthy introduction to what
continuum solvation models can do, in closing it feels a propos to comment
briefly on their limitations. These are very crude models. That is not
inconsistent with being useful models, but one should not demand too much of
something so simple. It is perhaps best to think of continuum solvation as
improved boundary conditions (as compared to vacuum boundary conditions) for
condensed-phase electronic structure calculations. The primary way in which
most users will encounter the crudeness of these models is in the fact that
there is considerable arbitrariness in construction of the solute cavity, for
which there is no “right” choice, although there are certainly plenty of wrong
ones. In particular, there is no “magical” cavity construction or scaling
factor for the van der Waals radii that will make these quantities universal.
Small tweaks that might provide better answers in one system may very well
degrade the accuracy in other cases. Isodensity cavity constructions, and
smooth interfaces or “soft” cavities, based on permittivity functionals
$\varepsilon[\rho(\mathbf{r})]$, seem like the least arbitrary choices.

With that in mind, differences between how $\mathcal{G}_{\text{elst}}$ is
computed amongst different PCMs seem quite minor in comparison to the overall
quality of these models, to the point where these differences can likely be
parameterized away, or are simply washed out, by minor changes in cavity
construction. The COSMO method, for example, performs well in comparison to
more exact formulations of the continuum electrostatics problem, even in low-
dielectric solvents.170 For the SS(V)PE approach, discretization of
$\hat{D}^{\dagger}\sigma(\mathbf{s})$ proves to be challenging,20, 185 leading
to a wrong conductor limit in some cases,20 nevertheless this model has
successfully been used as the basis of a composite model of solvation.251
Errors in the implementation of certain nonelectrostatic are completely
absorbed in the parameterization of other terms.251

In view of all this, there would seem to be little room to further improve
upon the electrostatic part of continuum solvation models. A corollary is that
efforts to make the solvation model fully consistent with correlated wave
function electronic structure methods, either in the ground or excited states,
seem misguided. The “zeroth-order” model, in which solvent-polarized MOs are
inserted into a post-SCF correlation calculation, recovers the most important
effects, and a ptSS-style correction for the fast polarization response
affords a simple-to-use model for what remains. Differences with respect to a
fully self-consistent model are likely considerably smaller than errors
introduced by continuum representation of the solvent itself. Keeping those
caveats in mind, the continuum solvation approach can be highly effective in
situations were vacuum boundary conditions are dubious, e.g., due to
significant charge rearrangement in a polar solvent (including redox
chemistry), or to modulate the energy levels of the frontier, solvent-exposed
orbitals at play in electronic spectroscopy.

Regarding calculation of solvation energies, which is arguably the primary
application of continuum solvation models in chemistry, the nonelectrostatic
contributions that are needed for favorable agreement with experiment can now
be modeled (in a variety of ways) with an accuracy of $<1$ kcal/mol for
charge-neutral solutes. For ions, the best methods approach the accuracy ($\pm
3$ kcal/mol)247 of the experimental data themselves.248, 251, 216 That said,
these models are trained on relatively small solutes (for which experimental
values of $\Delta_{\text{solv}}\mathcal{G}$ are available); it remains unclear
whether the aforementioned level of statistical accuracy is transferrable to
significantly larger solutes.

Regarding macromolecular solutes, it is clear that PCMs with linear-scaling
solvers ought to be seriously considered as replacements for biomolecular
electrostatics calculations based on finite-difference solution of the
Poisson-Boltzmann equation. The former approaches provide an exact solution to
the classical electrostatics problem,20, 19 up to controllable discretization
errors, and can be formulated in such a way that potential energy surfaces are
inherently continuous and smooth.17, 18 Indeed, it is this author’s opinion
that theorists should not accept as their starting point any approach that
does not intrinsically provide a smooth potential energy surface, as the
finite-difference approach clearly fails to do. The ability to explore the
potential energy surface, and thus to have a well-defined “model
chemistry”,536 is too important to sacrifice.

## Acknowledgements

This work was supported by National Science Foundation, grant nos. CHE-1665322
and CHE-1955282. It is a pleasure to acknowledge current and former group
members who have worked on these projects: Dr. Adrian Lange, Benjamin
Albrecht, Dr. Zhi-Qiang You, Dr. Marc Coons, and Suranjan Paul. Collaborations
with Dr. Jan-Michael Mewes (from the group of Prof. Andreas Dreuw) and Dr.
Christopher Stein (from the group of Prof. Martin Head-Gordon) are also
happily acknowledged.

## Conflict of Interest

J.M.H. serves on the board of directors of Q-Chem Inc.

## References

  * 1 D. Rinaldi and J.-L. Rivail, “Polarisabilités moléculaires et effet diélectrique de milieu à l’état liquide. Étude théorique de la molécule d’eau et de ses dimères”, Theor. Chem. Acc., 32, 57–70 (1973).
  * 2 J.-L. Rivail and D. Rinaldi, “A quantum chemical approach to dielectric solvent effects in molecular liquids”, J. Chem. Phys., 18, 233–242 (1976).
  * 3 S. Miertuš, E. Scrocco, and J. Tomasi, “Electrostatic interaction of a solute with a continuum. A direct utilization of ab initio molecular potentials for the prevision of solvent effects”, Chem. Phys., 55, 117–129 (1981).
  * 4 S. Miertuš and J. Tomasi, “Approximate evaluations of the electrostatic free energy and internal energy changes in solution processes”, Chem. Phys., 65, 239–245 (1982).
  * 5 R. Bonaccorsi, R. Cimiraglia, and J. Tomasi, “Ab initio evaluation of absorption and emission transitions for molecular solutes, including separate consideration of orientational and inductive solvent effects”, J. Comput. Chem., 4, 567–577 (1983).
  * 6 J. Tomasi, “Thirty years of continuum solvation chemistry: A review, and prospects for the near future”, Theor. Chem. Acc., 112, 184–203 (2004).
  * 7 P. Grochowski and J. Trylska, “Continuum molecular electrostatics, salt effects, and counterion binding—A review of the Poisson–Boltzmann theory and its modifications”, Biopolymers, 89, 93–113 (2008).
  * 8 J. R. Pliego Jr. and J. M. Riveros, “Hybrid discrete-continuum solvation methods”, WIREs Comput. Mol. Sci., 10, e1440:1–25 (2020).
  * 9 J. Tomasi and M. Persico, “Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent”, Chem. Rev., 94, 2027–2094 (1994).
  * 10 C. Amovilli, V. Barone, R. Cammi, E. Cancès, M. Cossi, B. Mennucci, C. S. Pomelli, and J. Tomasi, “Recent advances in the description of solvent effects with the polarizable continuum model”, Adv. Quantum Chem., 32, 227–261 (1999).
  * 11 J. Tomasi, B. Mennucci, R. Cammi, and M. Cossi, “Quantum mechanical models for reactions in solution”, in Computational Approaches to Biochemical Reactivity, G. Nâaray-Szabâo and A. Warshel, Eds., Vol. 19 of Understanding Chemical Reactivity; Springer: New York, 2002;  chapter 1, pages 1–102.
  * 12 R. Cammi, B. Mennucci, and J. Tomasi, “Computational modelling of the solvent effects on molecular properties: An overview of the polarizable continuum model (PCM) approach”, in Computational Chemistry: Reviews of Current Trends, J. Leszczynski, Ed., Vol. 8;  World Scientific: Singapore, 2003;  chapter 1, pages 1–79.
  * 13 J. Tomasi, B. Mennucci, and R. Cammi, “Quantum mechanical continuum solvation models”, Chem. Rev., 105, 2999–3093 (2005).
  * 14 J. Tomasi, “Selected features of the polarizable continuum model for the representation of solvation”, WIREs Comput. Mol. Sci., 1, 855–867 (2011).
  * 15 B. Mennucci, “Polarizable continuum model”, WIREs Comput. Mol. Sci., 2, 386–404 (2012).
  * 16 F. Lipparini and B. Mennucci, “Perspective: Polarizable continuum models for quantum-mechanical descriptions”, J. Chem. Phys., 144, 160901:1–9 (2016).
  * 17 A. W. Lange and J. M. Herbert, “Polarizable continuum reaction-field solvation models affording smooth potential energy surfaces”, J. Phys. Chem. Lett., 1, 556–561 (2010).
  * 18 A. W. Lange and J. M. Herbert, “A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: The switching/Gaussian approach”, J. Chem. Phys., 133, 244111:1–18 (2010).
  * 19 J. M. Herbert and A. W. Lange, “Polarizable continuum models for (bio)molecular electrostatics: Basic theory and recent developments for macromolecules and simulations”, in Many-Body Effects and Electrostatics in Biomolecules, Q. Cui, P. Ren, and M. Meuwly, Eds.;  CRC Press: Boca Raton, 2016;  chapter 11, pages 363–416.
  * 20 A. W. Lange and J. M. Herbert, “Symmetric versus asymmetric discretization of the integral equations in polarizable continuum solvation models”, Chem. Phys. Lett., 509, 77–87 (2011).
  * 21 A. W. Lange, J. M. Herbert, B. J. Albrecht, and Z.-Q. You, “Intrinsically smooth discretization of Connolly’s solvent-excluded molecular surface”, Mol. Phys., 118, e1644384:1–18 (2020).
  * 22 C. Cappelli, “Integrated QM/polarizable MM/continuum approaches to model chiroptical properties of strongly interacting solute–solvent systems”, Int. J. Quantum Chem., 116, 1532–1542 (2016).
  * 23 M. Caricato, “Coupled cluster theory with the polarizable continuum model of solvation”, Int. J. Quantum Chem., 119, e25674:1–19 (2019).
  * 24 O. Andreussi and G. Fisicaro, “Continuum embeddings in condensed-matter simulations”, Int. J. Quantum Chem., 119, e25725:1–17 (2019).
  * 25 J. Sadlej and M. Pecul, “Computational modelling of the solvent–solute effect on NMR parameters by a polarizable continuum model”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 125–144.
  * 26 V. Barone, P. Cimino, and M. Pavone, “EPR spectra of organic free radicals in solution from an integrated computational approach”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 145–166.
  * 27 C. Cappelli, “Continuum solvation approach to vibrational properties”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 167–179.
  * 28 M. Pecul and K. Ruud, “Solvent effects on natural optical activity”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 180–205.
  * 29 R. Cammi and B. Mennucci, “Macroscopic nonlinear optical properties from cavity models”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 238–251.
  * 30 H. Ågren and K. V. Mikkelsen, “Homogeneous and heterogeneous solvent models for nonlinear optical properties”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 282–299.
  * 31 R. Improta, “UV–visible absorption and emission energies in condensed phase by PCM/TD-DFT methods”, in Computational Strategies for Spectroscopy: From Small Molecules to Nano Systems, V. Barone, Ed.;  John Wiley & Sons: Hoboken, NJ, 1st ed., 2012;  chapter 1, pages 39–76.
  * 32 M. V. Basilevsky and D. F. Parsons, “An advanced continuum medium model for treating solvation effects: Nonlocal electrostatics with a cavity”, J. Chem. Phys., 105, 3734–3746 (1996).
  * 33 O. S. Jenkins and K. L. C. Hunt, “Nonlocal dielectric functions on the nanoscale: Screened forces from unscreened potentials”, J. Chem. Phys., 119, 8250–8256 (2003).
  * 34 M. V. Basilevsky and G. N. Chuev, “Nonlocal solvation theories”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 94–109.
  * 35 D. H. Whiffen, “Manual of symbols and terminology for physicochemical quantities and units”, Pure Appl. Chem., 51, 1–41 (1979).
  * 36 A. Rizzo, “Birefringences in liquids”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 252–264.
  * 37 A. Ferrarini, “Anisotropic fluids”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 265–281.
  * 38 R. K. Wangsness, Electromagnetic Fields, John Wiley & Sons: Hoboken, 2nd ed., 1986.
  * 39 M. P. Coons and J. M. Herbert, “Quantum chemistry in arbitrary dielectric environments: Theory and implementation of nonequilibrium Poisson boundary conditions and application to compute vertical ionization energies at the air/water interface”, J. Chem. Phys., 148, 222834:1–21 (2018). Erratum: J. Chem. Phys., 151, 189901:1–2 (2019).
  * 40 K. Patkowski, “Recent developments in symmetry-adapted perturbation theory”, WIREs Comput. Mol. Sci., 10, e1452:1–47 (2020).
  * 41 C. J. F. Böttcher, Theory of Electric Polarization, Vol. 1, Elsevier: Amsterdam, 2nd ed., 1976.
  * 42 L. D. Jacobson, C. F. Williams, and J. M. Herbert, “The static-exchange electron-water pseudopotential, in conjunction with a polarizable water model: A new Hamiltonian for hydrated-electron simulations”, J. Chem. Phys., 130, 124115:1–18 (2009).
  * 43 M. Caricato, G. Scalmani, and M. J. Frisch, “A Lagrangian formulation for continuum models”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 64–81.
  * 44 H. Nakamura, “Roles of electrostatic interaction in proteins”, Q. Rev. Biophys., 29, 1–90 (1996).
  * 45 E. Alexov, E. L. Mehler, N. Baker, A. M. Baptista, Y. Huang, F. Milletti, J. E. Nielsen, D. Farrell, T. Carstensen, M. H. M. Olsson, J. K. Shen, J. Warwicker, S. Williams, and J. M. Word, “Progress in the prediction of p${K}_{\text{a}}$ values in proteins”, Proteins, 79, 3260–3275 (2011).
  * 46 J. Antosiewicz, J. A. McCammon, and M. K. Gilson, “Prediction of pH-dependent properties of proteins”, J. Mol. Biol., 238, 415–436 (1994).
  * 47 E. Demchuk and R. C. Wade, “Improving the continuum dielectric approach to calculating p${K}_{\text{a}}$s of ionizable groups in proteins”, J. Phys. Chem., 100, 17373–17387 (1996).
  * 48 T. Grycuk, “Revision of the model system concept for the prediction of p${K}_{\text{a}}$’s in proteins”, J. Phys. Chem. B, 106, 1434–1445 (2002).
  * 49 A. Warshel and S. T. Russell, “Calculations of electrostatic interactions in biological systems and in solutions”, Q. Rev. Biophys., 17, 283–422 (1984).
  * 50 C. N. Schutz and A. Warshel, “What are the dielectric “constants” of proteins and how to validate electrostatic models?”, Proteins, 44, 400–417 (2001).
  * 51 A. Warshel, P. K. Sharma, M. Kato, and W. W. Parson, “Modeling electrostatic effects in proteins”, Biochim. Biophys. Acta, 1764, 1647–1676 (2006).
  * 52 L. Li, C. Li, Z. Zhang, and E. Alexov, “On the dielectric “constant” of proteins: Smooth dielectric function for macromolecular modeling and its implementation in DelPhi”, J. Chem. Theory Comput., 9, 2126–2136 (2013).
  * 53 K. A. Sharp and B. Honig, “Electrostatic interactions in macromolecules”, Annu. Rev. Biophys. Biophys. Chem., 19, 301–332 (1990).
  * 54 F. Fogolari, A. Brigo, and H. Molinari, “The Poisson–Boltzmann equation for biomolecular electrostatics: A tool for structural biology”, J. Mol. Recognit., 15, 377–392 (2002).
  * 55 N. A. Baker, “Biomolecular applications of Poisson-Boltzmann methods”, in Reviews in Computational Chemistry, K. Lipkowitz, R. Larter, and T. R. Cundari, Eds., Vol. 21;  John Wiley & Sons: Hoboken, 2005;  pages 349–379.
  * 56 W. M. Botello-Smith, Q. Cai, and R. Luo, “Biological applications of classical electrostatics methods”, J. Theor. Comput. Chem., 13, 1440008:1–25 (2014).
  * 57 M. J. Schnieders and J. W. Ponder, “Polarizable atomic multipole solutes in a generalized Kirkwood continuum”, J. Chem. Theory Comput., pages 2083–2097 (2007).
  * 58 M. J. Schnieders, N. A. Baker, P. Ren, and J. W. Ponder, “Polarizable atomic multipole solutes in a Poisson-Boltzmann continuum”, J. Chem. Phys., 126, 124114:1–21 (2007).
  * 59 K. A. Sharp and B. Honig, “Calculating total electrostatic energies with the nonlinear Poisson–Boltzmann equation”, J. Phys. Chem., 94, 7684–7692 (1990).
  * 60 M. Deserno and C. Holm, “Cell model and Poisson-Boltzmann theory: A brief introduction”, in Electrostatic Effects in Soft Matter and Biophysics, C. Holm, P. Kékicheff, and R. Podgornik, Eds., Vol. 46 of NATO Science Series; Springer Science+Business Media: Dordrecht, 2001;  pages 27–52.
  * 61 G. Lamm, “The Poisson-Boltzmann equation”, in Reviews in Computational Chemistry, K. B. Lipkowitz, R. Larter, T. R. Cundari, and D. B. Boyd, Eds., Vol. 19;  Wiley-VCH: New York, 2003;  chapter 4, pages 147–366.
  * 62 N. A. Baker, “Poisson–Boltzmann methods for biomolecular electrostatics”, Method. Enzymol., 383, 94–118 (2004).
  * 63 A. G. Moreira and R. R. Netz, “Field-theoretic approaches to classical charged systems”, in Electrostatic Effects in Soft Matter and Biophysics, C. Holm, P. Kékicheff, and R. Podgornik, Eds., Vol. 46 of NATO Science Series; Springer Science+Business Media: Dordrecht, 2001;  pages 367–408.
  * 64 C. J. Stein, J. M. Herbert, and M. Head-Gordon, “The Poisson–Boltzmann model for implicit solvation of electrolyte solutions: Quantum chemical implementation and assessment via Sechenov coefficients”, J. Chem. Phys., 151, 224111:1–14 (2019).
  * 65 H.-X. Zhou, “Macromolecular electrostatic energy within the nonlinear Poisson–Boltzmann equation”, J. Chem. Phys., 100, 3152–3162 (1993).
  * 66 F. Fogolari, P. Zuccato, G. Esposito, and P. Viglino, “Biomolecular electrostatics with the linearized Poisson-Boltzmann equation”, Biophys. J., 76, 1–16 (1999).
  * 67 P. Debye and E. Hückel, “On the theory of electrolytes. I. Freezing point depression and related phenomena”, in Collected Papers of Peter J. W. Debye; Interscience: New York, 1954;  pages 217–263.
  * 68 A. W. Lange and J. M. Herbert, “A simple polarizable continuum solvation model for electrolyte solutions”, J. Chem. Phys., 134, 204110:1–15 (2011).
  * 69 C. Wang, P. Ren, and R. Luo, “Ionic solution: What goes right and wrong with continuum solvation modeling”, J. Phys. Chem. B, 121, 11159–11179 (2017).
  * 70 V. Vlachy, “Ionic effects beyond Poisson-Boltzmann theory”, Annu. Rev. Phys. Chem., 50, 145–165 (1999).
  * 71 J. Dziedzic, A. Bhandari, L. Anton, C. Peng, J. Womack, M. Famili, D. Kramer, and C.-K. Skylaris, “Practical approach to large-scale electronic structure calculations in electrolyte solutions via continuum-embedded linear-scaling density functional theory”, J. Phys. Chem. C, 124, 7860–7872 (2020).
  * 72 M. Born, “Volumen und Hydratationswärme der Ionen”, Z. Phys., 1, 45–48 (1920).
  * 73 A. A. Rashin and B. Honig, “Reevaluation of the Born model of ion hydration”, J. Phys. Chem., 89, 5588–5593 (1985).
  * 74 L. Onsager, “Electric moments of molecules in liquids”, J. Am. Chem. Soc., 58, 1486–1493 (1936).
  * 75 R. P. Bell, “The electrostatic energy of dipole molecules in different media”, Trans. Faraday Soc., 27, 797–802 (1931).
  * 76 C. J. F. Böttcher, “The dielectric constant of dipole liquids”, Physica, 5, 635–639 (1938).
  * 77 J. B. Hasted, “Liquid water: Dielectric properties”, in Water: A Comprehensive Treatise, F. Franks, Ed., Vol. 1;  Plenum Press: New York, 1972;  pages 255–309.
  * 78 J. S. Høye and G. Stell, “Statistical mechanics of polar systems. II”, J. Chem. Phys., 64, 1952–1966 (1976).
  * 79 J. H. Hannay, “The Clausius–Mossotti equation: An alternative derivation”, Eur. J. Phys., 4, 141–143 (1983).
  * 80 J. G. Kirkwood, “Theory of solutions of molecules containing widely separated charges with special application to zwitterions”, J. Chem. Phys., 2, 351–361 (1934).
  * 81 J. G. Kirkwood and F. H. Westheimer, “The electrostatic influence of substituents on the dissociation constants of organic acids. I.”, J. Chem. Phys., 6, 506–512 (1938).
  * 82 F. H. Westheimer and J. G. Kirkwood, “The electrostatic influence of substituents on the dissociation constants of organic acids. II.”, J. Chem. Phys., 6, 513–517 (1938).
  * 83 W. B. Bonner, “The electrostatic energy of molecules in solution”, Trans. Faraday Soc., 47, 1143–1152 (1951).
  * 84 D. Rinaldi, M. F. Ruiz-Lopez, and J.-L. Rivail, “Ab initio SCF calculations on electrostatically solvated molecules using a deformable three axes ellipsoidal cavity”, J. Chem. Phys., 78, 834–838 (1983).
  * 85 K. V. Mikkelsen, E. Dalgaard, and P. Swanstrøm, “Electron-transfer reactions in solution. An ab initio approach”, J. Phys. Chem., 91, 3081–3092 (1987).
  * 86 K. V. Mikkelsen, H. Ågren, H. J. A. Jensen, and T. Helgaker, “A multiconfigurational self-consistent reaction-field method”, J. Chem. Phys., 89, 3086–3095 (1988).
  * 87 Y. Kong and J. W. Ponder, “Calculation of the reaction field due to off-center point multipoles”, J. Chem. Phys., 107, 481–492 (1997).
  * 88 K. V. Mikkelsen, P. Jørgensen, and H. J. A. Jensen, “A multiconfigurational self-consistent reaction field response method”, J. Chem. Phys., 100, 6597–6607 (1994).
  * 89 M. Medved’, Š. Budzák, W. Bartkowiak, and H. Reis, “Solvent effects on molecular electric properties”, in Handbook of Computational Chemistry, J. Leszczynski, A. Kaczmarek-Kedziera, T. Puzyn, M. G. Papadopoulos, H. Reis, and M. K. Shukla, Eds.;  Springer International Publishing: Switzerland, 2nd ed., 2017;  chapter 17, pages 741–794.
  * 90 I. Lotan and T. Head-Gordon, “An analytical electrostatic model for salt screened interactions between multiple proteins”, J. Chem. Theory Comput., 2, 541–555 (2006).
  * 91 C.-G. Zhan, J. Bentley, and D. M. Chipman, “Volume polarization in reaction field theory”, J. Chem. Phys., 108, 177–192 (1998).
  * 92 D. M. Chipman, “Charge penetration in dielectric models of solvation”, J. Chem. Phys., 106, 10194–10206 (1997).
  * 93 D. M. Chipman, “Simulation of volume polarization in reaction field theory”, J. Chem. Phys., 110, 8012–8018 (1999).
  * 94 D. M. Chipman, “Reaction field treatment of charge penetration”, J. Chem. Phys., 112, 5558–5565 (2000).
  * 95 D. M. Chipman, “Comparison of solvent reaction field representations”, Theor. Chem. Acc., 107, 80–89 (2002).
  * 96 D. M. Chipman, “New formulation and implementation for volume polarization in dielectric continuum theory”, J. Chem. Phys., 124, 224111:1–10 (2006).
  * 97 D. M. Chipman, “Energy correction to simulation of volume polarization in reaction field theory”, J. Chem. Phys., 116, 10129–10138 (2002).
  * 98 F. J. Luque, C. Curutchet, J. Muñoz-Muriedas, A. Bidon-Chanal, I. Soteras, A. Morreale, J. L. Gelpi, and M. Orozco, “Continuum solvation models: Dissecting the free energy of solvation”, Phys. Chem. Chem. Phys., 5, 3827–3836 (2003).
  * 99 D. Beglov and B. Roux, “Finite representation of an infinite bulk system: Solvent boundary potential for computer simulations”, J. Chem. Phys., 100, 9050–9063 (1994).
  * 100 I. G. Tironi, R. Sperb, P. E. Smith, and W. F. van Gunsteren, “A generalized reaction field method for molecular dynamics simulations”, J. Chem. Phys., 102, 5451–5459 (1995).
  * 101 W. Im, S. Bernèche, and B. Roux, “Generalized solvent boundary potential for computer simulations”, J. Chem. Phys., 114, 2924–2937 (2001).
  * 102 P. Schaefer, D. Riccardi, and Q. Cui, “Reliable treatment of electrostatics in combined QM/MM simulation of macromolecules”, J. Chem. Phys., 123, 014905:1–14 (2005).
  * 103 T. Benighaus and W. Thiel, “Efficiency and accuracy of the generalized solvent boundary potential for hybrid QM/MM simulations: Implementation for semiempirical Hamiltonians”, J. Chem. Theory Comput., 4, 1600–1609 (2008).
  * 104 T. Benighaus and W. Thiel, “A general boundary potential for hybrid QM/MM simulations of solvated biomolecular systems”, J. Chem. Theory Comput., 5, 3114–3128 (2009).
  * 105 T. Benighaus and W. Thiel, “Long-range electrostatic effects in QM/MM studies of enzymatic reactions: Application of the solvated macromolecule boundary potential”, J. Chem. Theory Comput., 7, 238–249 (2011).
  * 106 A. Aleksandrov and M. Field, “Efficient solvent boundary potential for hybrid potential simulations”, Phys. Chem. Chem. Phys., 13, 10503–10509 (2011).
  * 107 J. Zienau and Q. Cui, “Implementation of the solvent macromolecular boundary potential and application to model and realistic enzyme systems”, J. Phys. Chem. B, 116, 12522–12534 (2012).
  * 108 X. Lu and Q. Cui, “Charging free energy calculations using the generalized solvent boundary potential (GSBP) and periodic boundary condition: A comparative analysis using ion solvation and oxidation free energy in proteins”, J. Phys. Chem. B, 117, 2005–2018 (2013).
  * 109 N. Rega, G. Brancato, and V. Barone, “Non-periodic boundary conditions for ab initio molecular dynamics in condensed phase using localized basis functions”, Chem. Phys. Lett., 422, 367–371 (2006).
  * 110 G. Brancato, N. Rega, and V. Barone, “Reliable molecular simulations of solute-solvent systems with a minimum number of solvent shells”, J. Chem. Phys., 124, 214505:1–9 (2006).
  * 111 G. Brancato, N. Rega, and V. Barone, “A hybrid explicit/implicit solvation method for first-principle molecular dynamics simulations”, J. Chem. Phys., 128, 144501:1–10 (2008).
  * 112 G. Brancato, N. Rega, and V. Barone, “Molecular dynamics simulations in a ${N}p{T}$ ensemble using non-periodic boundary conditions”, Chem. Phys. Lett., 483, 177–181 (2009).
  * 113 J. M. J. Swanson, S. A. Adcock, and J. A. McCammon, “Optimized radii for Poisson–Boltzmann calculations with the AMBER force field”, J. Chem. Theory Comput., 1, 484–493 (2005).
  * 114 A. Bondi, “Van der Waals volumes and radii”, J. Phys. Chem., 68, 441–451 (1964).
  * 115 R. S. Rowland and R. Taylor, “Intermolecular nonbonded contact distances in organic crystal structures: Comparison with distances expected from van der Waals radii”, J. Phys. Chem., 100, 7384–7391 (1996).
  * 116 M. Manjeera, A. C. Chamberlin, R. Valero, C. J. Cramer, and D. G. Truhlar, “Consistent van der Waals radii for the whole main group”, J. Phys. Chem. A, 113, 5806–5812 (2009).
  * 117 R. Bonaccorsi, P. Palla, and J. Tomasi, “Conformational energy of glycine in aqueous solutions and relative stability of the zwitterionic and neutral forms. An ab initio study”, J. Am. Chem. Soc., 106, 1945–1950 (1984).
  * 118 D. H. Brookes and T. Head-Gordon, “Family of oxygen–oxygen radial distribution functions for water”, J. Phys. Chem. Lett., 6, 2938–2943 (2015).
  * 119 A. V. Onufriev and B. Aguilar, “Accuracy of continuum electrostatic calculations based on three common dielectric boundary definitions”, J. Theor. Comput. Chem., 13, 1440006:1–25 (2014).
  * 120 M. L. Connolly, “Solvent-accessible surfaces of proteins and nucleic acids”, Science, 221, 709–713 (1983).
  * 121 F. M. Richards, “Areas, volumes, packing, and protein structure”, Annu. Rev. Biophys. Bio., 6, 151–176 (1977).
  * 122 B. Lee and F. M. Richards, “The interpretation of protein structures: Estimation of static accessibility”, J. Mol. Biol., 55, 379–400 (1971).
  * 123 D.-S. Kim, C.-I. Won, and J. Bhak, “A proposal for the revision of molecular boundary typology”, J. Biomol. Struct. Dyn., 28, 277–287 (2010).
  * 124 C. S. Pomelli, “Cavity surfaces and their discretization”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 49–63.
  * 125 C.-G. Zhan and D. M. Chipman, “Cavity size in reaction field theory”, J. Chem. Phys., 109, 10543–10558 (1998).
  * 126 V. Barone, M. Cossi, and J. Tomasi, “A new definition of cavities for the computation of solvation free energies by the polarizable continuum model”, J. Chem. Phys., 107, 3210–3221 (1997).
  * 127 B. Ginovska, D. M. Camaioni, M. Dupuis, C. A. Schwerdtfeger, and Q. Gil, “Charge-dependent cavity radii for an accurate dielectric continuum model of solvation with emphasis on ions: Aqueous solutes with oxo, hydroxo, amino, methyl, chloro, bromo, and fluoro functionalities”, J. Phys. Chem. A, 112, 10604–10613 (2008).
  * 128 J. B. Foresman, T. A. Keith, K. B. Wiberg, J. Snoonian, and M. J. Frisch, “Solvent effects. 5. Influence of cavity shape, truncation of electrostatics, and electron correlation on ab initio reaction field calculations”, J. Phys. Chem., 100, 16098–16104 (1996).
  * 129 D. M. Chipman and M. Dupuis, “Implementation of solvent reaction fields for electronic structure”, Theor. Chem. Acc., 107, 90–102 (2002).
  * 130 F. Chen and D. M. Chipman, “Boundary element methods for dielectric cavity construction and integration”, J. Chem. Phys., 119, 10289–10297 (2003).
  * 131 P. C. do Couto and D. M. Chipman, “How does dielectric solvation affect the size of an ion?”, J. Phys. Chem. A, 114, 12788–12793 (2010).
  * 132 M. Holst and F. Saied, “Multigrid solution of the Poisson–Boltzmann equation”, J. Comput. Chem., 14, 105–113 (1993).
  * 133 M. J. Holst and F. Saied, “Numerical solution of the nonlinear Poisson–Boltzmann equation: Developing more robust and efficient methods”, J. Comput. Chem., 16, 337–364 (1995).
  * 134 M. Holst, N. Baker, and F. Wang, “Adaptive multilevel finite element solution of the Poisson–Boltzmann equation I. Algorithms and examples”, J. Comput. Chem., 21, 1319–1342 (2000). Erratum: J. Comput. Chem., 22, 45 (2001).
  * 135 B. Z. Lu, Y. C. Zhou, M. J. Holst, and J. A. McCammon, “Recent progress in numerical methods for the Poisson-Boltzmann equation in biophysical applications”, Commun. Comput. Phys., 3, 973–1009 (2008).
  * 136 J. Wang and R. Luo, “Assessment of linear finite-difference Poisson–Boltzmann solvers”, J. Comput. Chem., 31, 1689–1698 (2010).
  * 137 E.-H. Yap and T. Head-Gordon, “New and efficient Poisson–Boltzmann solver for interaction of multiple proteins”, J. Chem. Theory Comput., 6, 2214–2224 (2010).
  * 138 A. H. Boschitsch and M. O. Fenley, “A fast and robust Poisson–Boltzmann solver based on adaptive Cartesian grids”, J. Chem. Theory Comput., 7, 1524–1540 (2011).
  * 139 M. Holst, J. A. McCammon, Z. Yu, and Y. C. Zhou, “Adaptive finite element modeling techniques for the Poisson-Boltzmann equation”, Commun. Comput. Phys., 11, 179–214 (2012).
  * 140 C. Li, L. Li, M. Petukh, and E. Alexov, “Progress in developing Poisson-Boltzmann equation solvers”, Mol.-Based Math. Biol., 1, 42–62 (2013).
  * 141 W. Geng and R. Krasny, “A treecode-accelerated boundary integral Poisson–Boltzmann solver for electrostatics of solvated biomolecules”, J. Comput. Phys., 247, 62–78 (2013).
  * 142 I. Sakalli, J. Schöberl, and E. W. Knapp, “mFES: A robust molecular finite element solver for electrostatic energy computations”, J. Chem. Theory Comput., 10, 5095–5112 (2014).
  * 143 G. Fisicaro, L. Genovese, O. Andreussi, N. Marzari, and S. Goedecker, “A generalized Poisson and Poisson-Boltzmann solver for electrostatic environments”, J. Chem. Phys., 144, 014103:1–12 (2016).
  * 144 S. Ringe, H. Oberhofer, C. Hille, S. Matera, and K. Reuter, “Function-space-based solution scheme for the size-modified Poisson–Boltzmann equation in full-potential DFT”, J. Chem. Theory Comput., 12, 4052–4066 (2016).
  * 145 J. C. Womack, L. Anton, J. Dziedzic, P. J. Hasnip, M. I. J. Probert, and C.-K. Skylaris, “DL_MG: A parallel multigrid Poisson and Poisson–Boltzmann solver for electronic structure calculations in vacuum and solution”, J. Chem. Theory Comput., 14, 1412–1432 (2018).
  * 146 V. Luzhkov and A. Warshel, “Microscopic models for quantum mechanical calculations of chemical processes in solutions: LD/AMPAC and SCAAS/AMPAC calculations of solvation energies”, J. Comput. Chem., 13, 199–213 (1992).
  * 147 J. Florián and A. Warshel, “Langevin dipoles model for ab initio calculations of chemical processes in solution: Parameterization and application to hydration free energies of neutral and ionic solutes and conformational analysis in aqueous solution”, J. Phys. Chem. B, 101, 5583–5595 (1997).
  * 148 A. Papazyan and A. Warshel, “Continuum and dipole-lattice models of solvation”, J. Phys. Chem. B, 101, 11254–11264 (1997).
  * 149 J. Langlet, P. Claverie, J. Caillet, and A. Pullman, “Improvements of the continuum model. 1. Application to the calculation of the vaporization of thermodynamic quantities of nonassociated liquids”, J. Phys. Chem., 92, 1617–1631 (1988).
  * 150 M. Orozco and F. J. Luque, “Theoretical methods for the description of the solvent effect in biomolecular systems”, Chem. Rev., 100, 4187–4225 (2000). Erratum:Chem. Rev., 101, 203 (2001).
  * 151 R. Constanciel, “Theoretical basis of the empirical reaction field approximations through continuum model”, Theor. Chem. Acc., 69, 505–523 (1986).
  * 152 E. Cancès, “Integral equation approaches for continuum models”, in Continuum Solvation Models in Chemical Physics, B. Mennucci and R. Cammi, Eds.;  Wiley: Chichester, UK, 2007;  pages 29–48.
  * 153 M. Cossi, G. Scalmani, N. Rega, and V. Barone, “New developments in the polarizable continuum model for quantum mechanical and classical calculations on molecules in solution”, J. Chem. Phys., 117, 43–54 (2002).
  * 154 J. Tomasi, B. Mennucci, and E. Cancès, “The IEF version of the PCM solvation method: An overview of a new method addressed to study molecular solutes at the QM ab initio level”, J. Mol. Struct. (Theochem), 464, 211–226 (1999).
  * 155 E. Cancés, B. Mennucci, and J. Tomasi, “A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics”, J. Chem. Phys., 107, 3032–3041 (1997).
  * 156 B. Mennucci, E. Cancés, and J. Tomasi, “Evaluation of solvent effects in isotropic and anisotropic dielectrics and in ionic solutions with a unified integral equation method: Theoretical bases, computational implementation, and numerical applications”, J. Phys. Chem. B, 101, 10506–10517 (1997).
  * 157 E. Cancès and B. Mennucci, “New applications of integral equations methods for solvation continuum models: Ionic solutions and liquid crystals”, J. Math. Chem., 23, 309–326 (1998).
  * 158 E. Cancès and B. Mennucci, “Comment on “Reaction field treatment of charge penetration””, J. Chem. Phys., 114, 4744–4745 (2001).
  * 159 M. E. Davis and J. A. McCammon, “Solving the finite difference linearized Poisson-Boltzmann equation: A comparison of relaxation and conjugate gradient methods”, J. Comput. Chem., 10, 386–391 (1989).
  * 160 B. A. Luty, M. E. Davis, and J. A. McCammon, “Solving the finite-difference non-linear Poisson–Boltzmann equation”, J. Comput. Chem., 13, 1114–1118 (1992).
  * 161 J. Wang, Q. Cai, Y. Xiang, and R. Luo, “Reducing grid dependence in finite-difference Poisson–Boltzmann calculations”, J. Chem. Theory Comput., 8, 2741–2751 (2012).
  * 162 L. Xiao, Q. Cai, X. Ye, J. Wang, and R. Luo, “Electrostatic forces in the Poisson-Boltzmann systems”, J. Chem. Phys., 139, 094106:1–12 (2013).
  * 163 L. Xiao, C. Wang, and R. Luo, “Recent progress in adapting Poisson–Boltzmann methods to molecular simulations”, J. Theor. Comput. Chem., 13, 1430001:1–19 (2014).
  * 164 D. M. Chipman, “Solution of the linearized Poisson–Boltzmann equation”, J. Chem. Phys., 120, 5566–5575 (2004).
  * 165 F. Lipparini, B. Stamm, E. Cancès, Y. Maday, and B. Mennucci, “Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives”, J. Chem. Theory Comput., 9, 3637–3648 (2013).
  * 166 F. Lipparini, L. Lagardère, G. Scalmani, B. Stamm, E. Cancès, Y. Maday, J.-P. Piquemal, M. J. Frisch, and B. Mennucci, “Quantum calculations in solution for large to very large molecules: A new linear scaling QM/continuum approach”, J. Phys. Chem. Lett., 5, 953–958 (2014).
  * 167 S. Caprasecca, S. Jurinovich, L. Lagardère, B. Stamm, and F. Lipparini, “Achieving linear scaling in computational cost for a fully polarizable MM/continuum embedding”, J. Chem. Theory Comput., 11, 694–704 (2015).
  * 168 A. Klamt and G. Schüürmann, “COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient”, J. Chem. Soc., Perkin Trans. 2, pages 799–805 (1993).
  * 169 E. Cancès and B. Mennucci, “The escaped charge problem in solvation continuum models”, J. Chem. Phys., 115, 6130–6135 (2001).
  * 170 A. Klamt, C. Moya, and J. Palomar, “A comprehensive comparison of the IEFPCM and SS(V)PE continuum solvation methods with the COSMO approach”, J. Chem. Theory Comput., 11, 4220–4225 (2015).
  * 171 E. V. Stefanovich and T. N. Truong, “Optimized atomic radii for quantum dielectric continuum solvation models”, Chem. Phys. Lett., 244, 65–74 (1995).
  * 172 T. N. Truong and E. V. Stefanovich, “A new method for incorporating solvent effects into the classical, ab initio molecular orbital and density functional theory frameworks for arbitrary cavity shape”, Chem. Phys. Lett., 240, 253–260 (1995).
  * 173 T. N. Truong and E. V. Stefanovich, “Analytical first and second energy derivatives of the generalized conductorlike screening model for free energy of solvation”, J. Chem. Phys., 103, 3709–3717 (1995).
  * 174 T. N. Truong, U. N. Nguyen, and E. V. Stefanovich, “Generalized conductor-like screening model (GCOSMO) for solvation: An assessment of its accuracy and applicability”, Int. J. Quantum Chem. Symp., 60, 1615–1622 (1996).
  * 175 V. Barone and M. Cossi, “Quantum calculation of molecular energies and energy gradients in solution by a conductor solvent model”, J. Phys. Chem. A, 102, 1995–2001 (1998).
  * 176 M. Cossi, N. Rega, G. Scalmani, and V. Barone, “Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model”, J. Comput. Chem., 24, 669–681 (2003).
  * 177 R. Cammi and J. Tomasi, “Analytical derivatives for molecular solutes. I. Hartree–Fock energy first derivatives with respect to external parameters in the polarizable continuum model”, J. Chem. Phys., 100, 7495–7502 (1994).
  * 178 A. Klamt and V. Jonas, “Treatment of the outlying charge in continuum solvation models”, J. Chem. Phys., 105, 9972–9981 (1996).
  * 179 B. Mennucci and J. Tomasi, “Continuum solvation models: A new approach to the problem of solute’s charge distribution and cavity boundaries”, J. Chem. Phys., 106, 5151–5158 (1997).
  * 180 R. Cammi, M. Cossi, and J. Tomasi, “Analytical derivatives for molecular solutes. III. Hartree–Fock static polarizabilities in the polarizable continuum model”, J. Chem. Phys., 104, 4611–4620 (1996).
  * 181 J. Gauss, “Molecular properties”, in Modern Methods and Algorithms of Quantum Chemistry, J. Grotendorst, Ed., Vol. 3 of NIC Series; John von Neumann Institute for Computing: Jülich, 2nd ed., 2000;  pages 541–592.
  * 182 A. Rizzo, S. Coriani, and K. Ruud, “Response function theory computational approaches to linear and nonlinear optical spectroscopy”, in Computational Strategies for Spectroscopy: From Small Molecules to Nano Systems, V. Barone, Ed.;  John Wiley & Sons: Hoboken, NJ, 1st ed., 2012;  chapter 2, pages 77–136.
  * 183 T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen, and K. Ruud, “Recent advances in wave function-based methods of molecular-property calculations”, Chem. Rev., 112, 543–631 (2012).
  * 184 E. Jurrus, D. Engel, K. Star, K. Monson, J. Brandi, L. E. Felberg, D. H. Brookes, L. Wilson, J. Chen, K. Liles, M. Chun, P. Li, D. W. Gohara, T. Dolinsky, R. K. D. R. Koes, J. E. Nielsen, T. Head-Gordon, W. Geng, R. Krasny, G.-W. Wei, M. J. Holst, J. A. McCammon, and N. A. Baker, “Improvements to the APBS biomolecular solvation software suite”, Protein Sci., 27, 112–128 (2018).
  * 185 Z.-Q. You, J.-M. Mewes, A. Dreuw, and J. M. Herbert, “Comparison of the Marcus and Pekar partitions in the context of non-equilibrium, polarizable-continuum reaction-field solvation models”, J. Chem. Phys., 143, 204107:1–14 (2015).
  * 186 E. Silla, F. Villar, O. Nilsson, J. L. Pascual-Ahuir, and O. Tapia, “Molecular volumes and surfaces of biomacromolecules via GEPOL: A fast and efficient algorithm”, J. Mol. Graphics, 8, 168–172 (1990).
  * 187 J. L. Pascual-Ahuir and E. Silla, “GEPOL: An improved description of molecular surfaces. I. Building the spherical surface set”, J. Comput. Chem., 11, 1047–1060 (1990).
  * 188 E. Silla, I. Tuñon, and J. L. Pascual-Ahuir, “GEPOL: An improved description of molecular surfaces. II. Computing the molecular area and volume”, J. Comput. Chem., 12, 1077–1088 (1991).
  * 189 J. L. Pascual-Ahuir, E. Silla, and I. Tuñon, “GEPOL: An improved description of molecular surfaces. III. A new algorithm for the computation of a solvent-excluding surface”, J. Comput. Chem., 15, 1127–1138 (1994).
  * 190 D. A. Liotard, G. D. Hawkins, G. C. Lynch, C. J. Cramer, and D. G. Truhlar, “Improved methods for semiempirical solvation models”, J. Comput. Chem., 16, 422–440 (1995).
  * 191 M. Cossi, B. Mennucci, and R. Cammi, “Analytical first derivatives of molecular surfaces with respect to nuclear coordinates”, J. Comput. Chem., 17, 57–73 (1996).
  * 192 D. M. York and M. Karplus, “Smooth solvation potential based on the conductor-like screening model”, J. Phys. Chem. A, 103, 11060–11079 (1999).
  * 193 B. A. Gregersen and D. M. York, “High-order discretization schemes for biochemical applications of boundary element solvation and variational electrostatic projection methods”, J. Chem. Phys., 122, 194110:1–4 (2005).
  * 194 J. Khandogin, B. A. Gregersen, W. Thiel, and D. M. York, “Smooth solvation method for d-orbital semiemprical calculations of biological reactions. 1. Implementation”, J. Phys. Chem. B, 109, 9799–9809 (2005).
  * 195 C. W. Murray, N. C. Handy, and G. J. Laming, “Quadrature schemes for integrals of density functional theory”, Mol. Phys., 78, 997–1014 (1993).
  * 196 P. M. W. Gill, B. G. Johnson, and J. A. Pople, “A standard grid for density-functional calculations”, Chem. Phys. Lett., 209, 506–512 (1993).
  * 197 S.-H. Chien and P. M. W. Gill, “SG-0: A small standard grid for DFT quadrature on large systems”, J. Comput. Chem., 27, 730–739 (2006).
  * 198 S. Dasgupta and J. M. Herbert, “Standard grids for high-precision integration of modern density functionals: SG-2 and SG-3”, J. Comput. Chem., 38, 869–882 (2017).
  * 199 J. Liu and W. Liang, “Analytical second derivatives of excited-state energy within the time-dependent density functional theory coupled with a conductor-like polarizable continuum model”, J. Chem. Phys., 138, 024101:1–10 (2013).
  * 200 R. J. Wawak, K. D. Gibson, and H. A. Scheraga, “Gradient discontinuities in calculations involving molecular surface area”, J. Math. Chem., 15, 207–232 (1994).
  * 201 H. Li and J. H. Jensen, “Improving the efficiency and convergence of geometry optimization with the polarizable continuum model: New energy gradients and molecular surface tessellation”, J. Comput. Chem., 25, 1449–1462 (2004).
  * 202 P. Su and H. Li, “Continuous and smooth potential energy surface for conductor-like screening solvation model using fixed points with variable areas”, J. Chem. Phys., 130, 074109:1–13 (2009).
  * 203 A. I. Krylov and P. M. W. Gill, “Q-Chem: An engine for innovation”, WIREs Comput. Mol. Sci., 3, 317–326 (2013).
  * 204 G. Scalmani and M. J. Frisch, “Continuous surface charge polarizable continuum models of solvation. I. General formalism”, J. Chem. Phys., 132, 114110:1–15 (2010).
  * 205 M. Garcia-Ratès and F. Neese, “Effect of the solute cavity on the solvation energy and its derivatives within the framework of the Gaussian charge scheme”, J. Comput. Chem., 41, 922–939 (2020).
  * 206 H. Harbrecht and M. Randrianarivony, “Wavelent BEM on molecular surfaces: Parameterization and implementation”, Computing, 86, 1–22 (2009).
  * 207 V. Weijo, M. Randrianarivony, H. Harbrecht, and L. Frediani, “Wavelet formulation of the polarizable continuum model”, J. Comput. Chem., 31, 1469–1477 (2010).
  * 208 M. Bugeanu, R. Di Remigio, K. Mozgawa, S. S. Reine, H. Harbrecht, and L. Frediani, “Wavelet formulation of the polarizable continuum model. II. Use of piecewise bilinear boundary elements”, Phys. Chem. Chem. Phys., 17, 31566–31581 (2015).
  * 209 M. Bugeanu and H. Harbrecht, “Parametric representation of molecular surfaces”, Int. J. Quantum Chem., 119, e25695:1–12 (2019).
  * 210 C.-G. Zhan and D. M. Chipman, “Reaction field effects on nitrogen shielding”, J. Chem. Phys., 110, 1611–1622 (1999).
  * 211 H. Tjong and H.-X. Zhou, “On the dielectric boundary in Poisson–Boltzmann calculations”, J. Chem. Theory Comput., 4, 507–514 (2008).
  * 212 W. E. Lorensen and H. E. Cline, “Marching cubes: A high resolution 3D surface construction algorithm”, Comp. Graph., 21, 163–169 (1987).
  * 213 D. A. Rajon and W. E. Bolch, “Marching cubes algorithm: Review and trilinear interpolation adaptation for image-based dosimetric models”, Comput. Med. Imag. Grap., 27, 411–435 (2003).
  * 214 Z. Yu, M. P. Jacobson, and R. Friesner, “What role do surfaces play in GB models? A new-generation of surface-generalized Born model based on a novel Gaussian surface for biomolecules”, J. Comput. Chem., 27, 72–89 (2005).
  * 215 B. Zhou, M. Agarwal, and C. F. Wong, “Variable atomic radii for continuum-solvent electrostatics calculation”, J. Chem. Phys., 129, 014509:1–9 (2008).
  * 216 G. Fisicaro, L. Genovese, O. Andreussi, S. Mandai, N. N. Nair, N. Marzari, and S. Goedecker, “Soft-sphere continuum solvation in electronic-structure calculations”, J. Chem. Theory Comput., 13, 3829–3845 (2017).
  * 217 J.-L. Fattebert and F. Gygi, “Density functional theory for efficieint ab initio molecular dynamics simulations in solution”, J. Comput. Chem., 23, 662–666 (2002).
  * 218 J.-L. Fattebert and F. Gygi, “First-principles molecular dynamics simulations in a continuum solvent”, Int. J. Quantum Chem., 93, 139–147 (2003).
  * 219 D. A. Scherlis, J.-L. Fattebert, F. Gygi, M. Cococcioni, and N. Marzari, “A unified electrostatic and cavitation model for first-principles molecular dynamics in solution”, J. Chem. Phys., 124, 074103:1–12 (2006).
  * 220 J. Dziedzic, H. H. Helal, C.-K. Skylaris, A. A. Mostofi, and M. C. Payne, “Minimal parameter implicit solvent model for ab initio electronic-structure calculations”, Europhys. Lett., 95, 43001:1–6 (2011).
  * 221 O. Andreussi, I. Dabo, and N. Marzari, “Revised self-consistent continuum solvation in electronic-structure calculations”, J. Chem. Phys., 136, 064102:1–20 (2012).
  * 222 K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T. A. Arias, and R. G. Hennig, “Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways”, J. Chem. Phys., 140, 084106:1–8 (2014).
  * 223 V. M. Sánchez, M. Sued, and D. A. Scherlis, “First-principles molecular dynamics simulations at solid-liquid interfaces with a continuum solvent”, J. Chem. Phys., 131, 174108:1–9 (2009).
  * 224 F. Nattino, M. Truscott, N. Marzari, and O. Andreussi, “Continuum models of the electrochemical diffuse layer in electronic-structure calculations”, J. Chem. Phys., 150, 041722:1–17 (2019).
  * 225 O. Andreussi, N. G. Hörmann, F. Nattino, G. Fisicaro, S. Goedecker, and N. Marzari, “Solvent-aware interfaces in continuum solvation”, J. Chem. Theory Comput., 15, 1996–2009 (2019).
  * 226 R. Sundararaman and K. Schwarz, “Evaluating continuum solvation models for the electrode-electrolyte interface: Challenges and strategies for improvement”, J. Chem. Phys., 146, 084111:1–5 (2017).
  * 227 R. Sundararaman, K. Letchworth-Weaver, and K. A. Schwarz, “Improving accuracy of electrochemical capacitance and solvation energetics in first-principles calculations”, J. Chem. Phys., 148, 144105:1–7 (2018).
  * 228 K. Schwarz and R. Sundararaman, “The electrochemical interface in first-principles calculations”, Surf. Sci. Rep., 75, 100492:1–22 (2020).
  * 229 G. Bramley, M.-T. Nguyen, V.-A. Glezakou, R. Rousseau, and C.-K. Sylaris, “Reconciling work functions and adsorption enthalpies for implicit solvent models: A Pt(111)/water interface case study”, J. Chem. Theory Comput., 16, 2703–2715 (2020).
  * 230 A. Bhandari, L. Anton, J. Dziedzic, C. Peng, D. Kramer, and C.-K. Skylaris, “Electronic structure calculations in electrolyte solutions: Methods for neutralization of extended charged interfaces”, J. Chem. Phys., 153, 124101:1–12 (2020).
  * 231 M. E. Davis and J. A. McCammon, “Dielectric boundary smoothing in finite difference solutions of the Poisson equation: An approach to improve accuracy and convergence”, J. Comput. Chem., 12, 909–912 (1991).
  * 232 M. K. Gilson, M. E. Davis, B. A. Luty, and J. A. McCammon, “Computation of electrostatic forces on solvated molecules using the Poisson-Boltzmann equation”, J. Phys. Chem., 97, 3591–3600 (1993).
  * 233 W. Im, D. Beglov, and B. Roux, “Continuum solvation model: Computation of electrostatic forces from numerical solutions to the Poisson–Boltzmann equation”, Comput. Phys. Commun., 111, 59–75 (1998).
  * 234 B. Lu, D. Zhang, and J. A. McCammon, “Computation of electrostatic forces between solvated molecules determined by the Poisson–Boltzmann equation using a boundary element method”, J. Chem. Phys., 122, 214102:1–7 (2005).
  * 235 J. Wang, C. Tan, E. Chanco, and R. Luo, “Quantitative analysis of Poisson–Boltzmann implicit solvent in molecular dynamics”, Phys. Chem. Chem. Phys., 12, 1194–1202 (2010).
  * 236 Q. Cai, X. Ye, J. Wang, and R. Luo, “Dielectric boundary force in numerical Poisson–Boltzmann methods: Theory and numerical strategies”, Chem. Phys. Lett., 514, 368–373 (2011).
  * 237 N. Rega, M. Cossi, and V. Barone, “Towards linear scaling in continuum solvent models. A new iterative procedure for energies and geometry optimizations”, Chem. Phys. Lett., 293, 221–229 (1998).
  * 238 G. Scalmani, V. Barone, K. N. Kudin, C. S. Pomelli, G. E. Scuseria, and M. J. Frisch, “Achieving linear-scaling computational cost for the polarizable continuum model of solvation”, Theor. Chem. Acc., 111, 90–100 (2004).
  * 239 P. Gatto, F. Lipparini, and B. Stamm, “Computation of forces arising from the polarizable continuum model within the domain-decomposition paradigm”, J. Chem. Phys., 147, 224108:1–11 (2017).
  * 240 B. Stamm, L. Lagardère, G. Scalmani, P. Gatto, E. Cancès, J.-P. Piquemal, Y. Maday, B. Mennucci, and F. Lipparini, “How to make continuum solvation incredibly fast in a few simple steps: A practical guide to the domain decomposition paradigm for the conductor-like screening model”, Int. J. Quantum Chem., 119, e25669:1–15 (2019).
  * 241 M. Nottoli, B. Stamm, G. Scalmani, and F. Lipparini, “Quantum calculations in solution of energies, structures, and properties with a domain decomposition polarizable continuum model”, J. Chem. Theory Comput., 15, 6061–6073 (2019).
  * 242 L. Lagardère, L.-H. Jolly, F. Lipparini, F. Aviat, B. Stamm, Z. F. Jing, M. Harger, H. Torabifard, G. A. Cisneros, M. J. Schneiders, N. Gresh, Y. Maday, P. Y. Ren, J. W. Ponder, and J.-P. Piquemal, “Tinker-HP: A massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields”, Chem. Sci., 9, 956–972 (2018).
  * 243 F. Lipparini, L. Lagardère, C. Raynaud, B. Stamm, E. Cancès, B. Mennucci, M. Schnieders, P. Ren, Y. Maday, and J.-P. Piquemal, “Polarizable molecular dynamics in a polarizable continuum solvent”, J. Chem. Theory Comput., 11, 623–634 (2014).
  * 244 C. J. Cramer and D. G. Truhlar, “A universal approach to solvation modeling”, Acc. Chem. Res., 41, 760–768 (2008).
  * 245 C. J. Cramer and D. G. Truhlar, “SM$x$ continuum models for condensed phases”, in Trends and Perspectives in Modern Computational Science, G. Maroulis and T. E. Simos, Eds., Vol. 6 of Lecture Series on Computer and Computational Sciences; Brill/VSP: Leiden, 2006;  pages 112–140.
  * 246 J. D. Thompson, C. J. Cramer, and D. G. Truhlar, “New universal solvation model and comparison of the accuracy of the SM5.42R, SM5.43R, C-PCM, D-PCM, and IEF-PCM continuum solvation models for aqueous and organic solvation free energies and for vapor pressures”, J. Phys. Chem. A, 108, 6532–6542 (2004).
  * 247 C. P. Kelly, C. J. Cramer, and D. G. Truhlar, “SM6: A density functional theory continuum solvation model for calculating aqueous solvation free energies of neutrals, ions, and solute–water clusters”, J. Chem. Theory Comput., 1, 1133–1152 (2005).
  * 248 A. V. Marenich, C. J. Cramer, and D. G. Truhlar, “Generalized Born solvation model SM12”, J. Chem. Theory Comput., 9, 609–620 (2013).
  * 249 A. V. Marenich, C. J. Cramer, and D. G. Truhlar, “Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions”, J. Phys. Chem. B, 113, 6378–6396 (2009).
  * 250 A. V. Marenich, C. J. Cramer, and D. G. Truhlar, “Performance of SM6, SM8, and SMD on the SAMPL1 test set for the prediction of small-molecule solvation free energies”, J. Phys. Chem. B, 113, 4538–4543 (2009).
  * 251 Z.-Q. You and J. M. Herbert, “Reparameterization of an accurate, few-parameter implicit solvation model for quantum chemistry: Composite method for implicit representation of solvent, CMIRS v. 1.1”, J. Chem. Theory Comput., 12, 4338–4346 (2016).
  * 252 C. Dupont, O. Andreussi, and N. Marzari, “Self-consistent continuum solvation (SCCS): The case of charged systems”, J. Chem. Phys., 139, 214110:1–8 (2013).
  * 253 A. Klamt, B. Mennucci, J. Tomasi, V. Barone, C. Curutchet, M. Orozco, and F. J. Luque, “On the performance of continuum solvation methods. A comment on “Universal approaches to solvation modeling””, Acc. Chem. Res., 42, 489–492 (2009).
  * 254 L. Tomaník, E. Muchová, and P. Slavíček, “Solvation energies of ions with ensemble cluster-continuum approach”, Phys. Chem. Chem. Phys., 22, 22357–22368 (2020).
  * 255 R. Casasnovas, J. Ortega-Castro, J. Frau, J. Donoso, and F. Muñoz, “Theoretical p${K}_{\text{a}}$ calculations with continuum model solvents, alternative protocols to thermodynamic cycles”, Int. J. Quantum Chem., 114, 1350–1363 (2014).
  * 256 B. Thapa and H. B. Schlegel, “Calculations of p${K}_{\text{a}}$’s and redox potentials of nucleobases with explicit waters and polarizable continuum solvation”, J. Phys. Chem. A, 119, 5134–5144 (2015).
  * 257 B. Thapa and H. B. Schlegel, “Density functional theory calculation of p$K_{\text{a}}$’s of thiols in aqueous solution using explicit water molecules and the polarizable continuum model”, J. Phys. Chem. A, 120, 5726–5735 (2016).
  * 258 V. S. Bryantsev, M. S. Diallo, and W. A. Goddard III, “Calculation of solvation free energies of charged solutes using mixed cluster/continuum models”, J. Phys. Chem. B, 112, 9709–9719 (2008).
  * 259 D. Riccardi, H.-B. Guo, J. M. Parks, B. Gu, L. Liang, and J. C. Smith, “Cluster-continuum calculations of hydration free energies of anions and group 12 divalent cations”, J. Chem. Theory Comput., 9, 555–569 (2013).
  * 260 S. Dhillon and A. L. East, “Challenges in predicting $\Delta_{\text{rxn}}G$ in solution: Hydronium, hydroxide, and water autoionization”, Int. J. Quantum Chem., 118, e25703:1–7 (2018).
  * 261 D. H. Patel and A. L. L. East, “Semicontinuum (cluster-continuum) modeling of acid-catalyzed aqueous reactions: Alkene hydration”, J. Phys. Chem. A, 124, 9088–9104 (2020).
  * 262 I. Soteras, C. Curutchet, A. Bidon-Chanal, M. Orozco, and F. J. Luque, “Extension of the MST model to the IEF formalism: HF and B3LYP parameterizations”, J. Mol. Struct. (Theochem), 727, 29–40 (2005).
  * 263 C. Curutchet, M. Orozco, F. J. Luque, B. Mennucci, and J. Tomasi, “Dispersion and repulsion contributions to the solvation free energy: Comparison of quantum mechanical and classical approaches in the polarizable continuum model”, J. Comput. Chem., 27, 1769–1780 (2006).
  * 264 C. J. Cramer and D. G. Truhlar, “Reply to comment on “A universal approach to solvation modeling””, Acc. Chem. Res., 42, 493–497 (2009).
  * 265 A. V. Marenich, R. M. Olson, C. P. Kelly, C. J. Cramer, and D. G. Truhlar, “Self-consistent reaction field model for aqueous and nonaqueous solutions based on accurate polarized partial charges”, J. Chem. Theory Comput., 3, 2011–2033 (2007).
  * 266 A. C. Chamberlin, C. J. Cramer, and D. G. Truhlar, “Performance of SM8 on a test to predict small-molecule solvation free energies”, J. Phys. Chem. B, 112, 8651–8655 (2008).
  * 267 A. Klamt, “The COSMO and COSMO-RS solvation models”, WIREs Comput. Mol. Sci., 1, 699–709 (2011).
  * 268 A. Klamt, “The COSMO and COSMO-RS solvation models”, WIREs Comput. Mol. Sci., 8, e1338:1–11 (2018).
  * 269 D. Bashford and D. A. Case, “Generalized Born models of macromolecular solvation effects”, Annu. Rev. Phys. Chem., 51, 129–152 (2000).