11.2 Chemical Solvent Models

11.2.2 Polarizable Continuum Models

(June 30, 2021)

The multipolar expansion model is based on exact formulas for the solvation energy of a point multipole in a spherical cavity, 760 Mikkelsen K. V. et al.
J. Chem. Phys.
(1988), 89, pp. 3086.
Link
which is a crude approximation except (or perhaps even) for small molecules, and the Kirkwood-Onsager model has been largely superseded by the more general class of “apparent surface charge” SCRF solvation models, typically known as PCMs. 1093 Tomasi J., Mennucci B., Cammi R.
Chem. Rev.
(2005), 106, pp. 2999.
Link
These models improve upon the multipolar expansion method 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.

11.2.2.1 Formal Theory and Discussion of Different Models

The PCM literature has a long history 1093 Tomasi J., Mennucci B., Cammi R.
Chem. Rev.
(2005), 106, pp. 2999.
Link
and there are several different models in widespread use; connections between these models have not always been appreciated. 199 Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558.
Link
, 153 Cancès E., Mennucci B.
J. Chem. Phys.
(2001), 114, pp. 4744.
Link
, 200 Chipman D. M.
Theor. Chem. Acc.
(2002), 107, pp. 80.
Link
, 609 Lange A. W., Herbert J. M.
Chem. Phys. Lett.
(2011), 509, pp. 77.
Link
Chipman 199 Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558.
Link
, 200 Chipman D. M.
Theor. Chem. Acc.
(2002), 107, pp. 80.
Link
has shown how various PCMs can be formulated within a common theoretical framework; see Ref. Herbert:2021b for a review. The PCM takes the form of a set of linear equations,

𝐊𝐪=𝐑𝐯, (11.2)

in which the induced charges qi at the cavity surface discretization points [organized into a vector 𝐪 in Eq. (11.2)] are computed from the values vi of the solute’s electrostatic potential at those same discretization points. The form of the matrices 𝐊 and 𝐑 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 𝐊 Matrix 𝐑 Scalar fε
Refs.
COSMOa 556 𝐒 -fε𝟏 (ε-1)/(ε+1/2)
C-PCM 1105, 61 𝐒 -fε𝟏 (ε-1)/ε
IEF-PCM 199, 153 𝐒-(fε/2π)𝐃𝐀𝐒 -fε(𝟏-12π𝐃𝐀) (ε-1)/(ε+1)
SS(V)PE 199 𝐒-(fε/4π)(𝐃𝐀𝐒+𝐒𝐀𝐃) -fε(𝟏-12π𝐃𝐀) (ε-1)/(ε+1)
aAlso includes a charge renormalization correction; see Section 11.2.7.
Table 11.3: Definition of the matrices in Eq. (11.2) for the various PCMs that are available in Q-Chem. The matrix 𝐒 consists of Coulomb interactions between the cavity charges and 𝐃 is the discretized version of the matrix that generates the outward-pointing normal electric field vector. (See Refs. 200, 196, 446 for detailed definitions.) The matrix 𝐀 is diagonal and contains the surface areas of the cavity discretization elements, and 𝟏 is a unit matrix.

The oldest PCM is the so-called D-PCM model of Tomasi and coworkers, 758 Miertuš S., Scrocco E., Tomasi J.
Chem. Phys.
(1981), 55, pp. 117.
Link
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). 199 Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558.
Link
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 model, C-PCM. 61 Barone V., Cossi M.
J. Phys. Chem. A
(1998), 102, pp. 1995.
Link
, 226 Cossi M. et al.
J. Comput. Chem.
(2003), 24, pp. 669.
Link
Originally derived by Klamt and Schüürmann based on arguments invoking the conductor limit (ε), this model can also be derived as an approximation to more formally correct models. 198 Chipman D. M.
J. Chem. Phys.
(1999), 110, pp. 8012.
Link
, 608 Lange A. W., Herbert J. M.
J. Chem. Phys.
(2011), 134, pp. 204110.
Link
Over the years, the dielectric-dependent factor

fε=ε-1ε+x (11.3)

that appears in this model (see Table 11.3) has been used with different values of x. The value x=0 is typically used in C-PCM calculations and x=1/2 in COSMO calculations, although Klamt and co-workers later suggested using x=1/2 for neutral solutes and x=0 for ions. 555 Klamt A., Moya C., Palomar J.
J. Chem. Theory Comput.
(2015), 11, pp. 4220.
Link
The specific choice of fε is controllable via the $pcm input section that is described in Section 11.2.3.

Whereas from Table 11.3 the C-PCM and COSMO methods would appear to be the same up to a minor rescaling of the surface charge (i.e., up to the precise choice of fε), historically the term “COSMO" has been used by Klamt to mean a particular “dual-cavity” implementation of this model that makes it different from other PCMs. This construction is equivalent to the “outlying charge correction” that is discussed in Section 11.2.7, and was intended to account for the effects of the tail of the solute’s charge density that penetrates beyond the cavity surface. 554 Klamt A., Jonas V.
J. Chem. Phys.
(1996), 105, pp. 9972.
Link
Subsequent work cast considerable doubt on the theoretical justification for this correction, since the C-PCM/COSMO ansatz was shown to include already an implicit correction for outlying charge. 198 Chipman D. M.
J. Chem. Phys.
(1999), 110, pp. 8012.
Link
Further discussion of this construction and of COSMO is deferred to Section 11.2.7,

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. 199 Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558.
Link
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, 199 Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558.
Link
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.. 152 Cancès E., Mennucci B., Tomasi J.
J. Chem. Phys.
(1997), 107, pp. 3032.
Link
, 1094 Tomasi J., Mennucci B., Cancès E.
J. Mol. Struct. (Theochem)
(1999), 464, pp. 211.
Link
Some difference do arise when the integral equations are discretized to form finite-dimensional matrix equations, 609 Lange A. W., Herbert J. M.
Chem. Phys. Lett.
(2011), 509, pp. 77.
Link
and it should be noted from Table 11.3 that SS(V)PE uses a symmetrized form of the 𝐊 matrix as compared to IEF-PCM. The asymmetric IEF-PCM is the recommended approach, 609 Lange A. W., Herbert J. M.
Chem. Phys. Lett.
(2011), 509, pp. 77.
Link
although only the symmetrized version is available in the isodensity implementation of SS(V)PE that is discussed in Section 11.2.5. That said, differences between symmetry and asymmetric versions are only important in the case of van der Waals cavity surfaces; they are insignificant for the isodensity cavity construction.

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. 199 Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558.
Link
, 153 Cancès E., Mennucci B.
J. Chem. Phys.
(2001), 114, pp. 4744.
Link
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 ε, 199 Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558.
Link
, 609 Lange A. W., Herbert J. M.
Chem. Phys. Lett.
(2011), 509, pp. 77.
Link
which means that C-PCM must somehow include an implicit correction for volume polarization, even if this was not by design. 554 Klamt A., Jonas V.
J. Chem. Phys.
(1996), 105, pp. 9972.
Link
For ε50, numerical calculations reveal that there is essentially no difference between SS(V)PE and C-PCM results. 609 Lange A. W., Herbert J. M.
Chem. Phys. Lett.
(2011), 509, pp. 77.
Link
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.. 1093 Tomasi J., Mennucci B., Cammi R.
Chem. Rev.
(2005), 106, pp. 2999.
Link
For a briefer discussion of the connections between these models, see Refs. 200, 609, 446.

11.2.2.2 Cavity Construction and Discretization


Illustration of various solute cavity surface definitions for
PCMs.
Figure 11.1: Illustration of various solute cavity surface definitions for PCMs. 604 Lange A. W. et al.
Mol. Phys.
(2020), 118, pp. e1644384.
Link
The union of atomic van der Waals spheres (shown in gray) defines the van der Waals (vdW) surface, in black. Note that actual vdW radii from the literature are sometimes scaled in constructing the vdW surface. If a probe sphere (representing the assumed size of a solvent molecule) is rolled over the van der Waals surface, then its center point traces out the solvent accessible surface (SAS), shown in green; the SAS is equivalent to a vdW surface where the atomic radii are increases by the radius of the probe sphere. Finally, one can use the probe sphere to smooth out the sharp crevasses in the vdW surface using the re-entrant surface elements shown in red, resulting in the solvent-excluded surface (SES).

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. The radius for the sphere centered on atom A can be written in the form

RA=αvdWRvdW,A+Rprobe (11.4)

where RvdW,A is the vdW radius for atom A, taken for example from the set of vdW radii published by Bondi. 119 Bondi A.
J. Phys. Chem.
(1964), 68, pp. 441.
Link
Traditionally, the vdW radii RvdW,A that are extracted from crystallographic data are scaled by a factor αvdW=1.1–1.2. 118 Bonaccorsi R., Palla P., Tomasi J.
J. Am. Chem. Soc.
(1984), 106, pp. 1945.
Link
, 1095 Tomasi J., Persico M.
Chem. Rev.
(1994), 94, pp. 2027.
Link
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 the same value ought to be optimal in all cases. (The scaling factor defaults to αvdW=1.2 but can be modified by the user.)

An alternative to scaling the atomic radii is to add a certain fixed increment Rprobe to each, representing the approximate size of a solvent molecule, and leading to what is known as the solvent accessible surface (SAS). The choice Rprobe=1.4 Å is common for water and represents the approximate physical size of a water molecule, although values in the range Rprobe=0.2–0.5 Å often afford better solvation energies. In any case, if a nonzero value of Rprobe is used, then the scaling factor in Eq. (11.4) should be set to αvdW=1.0, since these two parameters are intended to model the same effect, namely, that a solvent molecule’s finite size prevents it from approaching all the way to the vdW radii of the solute.

Note from Fig. 11.1 that 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. 604 Lange A. W. et al.
Mol. Phys.
(2020), 118, pp. e1644384.
Link

Having chosen a model for the cavity surface, this surface is discretized using atom-centered Lebedev grids 627 Lebedev V. I.
Zh. Vychisl. Mat. Mat. Fix.
(1976), 16, pp. 293.
Link
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. 604 Lange A. W. et al.
Mol. Phys.
(2020), 118, pp. e1644384.
Link
) Surface charges qi are located at these grid points and the Lebedev quadrature weights can be used to define the surface area associated with each discretization point. 606 Lange A. W., Herbert J. M.
J. Phys. Chem. Lett.
(2010), 1, pp. 556.
Link

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. 606 Lange A. W., Herbert J. M.
J. Phys. Chem. Lett.
(2010), 1, pp. 556.
Link
It is also a fundamental hindrance to molecular dynamics calculations. 607 Lange A. W., Herbert J. M.
J. Chem. Phys.
(2010), 133, pp. 244111.
Link
Building upon earlier work by York and Karplus, 1223 York D. M., Karplus M.
J. Phys. Chem. A
(1999), 103, pp. 11060.
Link
Lange and Herbert 606 Lange A. W., Herbert J. M.
J. Phys. Chem. Lett.
(2010), 1, pp. 556.
Link
, 607 Lange A. W., Herbert J. M.
J. Chem. Phys.
(2010), 133, pp. 244111.
Link
, 604 Lange A. W. et al.
Mol. Phys.
(2020), 118, pp. e1644384.
Link
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.607, 446, Herbert:2021b 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. 606 Lange A. W., Herbert J. M.
J. Phys. Chem. Lett.
(2010), 1, pp. 556.
Link
, 607 Lange A. W., Herbert J. M.
J. Chem. Phys.
(2010), 133, pp. 244111.
Link
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). This represents a “zeroth-order” inclusion of solvent effects that captures the leading-order effect of continuum solvation on molecular properties. Given the crudeness of the model itself, more consistent inclusion of post-Hartree–Fock solvation effects is not expected to be important.

Researchers who use these PCMs are asked to cite Refs. 607 and 609, which provide the details of Q-Chem’s implementation, and Ref. 604 if the SES is used. We point the reader in particular to Refs. 607 and 1227, 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, and this should be done without resorting to software-specific keywords, the use of which has significantly muddled the literature on continuum electrostatics. 229 Cramer C. J., Truhlar D. G.
Acc. Chem. Res.
(2009), 42, pp. 493.
Link
As an example of good practice, the default cavity construction in Q-Chem is a vdW cavity using Bondi atomic radii, 119 Bondi A.
J. Phys. Chem.
(1964), 68, pp. 441.
Link
except that for hydrogen we use the modified radius of 1.1 Å, following a reassessment that judged Bondi’s original value of 1.2 Å for hydrogen to be too large. 949 Rowland R. S., Taylor R.
J. Phys. Chem.
(1996), 100, pp. 7384.
Link
Each of these radii RvdW,A in Eq. (11.4) is then scaled by a factor αvdW=1.2 for use in cavity construction. Radii for main-group elements that were not provided by Bondi are taken from Ref. 711. Absent details such as these, PCM calculations will be difficult to reproduce in other electronic structure programs.

11.2.2.3 Nonequilibrium Solvation for Vertical Excitation, Ionization and Emission

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 149 Cammi R., Mennucci B.
J. Chem. Phys.
(1999), 110, pp. 9877.
Link
, 224 Cossi M., Barone V.
J. Chem. Phys.
(2001), 115, pp. 4708.
Link
and the state-specific (SS) approach. 1095 Tomasi J., Persico M.
Chem. Rev.
(1994), 94, pp. 2027.
Link
, 150 Cammi R., Tomasi J.
Int. J. Quantum Chem. Symp.
(1995), 29, pp. 465.
Link
, 222 Cossi M., Barone V.
J. Phys. Chem. A
(2000), 104, pp. 10614.
Link
, 485 Improta R. et al.
J. Chem. Phys.
(2006), 125, pp. 054103.
Link
Both are implemented in Q-Chem, 1227 You Z.-Q. et al.
J. Chem. Phys.
(2015), 143, pp. 204104.
Link
,at the SCF level for vertical ionization and at the corresponding level (CIS, TDDFT or ADC, see Section 7.11.9) for vertical excitation. A brief introduction to these methods is given below, and users of the nonequilibrium PCM features are asked to cite Refs. 1227 and 754. 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.11.9. Users of this equilibrium-solvation PCM please cite and be referred to Ref. 753.

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 473 Hsu C.-P. et al.
J. Chem. Phys.
(2001), 114, pp. 3065.
Link

ω=𝑑𝐫𝑑𝐫𝑑𝐫′′𝑑𝐫′′′ρtr*(𝐫)(1|𝐫-𝐫|+gXC(𝐫,𝐫))×χ*(𝐫,𝐫′′,ω)(1|𝐫′′-𝐫′′′|+gXC(𝐫′′,𝐫′′′))ρtr(𝐫′′′), (11.5)

where χ is the charge density response function of the solvent and ρtr(𝐫) 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 (ω). Within a PCM, only classical Coulomb interactions are taken into account, and Eq. (11.5) becomes

ωPCM=𝑑𝐫𝑑𝐬ρtr*(𝐫)|𝐫-𝐬|𝑑𝐬𝑑𝐫𝒬(𝐬,𝐬,ε)ρtr(𝐫)|𝐬-𝐫|, (11.6)

where 𝒬 is PCM solvent response operator for a generic dielectric constant, ε. The integral of 𝒬 and the potential of the density ρ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 Ψi is obtained as a solution of a nonlinear Schrödinger equation

(H^vac+V^0slow+V^ifast)|Ψi=EiSS|Ψi (11.7)

that depends upon the solute’s charge distribution. Here H^vac is the usual Hamiltonian for the solute in vacuum and the reaction field operator 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, 485 Improta R. et al.
J. Chem. Phys.
(2006), 125, pp. 054103.
Link
, 729 Marenich A. V. et al.
Chem. Sci.
(2011), 2, pp. 2143.
Link
including both CIS and TDDFT.

In a linear dielectric medium, the solvent polarization is governed by the electric susceptibility, χ=[ε(ω)-1]/4π, where ε(ω) is the frequency-dependent permittivity. In case of very fast vertical transitions, the dielectric response is ruled by the optical dielectric constant, εopt=n2, 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: 222 Cossi M., Barone V.
J. Phys. Chem. A
(2000), 104, pp. 10614.
Link

𝐊εopt𝐪ifast,SS=𝐑εopt[𝐯i+𝐯(𝐪0slow)]. (11.8)

Here 𝐪fast,SS is the discretized fast surface charge. The dielectric constants in the matrices 𝐊 and 𝐑 (Section 11.2.2.1) are replaced with the optical dielectric constant, and 𝐯i is the potential of the solute’s excited state density, ρi. The quantity 𝐯(𝐪0slow) is the potential of the slow part of the apparent surface charges in the ground state, which are given by

𝐪0slow=(ε-εoptε-1)𝐪0. (11.9)

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.8) becomes

𝐊εopt𝐪ifast,LR=𝐑εopt𝐯(ρitr). (11.10)

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 Ψi is an eigenfunction of a different Hamiltonian, since Eq. (11.7) 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, 148 Cammi R. et al.
J. Chem. Phys.
(2005), 122, pp. 104513.
Link
, 155 Caricato M. et al.
J. Chem. Phys.
(2006), 124, pp. 124520.
Link
in what we have termed the “ptSS” method. 754 Mewes J.-M. et al.
J. Phys. Chem. A
(2015), 119, pp. 5446.
Link
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. 1227 You Z.-Q. et al.
J. Chem. Phys.
(2015), 143, pp. 204104.
Link
, 754 Mewes J.-M. et al.
J. Phys. Chem. A
(2015), 119, pp. 5446.
Link

The LR and SS approaches can also be used in the study of photon emission processes. 486 Improta R. et al.
J. Chem. Phys.
(2007), 127, pp. 074504.
Link
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.11.9.

11.2.2.4 Absorption and Emission Spectra Based on the Constrained Equilibrium Principle

For ultrafast processes in solution, such as electron transfer, photo-absorption/emission and photo-ionization, a continuum model should combine a proper nonequilibrium solvation theory to account for nonequilibrium solute–solvent interactions. In the traditional treatments, the nonequilibrium electrostatic solvation free energy was derived from the so-called reversible electric work integration along the path linking the initial equilibrium state (eq) and the intermediate nonequilibrium state (neq), 722 Marcus R. A.
J. Chem. Phys.
(1956), 24, pp. 966.
Link
, 723 Marcus R. A.
J. Chem. Phys.
(1956), 24, pp. 979.
Link
i.e.,

[ρ=0,Φ=0]𝜀[ρ1,Φ1eq]εopt[ρ2,Φ2neq] (11.11)

where Φ denotes the total electric potential including both the potential ψ due to the solute charge ρ in vacuum and polarization potential φ due to the medium. In order to deal with electron absorption and emission spectra in solution, the numerical expression of nonequilibrium solvation free energy which was established by intuitively collecting a series of energy terms from the interactions of solute charges and polarized charges, has been implemented using TDDFT with PCM model. 222 Cossi M., Barone V.
J. Phys. Chem. A
(2000), 104, pp. 10614.
Link
, 223 Cossi M., Barone V.
J. Chem. Phys.
(2000), 112, pp. 2427.
Link
It is easy to verify this numerical form can be achieved through the discretization of analytical expression of nonequilibrium solvation energy by traditional treatments. 665 Li X. Y.
Int. J. Quantum Chem.
(2015), 115, pp. 700.
Link
However, there exist a number of doubts on the overestimation of the solvent reorganization energy in ultrafast processes by this reversible electric work method. 109 Blackbourn R. L., Hupp J. T.
J. Phys. Chem.
(1990), 94, pp. 1788.
Link
, 376 Gomes P. J. S. et al.
J. Phys. Chem. A
(2010), 114, pp. 2778.
Link
It becomes clear now that there is no possibility to find a reversible pathway between the initial equilibrium state and the intermediate nonequilibrium state. Thus, the integrated electric work can not equal to the change of electrostatic free energy. 665 Li X. Y.
Int. J. Quantum Chem.
(2015), 115, pp. 700.
Link

Xiangyuan Li et al. 665 Li X. Y.
Int. J. Quantum Chem.
(2015), 115, pp. 700.
Link
, 1203 Wu H. W. et al.
Phys. Chem. Chem. Phys.
(2012), 14, pp. 5538.
Link
, 925 Ren H. S. et al.
J. Phys. Chem. A
(2013), 117, pp. 8017.
Link
established the new theory for nonequilibrium solvation by employing the constrained equilibrium principle using the following pathway

[ρ1,φ1eq]fast[ρ2,φ2neq]-λs[ρ2,φ2eq]+ρex,quasistaticC[ρ2+ρex,φ2neq] (11.12)

where C stands for the constrained equilibrium state which is constructed and mapped to the true nonequilibrium state by introducing the proper external charge ρex which is used to equilibrate the “residual” polarization potential, φ=φ2neq-φ2eq. In this way the solvent reorganization energy can be derived as

λs=-12VρexφdV (11.13)

Then the nonequilibrium solvation free energy is simply given by

F2neq=12Vρ2φ2eqdV+λs (11.14)

For more detailed descriptions of the gain of the external (constraining) charge ρex, or the equivalent constraining external electric field 𝐄ex, please refer to the review. 665 Li X. Y.
Int. J. Quantum Chem.
(2015), 115, pp. 700.
Link
Within the framework of continuum model, the discretization and numerical solution of Eq. (11.14) is expressed as 107 Bi T. J. et al.
Phys. Chem. Chem. Phys.
(2017), 19, pp. 32242.
Link

Fineq=mVi,mQi,mneq-12mQi,mneqD(ε)Qi,mneq (11.15)

where the subscript i denotes the ground (i = 1) or excited (i = 2) electronic state. The value Vi,m refers to the solute electrostatic potential at the mth tesserae. Qi,mneq is the apparent charge for the nonequilibrium state. D(ε) is the square matrix based on PCM versions (CPCM, IEFPCM, SSVPE, etc.).

The vertical excitation energy for absorption is defined as

hvab=G2neq-G1eq (11.16)

where G stands for total free energy of the solute in solution,

G1eq=E1+F1eq=E1+12mV1,mQ1,meq (11.17)

G1eq is calculated by the self-consistent reaction field (SCRF) method based on equilibrium ground-state reaction field and E1 means ground-state electronic energy of the solute. Based on the equilibrium ground-state reaction field, with a self-consistent state-specific (SS) method in the framework of TDDFT, G2neq is given by 107 Bi T. J. et al.
Phys. Chem. Chem. Phys.
(2017), 19, pp. 32242.
Link

G2neq=E1+ωab+mV1,mQ1,meq+mV2,m(Q2,mneq-Q1,meq)-12mQ2,mneqD(ε)Q2,mneq (11.18)

ωab is the excitation energy from TDDFT calculation.Alternatively, based on the nonequilibrium excited-state reaction field, G2neq is given by 106 Bi T. J. et al.
Phys. Chem. Chem. Phys.
(2018), 20, pp. 13178.
Link

G2,kneq=E1k+ωabk-mV1,mkQ2,mneq,k-12mQ2,mneq,kD(ε)Q2,mneq,k (11.19)

where k represents the kth iteration of the nonequilibrium excited-state reaction field at the ground-state geometry of solute. ωabk is the excitation energy from TDDFT calculation in the presence of the nonequilibrium excited-state reaction field. E1k stands for the ground-state electronic energy of solute at the kth iteration. Clearly, Eq. (11.19) is more physically-meaningful than Eq. (11.18).

Similarly, the vertical excitation energy for emission is given by 106 Bi T. J. et al.
Phys. Chem. Chem. Phys.
(2018), 20, pp. 13178.
Link

hvem=G2eq-G1neq (11.20)

where G2eq and G1neq represent the free energies of the equilibrium excited state and the nonequilibrium ground state at the excited-state equilibrium geometry, respectively, which can be expressed as 106 Bi T. J. et al.
Phys. Chem. Chem. Phys.
(2018), 20, pp. 13178.
Link

G2eq,k=ωemk+E1k-12mV2,mkQ2,meq,k+mV1,mkQ2,meq,k (11.21)
G1neq=E1′′+mV1,mQ1,mneq-12mQ1,mneqD(ε)Q1,mneq (11.22)

ωemk can be directly obtained by TDDFT calculation in the equilibrium excited-state reaction field at the kth iteration and E1k is the corresponding ground-state electronic energy of the solute. E1′′ is the ground-state electronic energy of the solute at the excited-state equilibrium geometry in the presence of the nonequilibrium ground-state reaction field.

The keyword TdNonEq is requested in the $pcm section. Refs. 665 should be cited if constrained equilibrium principle is employed to obtained the vertical absorption/emission energies in solution using the self-consistent SS-TDDFT method.