Clearly, the Kirkwood-Onsager model is inappropriate if the solute is very
non-spherical. Nowadays, a more general class of “apparent surface charge”
SCRF solvation models are much more popular, to the extent that the generic
term “polarizable continuum model” (PCM) is typically used to denote these
methods.^{Tomasi:2005} Apparent surface charge PCMs improve upon the
Kirkwood-Onsager model in two ways. Most importantly, they provide a much more
realistic description of molecular shape, typically by constructing the
“solute cavity” (*i.e.*, the interface between the atomistic region and the
dielectric continuum) from a union of atom-centered spheres, an aspect of the
model that is discussed in Section 11.2.2.2. In addition,
the exact electron density of the solute (rather than a multipole expansion) is
used to polarize the continuum. Electrostatic interactions between the solute
and the continuum manifest as an induced charge density on the cavity surface,
which is discretized into point charges for practical calculations. The
surface charges are determined based upon the solute’s electrostatic potential
at the cavity surface, hence the surface charges and the solute wave function
must be determined self-consistently.

The PCM literature has a long history^{Tomasi:2005} and there are several
different models in widespread use; connections between these models have not
always been
appreciated.^{Chipman:2000, Cances:2001a, Chipman:2002a, Lange:2011b}
Chipman^{Chipman:2000, Chipman:2002a} has shown how various PCMs can be
formulated within a common theoretical framework; see
Ref. Herbert:2016a for a pedagogical introduction. The PCM takes
the form of a set of linear equations,

$$\mathrm{\mathbf{K}\mathbf{q}}=\mathrm{\mathbf{R}\mathbf{v}},$$ | (11.2) |

in which the induced charges ${q}_{i}$ at the cavity surface discretization points (organized into a vector $\mathbf{q}$ in Eq. (11.2)) are computed from the values ${v}_{i}$ of the solute’s electrostatic potential at those same discretization points. The form of the matrices $\mathbf{K}$ and $\mathbf{R}$ depends upon the particular PCM in question. These matrices are given in Table 11.3 for the PCMs that are available in Q-Chem.

Model | Literature | Matrix $\mathbf{K}$ | Matrix $\mathbf{R}$ | Scalar ${f}_{\epsilon}$ |
---|---|---|---|---|

Refs. | ||||

COSMO | Klamt:1993 | $\mathbf{S}$ | $-{f}_{\epsilon}\mathrm{\U0001d7cf}$ | $(\epsilon -1)/(\epsilon +1/2)$ |

C-PCM | Truong:1995, Barone:1998 | $\mathbf{S}$ | $-{f}_{\epsilon}\mathrm{\U0001d7cf}$ | $(\epsilon -1)/\epsilon $ |

IEF-PCM | Chipman:2000, Cances:2001a | $\mathbf{S}-({f}_{\epsilon}/2\pi )\mathrm{\mathbf{D}\mathbf{A}\mathbf{S}}$ | $-{f}_{\epsilon}\left(\mathrm{\U0001d7cf}-\frac{1}{2\pi}\mathrm{\mathbf{D}\mathbf{A}}\right)$ | $(\epsilon -1)/(\epsilon +1)$ |

SS(V)PE | Chipman:2000, Chipman:2002b | $\mathbf{S}-({f}_{\epsilon}/4\pi )\left(\mathrm{\mathbf{D}\mathbf{A}\mathbf{S}}+{\mathrm{\mathbf{S}\mathbf{A}\mathbf{D}}}^{\u2020}\right)$ | $-{f}_{\epsilon}\left(\mathrm{\U0001d7cf}-\frac{1}{2\pi}\mathrm{\mathbf{D}\mathbf{A}}\right)$ | $(\epsilon -1)/(\epsilon +1)$ |

The oldest PCM is the so-called D-PCM model of Tomasi and
coworkers,^{Miertus:1981} but unlike the models listed in
Table 11.3, D-PCM requires explicit evaluation of the electric
field normal to the cavity surface, This is undesirable, as evaluation of the
electric field is both more expensive and more prone to numerical problems as
compared to evaluation of the electrostatic potential. Moreover, the
dependence on the electric field can be formally eliminated at the level of the
integral equation whose discretized form is given in
Eq. (11.2).^{Chipman:2000} As such, D-PCM is essentially
obsolete, and the PCMs available in Q-Chem require only the evaluation of the
electrostatic potential, not the electric field.

The simplest PCM that continues to enjoy widespread use is the Conductor-Like
Screening Model (COSMO) introduced by Klamt and
Schüürmann.^{Klamt:1993} Truong and Stefanovich^{Truong:1995}
later implemented the same model with a slightly different dielectric scaling
factor (${f}_{\epsilon}$ in Table 11.3), and called this
modification GCOSMO. The latter was implemented within the PCM formalism by
Barone and Cossi *et al.*,^{Barone:1998, Cossi:2003} who called the model
C-PCM (for “conductor-like” PCM). In each case, the dielectric
screening factor has the form

$${f}_{\epsilon}=\frac{\epsilon -1}{\epsilon +x},$$ | (11.3) |

where Klamt and Schüürmann proposed $x=1/2$ but $x=0$ was used in GCOSMO
and C-PCM. The latter value is the correct choice for a single charge in a
spherical cavity (*i.e.*, the Born ion model), although Klamt and coworkers
suggest that $x=1/2$ is a better compromise, given that the Kirkwood-Onsager
analytical result is $x=\mathrm{\ell}/(\mathrm{\ell}+1)$ for an $\mathrm{\ell}$th-order multipole centered
in a spherical cavity.^{Klamt:1993, Baldridge:1997} The distinction is
irrelevant in high-dielectric solvents; the $x=0$ and $x=1/2$ values of
${f}_{\epsilon}$ differ by only 0.6% for water at 25${}^{\circ}$C, for example.
Truong^{Truong:1995} argues that $x=0$ does a better job of preserving
Gauss’ Law in low-dielectric solvents, but more accurate solvation energies (at
least for neutral molecules, as compared to experiment) are sometimes obtained
using $x=1/2$ (Ref. Barone:1998). This result is likely highly
sensitive to cavity construction, and in any case, both versions are available in Q-Chem.

Whereas the original COSMO model introduced by Klamt and
Schüürmann^{Klamt:1993} corresponds to Eq. (11.2) with
$\mathbf{K}$ and $\mathbf{R}$ as defined in Table 11.3, Klamt
and coworkers later introduced a correction for outlying charge that goes
beyond Eq. (11.2).^{Klamt:1996, Baldridge:1997} Klamt now
consistently refers to this updated model as “COSMO”,^{Klamt:2011} and
we shall adopt this nomenclature as well. COSMO, with the outlying charge
correction, is available in Q-Chem and is described in
Section 11.2.7. In contrast, C-PCM consists entirely of
Eq. (11.2) with matrices $\mathbf{K}$ and $\mathbf{R}$ as defined in
Table 11.3, although it is possible to modify the dielectric
screening factor to use the $x=1/2$ value (as in COSMO) rather than the $x=0$
value. Additional non-electrostatic terms can be added at the user’s
discretion, as discussed below, but there is no explicit outlying charge
correction in C-PCM. These and other fine-tuning details for PCM jobs are
controllable via the *$pcm* input section that is described in
Section 11.2.3.

As compared to C-PCM, a more sophisticated treatment of continuum electrostatic
interactions is afforded by the “surface and simulation of volume polarization
for electrostatics” [SS(V)PE] approach.^{Chipman:2000} Formally speaking,
this model provides an exact treatment of the surface polarization
(*i.e.*, the surface charge induced by the solute charge that is contained within
the solute cavity, which induces a surface polarization owing to the
discontinuous change in dielectric constant across the cavity boundary) but
also an approximate treatment of the volume polarization (arising from
the aforementioned outlying charge). The “SS(V)PE” terminology is Chipman’s
notation,^{Chipman:2000} but this model is formally equivalent, at the
level of integral equations, to the “integral equation formalism” (IEF-PCM)
that was developed originally by Cancès *et al.*.^{Cances:1997, Tomasi:1999}
Some difference do arise when the integral equations are discretized to form
finite-dimensional matrix equations,^{Lange:2011b} and it should be noted
from Table 11.3 that SS(V)PE uses a symmetrized form of the
$\mathbf{K}$ matrix as compared to IEF-PCM. The *asymmetric* IEF-PCM is
the recommended approach,^{Lange:2011b} although only the symmetrized
version is available in the isodensity implementation of SS(V)PE that is
discussed in Section 11.2.5. As with the obsolete D-PCM
approach, the original version of IEF-PCM explicitly required evaluation of the
normal electric field at the cavity surface, but it was later shown that this
dependence could be eliminated to afford the version described in
Table 11.3.^{Chipman:2000, Cances:2001a} This version
requires only the electrostatic potential, and is thus preferred, and it is
this version that we designate as IEF-PCM. The C-PCM model becomes equivalent
to SS(V)PE in the limit
$\epsilon \to \mathrm{\infty}$,^{Chipman:2000, Lange:2011b} which means that
C-PCM must somehow include an implicit correction for volume
polarization, even if this was not by design.^{Klamt:1996} For
$\epsilon \gtrsim 50$, numerical calculations reveal that there is
essentially no difference between SS(V)PE and C-PCM results.^{Lange:2011b}
Since C-PCM is less computationally involved as compared to SS(V)PE, it is the
PCM of choice in high-dielectric solvents. The computational savings relative
to SS(V)PE may be particularly significant for large QM/MM/PCM
jobs.
For a more detailed discussion of the history of these models, see the lengthy
and comprehensive review by Tomasi *et al.*.^{Tomasi:2005} For a
briefer discussion of the connections between these models, see
Refs. Chipman:2002a, Lange:2011b, Herbert:2016a.

Construction of the cavity surface is a crucial aspect of PCMs, as computed
properties are quite sensitive to the details of the cavity construction. Most
cavity constructions are based on a union of atom-centered spheres (see
Fig. 11.1), but there are yet several different constructions
whose nomenclature is occasionally confused in the literature. Simplest and
most common is the van der Waals (vdW) surface consisting of a union of
atom-centered spheres. Traditionally,^{Bonaccorsi:1984, Tomasi:1994} and
by default in Q-Chem, the atomic radii are taken to be 1.2 times larger than
vdW radii extracted from crystallographic data, originally by
Bondi (and thus sometimes called “Bondi radii”).^{Bondi:1964} This 20%
augmentation is intended to mimic the fact that solvent molecules cannot
approach all the way to the vdW radius of the solute atoms, though it’s not
altogether clear that this is an optimal value. (The default scaling factor in
Q-Chem is 1.2 but can be modified by the user.) An alternative to scaling
the atomic radii is to add a certain fixed increment to each, representing the
approximate size of a solvent molecule (*e.g.*, 1.4 Å for water) and leading
to what is known as the solvent accessible surface (SAS). From
another point of view, the SAS represents the surface defined by the center of
a spherical solvent molecule as it rolls over the vdW surface, as suggested in
Fig. 11.1. Both the vdW surface and the SAS possess cusps
where the atomic spheres intersect, although these become less pronounced as
the atomic radii are scaled or augmented. These cusps are eliminated in what
is known as the solvent-accessible surface (SES), sometimes called the
Connolly surface or the “molecular surface". The SES uses the surface of the
probe sphere at points where it is simultaneously tangent to two or more atomic
spheres to define elements of a “re-entrant surface” that smoothly connects
the atomic (or “contact”) surface.^{Lange:2020}

Having chosen a model for the cavity surface, this surface is discretized using
atom-centered Lebedev grids^{Lebedev:1976, Lebedev:1977, Lebedev:1999} of
the same sort that are used to perform the numerical integrations in DFT.
(Discretization of the re-entrant facets of the SES is somewhat more
complicated but similar in spirit.^{Lange:2020}) Surface charges ${q}_{i}$
are located at these grid points and the Lebedev quadrature weights can be used
to define the surface area associated with each discretization point.^{Lange:2010a}

A long-standing (though not well-publicized) problem with the aforementioned
discretization procedure is that it fails to afford continuous potential
energy surfaces as the solute atoms are displaced, because certain surface grid
points may emerge from, or disappear within, the solute cavity, as the atomic
spheres that define the cavity are moved. This undesirable behavior can
inhibit convergence of geometry optimizations and, in certain cases, lead to
very large errors in vibrational frequency calculations.^{Lange:2010a} It
is also a fundamental hindrance to molecular dynamics
calculations.^{Lange:2010b} Building upon earlier work by York and
Karplus,^{York:1999} Lange and
Herbert^{Lange:2010a, Lange:2010b, Lange:2020} developed a general scheme
for implementing apparent surface charge PCMs in a manner that affords smooth
potential energy surfaces, even for *ab initio* molecular dynamics
simulations involving bond breaking.^{Lange:2010b, Herbert:2016a} Notably,
this approach is faithful to the properties of the underlying integral
equation theory on which the PCMs are based, in the sense that the smoothing
procedure does not significantly perturb solvation energies or cavity surface
areas.^{Lange:2010a, Lange:2010b} The smooth discretization procedure
combines a switching function with Gaussian blurring of the cavity surface
charge density, and is thus known as the “Switching/Gaussian” (SWIG)
implementation of the PCM.

Both single-point energies and analytic energy gradients are available for SWIG
PCMs, when the solute is described using molecular mechanics or an SCF
(Hartree-Fock or DFT) electronic structure model, except that for the SES
cavity model only single-point energies are available. Analytic Hessians are
available for the C-PCM model only. (As usual, vibrational frequencies for
other models will be computed, if requested, by finite difference of analytic
energy gradients.) Single-point energy calculations using correlated wave
functions can be performed in conjunction with these solvent models, in which
case the correlated wave function calculation will use Hartree-Fock molecular
orbitals that are polarized in the presence of the continuum dielectric solvent
(*i.e.*, there is no post-Hartree–Fock PCM correction).

Researchers who use these PCMs are asked to cite Refs. Lange:2010b, Lange:2011b, which provide the details of Q-Chem’s implementation, and Ref. Lange:2020 if the SES is used. (We point the reader in particular to Ref. Lange:2010b, which provides an assessment of the discretization errors that can be anticipated using various PCMs and Lebedev grids; default grid values in Q-Chem were established based on these tests.) When publishing results based on PCM calculations, it is essential to specify both the precise model that is used (see Table 11.3) as well as how the cavity was constructed.

For example, “Bondi radii multiplied by 1.2”, which is the Q-Chem default,
except for hydrogen, where the factor is reduced to 1.1,^{Rowland:1996} as
per usual. Radii for main-group elements that were not provided by Bondi are
taken from Ref. Mantina:2009. Absent details such as these,
PCM calculations will be difficult to reproduce in other electronic structure programs.

In vertical excitation or ionization, the solute undergoes a sudden change in
its charge distribution. Various microscopic motions of the solvent have
characteristic times to reach certain polarization response, and fast part of
the solvent response (electrons) can follow such a dynamic process while the
remaining degrees of freedom (nuclei) remain unchanged as in the initial state.
Such splitting of the solvent response gives rise to nonequilibrium solvation.
In the literature, two different approaches have been developed for describing
nonequilibrium solvent effects: the linear response (LR)
approach^{Cammi:1999, Cossi:2001} and the state-specific (SS)
approach.^{Tomasi:1994, Cammi:1995, Cossi:2000, Improta:2006} Both are
implemented in Q-Chem,^{You:2015},at the SCF level for vertical
ionization and at the corresponding level (CIS, TDDFT or ADC, see
Section 7.10.7) for vertical excitation. A brief introduction to
these methods is given below, and users of the nonequilibrium PCM features are
asked to cite Refs. You:2015 and Mewes:2015a.
State-specific solvent-field equilibration for long-lived excited states to
compute e.g. emission energies is implemented for the ADC-suite of methods as
described in section 7.10.7. Users of this equilibrium-solvation
PCM please cite and be referred to Ref. Mewes:2017.

The LR approach considers the solvation effects as a coupling between a pair of
transitions, one for solute and the other for solvent. The transition
frequencies when the interaction between the solute and solvent is turned on
may be determined by considering such an interaction as a perturbation. In the
framework of TDDFT, the solvent/solute interaction is given by^{Hsu:2001}

$$\begin{array}{cc}\hfill {\omega}^{\prime}=& \int \mathit{d}\mathbf{r}\int \mathit{d}{\mathbf{r}}^{\prime}\int \mathit{d}{\mathbf{r}}^{\prime \prime}\int \mathit{d}{\mathbf{r}}^{\prime \prime \prime}{\rho}^{\mathrm{tr}*}(\mathbf{r})\left(\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime}|}+{g}_{\mathrm{XC}}(\mathbf{r},{\mathbf{r}}^{\prime})\right)\hfill \\ & \times {\chi}^{*}({\mathbf{r}}^{\prime},{\mathbf{r}}^{\prime \prime},\omega )\left(\frac{1}{|{\mathbf{r}}^{\prime \prime}-{\mathbf{r}}^{\prime \prime \prime}|}+{g}_{\mathrm{XC}}({\mathbf{r}}^{\prime \prime},{\mathbf{r}}^{\prime \prime \prime})\right){\rho}^{\mathrm{tr}}({\mathbf{r}}^{\prime \prime \prime}),\hfill \end{array}$$ | (11.4) |

where $\chi $ is the charge density response function of the solvent and ${\rho}^{\mathrm{tr}}(\mathbf{r})$ is the solute’s transition density. This term accounts for a dynamical correction to the transition energy so that it is related to the response of the solvent to the charge density of the solute oscillating at the solute transition frequency ($\omega $). Within a PCM, only classical Coulomb interactions are taken into account, and Eq. (11.4) becomes

$$\begin{array}{cc}\hfill {\omega}_{\mathrm{PCM}}^{\prime}=& \int \mathit{d}\mathbf{r}\int \mathit{d}\mathbf{s}\frac{{\rho}^{\mathrm{tr}*}(\mathbf{r})}{|\mathbf{r}-\mathbf{s}|}\int \mathit{d}{\mathbf{s}}^{\prime}\int \mathit{d}{\mathbf{r}}^{\prime}\mathcal{Q}(\mathbf{s},{\mathbf{s}}^{\prime},\epsilon )\frac{{\rho}^{\mathrm{tr}}({\mathbf{r}}^{\prime})}{|{\mathbf{s}}^{\prime}-{\mathbf{r}}^{\prime}|},\hfill \end{array}$$ | (11.5) |

where $\mathcal{Q}$ is PCM solvent response operator for a generic dielectric constant, $\epsilon $. The integral of $\mathcal{Q}$ and the potential of the density ${\rho}^{\mathrm{tr}}$ gives the surface charge density for the solvent polarization.

The state-specific (SS) approach takes into account the capability of a part of the solvent degrees of freedom to respond instantaneously to changes in the solute wave function upon excitation. Such an effect is not accounted for in the LR approach. In SS, a generic solvated-solute excited state ${\mathrm{\Psi}}_{i}$ is obtained as a solution of a nonlinear Schrödinger equation

$$\left({\widehat{H}}^{\mathrm{vac}}+{\widehat{V}}_{0}^{\mathrm{slow}}+{\widehat{V}}_{i}^{\mathrm{fast}}\right)|{\mathrm{\Psi}}_{i}\u27e9={E}_{i}^{\mathrm{SS}}|{\mathrm{\Psi}}_{i}\u27e9$$ | (11.6) |

that depends upon the solute’s charge distribution. Here ${\widehat{H}}^{\mathrm{vac}}$ is
the usual Hamiltonian for the solute in vacuum and the reaction field operator
${\widehat{V}}_{i}$ generates the electrostatic potential of the apparent surface
charge density (Section 11.2.2.1), corresponding to
slow and fast polarization response.
The solute is polarized self-consistently with respect to the solvent’s
reaction field. In case of vertical ionization rather than excitation, both the
ionized and non-ionized states can be treated within a ground-state formalism.
For vertical excitations, self-consistent SS models have been developed for
various excited-state methods,^{Improta:2006, Marenich:2011} including both
CIS and TDDFT.

In a linear dielectric medium, the solvent polarization is governed by the
electric susceptibility, $\chi =[\epsilon (\omega )-1]/4\pi $, where
$\epsilon (\omega )$ is the frequency-dependent permittivity.
In case of very fast vertical transitions, the dielectric response is ruled by
the optical dielectric constant, ${\epsilon}_{\mathrm{opt}}={n}^{2}$, where $n$ is
the solvent’s index of refraction. In both LR and SS, the fast part of the
solvent’s degrees of freedom is in equilibrium with the solute density change.
Within PCM, the fast solvent polarization charges for the SS excited state $i$
can be obtained by solving the following equation:^{Cossi:2000}

$${\mathbf{K}}_{{\epsilon}_{\mathrm{opt}}}{\mathbf{q}}_{i}^{\mathrm{fast},\mathrm{SS}}={\mathbf{R}}_{{\epsilon}_{\mathrm{opt}}}\left[{\mathbf{v}}_{i}+\mathbf{v}({\mathbf{q}}_{0}^{\mathrm{slow}})\right].$$ | (11.7) |

Here ${\mathbf{q}}^{\mathrm{fast},\mathrm{SS}}$ is the discretized fast surface charge. The dielectric constants in the matrices $\mathbf{K}$ and $\mathbf{R}$ (Section 11.2.2.1) are replaced with the optical dielectric constant, and ${\mathbf{v}}_{i}$ is the potential of the solute’s excited state density, ${\rho}_{i}$. The quantity $\mathbf{v}({\mathbf{q}}_{0}^{\mathrm{slow}})$ is the potential of the slow part of the apparent surface charges in the ground state, which are given by

$${\mathbf{q}}_{0}^{\mathrm{slow}}=\left(\frac{\epsilon -{\epsilon}_{\mathrm{opt}}}{\epsilon -1}\right){\mathbf{q}}_{0}.$$ | (11.8) |

For LR-PCM, the solvent polarization is subjected to the first-order changes to the electron density (TDDFT linear density response), and thus Eq. (11.7) becomes

$${\mathbf{K}}_{{\epsilon}_{\mathrm{opt}}}{\mathbf{q}}_{i}^{\mathrm{fast},\mathrm{LR}}={\mathbf{R}}_{{\epsilon}_{\mathrm{opt}}}\mathbf{v}({\rho}_{i}^{\mathrm{tr}}).$$ | (11.9) |

The LR approach for CIS/TDDFT excitations and the self-consistent SS
method (using the ground-state SCF) for vertical ionizations are available in
Q-Chem. The self-consistent SS method for vertical excitations is not
available, because this method is problematic in the vicinity of (near-)
degeneracies between excited states, such as in the vicinity of a conical
intersection. The fundamental problem in the SS approach is that each wave
function ${\mathrm{\Psi}}_{i}$ is an eigenfunction of a different Hamiltonian, since
Eq. (11.6) depend upon the specific state of interest. To avoid the
ordering and the non-orthogonality problems, we compute the vertical
excitation energy using a first-order, perturbative approximation to the SS
approach,^{Cammi:2005, Caricato:2006} in what we have termed the “ptSS”
method.^{Mewes:2015a} The zeroth-order excited-state wave function can be
calculated using various excited-state methods (currently available for CIS and
TDDFT in Q-Chem) with solvent-relaxed molecular orbitals obtained from a
ground-state PCM calculation. As mentioned previously, LR and SS describe
different solvent relaxation features in nonequilibrium solvation. In the
perturbation scheme, we can calculate the LR contribution using the
zeroth-order transition density, in what we have called the "ptLR" approach.
The combination of ptSS and ptLR yields quantitatively good solvatochromatic
shifts in combination with TDDFT but not with the correlated variants of ADC,
for which the pure ptSS approach was shown to be superior.^{You:2015, Mewes:2015a}

The LR and SS approaches can also be used in the study of photon emission
processes.^{Improta:2007} An emission process can be treated as a vertical
excitation at a stationary point on the excited-state potential surface. The
basic requirement therefore is to prepare the solvent-relaxed geometry for the
excited-state of interest. TDDFT/C-PCM analytic gradients and Hessian are available.

Section 7.3.5 for computational details regarding
excited-state geometry optimization with PCM. An emission process is slightly
more complicated than the absorption case. Two scenarios are discussed in
literature, depending on the lifetime of an excited state in question. In the limiting case
of ultra-fast excited state decay, when only fast solvent degrees of freedom are
expected to be equilibrated with the excited-state density. In this limit, the
emission energy can be computed exactly in the same way as the vertical
excitation energy. In this case, excited state geometry optimization should be
performed in the nonequilibrium limit.
The other limit is that of long-lived excited state, *e.g.*, strongly
fluorescent species and phosphorescence. In the long-lived case, excited state
geometry optimization should be performed with the solvent equilibrium limit.
Thus, the excited state should be computed using an equilibrium LR or SS
approach, and the ground state is calculated using nonequilibrium
self-consistent SS approach. The latter approach is implemented for the
ADC-based methods as described in Section 7.10.7.