Ab initio quantum chemistry makes possible the study of gasphase molecular properties from first principles. In liquid solution, however, these properties may change significantly, especially in polar solvents. Although it is possible to model solvation effects by including explicit solvent molecules in the quantumchemical calculation (e.g. a supermolecular cluster calculation, averaged over different configurations of the molecules in the first solvation shell), such calculations are very computationally demanding. Furthermore, cluster calculations typically do not afford accurate solvation energies, owing to the importance of longrange electrostatic interactions. Accurate prediction of solvation free energies is, however, crucial for modeling of chemical reactions and ligand/receptor interactions in solution.
QChem contains several different implicit solvent models, which differ greatly in their level of sophistication and realism. These are generally known as selfconsistent reaction field (SCRF) models, because the continuum solvent establishes a reaction field (i.e., additional terms in the solute Hamiltonian) that depends upon the solute electron density, and must therefore be updated selfconsistently during the iterative convergence of the wave function. The simplest and oldest of these models that is available in QChem is the KirkwoodOnsager model [641, 642, 643], in which the solute molecule is placed inside of a spherical cavity and its electrostatic potential is represented in terms of a singlecenter multipole expansion. More sophisticated models, which use a moleculeshaped cavity and the full molecular electrostatic potential, include the conductorlike screening model (COSMO) [644] and the closelyrelated conductorlike PCM (CPCM) [645, 646, 647], along with the “surface and simulation of volume polarization for electrostatics” [SS(V)PE] model [648]. The latter is also known as the “integral equation formalism” (IEFPCM) [649, 650].
The CPCM and IEFPCM/SS(V)PE are examples of what are called “apparent surface charge” SCRF models, although the term polarizable continuum models (PCMs), as popularized by Tomasi and coworkers [651], is now used almost universally to refer to this class of solvation models. QChem employs a Switching/Gaussian or “SWIG” implementation of these PCMs [652, 653, 654, 655]. This approach resolves a longstanding—though littlepublicized—problem with standard PCMs, namely, that the boundaryelement methods used to discretize the solute/continuum interface may lead to discontinuities in the potential energy surface for the solute molecule. These discontinuities inhibit convergence of geometry optimizations, introduce serious artifacts in vibrational frequency calculations, and make ab initio molecular dynamics calculations virtually impossible [652, 653]. In contrast, QChem’s SWIG PCMs afford potential energy surfaces that are rigorously continuous and smooth. Unlike earlier attempts to obtain smooth PCMs, the SWIG approach largely preserves the properties of the underlying integralequation solvent models, so that solvation energies and molecular surface areas are hardly affected by the smoothing procedure.
Other solvent models available in QChem include the “Langevin dipoles” model [656, 657] as well as versions 8 and 12 of the SM models and the SMD model developed at the University of Minnesota [658, 659, 660]. SM8 and SM12 are based upon the generalized Born method for electrostatics, augmented with atomic surface tensions intended to capture nonelectrostatic effects (cavitation, dispersion, exchange repulsion, and changes in solvent structure). Empirical corrections of this sort are also available for the PCMs mentioned above, but within SM8 and SM12 these parameters have been optimized to reproduce experimental solvation energies. SMD combines IEFPCM and the cavitydispersionsolventstructure term.
Model 
Cavity 
Non 
Supported 
Derivatives 

Construction 
Discretization 
Electrostatic 
Basis 
Available? 

Terms? 
Sets 

KirkwoodOnsager 
spherical 
point charges 
no 
all 
none 
Langevin Dipoles 
atomic spheres 
dipoles in 
no 
all 
none 
(userdefinable) 
3d space 

CPCM 
atomic spheres 
point charges or 
user 
all 
1st and 
(userdefinable) 
smooth Gaussians 
specified 
2nd 

SS(V)PE/ 
atomic spheres 
point charges or 
user 
all 
1st 
IEFPCM 
(userdefinable) 
smooth Gaussians 
specified 

COSMO 
predefined 
point charges 
none 
all 
1st and 
atomic spheres 
2nd 

Isodensity SS(V)PE 
isodensity contour 
point charges 
none 
all 
none 
SM8 
predefined 
generalized 
automatic 
631G* 
1st 
atomic spheres 
Born 
631+G* 

631+G** 

SM12 
predefined 
generalized 
automatic 
all 
none 
atomic spheres 
Born 

SMD 
predefined 
point charges 
automatic 
all 
1st 
atomic spheres 
Table 11.1 summarizes the implicit solvent models that are available in QChem. Solvent models are invoked via the SOLVENT_METHOD keyword, as shown below. Additional details about each particular solvent model can be found in the sections that follow, and Table 11.1 provides a summary of the implicit solvent models that are available in QChem. In general, these methods are available for any SCF level of electronic structure theory, though in the case of SM8 only certain basis sets are supported. PostHartree–Fock calculations can be performed by first running an SCF + PCM job, in which case the correlated wave function will employ MOs and HartreeFock energy levels that are polarized by the solvent. Note that the jobcontrol format for specifying implicit solvent models changed significantly starting in QChem version 4.2.1. This change was made in an attempt to simply and unify the input notation for a large number of different models.
SOLVENT_METHOD
Sets the preferred solvent method.
TYPE:
STRING
DEFAULT:
0
OPTIONS:
0
Do not use a solvation model.
ONSAGER
Use the KirkwoodOnsager model (Section 11.2.1).
PCM
Use an apparent surface charge, polarizable continuum model
(Section 11.2.2).
ISOSVP
Use the isodensity implementation of the SS(V)PE model
(Section 11.2.5).
COSMO
Use COSMO (similar to CPCM but with an outlying charge
SM8
Use version 8 of the CramerTruhlar SM model (Section 11.2.7.1).
SM12
Use version 12 of the SM model (Section 11.2.7.2).
SMD
Use SMD (Section 11.2.8).
CHEM_SOL
Use the Langevin Dipoles model (Section 11.2.9).
RECOMMENDATION:
Consult the literature. PCM is a collective name for a family of models and additional input options may be required in this case, in order to fully specify the model. (See Section 11.2.2.) Several versions of SM12 are available as well, as discussed in Section 11.2.7.2.
Before going into detail about each of these models, a few potential points of confusion warrant mention, with regards to nomenclature. First, “PCM” refers to a family of models that includes CPCM and SS(V)PE/IEFPCM (the latter two being completely equivalent [650]). One or the other of these models can be selected by additional job control variables in a $pcm input section, as described in Section 11.2.2. COSMO is very similar to CPCM but includes a correction for that part of the solute’s electron density that penetrates beyond the cavity (the socalled “outlying charge”) [661, 662]. This is discussed in Section 11.2.6.
Two implementations of the SS(V)PE model are also available. The PCM implementation (which is requested by setting SOLVENT_METHOD = PCM) uses a solute cavity constructed from atomcentered spheres, as with most other PCMs. On the other hand, setting SOLVENT_METHOD = ISOSVP requests an SS(V)PE calculation in which the solute cavity is defined by an isocontour of the solute’s own electron density, as advocated by Chipman [648, 663, 664]. This is an appealing, oneparameter cavity construction, although it is unclear that this construction alone is superior in its accuracy to carefullyparameterized atomic radii [665], at least not without additional, nonelectrostatic terms included [666, 667] that are unfortunately not yet available in QChem’s implementation of the isodensity version of SS(V)PE. Moreover, analytic energy gradients are not available for the isodensity cavity construction, whereas they are available when the cavity is constructed from atomcentered spheres. One additional subtlety, which is discussed in detail in Ref. [654], is the fact that the PCM implementation of the equation for the SS(V)PE surface charges [Eq. eq:Kq=Rv] uses an asymmetric matrix. In contrast, Chipman’s isodensity implementation uses a symmetrized matrix. Although the symmetrized version is somewhat more computationally efficient when the number of surface charges is large, the asymmetric version is better justified, theoretically [654]. (This admittedly technical point is clarified in Section 11.2.2 and in particular in Table 11.2.)
Regarding the accuracy of these models for solvation free energies (), SM8 achieves subkcal/mol accuracy for neutral molecules, based on comparison to a large database of experimental values, although average errors for ions are more like 4 kcal/mol [668]. To achieve comparable accuracy with IEFPCM/SS(V)PE, nonelectrostatic terms must be included [669, 666, 667]. The SM12 model does not improve upon SM8 in any statistical sense [659], but does lift one important restriction on the level of electronic structure that can be combined with these models. Specifically, the Generalized Born model used in SM8 is based on a variant of Mullikenstyle atomic charges, and is therefore parameterized only for a few small basis sets, e.g., 631G*. SM12, on the other hand, uses a variety of charge schemes that are stable with respect to basisset expansion, and can therefore be combined with any level of electronic structure theory for the solute. Like IEFPCM, the SMD model is also applicable to any basis sets, and its accuracy is comparable to SM8 and SM12 [660]. Quantitative fluidphase thermodynamics can also be obtained using Klamt’s COSMORS approach [670, 671], where RS stands for “real solvent”. The COSMORS approach is not included in QChem and requires the COSMOtherm program, which is licensed separately through COSMOlogic [672], but QChem can write the input files that are need by COSMOtherm.
The following sections provide more details regarding theory and job control for the various implicit solvent models that are available in QChem. In addition, recent review articles are available for PCM methods [651], SM [668], and COSMO [673]. Formal relationships between various PCMs have been discussed in Refs. [663, 654].
The simplest implicit solvation model available in QChem is the KirkwoodOnsager model [641, 642, 643], wherein the solute is placed inside of a spherical cavity that is surrounded by a homogeneous dielectric medium. This model is characterized by two parameters: the cavity radius, , and the solvent dielectric constant, . The former is typically calculated according to
(11.1) 
where is the solute’s molar volume, usually obtained from experiment (molecular weight or density [674]), and is Avogadro’s number. It is also common to add 0.5 to the value of in Eq. (11.1) in order to account for the first solvation shell [675]. Alternatively, is sometimes selected as the maximum distance between the solute center of mass and the solute atoms, plus the relevant van der Waals radii. A third option is to set (the cavity diameter) equal to the largest solute–solvent internuclear distance, plus the the van der Waals radii of the relevant atoms. Unfortunately, solvation energies are typically quite sensitive to the choice of (and to the construction of the solute cavity, more generally).
Unlike older versions of the KirkwoodOnsager model, in which the solute’s electron distribution was described entirely in terms of its dipole moment, QChem’s version can use multipoles of arbitrarily high order, including the Born (monopole) term [676] for charged solutes, in order to describe the solute’s electrostatic potential. The solute–continuum electrostatic interaction energy is then computed using analytic expressions for the interaction of the point multipoles with a dielectric continuum.
Energies and analytic gradients for the KirkwoodOnsager solvent model are available for HartreeFock, DFT, and CCSD calculations. It is often advisable to perform a gasphase calculation of the solute molecule first, which can serve as the initial guess for a subsequent KirkwoodOnsager implicit solvent calculation.
The KirkwoodOnsager SCRF is requested by setting SOLVENT_METHOD = ONSAGER in the $rem section (along with normal job control variables for an energy or gradient calculation), and furthermore specifying several additional options in a $solvent input section, as described below. Of these, the keyword CavityRadius is required. The $rem variable CC_SAVEAMPL may save some time for CCSD calculations using the KirkwoodOnsager model.
NOTE: The following three job control variables belong only in the $solvent section. Do not place them in the $rem section. As with other parts of the QChem input file, this input section is not casesensitive.
CavityRadius
Sets the radius of the spherical solute cavity.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
No default.
OPTIONS:
Desired cavity radius, in ngstroms.
RECOMMENDATION:
Use Eq. eq1000.
Dielectric
Sets the dielectric constant of the solvent continuum.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
78.39
OPTIONS:
Use a (dimensionless) value of .
RECOMMENDATION:
As per required solvent; the default corresponds to water at 25C.
MultipoleOrder
Determines the order to which the multipole expansion of the solute charge density is carried out.
INPUT SECTION: $solvent
TYPE:
INTEGER
DEFAULT:
15
OPTIONS:
Include up to th order multipoles.
RECOMMENDATION:
Use the default. The multipole expansion is usually converged by order = 15.
Example 11.245 Onsager model applied at the HartreeFock level to HO in acetonitrile
$molecule 0 1 O 0.00000000 0.00000000 0.11722303 H 0.75908339 0.00000000 0.46889211 H 0.75908339 0.00000000 0.46889211 $end $rem method HF basis 631g** solvent_method Onsager $end $solvent CavityRadius 1.8 ! 1.8 Angstrom Solute Radius Dielectric 35.9 ! Acetonitrile MultipoleOrder 15 ! this is the default value $end
Example 11.246 CCSD/Onsager calculation applied to 1,2dichloroethane molecule
$comment 1,2dichloroethane GAUCHE Conformation $end $molecule 0 1 C 0.6541334418569877 0.3817051480045552 0.8808840579322241 C 0.6541334418569877 0.3817051480045552 0.8808840579322241 Cl 1.7322599856434779 0.0877596094659600 0.4630557359272908 H 1.1862455146007043 0.1665749506296433 1.7960750032785453 H 0.4889356972641761 1.4444403797631731 0.8058465784063975 Cl 1.7322599856434779 0.0877596094659600 0.4630557359272908 H 1.1862455146007043 0.1665749506296433 1.7960750032785453 H 0.4889356972641761 1.4444403797631731 0.8058465784063975 $end $rem METHOD CCSD BASIS 631g** N_FROZEN_CORE FC CC_SAVEAMPL 1 ! Save the CC amplitudes on disk SOLVENT_METHOD ONSAGER $end $solvent CavityRadius 3.65 ! in Angstroms Dielectric 8.93 ! methylene chloride $end
Example 11.247 KirkwoodOnsager SCRF applied to hydrogen fluoride in water, performing a gasphase calculation first.
$molecule 0 1 H 0.000000 0.000000 0.862674 F 0.000000 0.000000 0.043813 $end $rem METHOD HF BASIS 631G* $end @@@ $molecule 0 1 H 0.000000 0.000000 0.862674 F 0.000000 0.000000 0.043813 $end $rem JOBTYPE FORCE METHOD HF BASIS 631G* SOLVENT_METHOD ONSAGER SCF_GUESS READ ! read vacuum solution as a guess $end $solvent CavityRadius 2.5 $end
Clearly, the KirkwoodOnsager model is inappropriate if the solute is very nonspherical. 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 [651]. Apparent surface charge PCMs improve upon the KirkwoodOnsager model in two ways. Most importantly, they provide a much more realistic description of molecular shape, typically by constructing the solute cavity from a union of atomcentered spheres. 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 selfconsistently.
The PCM literature has a long history [651] and there are several different models in widespread use; connections between these models have not always been appreciated [648, 650, 663, 654]. Chipman [648, 663] has shown how various PCMs can be formulated within a common theoretical framework; see Ref. [655] for a pedagogical introduction. The PCM takes the form of a set of linear equations,
(11.2) 
in which the induced charges at the cavity surface discretization points [organized into a vector in Eq. eq:Kq=Rv] are computed from the values 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.2 for the PCMs that are available in QChem.
Model 
Literature 
Matrix 
Matrix 
Scalar 
Refs. 

COSMO 
[644] 



CPCM 




IEFPCM 




SS(V)PE 



The oldest PCM is the socalled DPCM model of Tomasi and coworkers [677], but unlike the models listed in Table 11.2, DPCM 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. eq:Kq=Rv [648]. As such, DPCM is essentially obsolete, and the PCMs available in QChem require only the evaluation of the electrostatic potential, not the electric field.
The simplest PCM that continues to enjoy widespread use is the ConductorLike Screening Model (COSMO) introduced by Klamt and Schüürmann [644]. Truong and Stefanovich [645] later implemented the same model with a slightly different dielectric scaling factor ( in Table 11.2), and called this modification GCOSMO. The latter was implemented within the PCM formalism by Barone and Cossi et al. [646, 647], who called the model CPCM (for “conductorlike” PCM). In each case, the dielectric screening factor has the form
(11.3) 
where Klamt and Schüürmann proposed but was used in GCOSMO and CPCM. 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 is a better compromise, given that the KirkwoodOnsager analytical result is for an thorder multipole centered in a spherical cavity [644, 662]. The distinction is irrelevant in highdielectric solvents; the and values of differ by only 0.6% for water at 25C, for example. Truong [645] argues that does a better job of preserving Gauss’ Law in lowdielectric solvents, but more accurate solvation energies (at least for neutral molecules, as compared to experiment) are sometimes obtained using [646]. This result is likely highly sensitive to cavity construction, and in any case, both versions are available in QChem.
Whereas the original COSMO model introduced by Klamt and Schüürmann [644] corresponds to Eq. eq:Kq=Rv with and as defined in Table 11.2, Klamt and coworkers later introduced a correction for outlying charge that goes beyond Eq. eq:Kq=Rv [661, 662]. Klamt now consistently refers to this updated model as “COSMO” [673], and we shall adopt this nomenclature as well. COSMO, with the outlying charge correction, is available in QChem and is described in Section 11.2.6. In contrast, CPCM consists entirely of Eq. eq:Kq=Rv with matrices and as defined in Table 11.2, although it is possible to modify the dielectric screening factor to use the value (as in COSMO) rather than the value. Additional nonelectrostatic terms can be added at the user’s discretion, as discussed below, but there is no explicit outlying charge correction in CPCM. These and other finetuning details for PCM jobs are controllable via the $pcm input section that is described in Section 11.2.3.
As compared to CPCM, a more sophisticated treatment of continuum electrostatic interactions is afforded by the “surface and simulation of volume polarization for electrostatics” [SS(V)PE] approach [648]. 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 [648], but this model is formally equivalent, at the level of integral equations, to the “integral equation formalism” (IEFPCM) that was developed originally by Cancès et al. [649, 678]. Some difference do arise when the integral equations are discretized to form finitedimensional matrix equations [654], however, and the reader will note from Table 11.2 that SS(V)PE uses a symmetrized form of the matrix as compared to IEFPCM. The asymmetric IEFPCM is the recommended approach [654], 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 DPCM approach, the original version of IEFPCM 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.2 [648, 650]. This version requires only the electrostatic potential, and is thus preferred, and it is this version that we designate as IEFPCM. The CPCM model becomes equivalent to SS(V)PE in the limit [648, 654], which means that CPCM must somehow include an implicit correction for volume polarization, even if this was not by design [661]. For , numerical calculations reveal that there is essentially no difference between SS(V)PE and CPCM results [654]. Since CPCM is less computationally involved as compared to SS(V)PE, it is the PCM of choice in highdielectric 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:2011a,Herbert:2016a.
Construction of the cavity surface is a crucial aspect of PCMs, and computed properties are quite sensitive to the details of the cavity construction. Traditionally [679, 680], and by default in QChem, solute cavities are constructed from a union of atomcentered spheres whose radii are 1.2 times the atomic van der Waals radii suggested by Bondi [681]. (This 20% augmentation is intended to mimic the fact that solvent molecules cannot approach all the way to the van der Waals radius of the solute atoms, but it is not clear that this is an optimal value. The default value in QChem is 1.2, but this value can be altered by the user.) Once the atomic radii are selected, the cavity surface is discretized using atomcentered Lebedev grids [148, 146, 149] of the same sort that are used to perform the numerical integrations in DFT. Surface charges are located at these grid points and the Lebedev quadrature weights can be used to define the surface area associated with each discretization point [652].
A longstanding (though not wellpublicized) problem with the aforementioned discretization procedure is that 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 [652]. It is also a fundamental hindrance to molecular dynamics calculations [653]. Recently, Lange and Herbert [652, 653] (building upon earlier work by York and Karplus [682]) 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 [653, 655]. 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 [652, 653]. 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 singlepoint energies and analytic energy gradients are available for SWIG PCMs, when the solute is described using molecular mechanics or an SCF (HartreeFock or DFT) electronic structure model. Analytic Hessians are available for the CPCM model only. (As usual, vibrational frequencies for other models will be computed, if requested, by finite difference of analytic energy gradients.) Singlepoint energy calculations using correlated wave functions can be performed in conjunction with these solvent models, in which case the correlated wave function calculation will utilize HartreeFock molecular orbitals that are polarized in the presence of the continuum dielectric solvent (i.e., there is no postHartree–Fock PCM correction).
Researchers who use these PCMs are asked to cite Refs. Lange:2010b,Lange:2011a, which provide the details of QChem’s implementation. (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 QChem 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.2) as well as how the cavity was constructed (e.g., “Bondi radii multiplied by 1.2”). Absent these details, 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 [683, 684] and the statespecific (SS) approach [680, 685, 686, 687]. Both are implemented in QChem [688], at the SCF level for vertical ionization and at the corresponding (CIS or TDDFT) level for vertical excitation. A brief introduction to these methods is given below, and users who utilize the nonequilibrium PCM features are asked to cite Refs. You:2015 and Mewes:2015a.
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 [689]
(11.4) 
where is the charge density response function of the solvent and 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. eq:LRterm becomes
(11.7) 
where is PCM solvent response operator for a generic dielectric constant, . The integral of and the potential of the density gives the surface charge density for the solvent polarization.
The statespecific (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 solvatedsolute excited state is obtained as a solution of a nonlinear Schrödinger equation
(11.9) 
that depends upon the solute’s charge distribution. Here is the usual Hamiltonian for the solute in vacuum and the reaction field operator 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 selfconsistently with respect to the solvent’s reaction field. In case of vertical ionization rather than excitation, both the ionized and nonionized states can be treated within a groundstate formalism. For vertical excitations, selfconsistent SS models [687, 690] have been developed for various excitedstate methods, including both CIS and TDDFT.
In a linear dielectric medium, the solvent polarization is governed by the electric susceptibility, , where is the frequencydependent permittivity, In case of very fast vertical transitions, the dielectric response is ruled by the optical dielectric constant, , where 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 can be obtained by solving the following equation [686]:
(11.10) 
Here 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 is the potential of the solute’s excited state density, . The quantity is the potential of the slow part of the apparent surface charges in the ground state, which are given by
(11.11) 
For LRPCM, the solvent polarization is subjected to the firstorder changes to the electron density (TDDFT linear density response), and thus Eq. eq:Kq=Rvfast becomes
(11.12) 
The LR approach for CIS/TDDFT excitations and the selfconsistent SS method (using the groundstate SCF) for vertical ionizations are available in QChem. The selfconsistent 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 is an eigenfunction of a different Hamiltonian, since Eq. eq:SSsdeq depend upon the specific state of interest. To avoid the ordering and the nonorthogonality problems, we compute the vertical excitation energy using a firstorder, perturbative approximation to the SS approach [691, 692], in what we have termed the "ptSS" method [440]. The zerothorder excitedstate wave function can be calculated using various excitedstate methods (currently available for CIS and TDDFT in QChem) with solventrelaxed molecular orbitals obtained from a groundstate 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 zerothorder transition density, in what we have called the "ptLR" approach. The combination of ptSS and ptLR yields quantitatively good solvatochromatic shifts, as compared to experimental results [688, 440].
The LR and SS approaches can also be used in the study of photon emission processes [693]. An emission process can be treated as a vertical excitation at a stationary point on the excitedstate potential surface. The basic requirement therefore is to prepare the solventrelaxed geometry for the excitedstate of interest. TDDFT/CPCM analytic gradients and Hessian are available in QChem; see Section 6.3.5 for computational details regarding excitedstate 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 ultrafast excited state decay, when only fast solvent degrees of freedom are expected to be equilibrated with the excitedstate 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 longlived excited state, e.g., strongly fluorescent species and phosphorescence. In the longlived 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 apporach, and the ground state is calculated using nonequilibrium selfconsistent SS approach.
A PCM calculation is requested by setting SOLVENT_METHOD = PCM in the $rem section. As mentioned above, there are a variety of different theoretical models that fall within the PCM family, so additional finetuning may be required, as described below.
Most PCM job control is accomplished via options specified in the $pcm input section, which allows the user to specify which flavor of PCM will be used, which algorithm will be used to solve the PCM equations, and other options. The format of the $pcm section is analogous to that of the $rem section:
$pcm <Keyword> <parameter/option> $end
NOTE: The following job control variables belong only in the $pcm section. Do not place them in the $rem section.
Theory
Specifies the which polarizable continuum model will be used.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
CPCM
OPTIONS:
CPCM
Conductorlike PCM with .
COSMO
Original conductorlike screening model with .
IEFPCM
IEFPCM with an asymmetric matrix.
SSVPE
SS(V)PE model, equivalent to IEFPCM with a symmetric matrix.
RECOMMENDATION:
The IEFPCM/SS(V)PE model is more sophisticated model than either CPCM or COSMO, and probably more appropriate for lowdielectric solvents, but it is also more computationally demanding. In highdielectric solvents there is little difference between these models. Note that the keyword COSMO in this context simply affects the dielectric screening factor ; to obtain the outlying charge correction suggested by Klamt [661, 662], one should use SOLVENT_METHOD = COSMO rather than SOLVENT_METHOD = PCM. (See Section 11.2.6.)
Method
Specifies which surface discretization method will be used.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
SWIG
OPTIONS:
SWIG
Switching/Gaussian method
ISWIG
“Improved” Switching/Gaussian method with an alternative switching function
Spherical
Use a single, fixed sphere for the cavity surface.
Fixed
Use discretization point charges instead of smooth Gaussians.
RECOMMENDATION:
Use of SWIG is recommended only because it is slightly more efficient than the switching function of ISWIG. On the other hand, ISWIG offers some conceptually more appealing features and may be superior in certain cases. Consult Refs. Lange:2010b,Lange:2011a for a discussion of these differences. The Fixed option uses the Variable Tesserae Number (VTN) algorithm of Li and Jensen [694], with Lebedev grid points. VTN uses point charges with no switching function or Gaussian blurring, and is therefore subject to discontinuities in geometry optimizations. It is not recommended, except to make contact with other calculations in the literature.
SwitchThresh
The threshold for discarding grid points on the cavity surface.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
8
OPTIONS:
Discard grid points when the switching function is less than .
RECOMMENDATION:
Use the default, which is found to avoid discontinuities within machine precision. Increasing reduces the cost of PCM calculations but can introduce discontinuities in the potential energy surface.
Construction of the solute cavity is an important part of the model and users should consult the literature in this capacity, especially with regard to the radii used for the atomic spheres. The default values provided in QChem correspond to the consensus choice that has emerged over several decades, namely, to use van der Waals radii scaled by a factor of 1.2. The most widelyused set of van der Waals radii are those determined from crystallographic data by Bondi [681], although the radius for hydrogen was later adjusted to 1.1 [695]. Bondi’s analysis was later extended to the whole main group [696], and this entire extended set of van der Waals radii is available in QChem. (For simplicity, we call these “Bondi radii” regardless of whether they come from Bondi’s original paper or the later work.) Alternatively, atomic radii from the Universal Force Field (UFF) are available [697]. The main appeal of UFF radii is that they are defined for all atoms of the periodic table, though the quality of these radii for PCM applications is unclear. Finally, the user may specify his or her own radii for cavity construction using a $van_der_waals input section, the format for which is described in Section 11.2.9. No scaling factor is applied to userdefined radii. Note that is allowed for a particular atomic radius, in which case the atom in question is not used to construct the cavity surface. This feature facilitates the construction of “united atom” cavities, in which the hydrogen atoms do not get their own spheres and the heavyatom radii are increased to compensate [665].
Scaling the atomic radii is a crude way to account for the fact that solvent molecules should not typically penetrate all the way to the van der Waals radii of the solute. Another way to account for this effect is to employ a “solventaccessible” cavity surface, which is constructed from the van der Waals surface by adding a certain value to each atomic radius. This value is the presumed radius of a solvent molecule. (The value 1.4 is often used for water.) The capability to specify a solvent “probe” radius to be added to the atomic radii of choice is available in QChem.
Radii
Specifies which set of atomic van der Waals radii will be used to define the solute cavity.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
BONDI
OPTIONS:
BONDI
Use the (extended) set of Bondi radii.
FF
Use LennardJones radii from a molecular mechanics force field.
UFF
Use radii form the Universal Force Field.
READ
Read the atomic radii from a $van_der_waals input section.
RECOMMENDATION:
Bondi radii are widely used. The FF option requires the user to specify an MM force field using the FORCE_FIELD $rem variable, and also to define the atom types in the $molecule section (see Section 11.3). This is not required for UFF radii.
vdwScale
Scaling factor for the atomic van der Waals radii used to define the solute cavity.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
1.2
OPTIONS:
Use a scaling factor of .
RECOMMENDATION:
The default value is widely used in PCM calculations, although a value of 1.0 might be appropriate if using a solventaccessible surface.
SASrad
Form a “solvent accessible” surface with the given solvent probe radius.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
None
OPTIONS:
Use a solvent probe radius of , in .
RECOMMENDATION:
The solvent probe radius is added to the scaled van der Waals radii of the solute atoms. A common solvent probe radius for water is 1.4 , but the user should consult the literature regarding the use of solventaccessible surfaces.
Historically, discretization of the cavity surface has involved “tessellation” methods that divide the cavity surface area into finite polygonal “tesserae”. (The GEPOL algorithm [698] is perhaps the most widelyused tessellation scheme.) Tessellation methods, however, suffer not only from discontinuities in the cavity surface area and solvation energy as a function of the nuclear coordinates, but in addition they lead to analytic energy gradients that are complicated to derive and implement. To avoid these problems, QChem’s SWIG PCM implementation [652, 653, 654] uses Lebedev grids to discretize the atomic spheres. These are atomcentered grids with icosahedral symmetry, and may consist of anywhere from 26 to 5294 grid points per atomic sphere. The default values used by QChem were selected based on extensive numerical tests [653, 654]. The default for MM atoms (in MM/PCM or QM/MM/PCM jobs) is Lebedev points per atomic sphere, whereas the default for QM atoms is . (This represents a change relative to QChem versions earlier than 4.2.1, where the default for QM atoms was .) These default values exhibit good rotational invariance and absolute solvation energies that, in most cases, lie within 0.5–1.0 kcal/mol of the limit [653], although exceptions (especially where charged solutes are involved) can be found [654]. The acceptable values for the number of Lebedev points per sphere are , 50, 110, 194, 302, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, or 5294.
HPoints
The number of Lebedev grid points to be placed on H atoms in the QM system.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
302
OPTIONS:
Acceptable values are listed above.
RECOMMENDATION:
Use the default for geometry optimizations. For absolute solvation energies, the user may want to examine convergence with respect to .
HeavyPoints
The number of Lebedev grid points to be placed nonhydrogen atoms in the QM system.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
302
OPTIONS:
Acceptable values are listed above.
RECOMMENDATION:
Use the default for geometry optimizations. For absolute solvation energies, the user may want to examine convergence with respect to .
MMHPoints
The number of Lebedev grid points to be placed on H atoms in the MM subsystem.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
110
OPTIONS:
Acceptable values are listed above.
RECOMMENDATION:
Use the default for geometry optimizations. For absolute solvation energies, the user may want to examine convergence with respect to . This option applies only to MM/PCM or QM/MM/PCM calculations.
MMHeavyPoints
The number of Lebedev grid points to be placed on nonhydrogen atoms in the MM subsystem.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
110
OPTIONS:
Acceptable values are listed above.
RECOMMENDATION:
Use the default for geometry optimizations. For absolute solvation energies, the user may want to examine convergence with respect to . This option applies only to MM/PCM or QM/MM/PCM calculations.
Especially for complicated molecules, the user may want to visualize the cavity surface. This can be accomplished by setting PrintLevel , which will trigger the generation of several “.PQR” files that describe the cavity surface. (These are written to the QChem output file.) The .PQR format is similar to the common .PDB (Protein Data Bank) format, but also contains charge and radius information for each atom. One of the output .PQR files contains the charges computed in the PCM calculation and radii (in ) that are half of the square root of the surface area represented by each surface grid point. Thus, in examining this representation of the surface, larger discretization points are associated with larger surface areas. A second .PQR file contains the solute’s electrostatic potential (in atomic units), in place of the charge information, and uses uniform radii for the grid points. These .PQR files can be visualized using various thirdparty software, including the freelyavailable Visual Molecular Dynamics (VMD) program [550, 551], which is particularly useful for coloring the .PQR surface grid points according to their charge, and sizing them according to their contribution to the molecular surface area. (Examples of such visualizations can be found in Ref. Lange:2010a.)
PrintLevel
Controls the printing level during PCM calculations.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
Prints PCM energy and basic surface grid information. Minimal additional printing.
1
Level 0 plus PCM solutesolvent interaction energy components and Gauss’ Law error.
2
Level 1 plus surface grid switching parameters and a .PQR file for visualization of
the cavity surface apparent surface charges.
3
Level 2 plus a .PQR file for visualization of the electrostatic potential at the surface
grid created by the converged solute.
4
Level 3 plus additional surface grid information, electrostatic potential and apparent
surface charges on each SCF cycle.
5
Level 4 plus extensive debugging information.
RECOMMENDATION:
Use the default unless further information is desired.
Finally, note that setting Method to Spherical in the $pcm input selection requests the construction of a solute cavity consisting of a single, fixed sphere. This is generally not recommended but is occasionally useful for making contact with the results of Born models in the literature, or the KirkwoodOnsager model discussed in Section 11.2.1. In this case, the cavity radius and its center must also be specified in the $pcm section. The keyword HeavyPoints controls the number of Lebedev grid points used to discretize the surface.
CavityRadius
Specifies the solute cavity radius.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
None
OPTIONS:
Use a radius of , in ngstroms.
RECOMMENDATION:
None.
CavityCenter
Specifies the center of the spherical solute cavity.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
0.0 0.0 0.0
OPTIONS:
Coordinates of the cavity center, in ngstroms.
RECOMMENDATION:
The format is CavityCenter followed by three floatingpoint values, delineated by spaces. Uses the same coordinate system as the $molecule section.
The following example shows a very basic PCM job. The solvent dielectric is specified in the $solvent section, which is described below.
Example 11.248 A basic example of using the PCMs: optimization of trifluoroethanol in water.
$rem JOBTYPE OPT BASIS 631G* METHOD B3LYP SOLVENT_METHOD PCM $end $pcm Theory CPCM Method SWIG Solver Inversion HeavyPoints 194 HPoints 194 Radii Bondi vdwScale 1.2 $end $solvent Dielectric 78.39 $end $molecule 0 1 C 0.245826 0.351674 0.019873 C 0.244003 0.376569 1.241371 O 0.862012 0.527016 2.143243 F 0.776783 0.909300 0.666009 F 0.858739 0.511576 0.827287 F 1.108290 1.303001 0.339419 H 0.587975 0.878499 1.736246 H 0.963047 1.147195 0.961639 H 0.191283 1.098089 2.489052 $end
The next example uses a single spherical cavity and should be compared to the KirkwoodOnsager job, Example 11.11.1 on page *.
Example 11.249 PCM with a single spherical cavity, applied to HO in acetonitrile
$molecule 0 1 O 0.00000000 0.00000000 0.11722303 H 0.75908339 0.00000000 0.46889211 H 0.75908339 0.00000000 0.46889211 $end $rem method HF basis 631g** solvent_method pcm $end $pcm method spherical ! single spherical cavity with 590 discretization points HeavyPoints 590 CavityRadius 1.8 ! Solute Radius, in Angstrom CavityCenter 0.0 0.0 0.0 ! Will be at center of Standard Nuclear Orientation Theory SSVPE $end $solvent Dielectric 35.9 ! Acetonitrile $end
Finally, we consider an example of a unitedatom cavity. Note that a userdefined van der Waals radius is supplied only for carbon, so the hydrogen radius is taken to be zero and thus the hydrogen atoms are not used to construct the cavity surface. (As mentioned above, the format for the $van_der_waals input section is discussion in Section 11.2.9).
Example 11.250 Unitedatom cavity construction for ethylene.
$comment Benzene (in benzene), with a unitedatom cavity construction R = 2.28 A for carbon, R = 0 for hydrogen $end $molecule 0 1 C 1.3862000000 0.0000000000 0.0000000000 C 0.6931000000 1.2004844147 0.0000000000 C 0.6931000000 1.2004844147 0.0000000000 C 1.3862000000 0.0000000000 0.0000000000 C 0.6931000000 1.2004844147 0.0000000000 C 0.6931000000 1.2004844147 0.0000000000 H 2.4618000000 0.0000000000 0.0000000000 H 1.2309000000 2.1319813390 0.0000000000 H 1.2309000000 2.1319813390 0.0000000000 H 2.4618000000 0.0000000000 0.0000000000 H 1.2309000000 2.1319813390 0.0000000000 H 1.2309000000 2.1319813390 0.0000000000 $end $rem exchange hf basis 631G* solvent_method pcm $end $pcm theory iefpcm ! this is a synonym for ssvpe method swig printlevel 1 radii read $end $solvent dielectric 2.27 $end $van_der_waals 1 6 2.28 $end
The solvent for PCM calculations is specified using the $solvent section, as documented below. In addition, the $solvent section can be used to incorporate nonelectrostatic interaction terms into the solvation energy. (The Theory keyword in the $pcm section specifies only how the electrostatic interactions are handled.) The general form of the $solvent input section is shown below. The $solvent section was used above to specify parameters for the KirkwoodOnsager SCRF model, and will be used again below to specify the solvent for SM calculations (Section 11.2.7); in each case, the particular options that can be listed in the $solvent section depend upon the value of the $rem variable SOLVENT_METHOD.
$solvent NonEls <Option> NSolventAtoms <Number unique of solvent atoms> SolventAtom <Number1> <Number2> <Number3> <SASrad> SolventAtom <Number1> <Number2> <Number3> <SASrad> . . . <Keyword> <parameter/option> . . . $end
The keyword SolventAtom requires multiple parameters, whereas all other keywords require only a single parameter. In addition to any (optional) nonelectrostatic parameters, the $solvent section is also used to specify the solvent’s dielectric constant. If nonelectrostatic interactions are ignored, then this is the only keyword that is necessary in the $solvent section. For nonequilibrium TDDFT/CPCM calculations (Section 11.2.2.3), the optical dielectric constant should be specified in the $solvent section as well.
Dielectric
The (static) dielectric constant of the PCM solvent.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
78.39
OPTIONS:
Use a dielectric constant of .
RECOMMENDATION:
The static (i.e., zerofrequency) dielectric constant is what is usually called “the” dielectric constant. The default corresponds to water at 25C.
OpticalDielectric
The optical dielectric constant of the PCM solvent.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
1.78
OPTIONS:
Use an optical dielectric constant of .
RECOMMENDATION:
The default corresponds to water at 25C. Note that , where is the solvent’s index of refraction.
The nonelectrostatic interactions currently available in QChem are based on the work of Cossi et al. [699], and are computed outside of the SCF procedure used to determine the electrostatic interactions. The nonelectrostatic energy is highly dependent on the input parameters and can be extremely sensitive to the radii chosen to define the solute cavity. Accordingly, the inclusion of nonelectrostatic interactions is highly empirical and probably needs to be considered on a casebycase basis. Following Ref. Cossi:1996, the cavitation energy is computed using the same solute cavity that is used to compute the electrostatic energy, whereas the dispersion/repulsion energy is computed using a solventaccessible surface.
The following keywords (in the $solvent section) are used to define nonelectrostatic parameters for PCM calculations.
NonEls
Specifies what type of nonelectrostatic contributions to include.
INPUT SECTION: $solvent
TYPE:
STRING
DEFAULT:
None
OPTIONS:
Cav
Cavitation energy
Buck
Buckingham dispersion and repulsion energy from atomic number
LJ
LennardJones dispersion and repulsion energy from force field
BuckCav
Buck + Cav
LJCav
LJ + Cav
RECOMMENDATION:
A very limited set of parameters for the Buckingham potential is available at present.
NSolventAtoms
The number of different types of atoms.
INPUT SECTION: $solvent
TYPE:
INTEGER
DEFAULT:
None
OPTIONS:
Specifies that there are different types of atoms.
RECOMMENDATION:
This keyword is necessary when NonEls = Buck, LJ, BuckCav, or LJCav. Methanol (), for example, has three types of atoms (C, H, and O).
SolventAtom
Specifies a unique solvent atom.
INPUT SECTION: $solvent
TYPE:
Various
DEFAULT:
None.
OPTIONS:
Input (TYPE)
Description
Number1 (INTEGER):
The atomic number of the atom
Number2 (INTEGER):
How many of this atom are in a solvent molecule
Number3 (INTEGER):
Force field atom type
SASrad (FLOAT):
Probe radius (in ) for defining the solvent accessible surface
RECOMMENDATION:
If not using LJ or LJCav, Number3 should be set to 0. The SolventAtom keyword is necessary when NonEls = Buck, LJ, BuckCav, or LJCav.
Temperature
Specifies the solvent temperature.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
300.0
OPTIONS:
Use a temperature of , in Kelvin.
RECOMMENDATION:
Used only for the cavitation energy.
Pressure
Specifies the solvent pressure.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
1.0
OPTIONS:
Use a pressure of , in bar.
RECOMMENDATION:
Used only for the cavitation energy.
SolventRho
Specifies the solvent number density
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
Determined for water, based on temperature.
OPTIONS:
Use a density of , in molecules/.
RECOMMENDATION:
Used only for the cavitation energy.
SolventRadius
The radius of a solvent molecule of the PCM solvent.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
None
OPTIONS:
Use a radius of , in .
RECOMMENDATION:
Used only for the cavitation energy.
The following example illustrates the use of the nonelectrostatic interactions.
Example 11.251 Optimization of trifluoroethanol in water using both electrostatic and nonelectrostatic PCM interactions. OPLSAA parameters are used in the LennardJones potential for dispersion and repulsion.
$rem JOBTYPE OPT BASIS 631G* METHOD B3LYP SOLVENT_METHOD PCM FORCE_FIELD OPLSAA $end $pcm Theory CPCM Method SWIG Solver Inversion HeavyPoints 194 HPoints 194 Radii Bondi vdwScale 1.2 $end $solvent NonEls LJCav NSolventAtoms 2 SolventAtom 8 1 186 1.30 SolventAtom 1 2 187 0.01 SolventRadius 1.35 Temperature 298.15 Pressure 1.0 SolventRho 0.03333 Dielectric 78.39 $end $molecule 0 1 C 0.245826 0.351674 0.019873 23 C 0.244003 0.376569 1.241371 22 O 0.862012 0.527016 2.143243 24 F 0.776783 0.909300 0.666009 26 F 0.858739 0.511576 0.827287 26 F 1.108290 1.303001 0.339419 26 H 0.587975 0.878499 1.736246 27 H 0.963047 1.147195 0.961639 27 H 0.191283 1.098089 2.489052 25 $end
The OpticalDielectric keyword in $solvent is always needed. The LR energy is automatically calculated while the CIS/TDDFT calculations are performed with PCM, but it is turned off while the perturbation scheme is employed.
ChargeSeparation
Partition fast and slow charges in solvent equilibrium state
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
No default.
OPTIONS:
Marcus
Do slowfast charge separation in the ground state.
Excited
Do slowfast charge separation in an excitedstate.
RECOMMENDATION:
Charge separation is used in conjunction with the StateSpecific keyword in $pcm.
StateSpecific
Specifies which the statespecific method will be used.
INPUT SECTION: $pcm
TYPE:
Various
DEFAULT:
No default.
OPTIONS:
Ground
Run selfconsistent SS method in the groundstate with a given slow polarization charges.
Perturb
Perform ptSS and ptLR for vertical excitations.
The th excitedstate used for charge separation (for emission).
RECOMMENDATION:
NoneqGrad
Control whether perform excited state geometry optimization in equilibrium or nonequilibrium.
INPUT SECTION: $pcm
TYPE:
NONE
DEFAULT:
No default.
OPTIONS:
RECOMMENDATION:
Specify it for nonequilibrium optimization otherwise equilibrium geometry optimization will be performed.
Example 11.252 LRTDDFT/CPCM lowlying vertical excitation energy
$molecule 0 1 C 0 0 0.0 O 0 0 1.21 $end $rem EXCHANGE B3LP CIS_N_ROOTS 10 CIS_SINGLETS TRUE CIS_TRIPLETS TRUE RPA TRUE BASIS 631+G* SOLVENT_METHOD PCM $end $pcm Theory CPCM Method SWIG Solver Inversion Radii Bondi $end $solvent Dielectric 78.39 OpticalDielectric 1.777849 $end
Example 11.253 PCM solvation effects on the vertical excitation energies of planar DMABN using the ptSS and ptLR methods.
$molecule 0 1 C 0.0000466670 0.0003989920 1.9049533013 C 1.2100273911 0.0003798174 1.1860514263 C 1.2146402411 0.0000659664 0.1945152488 C 0.0001642918 0.0006163800 0.9338323670 C 1.2143493593 0.0015573129 0.1946877858 C 1.2097539311 0.0018468783 1.1857751832 H 2.1519498961 0.0013779318 1.7220180544 H 2.1643711416 0.0004818107 0.7096402413 H 2.1640827320 0.0020086712 0.7097815919 H 2.1517638721 0.0022872146 1.7216153743 C 0.0002270265 0.0010614144 3.3253023786 N 0.0004752169 0.0024050570 4.4843211392 N 0.0000539889 0.0001561808 2.2973728307 C 1.2586568849 0.0012846735 3.0369947926 H 1.0410420397 0.0016150528 4.1023766685 H 1.8608970185 0.8856472814 2.8111179665 H 1.8592477277 0.8891330458 2.8102379603 C 1.2585638668 0.0006605040 3.0372852058 H 1.8606513946 0.8862082585 2.8107557102 H 1.8593622362 0.8886046730 2.8114618382 H 1.0406646935 0.0000970072 4.1026090234 $end $rem JOBTYPE SP EXCHANGE LRCwPBEPBE OMEGA 260 BASIS 631G* CIS_N_ROOTS 10 RPA 2 CIS_SINGLETS 1 CIS_TRIPLETS 0 CIS_DYNAMIC_MEM TRUE CIS_RELAXED_DENSITY TRUE USE_NEW_FUNCTIONAL TRUE SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm Theory IEFPCM ChargeSeparation Marcus StateSpecific Perturb $end $solvent Dielectric 35.688000 ! Acetonitrile OpticalDielectric 1.806874 $end
Example 11.254 Aqueous phenol ionization using statespecific nonequilibrium PCM
$molecule 0 1 C 0.189057 1.215927 0.000922 H 0.709319 2.157526 0.001587 C 1.194584 1.155381 0.000067 H 1.762373 2.070036 0.000230 C 1.848872 0.069673 0.000936 H 2.923593 0.111621 0.001593 C 1.103041 1.238842 0.001235 H 1.595604 2.196052 0.002078 C 0.283047 1.185547 0.000344 H 0.862269 2.095160 0.000376 C 0.929565 0.042566 0.000765 O 2.287040 0.159171 0.001759 H 2.663814 0.725029 0.001075 $end $rem JOBTYPE SP EXCHANGE wPBE CORRELATION PBE LRC_DFT 1 OMEGA 300 BASIS 631+G* SCF_CONVERGENCE 8 SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm ChargeSeparation Marcus $end $solvent Dielectric 78.39 OpticalDielectric 1.777849 $end @@@ $molecule 1 2 READ $end $rem JOBTYPE SP EXCHANGE wPBE CORRELATION PBE LRC_DFT 1 OMEGA 300 BASIS 631+G* SCF_CONVERGENCE 8 SOLVENT_METHOD PCM PCM_PRINT 1 SCF_GUESS READ $end $pcm StateSpecific Ground $end $solvent Dielectric 78.39 OpticalDielectric 1.777849 $end
Example 11.255 PCM solvation effects on the emission energy of twisted DMABN in acetonitrile.
$molecule 0 1 C 0.0000020664 0.0000030657 1.8944472102 C 1.2301788386 0.0000034932 1.1607654421 C 1.2448623254 0.0000019224 0.2003601210 C 0.0000035154 0.0000130458 0.9179014646 C 1.2448450017 0.0000007858 0.2003542615 C 1.2301606293 0.0000009019 1.1607531534 H 2.1711288221 0.0000135363 1.7040640230 H 2.1806701315 0.0000099280 0.7501405321 H 2.1806573730 0.0000061271 0.7501203529 H 2.1711120141 0.0000069286 1.7040377372 C 0.0000236776 0.0000026392 3.2935096257 N 0.0000117329 0.0000010959 4.4594155390 N 0.0000087428 0.0000059515 2.3091936565 C 0.0000010176 1.2538813765 3.0212357565 H 0.0000057176 1.0990763913 4.0986604824 H 0.8829081343 1.8189750963 2.7058150336 H 0.8829113886 1.8189675327 2.7058056136 C 0.0000093023 1.2538863802 3.0212468120 H 0.8829172783 1.8189662532 2.7058494734 H 0.8829018393 1.8189927074 2.7058033199 H 0.0000385216 1.0990713095 4.0986701147 $end $rem JOBTYPE SP EXCHANGE LRCwPBEPBE OMEGA 280 BASIS 631G* CIS_N_ROOTS 10 RPA 2 CIS_SINGLETS 1 CIS_TRIPLETS 0 CIS_DYNAMIC_MEM TRUE CIS_RELAXED_DENSITY TRUE USE_NEW_FUNCTIONAL TRUE SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm ChargeSeparation Excited StateSpecific 1 $end $solvent Dielectric 35.688000 ! Acetonitrile OpticalDielectric 1.806874 $end @@@ $molecule READ $end $rem JOBTYPE SP EXCHANGE LRCwPBEPBE OMEGA 280 BASIS 631G* SCF_GUESS READ SCF_CONVERGENCE 8 SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm StateSpecific Ground $end $solvent Dielectric 35.688000 ! Acetonitrile OpticalDielectric 1.806874 $end
Recall that PCM electrostatics calculations require the solution of the set of linear equations given in Eq. eq:Kq=Rv, to determine the vector of apparent surface charges. The precise forms of the matrices and depend upon the particular PCM (Table 11.2), but in any case they have dimension , where is the number of Lebedev grid points used to discretize the cavity surface. Construction of the matrix affords a numerically exact solution to Eq. eq:Kq=Rv, whose cost scales as in CPU time and in memory. This cost is exacerbated by smooth PCMs, which discard fewer interior grid points so that tends to be larger, for a given solute, as compared to traditional discretization schemes. For QM solutes, the cost of inverting is usually negligible relative to the cost of the electronic structure calculation, but for the large values of that are encountered in MM/PCM or QM/MM/PCM jobs, the cost of inverting is often prohibitively expensive.
To avoid this bottleneck, Lange and Herbert [655] have developed an iterative conjugate gradient (CG) solver for Eq. eq:Kq=Rv whose cost scales as in CPU time and in memory. A number of other costsaving options are available, including efficient preconditioners and matrix factorizations that speed up convergence of the CG iterations, and a fast multipole algorithm for computing the electrostatic interactions [700]. Together, these features lend themselves to a solution of Eq. eq:Kq=Rv whose cost scales as in both memory and CPU time, for sufficiently large systems [655]. Currently, these options are available only for CPCM, not for SS(V)PE/IEFPCM.
Listed below are job control variables for the CG solver, which should be specified within the $pcm input section. Researchers who use this feature are asked to cite the original SWIG PCM references [653, 654] as well as the reference for the CG solver [655].
Solver
Specifies the algorithm used to solve the PCM equations.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
INVERSION
OPTIONS:
INVERSION
Direct matrix inversion
CG
Iterative conjugate gradient
RECOMMENDATION:
Matrix inversion is faster for small solutes because it needs to be performed only once in a singlepoint calculation. However, the CG solver (which must be applied at each SCF iteration) is recommended for large MM/PCM or QM/MM/PCM calculations.
CGThresh
The threshold for convergence of the conjugate gradient solver.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
6
OPTIONS:
Conjugate gradient converges when the maximum residual is less than .
RECOMMENDATION:
The default typically affords PCM energies on par with the precision of matrix inversion for small systems. For systems that have difficulty with SCF convergence, one should increase or try the matrix inversion solver. For wellbehaved or very large systems, a smaller might be permissible.
DComp
Controls decomposition of matrices to reduce the matrix norm for the CG Solver.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
1
OPTIONS:
0
Turns off matrix decomposition
1
Turns on matrix decomposition
3
Option 1 plus only stores upper half of matrix and enhances gradient evaluation
RECOMMENDATION:
None
PreCond
Controls the use of the preconditioner for the CG solver.
INPUT SECTION: $pcm
TYPE:
None
DEFAULT:
Off
OPTIONS:
No options. Specify the keyword to enable preconditioning.
RECOMMENDATION:
A Jacobi blockdiagonal preconditioner is applied during the conjugate gradient algorithm to improve the rate of convergence. This reduces the number of CG iterations, at the expense of some overhead. Preconditioning is generally recommended for large systems.
NoMatrix
Specifies whether PCM matrices should be explicitly constructed and stored.
INPUT SECTION: $pcm
TYPE:
None
DEFAULT:
Off
OPTIONS:
No options. Specify the keyword to avoid explicit construction of PCM matrices.
RECOMMENDATION:
Storing the PCM matrices requires memory. If this is prohibitive, the NoMatrix option forgoes explicit construction of the PCM matrices, and instead constructs the matrix elements as needed, reducing the memory requirement to at the expense of additional computation.
UseMultipole
Controls the use of the adaptive fast multipole method in the CG solver.
INPUT SECTION: $pcm
TYPE:
None
DEFAULT:
Off
OPTIONS:
No options. Specify the keyword in order to enable the fast multipole method.
RECOMMENDATION:
The fast multipole approach formally reduces the CPU time to , but is only beneficial for spatially extended systems with several thousand cavity grid points. Requires the use of NoMatrix.
MultipoleOrder
Specifies the highest multipole order to use in the FMM.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
4
OPTIONS:
The highest order multipole in the multipole expansion.
RECOMMENDATION:
Increasing the multipole order improves accuracy but also adds more computational expense. The default yields satisfactory performance in common QM/MM/PCM applications.
Theta
The multipole acceptance criterion.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
0.6
OPTIONS:
A number between zero and one.
RECOMMENDATION:
The default is recommended for general usage. This variable determines when the use of a multipole expansion is valid. For a given grid point and box center in the FMM, a multipole expansion is accepted when Theta, where is the distance from the grid point to the box center and is the radius of the box. Setting Theta to one will accept all multipole expansions, whereas setting it to zero will accept none. If not accepted, the grid point’s interaction with each point inside the box is computed explicitly. A low Theta is more accurate but also more expensive than a higher Theta.
NBox
The FMM boxing threshold.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
100
OPTIONS:
The maximum number of grid points for a leaf box.
RECOMMENDATION:
The default is recommended. This option is for advanced users only. The adaptive FMM boxing algorithm divides space into smaller and smaller boxes until each box has less than or equal to NBox grid points. Modification of the threshold can lead to speedup or slowdown depending on the molecular system and other FMM variables.
A sample input file for the linearscaling QM/MM/PCM methodology can be found in the $QC/samples directory, under the name QMMMPCM_crambin.in. This sample involves a QM/MM description of a protein (crambin) in which a single tyrosine side chain is taken to be the QM region. The entire protein is immersed in a dielectric using CPCM with SWIG discretization.
As discussed above, results obtained various types of PCMs are quite sensitive to the details of the cavity construction. QChem’s implementation of PCMs, using Lebedev grids, simplifies this construction somewhat, but leaves the radii of the atomic spheres as empirical parameters (albeit ones for which widelyused default values are provided). An alternative implementation of the SS(V)PE solvation model is also available [664], which attempts to further eliminate empiricism associated with cavity construction by taking the cavity surface to be a specified isocontour of the solute’s electron density. (We call this the isodensity implementation of SS(V)PE in Table 11.2, and it is based on Chipman’s “symmetrized” form of the matrix [664, 654].) In this case, the cavity surface is discretized by projecting a singlecenter Lebedev grid onto the isocontour surface. Unlike the PCM implementation discussed in Section 11.2.2, for which pointgroup symmetry is disabled, this implementation of SS(V)PE supports full symmetry for all Abelian point groups. The larger and/or the less spherical the solute molecule is, the more points are needed to get satisfactory precision in the results. Further experience will be required to develop detailed recommendations for this parameter. Values as small as 110 points are usually sufficient for diatomic or triatomic molecules. The default value of 1202 points is adequate to converge the energy within 0.1 kcal/mol for solutes the size of monosubstituted benzenes.
Although cavitation, dispersion, and and specific hydrogenbonding corrections for use with the isodensity SS(V)PE model have been reported in the literature [666, 667], these are not yet available in QChem. Energy gradients are also not available for this implementation of SS(V)PE, although they are available for the implementation described in Section 11.2.2 in which the cavity is constructed from atomcentered spheres. As with the PCMs discussed in that section, the solute may be described using HartreeFock theory or DFT; postHartree–Fock correlated wave functions can also take advantage of molecular orbitals that are polarized using SS(V)PE. Researchers who use the isodensity SS(V)PE feature are asked to cite Ref. Chipman:2000.
An isodensity SS(V)PE calculation is requested by setting SOLVENT_METHOD = ISOSVP in the $rem section, in addition to normal job control variables for a singlepoint energy calculation. Whereas the other solvation models described in this chapter use specialized input sections (e.g., $pcm) in lieu of a slew of $rem variables, the isodensity SS(V)PE code is an interface between QChem and a standalone FORTRAN code written by Chipman [664], and some $rem variables are used for job control of isodensity SS(V)PE calculations. These are listed below.
SVP_MEMORY
Specifies the amount of memory for use by the solvation module.
TYPE:
INTEGER
DEFAULT:
125
OPTIONS:
corresponds to the amount of memory in MB.
RECOMMENDATION:
The default should be fine for medium size molecules with the default Lebedev grid, only increase if needed.
SVP_PATH
Specifies whether to run a gas phase computation prior to performing the solvation procedure.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
runs a gasphase calculation and after
convergence runs the SS(V)PE computation.
1
does not run a gasphase calculation.
RECOMMENDATION:
Running the gasphase calculation provides a good guess to start the solvation stage and provides a more complete set of solvated properties.
SVP_CHARGE_CONV
Determines the convergence value for the charges on the cavity. When the change in charges fall below this value, if the electron density is converged, then the calculation is considered converged.
TYPE:
INTEGER
DEFAULT:
7
OPTIONS:
Convergence threshold set to .
RECOMMENDATION:
The default value unless convergence problems arise.
SVP_CAVITY_CONV
Determines the convergence value of the iterative isodensity cavity procedure.
TYPE:
INTEGER
DEFAULT:
10
OPTIONS:
Convergence threshold set to .
RECOMMENDATION:
The default value unless convergence problems arise.
SVP_GUESS
Specifies how and if the solvation module will use a given guess for the charges and cavity points.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
No guessing.
1
Read a guess from a previous QChem solvation computation.
2
Use a guess specified by the $svpirf section from the input
RECOMMENDATION:
It is helpful to also set SCF_GUESS to READ when using a guess from a previous QChem run.
This last $rem variable requires specification of a $svpirf input section, the format for which is the following:
$svpirf <# point> <x point> <y point> <z point> <charge> <grid weight> <# point> <x normal> <y normal> <z normal> $end
More refined control over SS(V)PE jobs is obtained using a $svp input section. These are read directly by Chipman’s SS(V)PE solvation module and and therefore must be specified in the context of a FORTRAN namelist. The format is as follows:
$svp <KEYWORD>=<VALUE>, <KEYWORD>=<VALUE>,... <KEYWORD>=<VALUE> $end
For example, the section may look like this:
$svp RHOISO=0.001, DIELST=78.39, NPTLEB=110 $end
The following keywords are supported in the $svp section:
DIELST
The static dielectric constant.
TYPE:
FLOAT
DEFAULT:
78.39
OPTIONS:
real number specifying the constant.
RECOMMENDATION:
The default value 78.39 is appropriate for water solvent.
ISHAPE
A flag to set the shape of the cavity surface.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
use the electronic isodensity surface.
1
use a spherical cavity surface.
RECOMMENDATION:
Use the default surface.
RHOISO
Value of the electronic isodensity contour used to specify the cavity surface. (Only relevant for ISHAPE = 0).
TYPE:
FLOAT
DEFAULT:
0.001
OPTIONS:
Real number specifying the density in electrons/bohr.
RECOMMENDATION:
The default value is optimal for most situations. Increasing the value produces a smaller cavity which ordinarily increases the magnitude of the solvation energy.
RADSPH
Sphere radius used to specify the cavity surface (Only relevant for ISHAPE=1).
TYPE:
FLOAT
DEFAULT:
Half the distance between the outermost atoms plus 1.4 ngstroms.
OPTIONS:
Real number specifying the radius in bohr (if positive) or in ngstroms (if negative).
RECOMMENDATION:
Make sure that the cavity radius is larger than the length of the molecule.
INTCAV
A flag to select the surface integration method.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
Single center Lebedev integration.
1
Single center spherical polar integration.
RECOMMENDATION:
The Lebedev integration is by far the more efficient.
NPTLEB
The number of points used in the Lebedev grid for the singlecenter surface integration. (Only relevant if INTCAV=0).
TYPE:
INTEGER
DEFAULT:
1202
OPTIONS:
Valid choices are:
6, 18, 26, 38, 50, 86, 110, 146, 170, 194, 302, 350, 434, 590, 770,
974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334,
4802, or 5294.
RECOMMENDATION:
The default value has been found adequate to obtain the energy to within 0.1 kcal/mol for solutes the size of monosubstituted benzenes.
NPTTHE, NPTPHI
The number of (,) points used for singlecentered surface integration (relevant only if INTCAV=1).
TYPE:
INTEGER
DEFAULT:
8,16
OPTIONS:
, specifying the number of points.
RECOMMENDATION:
These should be multiples of 2 and 4 respectively, to provide symmetry sufficient for all Abelian point groups. Defaults are too small for all but the tiniest and simplest solutes.
LINEQ
Flag to select the method for solving the linear equations that determine the apparent point charges on the cavity surface.
TYPE:
INTEGER
DEFAULT:
1
OPTIONS:
0
use LU decomposition in memory if space permits, else switch to LINEQ=2
1
use conjugate gradient iterations in memory if space permits, else use LINEQ=2
2
use conjugate gradient iterations with the system matrix stored externally on disk.
RECOMMENDATION:
The default should be sufficient in most cases.
CVGLIN
Convergence criterion for solving linear equations by the conjugate gradient iterative method (relevant if LINEQ=1 or 2).
TYPE:
FLOAT
DEFAULT:
1.0E7
OPTIONS:
Real number specifying the actual criterion.
RECOMMENDATION:
The default value should be used unless convergence problems arise.
Note that the singlecenter surface integration approach that is used to find the isodensity surface may fail for certain very nonspherical solute molecules. The program will automatically check for this, aborting with a warning message if necessary. The singlecenter approach succeeds only for what is called a “star surface”, meaning that an observer sitting at the center has an unobstructed view of the entire surface. Said another way, for a star surface any ray emanating out from the center will pass through the surface only once. Some cases of failure may be fixed by simply moving to a new center with the ITRNGR parameter described below. But some surfaces are inherently nonstar surfaces and cannot be treated with this program until more sophisticated surface integration approaches are developed and implemented.
ITRNGR
Translation of the cavity surface integration grid.
TYPE:
INTEGER
DEFAULT:
2
OPTIONS:
0
No translation (i.e., center of the cavity at the origin
of the atomic coordinate system)
1
Translate to the center of nuclear mass.
2
Translate to the center of nuclear charge.
3
Translate to the midpoint of the outermost atoms.
4
Translate to midpoint of the outermost nonhydrogen atoms.
5
Translate to userspecified coordinates in Bohr.
6
Translate to userspecified coordinates in ngstroms.
RECOMMENDATION:
The default value is recommended unless the singlecenter integrations procedure fails.
TRANX, TRANY, TRANZ
, , and value of userspecified translation (only relevant if ITRNGR is set to 5 or 6
TYPE:
FLOAT
DEFAULT:
0, 0, 0
OPTIONS:
, , and relative to the origin in the appropriate units.
RECOMMENDATION:
None.
IROTGR
Rotation of the cavity surface integration grid.
TYPE:
INTEGER
DEFAULT:
2
OPTIONS:
0
No rotation.
1
Rotate initial axes of the integration grid to coincide
with principal moments of nuclear inertia (relevant if ITRNGR=1)
2
Rotate initial axes of integration grid to coincide with
principal moments of nuclear charge (relevant if ITRNGR=2)
3
Rotate initial axes of the integration grid through userspecified
Euler angles as defined by Wilson, Decius, and Cross.
RECOMMENDATION:
The default is recommended unless the knowledgeable user has good reason otherwise.
ROTTHE ROTPHI ROTCHI
Euler angles (, , ) in degrees for userspecified rotation of the cavity surface. (relevant if IROTGR=3)
TYPE:
FLOAT
DEFAULT:
0,0,0
OPTIONS:
, , in degrees
RECOMMENDATION:
None.
IOPPRD
Specifies the choice of system operator form.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
Symmetric form.
1
Nonsymmetric form.
RECOMMENDATION:
The default uses more memory but is generally more efficient, we recommend its use unless there is shortage of memory available.
By default, QChem will check the validity of the singlecenter expansion by searching for the isodensity surface in two different ways: first, working inwards from a large distance, and next by working outwards from the origin. If the same result is obtained (within tolerances) using both procedures, then the cavity is accepted. If the two results do not agree, then the program exits with an error message indicating that the inner isodensity surface is found to be too far from the outer isodensity surface.
Some molecules, for example C, can have a hole in the middle. Such molecules have two different “legal” isodensity surfaces, a small inner one inside the “hole”, and a large outer one that is the desired surface for solvation. In such cases, the cavity check described in the preceding paragraph causes the program to exit. To avoid this, one can consider turning off the cavity check that works out from the origin, leaving only the outer cavity determined by working in from large distances.
ICVICK
Specifies whether to perform cavity check
TYPE:
INTEGER
DEFAULT:
1
OPTIONS:
0
no cavity check, use only the outer cavity
1
cavity check, generating both the inner and outer cavities and compare.
RECOMMENDATION:
Consider turning off cavity check only if the molecule has a hole and if a star (outer) surface is expected.
According to Table 11.2, COSMO and CPCM appear to differ only in the dielectric screening factor, in Eq. eq:f_epsilon. Indeed, surface charges in either model are computed according to
(11.13) 
and as discussed in Section 11.2.3 the user has the option to choose either the original value suggested by Klamt [644, 661], , or else as in, e.g., Refs. [645, 647, 654]. More importantly, however, COSMO differs from CPCM in that the former includes a correction for outlying charge that goes beyond Eq. eq:CPCM, whereas CPCM consists of nothing more than induced surface charges computed (selfconsistently) according to Eq. eq:CPCM
Upon solution of Eq. eq:CPCM, the outlying charge correction in COSMO [661, 662] is obtained by first defining a larger cavity that is likely to contain essentially all of the solute’s electron density; in practice, this typically means using atomic radii of 1.95 , where denotes the original atomic van der Waals radius that was used to compute . (Note that unlike the PCMs described in Sections 11.2.2 and 11.2.3, where the atomic radii have default values but a high degree of usercontrollability is allowed, the COSMO atomic radii are parameterized for this model and are fixed.) A new set of charges, , is then computed on this larger cavity surface, and the charges on the original cavity surface are adjusted to new values, . Finally, a corrected electrostatic potential on the original surface is computed according to . It is this potential that is used to compute the solute–continuum electrostatic interaction (polarization energy), . (For comparison, when the CPCM approach described in Section 11.2.2 is used, the electrostatic polarization energy is , computed using the original surface charges and surface electrostatic potential .) With this outlying charge correction, QChem’s implementation of COSMO resembles the one in Turbomole [701].
A COSMO calculation is requested by setting SOLVENT_METHOD = COSMO in the $rem section, in addition to normal job control variables. The keyword Dielectric in the $solvent section is used to set the solvent’s static dielectric constant, as described above for other solvation models. COSMO calculations can also be used as a starting point for “COSMORS” calculations [702, 671], where RS stands for “real solvent”. The COSMORS approach is not included in QChem and requires the COSMOtherm program, which is licensed separately through COSMOlogic [672]. QChem users who are interested in COSMOtherm can request special versions of QChem for the generation of surface files that are needed by COSMOtherm.
The SM models developed by Cramer, Truhlar, and coworkers at the University of Minnesota are a class of implicit solvation models that are designed to be “universal” in the sense that they can be applied to any solvent for which a small of descriptors is known [668]. (Note that the in SM is simply a version number; versions SM8 and SM12 are available in QChem.) In particular, the solvent descriptors are:
dielectric constant
refractive index
bulk surface tension
acidity on the Abraham scale
basicity on the Abraham scale
carbon aromaticity, which equals the fraction of nonhydrogenic solvent atoms that are aromatic carbon atoms
electronegative halogenicity, which equals the fraction of nonhydrogenic solvent atoms that are F, Cl, or Br).
These models consist of a generalized Born treatment of continuum electrostatic interactions, along with nonelectrostatic interactions that are parameterized in terms of atomic surface tensions. The nonelectrostatic interactions include cavitation, dispersion, and changes in solvent structure, and the treatment of these nonelectrostatic effects is crucial to obtaining accurate (free) energies of solvation.
The SM8 model is described in detail in Ref. Marenich:2007. It may be employed in conjunction with density functional theory (with any density functional available in QChem) or with HartreeFock theory, but is intended for use only with the 631G*, 631+G*, and 631+G** basis sets, for reasons discussed below.
Bulk (continuum) electrostatic interactions in SM8 are described in terms of a generalized Born (GB) SCRF, using a solute cavity constructed from atomcentered spheres. For the atoms H, C, N, O, F, Si, P, S, Cl, and Br, atomic radii have been specifically optimized for use with SM8, whereas for other atoms the Bondi radius is used [681], or else a value of 2.0 for atoms not included in Bondi’s paper. Geometrydependent radii are computed from these “intrinsic” Coulomb radii via a descreening approximation [658].
In addition to GB electrostatics, there are several other contributions to the SM8 standardstate free energy of solvation. The first of these is called the electronicnuclearpolarization (ENP) energy, or simply the electronic polarization (EP) energy if the solute geometry is assumed to be identical in the gas and solution phases. Another contribution to the free energy of solvation comes from shortrange interactions with solvent molecules in the first solvation shell, and is sometimes called the cavitation/dispersion/solventstructure (CDS) term. The CDS contribution to the solvation energy is a sum of terms that are each proportional (with geometrydependent proportionality constants called atomic surface tensions) to the solventaccessible surface areas (SASAs) of the individual solute atoms. The SASA of the solute molecule is the area of a surface generated by the center of a spherical effective solvent molecule rolling on the van der Waals surface of the solute molecule, as in the solventaccessible surface that was mentioned in Section 11.2.2. The SASA is computed using the Analytic Surface Area (ASA) algorithm of Ref. [703] and Bondi’s values for the van der Waals radii [681], or else a value of 2.0 if no Bondi radius is available. (Note that, as in the case of nonelectrostatic interactions in PCMs, this means that a different molecular surface is used for the bulk electrostatics as compared to the nonelectrostatic interactions.) The solvent probe radius used to generate the SASAs is set to 0.40 for all solvents. Note that the solventstructure part of the CDS term includes many aspects of solvent structure that are not described by bulk electrostatics, for example, hydrogen bonding, exchange repulsion, and the variation of the effective dielectric constant in the first solvation shell, relative to its bulk value. The semiempirical nature of the CDS term also makes up for errors due to () assuming fixed and modeldependent values of the intrinsic Coulomb radii, and () any systematic errors in the description of the solute–solvent electrostatic interactions using the GB approximation.
The final component of the SM8 solvation free energy is the concentration component. This is zero if the standardstate concentration of the solute is the same in the gas and solution phases (e.g., if it is 1 mole/liter in the gas phase as well as in the solution). Otherwise, this correction can be computed using ideal gas formulas, as discussed below.
SM8 does not require the user to assign molecular mechanics atom types to atoms or groups; all atomic surface tensions in the theory are unique and continuous functions of the solute geometry, defined by the model and calculated internally within QChem. In principle, SM8 can be used with any level of electronic structure theory so long as accurate partial charges can be computed, but QChem’s implementation of SM8 specifically uses selfconsistently polarized Charge Model 4 (CM4) class IV charges [704] CM4 charges are obtained from Löwdin population analysis charges, via a mapping whose parameters depend on the basis set (and only on the basis set, not on the density functional or anything else). Basis sets supported for SM8 calculations in QChem are:
631G*
631+G*
631+G**
The charge mapping parameters are given in Ref. Kelly:2005. Other basis sets should not be used in SM8 calculations.
The SM8 solvation free energy is output at K for a standardstate concentration of 1 M in both the gas and solution phase. However, solvation free energies in the literature are often tabulated using a standard state of atm for the gase. To convert 1 Mto1 M solvation free energies at 298 K to a standard state consisting of atm for the gas and a 1 M concentration in solution, add kcal/mol to the computed solvation free energy.
Solutionphase geometry optimizations can be carried out, but basis sets that use spherical harmonic functions, or angular momentum higher than (, , etc..) are not supported. Since, by definition, the 631G*, 631+G*, and 631+G** basis sets have Cartesian shells, they are examples of basis sets that may be used for geometry optimization with SM8. Solutionphase Hessian calculations can be carried out by numerical differentiation of analytical energy gradients or by double differentiation of energies, although the former procedure is both more stable and more economical. The analytic gradients of SM8 are based on the analytical derivatives of the polarization free energy and the analytical derivatives of the CDS terms derived in Ref. Zhu:1999.
An SM8 calculation is requested by setting SOLVENT_METHOD = SM8 in the $rem section, along with other jobcontrol variables appropriate for a HartreeFock or DFT calculation. Additional variables for SM calculations should be listed in a $smx input section; for SM8, the only additional variable that is required is the name of the solvent:
Solvent
Sets the SM solvent
INPUT SECTION: $smx
TYPE:
STRING
DEFAULT:
water
OPTIONS:
Any name from the list of solvents given below.
RECOMMENDATION:
NONE
1,1,1trichloroethane 
bromoethane 
ethylbenzoate 
1,1,2trichloroethane 
bromooctane 
ethylethanoate 
1,1dichloroethane 
butanal 
ethylmethanoate 
1,2,4trimethylbenzene 
butanoicacid 
ethylphenylketone 
1,4dioxane 
butanone 
ethylpropanoate 
1bromo2methylpropane 
butanonitrile 
ethylbutanoate 
1bromopentane 
butylethanoate 
ethylcyclohexane 
1bromopropane 
butylamine 
ethylformamide 
1butanol 
butylbenzene 
xylene 
1chloropentane 
carbon disulfide 
heptane 
1chloropropane 
carbon tetrachloride 
hexadecane 
1decanol 
chlorobenzene 
hexane 
1fluorooctane 
chlorotoluene 
nitrobenzene 
1heptanol 
1,2dimethylcyclohexane 
nitroethane 
1hexanol 
decalin 
nitromethane 
1hexene 
cyclohexane 
methylaniline 
1hexyne 
cyclohexanone 
nonane 
1iodobutane 
cyclopentane 
octane 
1iodopentene 
cyclopentanol 
pentane 
1iodopropane 
cyclopentanone 
chlorotoluene 
1nitropropane 
decane 
cresol 
1nonanol 
dibromomethane 
dichlorobenzene 
1octanol 
dibutyl ether 
nitrotoluene 
1pentanol 
dichloromethane 
xylene 
1pentene 
diethyl ether 
pentadecane 
1pentyne 
diethylsulfide 
pentanal 
1propanol 
diethylamine 
pentanoic acid 
2,2,20trifluoroethanol 
diiodomethane 
pentylethanoate 
2,2,4trimethylpentane 
dimethyldisulfide 
pentylamine 
2,4dimethylpentane 
dimethylacetamide 
perfluorobenzene 
2,4dimethylpyridine 
dimethylformamide 
phenyl ether 
2,6dimethylpyridine 
dimethylpyridine 
propanal 
2bromopropane 
dimethyl sulfoxide 
propanoic acid 
2chlorobutane 
dipropylamine 
propanonitrile 
2heptanone 
dodecane 
propylethanoate 
2hexanone 
1,2dichloroethene 
propylamine 
2methylpentane 
2pentene 
xylene 
2methylpyridine 
ethanethiol 
pyridine 
2nitropropane 
ethanol 
pyrrolidine 
2octanone 
ethylethanoate 
butanol 
2pentanone 
ethylmethanoate 
butanol 
2propanol 
ethylphenyl ether 
butylbenzene 
2propen1ol 
ethylbenzene 
tetrachloroethene 
3methylpyridine 
ethylene glycol 
tetrahydrofuran 
3pentanone 
fluorobenzene 
tetrahyrothiophenedioxide 
4heptanone 
formamide 
tetralin 
4methyl2pentanone 
formic acid 
thiophene 
4methylpyridine 
hexadecyliodide 
thiophenol 
5nonanone 
hexanoic acid 
toluene 
acetic acid 
iodobenzene 
decalin 
acetone 
iodoethane 
tribromomethane 
acetonitrile 
iodomethane 
tributylphosphate 
aniline 
isobutanol 
trichloroethene 
anisole 
isopropyl ether 
trichloromethane 
benzaldehyde 
isopropylbenzene 
triethylamine 
benzene 
isopropyltoluene 
undecane 
benzonitrile 
cresol 
water 
benzyl alcohol 
mesitylene 
1,2dichloroethene 
bromobenzene 
methanol 
other 
The final choice, Solvent = “other”, requires an additional freeformat file called “solvent_data” that should contain the floatpoint values of the following solvent descriptors:
Dielec 
dielectric constant, , of the solvent 
SolN 
index of refraction at optical frequencies at 293 K, 
SolA 
Abraham’s hydrogen bond acidity, 
SolB 
Abraham’s hydrogen bond basicity, 
SolG 
(default is 0.0), where is the macroscopic surface tension at air/solvent 
interface at 298 K, and is 1 cal mol (1 dyne/cm = 1.43932 cal mol ) 

SolC 
aromaticity, : the fraction of nonhydrogenic solvent atoms that are aromatic 
carbon atoms 

SolH 
electronegative “halogenicity”, : the fraction of nonhydrogenic solvent atoms that are 
F, Cl or Br 
For a desired solvent, these values can be derived from experiment or from interpolation or extrapolation of data available for other solvents. Solvent parameters for common organic solvents are tabulated in the Minnesota Solvent Descriptor Database. The latest version of this database is available at:
The SM8 test suite contains the following representative examples:
singlepoint solvation energy and analytical gradient calculation for 2,2dichloroethenyl dimethyl phosphate in water at the M062X/631G* level;
singlepoint solvation energy calculation for 2,2dichloroethenyl dimethyl phosphate in benzene at the M062X/631G* level;
singlepoint solvation energy calculation for 2,2dichloroethenyl dimethyl phosphate in ethanol at the M062X/631G* level;
singlepoint solvation energy calculation for 5fluorouracil in water at the M06/631+G* level;
singlepoint solvation energy calculation for 5fluorouracil in octanol at the M06L/631+G* level;
singlepoint solvation energy and analytical gradient calculation for 5fluorouracil in fluorobenzene at the M06HF/631+G** level;
geometry optimization for protonated methanol in water at the B3LYP/631G* level;
finitedifference frequency (with analytical gradient) calculation for protonated methanol in water at the B3LYP/631G* level.
Users who wish to calculate solubilities can calculate them from the free energies of solvation by the method described in Ref. Thompson:2003. The present model can also be used with confidence to calculate partition coefficients (e.g., Henry’s Law constants, octanol/water partition coefficients, etc..) by the method described in Ref. Cramer:2001.
The user should note that the free energies of solvation calculated by the SM8 model in the current version of QChem are all what may be called equilibrium free energies of solvation. The nonequilibrium algorithm required for vertical excitation energies [708] is not yet available in QChem.
The SM12 model [659] is also available in QChem. Similar to SM8, it employs (a) the generalized Born approximation for the bulk electrostatic contribution to the free energy of solvation, and (b) the same formulas (with reoptimized parameters) for CDS contributions. SM12 holds several advantages over SM8, and perhaps foremost among these is that it uses CM5 charges [523], which are based on Hirshfeld population analysis, or else charges derived from the electrostatic potential [709, 519], for the bulk electrostatics term. These charges are stable with respect to extension of the basis set, and thus SM12 can be used with larger basis sets whereas SM8 is limited to 631G*, 631+G*, and 631+G**, due to instabilities in the Löwdin charges in larger basis sets. In addition, SM12 is parameterized using a more diverse training set as compared to SM8, and is defined for the entire periodic table. However, the SM12 analytic gradient is not available in QChem at present.
An SM12 calculation is requested by setting SOLVENT_METHOD = SM12 in the $rem section. The manner in which the electrostatic term is computed is controlled by the Charges keyword in the $smx input section.
Charges
Sets the type of atomic charges for the SM12 electrostatic term.
INPUT SECTION: $smx
TYPE:
STRING
DEFAULT:
CM5
OPTIONS:
CM5
Charge Model 5 charges [523]
MK
MerzSinghKollman charges [709]
CHELPG
ChElPG charges [519]
RECOMMENDATION:
None. MerzSinghKollman and ChElPG charges are fit to reproduce the molecular electrostatic potential on the van der Waals surface or on a cubic grid, respectively, whereas CM5 is an empirical model based on Hirshfeld population analysis.
Example 11.256 SM12CM5 calculation of the solvation free energy of water in the 1octanol solvent.
$molecule 0 1 O 0.000000 0.125787 0.000000 H 0.758502 0.503148 0.000000 H 0.758502 0.503148 0.000000 $end $rem jobtype sp scf_guess core method b3lyp basis 631G* sym_ignore true solvent_method sm12 $end $smx solvent 1octanol charges chelpg $end
The SMD model [660] is also available in QChem. Within this model, the electrostatic contribution to the free energy solvation is described via the IEFPCM model, where the CDS contributions follow the formulas as SM8 and SM12 with the parameters reoptimized to be compatible with the IEFPCM electrostatics.
An SM12 energy or gradient calculation is requested by setting SOLVENT_METHOD = SMD in the $rem section. While QChem users can vary the parameters for the IEFPCM part of the SMD calculation, this should be done with caution because a modified IEFPCM electrostatics might be less compatible with CDS parameters and thus lead to less accurate results.
Example 11.257 SMD force calculation for tetrahydrofuran (THF) in the pentane solvent.
$molecule 0 1 C 0.3616581399 0.9869677370 0.2223661594 C 1.3310981828 0.1445973123 0.1083638467 O 0.5925741631 1.3541832621 0.0367381550 C 0.7980898779 1.0708993132 0.1365091554 C 0.9646828335 0.3961546963 0.2563198455 H 0.6256761390 1.9258626049 0.2670118459 H 0.3332291383 1.1581017321 1.3027531705 H 1.6975291807 0.0685183162 1.1404488482 H 2.1934121687 0.1816203090 0.5621291600 H 1.1301998060 1.2383993322 1.1698391042 H 1.3485248249 1.7543183024 0.5146978297 H 1.0506138962 0.4896466678 1.3431518473 H 1.8430658739 0.8558027094 0.1996591590 $end $rem jobtype force method b3lyp basis 631G* solvent_method smd $end $smx solvent pentane $end
QChem provides the option to calculate molecular properties in aqueous solution and the magnitudes of the hydration free energies by the Langevin dipoles (LD) solvation model developed by Jan Florián and Arieh Warshel [656, 657] at the University of Southern California. In this model, a solute molecule is surrounded by a sphere of point dipoles, with centers on a cubic lattice. Each of these “Langevin” dipoles changes its size and orientation in the electrostatic field of the solute and the other Langevin dipoles. The electrostatic field from the solute is determined rigorously by the integration of its charge density, whereas for dipole–dipole interactions, a 12 cutoff is used. The QChem/ChemSol 1.0 implementation of the LD model is fully selfconsistent in that the molecular quantum mechanical calculation takes into account solute–solvent interactions. Further details on the implementation and parameterization of this model can be found in the original literature [656, 657].
The results of ChemSol calculations are printed in the standard output file. Below is a part of the output for a calculation on the methoxide anion (corresponding to the sample input given later on, and the sample file in the $QC/samples directory).
Energy Component 
Value / kcal mol 
LD Electrostatic energy 

Hydrophobic energy 
0.28 
van der Waals energy 

Bulk correction 

Solvation free energy, (ILD) 

The total hydration free energy, (ILD) is calculated as a sum of several contributions. Note that the electrostatic part of is calculated by using the linearresponse approximation [656] and contains contributions from the polarization of the solute charge distribution due to its interaction with the solvent. This results from the selfconsistent implementation of the Langevin dipoles model within QChem.
To perform an LD calculation in QChem, specify normal jobcontrol variables for a HartreeFock or DFT calculation, and set SOLVENT_METHOD = CHEM_SOL in the $rem section. Additional finetuning is accomplished using a set of keywords in a $chem_sol input section. The remainder of this section summarizes these keywords.
EField
Determines how the solute charge distribution is approximated in evaluating the electrostatic field of the solute.
INPUT SECTION: $chemsol
TYPE:
INTEGER
DEFAULT:
1
OPTIONS:
1
Exact solute charge distribution is used.
0
Solute charge distribution is approximated by Mulliken atomic charges.
RECOMMENDATION:
None. The Mullikenbased procedure is faster but less rigorous.
NGrids
Sets the number of grids used to calculate the average hydration free energy.
INPUT SECTION: $chemsol
TYPE:
INTEGER
DEFAULT:
5
will be averaged over 5 different grids.
OPTIONS:
Use different grids.
RECOMMENDATION:
None. The maximum allowed value of is 20.
Controls printing in the ChemSol part of the QChem output file.
INPUT SECTION: $chemsol
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
Limited printout
1
Full printout
RECOMMENDATION:
None.
Accurate calculations of hydration free energies require a judicious choice of the solute–solvent boundary in terms of atomtype dependent parameters. The default atomic van der Waals radii available in QChem were chosen to provide reasonable hydration free energies for most solutes and basis sets. These parameters basically coincide with the ChemSol 2.0 radii given in Ref. Florian:1999. The only difference between the QChem and ChemSol 2.0 atomic radii stems from the fact that QChem parameter set uses radii for carbon and oxygen that are independent of the atom’s hybridization state. Userdefined atomic radii can be specified by declaring the option ReadRadii in the $chem_sol input section, and then placing the radii in the $van_der_waals section. Two different (and mutually exclusive) formats can be used, as shown below.
$van_der_waals 1 atomic_number vdW_radius ... $end $van_der_waals 2 sequential_atom_number vdW_radius ... $end
The purpose of the second format is to permit the user to customize the radius of specific atoms, in the order that they appear in the $molecule section, rather than simply by atomic numbers as in format 1. The radii of atoms that are not listed in the $van_der_waals input will be assigned default values. The atomic radii that were used in the calculation are printed in the ChemSol part of the output file in the column denoted rp. All radii should be given in ngstroms.
Example 11.258 A Langevin dipoles calculation on the methoxide anion. A customized value is specified for the radius of the C atom.
$molecule 1 1 C 0.0000 0.0000 0.5274 O 0.0000 0.0000 0.7831 H 0.0000 1.0140 1.0335 H 0.8782 0.5070 1.0335 H 0.8782 0.5070 1.0335 $end $rem METHOD hf BASIS 631G SCF_CONVERGENCE 6 SOLVENT_METHOD Chem_Sol $end $chemsol ReadRadii $end $van_der_waals 2 1 2.5 $end