• Previous: Introduction
• Up: Molecules in Complex Environments: Solvent Models, QM/MM and QM/EFP Features, Density Embedding
• Next: Stand-Alone QM/MM Calculations

11.2 Chemical Solvent Models

Ab initio quantum chemistry makes possible the study of gas-phase 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 quantum-chemical calculation (e.g. a super-molecular 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 long-range electrostatic interactions. Accurate prediction of solvation free energies is, however, crucial for modeling of chemical reactions and ligand/receptor interactions in solution.

Q-Chem contains several different implicit solvent models, which differ greatly in their level of sophistication. These are generally known as self-consistent reaction field (SCRF) models, because the continuum solvent establishes a “reaction field” (additional terms in the solute Hamiltonian) that depends upon the solute electron density, and must therefore be updated self-consistently during the iterative convergence of the wave function. The simplest and oldest of these models that is available in Q-Chem is the Kirkwood-Onsager model [722, 723, 724], in which the solute molecule is placed inside of a spherical cavity and its electrostatic potential is represented in terms of a single-center multipole expansion. More sophisticated models, which use a molecule-shaped cavity and the full molecular electrostatic potential, include the conductor-like screening model (COSMO) [725] and the closely related conductor-like PCM (C-PCM) [726, 727, 728], along with the “surface and simulation of volume polarization for electrostatics” [SS(V)PE] model [729]. The latter is also known as the “integral equation formalism” (IEF-PCM) [730, 731].

The C-PCM and IEF-PCM/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 co-workers [732], is now used almost universally to refer to this class of solvation models. Q-Chem employs a Switching/Gaussian or “SWIG” implementation of these PCMs [733, 734, 735, 736, 737]. This approach resolves a long-standing—though little-publicized—problem with standard PCMs, namely, that the boundary-element 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 [733, 734]. In contrast, Q-Chem’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 integral-equation solvent models, so that solvation energies and molecular surface areas are hardly affected by the smoothing procedure.

Other solvent models available in Q-Chem include the “Langevin dipoles” model [738, 739] as well as versions 8 and 12 of the SM models and the SMD model developed at the University of Minnesota [740, 741, 742]. SM8 and SM12 are based upon the generalized Born method for electrostatics, augmented with atomic surface tensions intended to capture non-electrostatic 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 (where the “D” is for “density") combines IEF-PCM with the non-electrostatic corrections, but because the electrostatics is based on the density rather than atomic point charges, it is supported for arbitrary basis sets whereas SM8 and SM12 are not.

Table 11.1 summarizes the implicit solvent models that are available in Q-Chem. 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. 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. Post-Hartree–Fock calculations can be performed by first running an SCF + PCM job, in which case the correlated wave function will employ MOs and Hartree-Fock energy levels that are polarized by the solvent.

Note: The job-control format for specifying implicit solvent models changed significantly starting in Q-Chem 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 Kirkwood-Onsager model (Section 11.2.1).
PCM	Use an apparent surface charge, polarizable continuum model
	(Section 11.2.2).
ISOSVP	Use the iso-density implementation of the SS(V)PE model
	(Section 11.2.5).
COSMO	Use COSMO (similar to C-PCM but with an outlying charge
	correction [743, 744]; see Section 11.2.7).
SM8	Use version 8 of the Cramer-Truhlar SM model (Section 11.2.8.1).
SM12	Use version 12 of the SM model (Section 11.2.8.2).
SMD	Use SMD (Section 11.2.8.3).
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.8.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 C-PCM and SS(V)PE/IEF-PCM (the latter two being completely equivalent [731]). 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 C-PCM but includes a correction for that part of the solute’s electron density that penetrates beyond the cavity (the so-called “outlying charge”) [743, 744]. This is discussed in Section 11.2.7.

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 atom-centered 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 [729, 745, 746]. This is an appealing, one-parameter cavity construction, although it is unclear that this construction alone is superior in its accuracy to carefully-parameterized atomic radii [747], at least not without additional, non-electrostatic terms included [748, 749, 750, 751] that are available in Q-Chem’s implementation of the isodensity version of SS(V)PE (Section 11.2.6). Moreover, analytic energy gradients are not available for the isodensity cavity construction, whereas they are available when the cavity is constructed from atom-centered spheres. One additional subtlety, which is discussed in detail in Ref. Lange:2011a, 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 [735]. (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 sub-kcal/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 [752]. To achieve comparable accuracy with IEF-PCM/SS(V)PE, non-electrostatic terms must be included [753, 748, 750]. The SM12 model does not improve upon SM8 in any statistical sense [741], 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 Mulliken-style atomic charges, and is therefore parameterized only for a few small basis sets, e.g., 6-31G*. SM12, on the other hand, uses a variety of charge schemes that are stable with respect to basis-set expansion, and can therefore be combined with any level of electronic structure theory for the solute. Like IEF-PCM, the SMD model is also applicable to any basis sets, and its accuracy is comparable to SM8 and SM12 [742]. Quantitative fluid-phase thermodynamics can also be obtained using Klamt’s COSMO-RS approach [754, 755], where RS stands for “real solvent”. The COSMO-RS approach is not  included in Q-Chem and requires the COSMOtherm program, which is licensed separately through COSMOlogic [756], but Q-Chem 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 Q-Chem. In addition, recent review articles are available for PCM methods [732], SM [752], and COSMO [757]. Formal relationships between various PCMs have been discussed in Refs. Chipman:2002a,Lange:2011a.

11.2.1 Kirkwood-Onsager Model

The simplest implicit solvation model available in Q-Chem is the Kirkwood-Onsager model [722, 723, 724], 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 [758]), 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 [759]. 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 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 Kirkwood-Onsager model, in which the solute’s electron distribution was described entirely in terms of its dipole moment, Q-Chem’s version can use multipoles of arbitrarily high order, including the Born (monopole) term [760] 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 Kirkwood-Onsager solvent model are available for Hartree-Fock, DFT, and CCSD calculations. It is often advisable to perform a gas-phase calculation of the solute molecule first, which can serve as the initial guess for a subsequent Kirkwood-Onsager implicit solvent calculation.

The Kirkwood-Onsager 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 Kirkwood-Onsager model.

Note: SCRF and CCSD combo works only in CCMAN (with CCMAN2 = FALSE).

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 Q-Chem input file, this input section is not case-sensitive.

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.268  Onsager model applied at the Hartree-Fock 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             6-31g**
   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.269  CCSD/Onsager calculation applied to 1,2-dichloroethane molecule

$comment
   1,2-dichloroethane GAUCHE Conformation
$end

$molecule
   0 1
   C     0.654133  -0.381705   0.880884
   C    -0.654133   0.381705   0.880884
   Cl    1.732259   0.087759  -0.463055
   H     1.186245  -0.166574   1.796075
   H     0.488935  -1.444440   0.805846
   Cl   -1.732259  -0.087759  -0.463055
   H    -1.186245   0.166574   1.796075
   H    -0.488935   1.444440   0.805846
$end

$rem
   METHOD           CCSD
   BASIS            6-31g**
   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.270  Kirkwood-Onsager SCRF applied to hydrogen fluoride in water, performing a gas-phase calculation first.

$molecule
   0 1
   H       0.000000     0.000000    -0.862674
   F       0.000000     0.000000     0.043813
 $end

$rem
   METHOD          HF      
   BASIS           6-31G*
$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             6-31G*
   SOLVENT_METHOD    ONSAGER
   SCF_GUESS         READ  ! read vacuum solution as a guess
$end

$solvent
   CavityRadius    2.5
$end

11.2.2 Polarizable Continuum Models

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

11.2.2.1 Formal Theory and Discussion of Different Models

The PCM literature has a long history [732] and there are several different models in widespread use; connections between these models have not always been appreciated [729, 731, 745, 735]. Chipman [729, 745] has shown how various PCMs can be formulated within a common theoretical framework; see Ref. Herbert:2016a for a pedagogical introduction. The PCM takes the form of a set of linear equations,

(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 Q-Chem.

Model	Literature	Matrix	Matrix	Scalar
	Refs.
COSMO	Klamt:1993
C-PCM	Truong:1995,Barone:1998
IEF-PCM	Chipman:2000,Cances:2001a
SS(V)PE	Chipman:2000,Chipman:2002b

Table 11.2: Definition of the matrices in Eq. eq:Kq=Rv for the various PCMs that are available in Q-Chem. The matrix consists of Coulomb interactions between the cavity charges and is the discretized version of the matrix that generates the outward-pointing normal electric field vector. (See Refs. Chipman:2002a,Chipman:2002b,Herbert:2016a for detailed definitions.) The matrix is diagonal and contains the surface areas of the cavity discretization elements, and is a unit matrix. At the level of Eq. eq:Kq=Rv, COSMO and C-PCM differ only in the dielectric screening factor , although COSMO includes an additional outlying charge correction that goes beyond Eq. eq:Kq=Rv [743, 744].

The oldest PCM is the so-called D-PCM model of Tomasi and co-workers [761], but unlike the models listed in Table 11.2, D-PCM requires explicit evaluation of the electric field normal to the cavity surface, This is undesirable, as evaluation of the electric field is both more expensive and more prone to numerical problems as compared to evaluation of the electrostatic potential. Moreover, the dependence on the electric field can be formally eliminated at the level of the integral equation whose discretized form is given in Eq. eq:Kq=Rv [729]. As such, D-PCM is essentially obsolete, and the PCMs available in Q-Chem require only the evaluation of the electrostatic potential, not the electric field.

The simplest PCM that continues to enjoy widespread use is the Conductor-Like Screening Model (COSMO) introduced by Klamt and Schüürmann [725]. Truong and Stefanovich [726] 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. [727, 728], who called the model C-PCM (for “conductor-like” 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 C-PCM. The latter value is the correct choice for a single charge in a spherical cavity (i.e., the Born ion model), although Klamt and co-workers suggest that is a better compromise, given that the Kirkwood-Onsager analytical result is for an th-order multipole centered in a spherical cavity [725, 744]. The distinction is irrelevant in high-dielectric solvents; the and values of differ by only 0.6% for water at 25C, for example. Truong [726] argues that does a better job of preserving Gauss’ Law in low-dielectric solvents, but more accurate solvation energies (at least for neutral molecules, as compared to experiment) are sometimes obtained using (Ref. Barone:1998). This result is likely highly sensitive to cavity construction, and in any case, both versions are available in Q-Chem.

Whereas the original COSMO model introduced by Klamt and Schüürmann [725] corresponds to Eq. eq:Kq=Rv with and as defined in Table 11.2, Klamt and co-workers later introduced a correction for outlying charge that goes beyond Eq. eq:Kq=Rv [743, 744]. Klamt now consistently refers to this updated model as “COSMO” [757], and we shall adopt this nomenclature as well. COSMO, with the outlying charge correction, is available in Q-Chem and is described in Section 11.2.7. In contrast, C-PCM consists entirely of Eq. 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 non-electrostatic terms can be added at the user’s discretion, as discussed below, but there is no explicit outlying charge correction in C-PCM. These and other fine-tuning details for PCM jobs are controllable via the $pcm input section that is described in Section 11.2.3.

As compared to C-PCM, a more sophisticated treatment of continuum electrostatic interactions is afforded by the “surface and simulation of volume polarization for electrostatics” [SS(V)PE] approach [729]. 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 [729], but this model is formally equivalent, at the level of integral equations, to the “integral equation formalism” (IEF-PCM) that was developed originally by Cancès et al. [730, 762]. Some difference do arise when the integral equations are discretized to form finite-dimensional matrix equations [735], and it should be noted from Table 11.2 that SS(V)PE uses a symmetrized form of the matrix as compared to IEF-PCM. The asymmetric  IEF-PCM is the recommended approach [735], although only the symmetrized version is available in the isodensity implementation of SS(V)PE that is discussed in Section 11.2.5. As with the obsolete D-PCM approach, the original version of IEF-PCM explicitly required evaluation of the normal electric field at the cavity surface, but it was later shown that this dependence could be eliminated to afford the version described in Table 11.2 [729, 731]. This version requires only the electrostatic potential, and is thus preferred, and it is this version that we designate as IEF-PCM. The C-PCM model becomes equivalent to SS(V)PE in the limit  [729, 735], which means that C-PCM must somehow include an implicit  correction for volume polarization, even if this was not by design [743]. For , numerical calculations reveal that there is essentially no difference between SS(V)PE and C-PCM results [735]. Since C-PCM is less computationally involved as compared to SS(V)PE, it is the PCM of choice in high-dielectric solvents. The computational savings relative to SS(V)PE may be particularly significant for large QM/MM/PCM jobs. For a more detailed discussion of the history of these models, see the lengthy and comprehensive review by Tomasi et al. [732]. For a briefer discussion of the connections between these models, see Refs. Chipman:2002a,Lange:2011a,Herbert:2016a.

11.2.2.2 Cavity Construction and Discretization

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

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

Having chosen a model for the cavity surface, this surface is discretized using atom-centered Lebedev grids [185, 183, 186] of the same sort that are used to perform the numerical integrations in DFT. (Discretization of the re-entrant facets of the SES is somewhat more complicated but similar in spirit [737].) 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 [733].

A long-standing (though not well-publicized) problem with the aforementioned discretization procedure is that it fails to afford continuous potential energy surfaces as the solute atoms are displaced, because certain surface grid points may emerge from, or disappear within, the solute cavity, as the atomic spheres that define the cavity are moved. This undesirable behavior can inhibit convergence of geometry optimizations and, in certain cases, lead to very large errors in vibrational frequency calculations [733]. It is also a fundamental hindrance to molecular dynamics calculations [734]. Building upon earlier work by York and Karplus [766], Lange and Herbert [733, 734, 737] 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 [734, 736]. 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 [733, 734]. The smooth discretization procedure combines a switching function with Gaussian blurring of the cavity surface charge density, and is thus known as the “Switching/Gaussian” (SWIG) implementation of the PCM.

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

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

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

11.2.2.3 Non-Equilibrium Solvation for Vertical Excitation, Ionization and Emission

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

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

(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:LR-term 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 state-specific (SS) approach takes into account the capability of a part of the solvent degrees of freedom to respond instantaneously to changes in the solute wave function upon excitation. Such an effect is not accounted for in the LR approach. In SS, a generic solvated-solute excited state 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 self-consistently with respect to the solvent’s reaction field. In case of vertical ionization rather than excitation, both the ionized and non-ionized states can be treated within a ground-state formalism. For vertical excitations, self-consistent SS models [773, 776] have been developed for various excited-state methods, including both CIS and TDDFT.

In a linear dielectric medium, the solvent polarization is governed by the electric susceptibility, , where is the frequency-dependent 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 [772]:

(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 LR-PCM, the solvent polarization is subjected to the first-order changes to the electron density (TDDFT linear density response), and thus Eq. eq:Kq=Rv-fast becomes

(11.12)

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

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

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

11.2.3 PCM Job Control

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 fine-tuning may be required, as described below.

11.2.3.1 $pcm section

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	Conductor-like PCM with .
COSMO	Original conductor-like screening model with .
IEFPCM	IEF-PCM with an asymmetric matrix.
SSVPE	SS(V)PE model, equivalent to IEF-PCM with a symmetric matrix.

RECOMMENDATION:

The IEF-PCM/SS(V)PE model is more sophisticated model than either C-PCM or COSMO, and probably more appropriate for low-dielectric solvents, but it is also more computationally demanding. In high-dielectric 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 [743, 744], one should use SOLVENT_METHOD = COSMO rather than SOLVENT_METHOD = PCM. (See Section 11.2.7.)

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 [780], 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

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 Q-Chem correspond to the consensus choice that has emerged over several decades, namely, to use vdW radii scaled by a factor of 1.2. The most widely-used set of vdW radii are those determined from crystallographic data by Bondi [765], although the radius for hydrogen was later adjusted to 1.1  [767], and radii for those main-group elements not addressed by Bondi were provided later [768]. This extended set of vdW is used by default in Q-Chem, and 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 [781]. 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 user-defined 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 [747], in which the hydrogen atoms do not get their own spheres and the heavy-atom radii are increased to compensate Finally, since the solvent molecules should not be able to penetrate all the way to the atomic vdW radii of the solute, it is traditional either to scale the atomic radii (vdW surface construction) or else to augment them with an assumed radius of a spherical solvent molecule (SAS construction), but not both.

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 Lennard-Jones 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 solvent-accessible surface.

SASradius

Form a “solvent accessible” surface with the given solvent probe radius.

INPUT SECTION: $pcm
TYPE:

FLOAT

DEFAULT:

0.0

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 solvent-accessible surfaces.

SurfaceType

Selects the solute cavity surface construction.

INPUT SECTION: $pcm
TYPE:

STRING

DEFAULT:

VDW_SAS

OPTIONS:

VDW_SAS	van der Waals or solvent-accessible surface
SES	solvent-excluded surface

RECOMMENDATION:

The vdW surface and the SAS are each comprised simply of atomic spheres and thus share a common option; the only difference is the specification of a solvent probe radius, SASradius. For a true vdW surface, the probe radius should be zero (which is the default), whereas for the SAS the atomic radii are traditionally not scaled, hence vdwScale should be set to zero (which is not the default). For the SES, only SWIG discretization is available, but this can be used with any set of (scaled or unscaled) atomic radii, or with radii that are augmented by SASradius.

Historically, discretization of the cavity surface has involved “tessellation” methods that divide the cavity surface area into finite polygonal “tesserae”. (The GEPOL algorithm [782] is perhaps the most widely-used 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, Q-Chem’s SWIG PCM implementation [733, 734, 735] uses Lebedev grids to discretize the atomic spheres. These are atom-centered grids with icosahedral symmetry, and may consist of anywhere from 26 to 5294 grid points per atomic sphere. The default values used by Q-Chem were selected based on extensive numerical tests [734, 735]. 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 Q-Chem 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 [734], although exceptions (especially where charged solutes are involved) can be found [735].

Note: 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 non-hydrogen 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 non-hydrogen 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 Q-Chem 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 third-party software, including the freely-available Visual Molecular Dynamics (VMD) program [621, 622], 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 solute-solvent 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 Kirkwood-Onsager 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 floating-point values, delineated by spaces. Uses the same coordinate system as the $molecule section.

11.2.3.2 Examples

The following example shows a very basic PCM job. The solvent dielectric is specified in the $solvent section, which is described below.

Example 11.271  A basic example of using the PCMs: optimization of trifluoroethanol in water.

$rem
   JOBTYPE         OPT
   BASIS           6-31G*
   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 Kirkwood-Onsager job, Example 11.11.1 on page *.

Example 11.272  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             6-31g**
   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 united-atom cavity. Note that a user-defined 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.273  United-atom cavity construction for ethylene.

$comment
   Benzene (in benzene), with a united-atom cavity construction
   R = 2.28 A for carbon, R = 0 for hydrogen
$end

$molecule
   0 1
   C    1.38620    0.000000   0.000000
   C    0.69310    1.200484   0.000000
   C   -0.69310    1.200484   0.000000
   C   -1.38620    0.000000   0.000000
   C   -0.69310   -1.200484   0.000000
   C    0.69310   -1.200484   0.000000
   H    2.46180    0.000000   0.000000
   H    1.23090    2.131981   0.000000
   H   -1.23090    2.131981   0.000000
   H   -2.46180    0.000000   0.000000
   H   -1.23090   -2.131981   0.000000
   H    1.23090   -2.131981   0.000000
$end

$rem
   EXCHANGE        hf
   BASIS           6-31G*
   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
  1     0.00
$end

11.2.3.3 $solvent section

The solvent for PCM calculations is specified using the $solvent section, as documented below. In addition, the $solvent section can be used to incorporate non-electrostatic 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 Kirkwood-Onsager SCRF model, and will be used again below to specify the solvent for SM calculations (Section 11.2.8); 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> <SASradius>
   SolventAtom <Number1> <Number2> <Number3> <SASradius>
   . . .
   <Keyword>  <parameter/option>
   . . .
$end

The keyword SolventAtom requires multiple parameters, whereas all other keywords require only a single parameter. In addition to any (optional) non-electrostatic parameters, the $solvent section is also used to specify the solvent’s dielectric constant. If non-electrostatic interactions are ignored, then this is the only keyword that is necessary in the $solvent section. For non-equilibrium TDDFT/C-PCM 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., zero-frequency) 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 non-electrostatic interactions currently available in Q-Chem are based on the work of Cossi et al. [783], and are computed outside of the SCF procedure used to determine the electrostatic interactions. The non-electrostatic 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 non-electrostatic interactions is highly empirical and probably needs to be considered on a case-by-case 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 solvent-accessible surface.

The following keywords (in the $solvent section) are used to define non-electrostatic parameters for PCM calculations.

NonEls

Specifies what type of non-electrostatic contributions to include.

INPUT SECTION: $solvent
TYPE:

STRING

DEFAULT:

None

OPTIONS:

Cav	Cavitation energy
Buck	Buckingham dispersion and repulsion energy from atomic number
LJ	Lennard-Jones 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
SASradius (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 non-electrostatic interactions.

Example 11.274  Optimization of trifluoroethanol in water using both electrostatic and non-electrostatic PCM interactions. OPLSAA parameters are used in the Lennard-Jones potential for dispersion and repulsion.

$rem
   JOBTYPE        OPT
   BASIS          6-31G*
   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

11.2.3.4 Job control and Examples for Non-Equilibrium Solvation

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 slow-fast charge separation in the ground state.
Excited	Do slow-fast charge separation in an excited-state.

RECOMMENDATION:

Charge separation is used in conjunction with the StateSpecific keyword in $pcm.

StateSpecific

Specifies which the state-specific method will be used.

INPUT SECTION: $pcm
TYPE:

Various

DEFAULT:

No default.

OPTIONS:

Ground	Run self-consistent SS method in the ground-state with a given slow polarization charges.
Perturb	Perform ptSS and ptLR for vertical excitations.
	The th excited-state used for charge separation (for emission).

RECOMMENDATION:

NoneqGrad

Control whether perform excited state geometry optimization in equilibrium or non-equilibrium.

INPUT SECTION: $pcm
TYPE:

NONE

DEFAULT:

No default.

OPTIONS:

RECOMMENDATION:

Specify it for non-equilibrium optimization otherwise equilibrium geometry optimization will be performed.

Example 11.275  LR-TDDFT/C-PCM low-lying 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               6-31+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.276  PCM solvation effects on the vertical excitation energies of planar DMABN using the ptSS and ptLR methods.

$molecule
   0  1
   C     0.000046   -0.000398    1.904953
   C     1.210027    0.000379    1.186051
   C     1.214640   -0.000065   -0.194515
   C     0.000164   -0.000616   -0.933832
   C    -1.214349   -0.001557   -0.194687
   C    -1.209753   -0.001846    1.185775
   H     2.151949    0.001377    1.722018
   H     2.164371    0.000481   -0.709640
   H    -2.164082   -0.002008   -0.709781
   H    -2.151763   -0.002287    1.721615
   C    -0.000227    0.001061    3.325302
   N    -0.000475    0.002405    4.484321
   N     0.000053   -0.000156   -2.297372
   C    -1.258656    0.001284   -3.036994
   H    -1.041042    0.001615   -4.102376
   H    -1.860897   -0.885647   -2.811117
   H    -1.859247    0.889133   -2.810237
   C     1.258563   -0.000660   -3.037285
   H     1.860651    0.886208   -2.810755
   H     1.859362   -0.888604   -2.811461
   H     1.040664   -0.000097   -4.102609
$end

$rem
   EXCHANGE              LRC-wPBEPBE
   OMEGA                 260
   BASIS                 6-31G*
   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.277  Aqueous phenol ionization using state-specific non-equilibrium 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
  EXCHANGE              wPBE
  CORRELATION           PBE
  LRC_DFT               1
  OMEGA                 300
  BASIS                 6-31+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
   EXCHANGE              wPBE
   CORRELATION           PBE
   LRC_DFT               1
   OMEGA                 300
   BASIS                 6-31+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.278  PCM solvation effects on the emission energy of twisted DMABN in acetonitrile.

$molecule
   0  1
   C    -0.000002   -0.000003    1.894447
   C     1.230178    0.000003    1.160765
   C     1.244862    0.000001   -0.200360
   C     0.000003   -0.000013   -0.917901
   C    -1.244845   -0.000000   -0.200354
   C    -1.230160    0.000000    1.160753
   H     2.171128    0.000013    1.704064
   H     2.180670    0.000009   -0.750140
   H    -2.180657    0.000006   -0.750120
   H    -2.171112    0.000006    1.704037
   C    -0.000023   -0.000002    3.293509
   N    -0.000011   -0.000001    4.459415
   N    -0.000008   -0.000005   -2.309193
   C     0.000001    1.253881   -3.021235
   H     0.000005    1.099076   -4.098660
   H    -0.882908    1.818975   -2.705815
   H     0.882911    1.818967   -2.705805
   C    -0.000009   -1.253886   -3.021246
   H     0.882917   -1.818966   -2.705849
   H    -0.882901   -1.818992   -2.705803
   H    -0.000038   -1.099071   -4.098670
$end

$rem
   EXCHANGE              LRC-wPBEPBE
   OMEGA                 280
   BASIS                 6-31G*
   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
   EXCHANGE              LRC-wPBEPBE
   OMEGA                 280
   BASIS                 6-31G*
   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

11.2.4 Linear-Scaling QM/MM/PCM Calculations

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 [736] 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 cost-saving options are available, including efficient pre-conditioners and matrix factorizations that speed up convergence of the CG iterations, and a fast multipole algorithm for computing the electrostatic interactions [784]. 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 [736]. Currently, these options are available only for C-PCM, not for SS(V)PE/IEF-PCM.

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 [734, 735] as well as the reference for the CG solver [736].

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 single-point 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 well-behaved 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 pre-conditioner for the CG solver.

INPUT SECTION: $pcm
TYPE:

None

DEFAULT:

Off

OPTIONS:

No options. Specify the keyword to enable pre-conditioning.

RECOMMENDATION:

A Jacobi block-diagonal pre-conditioner 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. Pre-conditioning 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 linear-scaling 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 C-PCM with SWIG discretization.

11.2.5 Iso-Density Implementation of SS(V)PE

11.2.5.1 Basic Job Control

As discussed above, results obtained various types of PCMs are quite sensitive to the details of the cavity construction. Q-Chem’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 widely-used default values are provided). An alternative implementation of the SS(V)PE solvation model is also available [746], which attempts to further eliminate empiricism associated with cavity construction by taking the cavity surface to be a specified iso-contour 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 [746, 735].) In this case, the cavity surface is discretized by projecting a single-center Lebedev grid onto the iso-contour surface. Unlike the PCM implementation discussed in Section 11.2.2, for which point-group 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 mono-substituted benzenes.

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 atom-centered spheres. As with the PCMs discussed in that section, the solute may be described using Hartree-Fock theory or DFT; post-Hartree–Fock correlated wave functions can also take advantage of molecular orbitals that are polarized using SS(V)PE. Researchers who use the iso-density SS(V)PE feature are asked to cite Ref. Chipman:2000.

In related work, Pomogaeva and Chipman[748, 749, 750, 751] recently introduced a “composite method for implicit representation of solvent” (CMIRS) that is based on SS(V)PE electrostatics but adds non-electrostatic terms. This model is available in Q-Chem [785] and is discussed in Section 11.2.6. In its current implementation, CMIRS requires an iso-density SS(V)PE calculation, However, the current implementation computes the non-electrostatic interactions using the cavity and the solute’s charge density generated from the iso-density SS(V)PE. To use the CMIRS model, an iso-density SS(V)PE calculation must be requested (as described below), and the IDEFESR must be set to 1 in the $svp input section (see Section 11.2.5.2).

An iso-density SS(V)PE calculation is requested by setting SOLVENT_METHOD = ISOSVP in the $rem section, in addition to normal job control variables for a single-point 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 iso-density SS(V)PE code is an interface between Q-Chem and a code written by Chipman [746], so 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 gas-phase calculation and after
	convergence runs the SS(V)PE computation.
1	does not run a gas-phase calculation.

RECOMMENDATION:

Running the gas-phase 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 iso-density 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 Q-Chem 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 Q-Chem 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

11.2.5.2 The $svp Input Section

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 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.

INPUT SECTION: $svp
TYPE:

FLOAT

DEFAULT:

78.39

OPTIONS:

real number specifying the constant.

RECOMMENDATION:

The default value 78.39 is appropriate for water solvent.

IDEFESR

Specifies whether to request a CMIRS calculation.

INPUT SECTION: $svp
TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

0	do not invoke a CMIRS calculation.
1	do invoke a CMIRS calculation.

RECOMMENDATION:

ISHAPE

A flag to set the shape of the cavity surface.

INPUT SECTION: $svp
TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

0	use the electronic iso-density surface.
1	use a spherical cavity surface.

RECOMMENDATION:

Use the default surface.

RHOISO

Value of the electronic iso-density contour used to specify the cavity surface. (Only relevant for ISHAPE = 0.)

INPUT SECTION: $svp
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.)

INPUT SECTION: $svp
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.

INPUT SECTION: $svp
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 single-center surface integration. (Only relevant if INTCAV = 0.)

INPUT SECTION: $svp
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 mono-substituted benzenes.

NPTTHE, NPTPHI

The number of (,) points used for single-centered surface integration (relevant only if INTCAV = 1.)

INPUT SECTION: $svp
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.

INPUT SECTION: $svp
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.)

INPUT SECTION: $svp
TYPE:

FLOAT

DEFAULT:

1.0E-7

OPTIONS:

Real number specifying the actual criterion.

RECOMMENDATION:

The default value should be used unless convergence problems arise.

Note that the single-center surface integration approach that is used to find the iso-density surface may fail for certain very non-spherical solute molecules. The program will automatically check for this, aborting with a warning message if necessary. The single-center 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 non-star 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.

INPUT SECTION: $svp
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 non-hydrogen atoms.
5	Translate to user-specified coordinates in Bohr.
6	Translate to user-specified coordinates in ngstroms.

RECOMMENDATION:

The default value is recommended unless the single-center integrations procedure fails.

TRANX, TRANY, TRANZ

, , and value of user-specified translation (only relevant if ITRNGR is set to 5 or 6).

INPUT SECTION: $svp
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.

INPUT SECTION: $svp
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 user-specified
	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 user-specified rotation of the cavity surface (relevant if IROTGR = 3).

INPUT SECTION: $svp
TYPE:

FLOAT

DEFAULT:

0,0,0

OPTIONS:

, , in degrees

RECOMMENDATION:

None.

IOPPRD

Specifies the choice of system operator form.

INPUT SECTION: $svp
TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

0	Symmetric form.
1	Non-symmetric 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, Q-Chem will check the validity of the single-center expansion by searching for the iso-density 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 iso-density surface is found to be too far from the outer iso-density surface.

Some molecules, for example C, can have a hole in the middle. Such molecules have two different “legal” iso-density 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

INPUT SECTION: $svp
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.

11.2.6 Composite Method for Implicit Representation of Solvent (CMIRS)

Whereas PCMs, including sophisticated ones like SS(V)PE and IEF-PCM, account for long-range solute–solvent interactions, an accurate model for free energies of solvation must also include a treatment of short-range, non-electrostatic interactions. Various models decompose these interactions in different ways, but usually the non-electrostatic terms attempt to model all or most of the following: solute–solvent dispersion (van der Waals) interactions, Pauli (exchange) repulsion between solute and solvent, the work associated with forming the solute cavity within the dielectric medium (the so-called “cavitation energy”), hydrogen-bonding and other specific interactions due to the molecular structure of the solvent, and changes in the structure (and therefore the entropy) of the neat solvent upon introduction of the solute

Pomogaeva and Chipman[748, 749, 750, 751] have introduced an implicit solvation model that attempts to model these non-electrostatic interactions alongside a PCM-style treatment of the bulk electrostatics. They call this approach the “composite method for implicit representation of solvent” (CMIRS), and it consists first of a self-consistent treatment of solute–continuum electrostatics using the SS(V)PE model (Section 11.2.5). To this electrostatics calculation, CMIRS adds a solute–solvent dispersion term that is modeled upon the non-local VV09 van der Waals dispersion density functional [34], a Pauli repulsion contribution that depends upon the tail of the solute’s electron density that extends beyond the solute cavity, and a hydrogen-bonding correction based on the maximum and minimum values of the normal component of the electric field generated by the solute at the cavity surface. The Gibbs free energy of solvation is thus modeled as

(11.13)

where is the continuum electrostatics contribution from the SS(V)PE model, which is based on a solute cavity defined as an isocontour of the solute’s charge density. The second term contains the short-range dispersion, exchange, and "field-extremum short-range" (DEFESR) interactions:

(11.15)

These terms are evaluated only once, using a converged charge density for the solute from a SS(V)PE calculation.

The CMIRS approach was implemented in Q-Chem by Z.-Q. You and J. M. Herbert [785]. In the course of this work, a serious error was discovered in the the original implementation of by Pomogaeva and Chipman, in the gamess program. Although reparameterization of a corrected version of the model leads to only small changes in the overall error statistics across a large database of experimental free energies of solvation, the apportionment between energy components in Eq. eq:DEFESR-free-energy changes significantly [785]. By request of Dan Chipman, the Q-Chem implementation is termed “CMIRS v. 1.1", reserving v. 1.0 for the original gamess implementation of Pomogaeva and Chipman.

The CMIRS model is independently parametrized for each solvent of interest, but uses no more than five empirical parameters per solvent. It is presently available in Q-Chem for water, acetonitrile, dimethyl sulfoxide, benzene, and cyclohexane. Error statistics for compare very favorably to those of the SM models that are described in Section 11.2.8, e.g., mean unsigned errors  kcal/mol in benzene and cyclohexane and  kcal/mol in water. The latter statistic includes challenging ionic solutes; errors for charge-neutral aqueous solutes are smaller still [785].

Solvent	SolvRho	A	B	C	D	Gamma
Benzene
Cyclohexane
DMSO					-162.07
					-0.33914
Water					-21.405

Table 11.3: Optimized CMIRS parameters (from Ref. You:2016) using RHOISO = 0.0010 a.u. and, for ions, the proton reference value from Ref. Kelly:2007.

Solvent	SolvRho	A	B	C	D	Gamma
Benzene
Cyclohexane
DMSO					-817.93
					-0.44492
Water					-30.503

Table 11.4: Optimized CMIRS parameters (from Ref. You:2016) using RHOISO = 0.0005 a.u. and, for ions, the proton reference value from Ref. Kelly:2007.

The current implementation of CMIRS in Q-Chem computes the electrostatic energy using Chipman’s iso-density SS(V)PE module. The resulting iso-density cavity and the solute charge density are then employed in the calculation of the DEFESR interactions. To request a CMIRS calculation, users must set IDEFESR = in an iso-density SS(V)PE calculation (see Section 11.2.5.2). The solvent-dependent empirical parameters A, B, C, D, Gamma in the CMIRS model need to be specified in the $pcm_nonels section. Three additional parameters are also required. One is the damping parameter Delta in the dispersion equation. We recommend the parameter fixed at a.u. (about 3.7 ), an optimized value that only considers dispersion at intermolecular distances larger than van der Waals contact distance. The second is solvent’s average electron density SolvRho. The last one is the number of Gauss–Laguerre points GauLag_N for the integration over the solvent region in the exchange equation. We recommend grid points for efficient integration with accuracy. Optimized parameters for the supported solvents are listed in Tables 11.3 and 11.4 for two different values of  [785].

Example 11.279  CMIRS calculation for methane in the water solvent.

$rem
   EXCHANGE        = b3lyp5
   BASIS           = 6-31+g*
   SCF_CONVERGENCE = 8
   MEM_TOTAL       = 4000
   MEM_STATIC      = 400
   SYM_IGNORE      = true
   XC_GRID         = 000096001202
$end

$molecule
   0  1
   C     0.000000     0.000005     0.000000
   H     0.748558     0.801458     0.000000
   H     0.506094    -0.972902     0.000000
   H    -0.627326     0.085706     0.895425
   H    -0.627326     0.085706    -0.895425
$end

@@@

$rem
   EXCHANGE        = b3lyp5
   BASIS           = 6-31+g*
   SCF_CONVERGENCE = 8
   MAX_SCF_CYCLES  = 100
   SOLVENT_METHOD  = isosvp
   SCF_GUESS       = read
   PCM_PRINT       = 1
   MEM_TOTAL       = 4000
   MEM_STATIC      = 400
   SVP_MEMORY      = 1000
   SYM_IGNORE      = true
   XC_GRID         = 000096001202
$end

$molecule
   read
$end

$svp
   RHOISO=0.001, DIELST=78.36, NPTLEB=1202,
   ITRNGR=2, IROTGR=2, IPNRF=1, IDEFESR=1
$end

$pcm_nonels
   A         -0.006736
   B          0.032698
   C      -1249.6
   D        -21.405
   Delta    7.0
   Gamma    3.7
   SolvRho  0.05
   GauLag_N 40
$end

11.2.7 COSMO

According to Table 11.2, COSMO and C-PCM 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.17)

and as discussed in Section 11.2.3 the user has the option to choose either the original value suggested by Klamt [725, 743], , or else as in, e.g., Refs. Truong:1995,Cossi:2003,Lange:2011a. More importantly, however, COSMO differs from C-PCM in that the former includes a correction for outlying charge that goes beyond Eq. eq:C-PCM, whereas C-PCM consists of nothing more than induced surface charges computed (self-consistently) according to Eq. eq:C-PCM

Upon solution of Eq. eq:C-PCM, the outlying charge correction in COSMO [743, 744] 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 user-controllability 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 C-PCM 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, Q-Chem’s implementation of COSMO resembles the one in Turbomole [787].

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 “COSMO-RS” calculations [788, 755], where RS stands for “real solvent”. The COSMO-RS approach is not  included in Q-Chem and requires the COSMOtherm program, which is licensed separately through COSMOlogic [756]. Q-Chem users who are interested in COSMOtherm can request special versions of Q-Chem for the generation of -surface files that are needed by COSMOtherm.

11.2.8 SM8, SM12, and SMD Models

The SM models developed by Cramer, Truhlar, and co-workers 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 [752]. (Note that the in SM is simply a version number; SM8, SM12, and SMD are available in Q-Chem.) 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 non-hydrogenic solvent atoms that are aromatic carbon atoms

• electronegative halogenicity, which equals the fraction of non-hydrogenic solvent atoms that are F, Cl, or Br).

These models consist of a generalized Born treatment of continuum electrostatic interactions, along with non-electrostatic interactions that are parameterized in terms of atomic surface tensions. The non-electrostatic interactions include cavitation, dispersion, and changes in solvent structure, and the treatment of these non-electrostatic effects is crucial to obtaining accurate (free) energies of solvation.

11.2.8.1 The SM8 Model

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 Q-Chem) or with Hartree-Fock theory, but is intended for use only with the 6-31G*, 6-31+G*, and 6-31+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 atom-centered 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 [765], or else a value of 2.0  for atoms not included in Bondi’s paper. Geometry-dependent radii are computed from these “intrinsic” Coulomb radii via a de-screening approximation [740].

In addition to GB electrostatics, there are several other contributions to the SM8 standard-state free energy of solvation. The first of these is called the electronic-nuclear-polarization (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 short-range interactions with solvent molecules in the first solvation shell, and is sometimes called the cavitation/dispersion/solvent-structure (CDS) term. The CDS contribution to the solvation energy is a sum of terms that are each proportional (with geometry-dependent proportionality constants called atomic surface tensions) to the solvent-accessible 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 solvent-accessible surface that was mentioned in Section 11.2.2. The SASA is computed using the Analytic Surface Area (ASA) algorithm of Ref. Liotard:1995 and Bondi’s values for the van der Waals radii [765], or else a value of 2.0  if no Bondi radius is available. (Note that, as in the case of non-electrostatic interactions in PCMs, this means that a different molecular surface is used for the bulk electrostatics as compared to the non-electrostatic interactions.) The solvent probe radius used to generate the SASAs is set to 0.40  for all solvents. Note that the solvent-structure 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 semi-empirical nature of the CDS term also makes up for errors due to () assuming fixed and model-dependent 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 standard-state 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 Q-Chem. In principle, SM8 can be used with any level of electronic structure theory so long as accurate partial charges can be computed, but Q-Chem’s implementation of SM8 specifically uses self-consistently polarized Charge Model 4 (CM4) class IV charges [790] 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 Q-Chem are:

• 6-31G*

• 6-31+G*

• 6-31+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 standard-state 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 M-to1 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.

Solution-phase 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 6-31G*, 6-31+G*, and 6-31+G** basis sets have Cartesian shells, they are examples of basis sets that may be used for geometry optimization with SM8. Solution-phase 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 job-control variables appropriate for a Hartree-Fock 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

The final choice, Solvent = “other”, requires an additional free-format file called “solvent_data” that should contain the float-point values of the following solvent descriptors:

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:

• single-point solvation energy and analytical gradient calculation for 2,2-dichloroethenyl dimethyl phosphate in water at the M06-2X/6-31G* level;

• single-point solvation energy calculation for 2,2-dichloroethenyl dimethyl phosphate in benzene at the M06-2X/6-31G* level;

• single-point solvation energy calculation for 2,2-dichloroethenyl dimethyl phosphate in ethanol at the M06-2X/6-31G* level;

• single-point solvation energy calculation for 5-fluorouracil in water at the M06/6-31+G* level;

• single-point solvation energy calculation for 5-fluorouracil in octanol at the M06-L/6-31+G* level;

• single-point solvation energy and analytical gradient calculation for 5-fluorouracil in fluorobenzene at the M06-HF/6-31+G** level;

• geometry optimization for protonated methanol in water at the B3LYP/6-31G* level;

• finite-difference frequency (with analytical gradient) calculation for protonated methanol in water at the B3LYP/6-31G* 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 Q-Chem are all what may be called equilibrium free energies of solvation. The non-equilibrium algorithm required for vertical excitation energies [794] is not yet available in Q-Chem.

11.2.8.2 The SM12 Model

The SM12 model [741] is also available in Q-Chem. 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 re-optimized parameters) for CDS contributions. SM12 holds several advantages over SM8, and perhaps foremost among these is that it uses CM5 charges [594], which are based on Hirshfeld population analysis, or else charges derived from the electrostatic potential [795, 590], 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 6-31G*, 6-31+G*, and 6-31+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 Q-Chem 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 [594]
MK	Merz-Singh-Kollman charges [795]
CHELPG	ChElPG charges [590]

RECOMMENDATION:

None. Merz-Singh-Kollman 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.280  SM12CM5 calculation of the solvation free energy of water in the 1-octanol 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
   SCF_GUESS        core
   METHOD           b3lyp
   BASIS            6-31G*
   SYM_IGNORE       true
   SOLVENT_METHOD   sm12
$end

$smx
   solvent 1octanol
   charges chelpg
$end

11.2.8.3 The SMD Model

The SMD model [742] is also available in Q-Chem. Within this model, the electrostatic contribution to the free energy solvation is described via the IEF-PCM model, where the CDS contributions follow the formulas as SM8 and SM12 with the parameters re-optimized to be compatible with the IEF-PCM electrostatics. Relative to SM8 or SM12, where the electrostatic interactions are defined in terms of atomic point charges that are sensitive to the choice of basis set (and therefore only certain basis sets are supported for use with these models), SMD can be used with any basis set.

An SMD energy or gradient calculation is requested by setting SOLVENT_METHOD = SMD in the $rem section. While Q-Chem users can vary the parameters for the IEF-PCM part of the SMD calculation, this should be done with caution because a modified IEF-PCM electrostatics might be less compatible with CDS parameters and thus lead to less accurate results.

Example 11.281  SMD force calculation for tetrahydrofuran (THF) in the pentane solvent.

$molecule
   0 1
   C    -0.361658    -0.986967     0.222366
   C    -1.331098     0.144597    -0.108363
   O    -0.592574     1.354183     0.036738
   C     0.798089     1.070899     0.136509
   C     0.964682    -0.396154    -0.256319
   H    -0.625676    -1.925862    -0.267011
   H    -0.333229    -1.158101     1.302753
   H    -1.697529     0.068518    -1.140448
   H    -2.193412     0.181620     0.562129
   H     1.130199     1.238399     1.169839
   H     1.348524     1.754318    -0.514697
   H     1.050613    -0.489646    -1.343151
   H     1.843065    -0.855802     0.199659
$end

$rem
   JOBTYPE          force
   METHOD           b3lyp
   BASIS            6-31G*
   SOLVENT_METHOD   smd 
$end

$smx
   solvent pentane
$end

11.2.9 Langevin Dipoles Model

Q-Chem 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 [738, 739] 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 Q-Chem/ChemSol 1.0 implementation of the LD model is fully self-consistent 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 [738, 739].

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

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 linear-response approximation [738] and contains contributions from the polarization of the solute charge distribution due to its interaction with the solvent. This results from the self-consistent implementation of the Langevin dipoles model within Q-Chem.

To perform an LD calculation in Q-Chem, specify normal job-control variables for a Hartree-Fock or DFT calculation, and set SOLVENT_METHOD = CHEM_SOL in the $rem section. Additional fine-tuning 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 Mulliken-based 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.

Print

Controls printing in the ChemSol part of the Q-Chem 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 atom-type dependent parameters. The default atomic van der Waals radii available in Q-Chem 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 Q-Chem and ChemSol 2.0 atomic radii stems from the fact that Q-Chem parameter set uses radii for carbon and oxygen that are independent of the atom’s hybridization state. User-defined 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.282  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             6-31G
   SCF_CONVERGENCE   6
   SOLVENT_METHOD    Chem_Sol
$end

$chemsol
ReadRadii
$end

$van_der_waals
2
1 2.5
$end

11.2.10 Poisson Equation Boundary Conditions

11.2.10.1 Theory

Each of the implicit solvation model described above is designed for isotropic, bulk solvation, which can be qualitatively described by a scalar dielectric constant, . For an anisotropic  environment, such as a liquid/vapor interface, a more general approach is to employ a spatially-varying dielectric function , and to then solve Poisson’s equation for the total electrostatic potential, in three-dimensional space:

(11.18)

The quantity on the right is the solute’s charge density arising from its electronic and nuclear components,

(11.19)

The total electrostatic potential on the left side of Eq. eq:Poisson is comprised of the solute potential as well as the induced polarization potential:

(11.20)

To compute , the electronic electrostatic potential in vacuum is first evaluated on a Cartesian grid via

(11.21)

where P is the one-electron density matrix, is an atom-centered Gaussian basis function, and the subscript signifies that the integral is performed over the Cartesian grid coordinates. is then computed as

(11.22)

where the Laplacian operator is discretized using an order eight finite difference scheme. is constructed classically: for nucleus , is added to the nearest Cartesian grid point, where is the nuclear charge of and is the volume element of the Cartesian grid. This is equivalent to smearing the nuclear charge over one voxel of the Cartesian grid. Formally,

(11.23)

where is the position vector of the Cartesian grid point nearest nucleus . The exact nuclear electrostatic potential in vacuum is computed as

(11.24)

The Dirac delta function in Eq. eq:rho_nuc has dimensions of volume and incorporates the factor of .

Having computed the electronic and nuclear charge densities and electrostatic potentials in vacuum, which exactly solves Eq. eq:Poisson for , the induced polarization charge density and potential are iteratively solved following closely the procedure outlined in Refs. Fisicaro:2016 and Andreussi:2012膔. Eq. eq:Poisson is rewritten to appear as Poisson’s equation in vacuum with an additional source charge density:

(11.25)

where the polarization charge density is

(11.28)

The iterative charge density takes the form

(11.29)

and is nonzero only in regions where the dielectric function is inhomogenous. The second term on the right side of Eq. eq:rho_pol is a scaled solute charge density that is exactly zero in regions of vacuum. is iterated to convergence by performing updates using Eq. eq:rho_iter1. Note that for iteration , in Eq. eq:rho_iter1 is the solute electrostatic potential. Each time is updated, a new polarization charge density is generated and Eq. eq:Poisson2 is solved to obtain a new total electrostatic potential. The updates are continued until the Euclidean norm of the residual vector falls below a threshold, taken to be 10 a.u. by default. To stabilize the updates of , a damping procedure is performed as follows:

(11.30)

where as in Refs. Fisicaro:2016 and Andreussi:2012膔.

Once the total charge density and electrostatic potential are iterated to self-consistently, the free energy of solvation is computed as

(11.33)

where

(11.34)

Since and are computed self-consistently at each SCF iteration, the polarization effects must be incorporated into the Kohn-Sham Hamiltonian. This is accomplished by adding the correction term to the Fock matrix, equivalent to

(11.35)

in matrix form. For additional details and a full description of the algorithm, consult Ref. Coons:2017 or the Supporting Information of Ref. Coons:2016.

11.2.10.2 Non-Equilibrium Solvation for Vertical Ionization

Non-equilibrium solvation within the framework of PCMs is discussed in Section 11.2.2.3 and references therein, and that methodology has been adapted for vertical ionization processes in conjunction with Poisson equation boundary conditions. A state-specific approach has been implemented, which requires the solution of a nonlinear Schrödinger equation given by Eq. eq:SS-sdeq. For the reference state, corresponding to in Eq. eq:SS-sdeq, this is equivalent to solving Eq. eq:Poisson for the induced polarization potential and computing the solvation free energy and Fock matrix correction according to Eqs. eq:Gpol and eq:fock, respectively. For the ionized state, , the solvent polarization must be partitioned into fast and slow components. Within the Marcus partitioning scheme, Eq. eq:Poisson2 becomes[798]:

(11.36)

where the fast polarization charge density, , of the ionized state is given by

(11.37)

and the fast polarization potential, , is

(11.38)

The superscript notation indicates the component of solvent polarization, and the additional subscript index indicates from which state (reference or ionized) the quantity derives. In Eq. eq:rho_pol_neq1, is the optical dielectric of the solvent, and is computed using Eq. eq:rho_iter2 with and . For the Marcus partitioning scheme, the slow component of the induced solvent polarization charge density is computed as the three-dimensional analog of Eq. slow_q, which utilizes the two-dimensional apparent surface charges.

11.2.10.3 Poisson Equation Job Control

A calculation utilizing the Poisson equation boundary conditions is initialized by setting SOLVE_PEQ = TRUE in the $rem section. For computational efficiency, solvent effects from the Poisson equation solver are not computed until the error of the vacuum SCF calculation falls below a certain threshold, taken to be  a.u. This is controlled by the user through the PEQ_SWITCH keyword, which is also set in the $rem section.

Note: When SOLVE_PEQ = TRUE, the keywords SYM_IGNORE and NO_REORIENT are also set to TRUE.

SOLVE_PEQ

Perform a solvation free energy calculation on a Cartesian grid utilizing Poisson equation boundary conditions.

TYPE:

STRING

DEFAULT:

False

OPTIONS:

TRUE/FALSE

RECOMMENDATION:

None.

PEQ_SWITCH

Inclusion of solvent effects begins when the SCF error falls below 10.

TYPE:

INTEGER

DEFAULT:

3

OPTIONS:

Corresponding to 10

RECOMMENDATION:

Use the default unless solvent effects need to be incorporated at earlier SCF cycles.

11.2.10.4 $peqs_grid and $peqs sections

Job control for the Poisson equation solver routines is set up in two sections: $peqs_grid and $peqs. The $peqs_grid section is where the user will specify the Cartesian grid utilized in the Poisson equation routines as follows:

$peqs_grid
   DimX <Nx> <Xmin> <Xmax>
   DimY <Ny> <Ymin> <Ymax>
   DimZ <Nz> <Zmin> <Zmax>
$end

Note: An order four multi-grid is implemented to speed up convergence of the conjugate gradient algorithm employed to solve Eq. eq:Poisson2. Therefore an odd number of grid points must be used with the additional constraint that the quantities Nx-1, Ny-1, and Nz-1 must all be divisible by 8.

The values in the last two columns, specifying the minimum and maximum values of the Cartesian grid for each coordinate, are assumed to be in ngstroms. The length of the grid along each axis is computed as, for example, Xmax - Xmin, and the grid spacings are defined as . The volume element of the Cartesian grid, , surrounding the discussion of Eq. eq:rho_nuc, is .

The $peqs section controls features such as the type of cavity construction of the solute and other required parameters for interfacial or non-equilibrium solvation. The format of the $peqs section follows that of the $rem section:

$peqs
   <Keyword>  <parameter/option>
$end

PEQMaxIter

Sets the maximum number of iterations used in the conjugate gradient solver.

INPUT SECTION: $peqs
TYPE:

INTEGER

DEFAULT:

50

OPTIONS:

User-defined.

RECOMMENDATION:

Use the default unless the system is slowly converging.

PEQSolverThresh

The electrostatic potential is considered converged when the error falls below 10.

INPUT SECTION: $peqs
TYPE:

INTEGER

DEFAULT:

5

OPTIONS:

Corresponding to 10

RECOMMENDATION:

Use the default unless a tighter convergence criterion is desired at a higher computational cost.

PolarIterScale

A mixing parameter used to stabilize the iterative method when computing the polarization charge density.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

0.6

OPTIONS:

Desired value of the mixing parameter (unit-less).

RECOMMENDATION:

Use the default, which was tested and used in Refs. Fisicaro:2016 and Andreussi:2012, unless the system is difficult to converge.

BatchSize

Improves the efficiency of obtaining the electronic electrostatic potential on the grid by computing atomic orbital integrals in batches (i.e. performing the necessary integrals on small parts of the Cartesian grid instead of in its entirety).

INPUT SECTION: $peqs
TYPE:

INTEGER

DEFAULT:

5000

OPTIONS:

Corresponding to a batch size of grid points.

RECOMMENDATION:

For large grids ( grid points), a batch size is recommended.

SolventDielectric

Sets the dielectric value of the solvent.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

78.39

OPTIONS:

Desired value of (unit-less).

RECOMMENDATION:

Consult the literature to determine the value corresponding to the desired solvent and temperature. The default value corresponds to water at 25C.

OpticalDielectric

Sets the optical dielectric value of the solvent.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

1.7778

OPTIONS:

Desired value of (unit-less).

RECOMMENDATION:

Consult the literature to determine the value corresponding to the desired solvent and temperature. Note that = where is the solvent index of refraction. The default value corresponds to water at 25C.

SoluteCavity

Specifies the type of solute cavity, which determines the form of .

INPUT SECTION: $peqs
TYPE:

STRING

DEFAULT:

Vacuum

OPTIONS:

Vacuum	Perform calculation in vacuum ().
RigidVDW	Create a cavity using spherically symmetric error functions.
Spherical	Create a spherical cavity around the solute.

RECOMMENDATION:

None.

In conventional PCMs, the solute cavity is a rigid, two-dimensional surface constructed from atomic spheres, which are generally taken to be the van der Waals radius () scaled by a factor of 1.2 (see discussion in Section 11.2.3). The cavity is rigid in the sense that once the atomic coordinates are specified, it remains unchanged during the SCF cycles. Furthermore, there is an abrupt and discontinuous change in the dielectric at the solute/solvent boundary. A three-dimensional analog of this cavity, using continuous and differentiable spherically symmetric error functions, has been implemented for the Poisson equation solver, and follows the procedure in Ref. Fisicaro:2016:

(11.39)

where

(11.40)

Eq. eq:rigidvdw smoothly interpolates the dielectric from vacuum to solvent () over a length scale of , which is controlled with the RigidScale keyword. In Eqs. eq:rigidvdw and eq:h_func, is the vector of Cartesian grid coordinates, specifies the atomic coordinates of the solute, and sets the center of the error function for atom . By default, is set to be scaled by a factor of 1.2, but this can be controlled by the user with the VDWType, VDWScale, and VDWShift keywords described below.

RigidScale

Sets the length scale on which the error function employed in the RigidVDW solute cavity construction transitions the dielectric from vacuum to solvent.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

0.265 ngstroms

OPTIONS:

Specifies the desired scaling factor in ngstroms for the error functions.

RECOMMENDATION:

Use the default value, which was tuned in Ref. Fisicaro:2016 and yielded 1.0 kcal/mol error for a variety of polar and non-polar molecules, or consult the literature.

VDWType

Specifies details for the RigidVDW cavity construction.

INPUT SECTION: $peqs
TYPE:

STRING

DEFAULT:

Scaled

OPTIONS:

Unscaled
Scaled
Shifted

RECOMMENDATION:

None. The values of and are set with the VDWScale and VDWShift keywords.

VDWScale

Sets the empirical scale factor applied to the atomic van der Waals radius.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

1.2

OPTIONS:

Specifies the desired scaling factor (unit-less) of the atomic van der Waals radius.

RECOMMENDATION:

Use the default value.

VDWShift

Adjusts the center of the spherically symmetric error functions when constructing the RigidVDW cavity.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

0.0 ngstroms

OPTIONS:

Specifies the desired shift in ngstroms.

RECOMMENDATION:

None. If VDWType is set to Shifted, the van der Waals scale factor is set to 1.0 by default. This can be adjusted by the user with the VDWScale keyword.

Setting the SoluteCavity keyword to Spherical requests the construction of a spherical cavity around the solute, and the dielectric is smoothly interpolated from vacuum to solvent using a hyperbolic tangent function as in Ref Coons:2016:

(11.41)

is where the dielectric assumes the value . The value of is taken to be the sum of the sphere radius, , and half of the interpolation length, : . The sphere radius and interpolation length are controlled with the SphereRadius and InterpolLength keywords described below. The parameter controls the sharpness of the switching process and is approximated as .

SphereRadius

Sets the radius of the spherical solute cavity.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

No default.

OPTIONS:

Desired spherical cavity radius in ngstroms.

RECOMMENDATION:

None. See the Supporting Information of Ref. Coons:2016 for more information.

InterpolLength

Sets the length scale on which the dielectric is smoothly interpolated from vacuum to solvent.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

No default.

OPTIONS:

Desired interpolation length in ngstroms.

RECOMMENDATION:

None.

InterpolScale

For a given interpolation length (), InterpolScale () sets the sharpness of the dielectric transition from vacuum to solvent.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

= 4.0

OPTIONS:

Desired interpolation scale factor in ngstroms.

RECOMMENDATION:

Use the default unless a broader (smaller ) or narrower (larger ) transition is desired.

Interface

Perform a solvation calculation at a solvent/vacuum interface.

INPUT SECTION: $peqs
TYPE:

STRING

DEFAULT:

False

OPTIONS:

True	Modify the dielectric function to simulate an interface between solvent and vacuum.
False	Perform a solvation calculation in bulk solvent.

RECOMMENDATION:

The user will also need to specify the length scale on which the dielectric is smoothly interpolated from the bulk solvent value to vacuum, and the location of the Gibbs dividing surface (see below).

By setting Interface = True, the dielectric function on the Cartesian grid is further modified to mimic a solvent/vacuum interface. First, the solute cavity is constructed as specified by the SoluteCavity keyword as discussed above. Then, the dielectric is smoothly interpolated in the direction across the Gibbs dividing surface () using the following hyperbolic tangent switching function:

(11.42)

Eq. eq:hyperbolic2 will interpolate the dielectric function over the interface length, , centered at the Gibbs dividing surface, . Both of these values are controlled by the keywords InterfaceLength and GibbsDS, respectively. Similar to used in the spherical cavity construction, the parameter controls the sharpness of the interpolation across the Gibbs dividing surface and is approximated as .

InterfaceLength

Sets the length scale over which the dielectric function is smoothly transitioned from bulk solvent to vacuum in the -direction

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

None

OPTIONS:

Desired interface length in ngstroms.

RECOMMENDATION:

See the Supporting Information of Ref. Coons:2016 for a full description of the solvent/vacuum interface construction.

GibbsDS

Sets the location of the Gibbs dividing surface in the direction for an interface.

INPUT SECTION: $peqs
TYPE:

FLOAT

DEFAULT:

None

OPTIONS:

Desired location of the Gibbs dividing surface in ngstroms.

RECOMMENDATION:

Consult the literature. One such way to determine this value is to compute a density profile of the solvent as a function of and set this location to be where the solvent density has decreased to 50% of the bulk value. Usually, , where is the interface length described above.

NonequilJob

Obtain the non-equilibrium free energy of solvation for a vertical ionization process.

INPUT SECTION: $peqs
TYPE:

STRING

DEFAULT:

False

OPTIONS:

True	Compute the non-equilibrium free energy for a vertical ionization process.
False	Compute the equilibrium solvation free energy.

RECOMMENDATION:

None.

NonequilPartition

Specifies the scheme used to partition the fast and slow charges of the solvent.

INPUT SECTION: $peqs
TYPE:

STRING

DEFAULT:

Marcus

OPTIONS:

Marcus	Employ the Marcus partitioning scheme.
Pekar	Employ the Pekar partitioning scheme.

RECOMMENDATION:

Use the default. The non-equilibrium free energies are partition independent[774] and the Pekar scheme is computationally more expensive than the Marcus scheme.

NonequilState

Specifies the state of interest for a non-equilibrium vertical ionization.

INPUT SECTION: $peqs
TYPE:

STRING

DEFAULT:

Reference

OPTIONS:

Reference	The reference state from which to vertically remove an electron.
Ionized	The vertically ionized state.

RECOMMENDATION:

None.

11.2.10.5 Examples

The following example illustrates a free energy solvation calculation of a water molecule immersed in water. The Cartesian grid is cubic with a side length of 20.0 ngstroms and 153 grid points along each axis. The solute cavity is chosen to be RigidVDW centered on the scaled van der Waals radii. The dielectric is not set here, but takes the default value of 78.39.

Example 11.283  Free energy of solvation of water in water.

$rem
   SCF_CONVERGENCE   5
   THRESH            14
   EXCHANGE          omegab97x-v
   BASIS             6-311G*
   SOLVE_PEQ         true
   PEQ_SWITCH        3
   MEM_TOTAL         75000
   MEM_STATIC        8000
$end

$molecule
   0 1
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$peqs
   SoluteCavity  RigidVDW
   VDWType       Scaled
   VDWScale      1.2
   BatchSize     1193859
$end

$peqs_grid
   DimX 153 -10.00 10.00
   DimY 153 -10.00 10.00
   DimZ 153 -10.00 10.00
$end

The next example illustrates a free energy solvation calculation of a water molecule at a water/vacuum interface with the Gibbs dividing surface placed at 0.85 ngstroms, and the length of the interface is set to be 1.5 ngstroms. The Cartesian grid, solute cavity, and solvent dielectric are the same as in the previous example.

Example 11.284  Free energy of solvation of water at a water/vacuum interface.

$rem
   SCF_CONVERGENCE   5
   THRESH            14
   EXCHANGE          omegab97x-v
   BASIS             6-311G*
   SOLVE_PEQ         true
   PEQ_SWITCH        3
   MEM_TOTAL         75000
   MEM_STATIC        8000
$end

$molecule
   0 1
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$peqs
   SoluteCavity      RigidVDW
   VDWType           Scaled
   VDWScale          1.2
   Interface         true
   InterfaceLength   1.5
   GibbsDS           0.85
   BatchSize         1193859
$end

$peqs_grid
   DimX 153 -10.00 10.00
   DimY 153 -10.00 10.00
   DimZ 153 -10.00 10.00
$end

The final example illustrates a non-equilibrium (NONEQUILJOB = TRUE) solvation calculation for the vertical ionization of HO in bulk water. This is a compound job that will first calculate the equilibrium solvation free energy of the anionic state (NONEQUILSTATE = REFERENCE), and then compute the non-equilibrium correction for the ionized state (NONEQUILSTATE = IONIZED). The two jobs are separated by ”@@@". The Cartesian grid and solvent dielectric are the same as the previous examples, but the solute cavity is chosen to be Spherical with a radius of 1.884503499 . The Marcus scheme (NONEQUILPARTITION = MARCUS) is employed to partition the fast and slow charges of the solvent.

Example 11.285  Non-equilibrium free energy of solvation for the vertical ionization of HO.

$rem
   SCF_CONVERGENCE   5
   THRESH            14
   EXCHANGE          omegab97x-v
   BASIS             6-311G*
   SOLVE_PEQ         true
   PEQ_SWITCH        3
   MEM_TOTAL         75000
   MEM_STATIC        8000
$end

$molecule
   -1 2
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$peqs
   SoluteCavity        Spherical
   SphereRadius        1.884503499
   InterpolLength      1.0
   NonequilJob         true
   NonequilState       Reference
   NonequilPartition   Marcus
   BatchSize           1193859
$end

$peqs_grid
   DimX 153 -10.00 10.00
   DimY 153 -10.00 10.00
   DimZ 153 -10.00 10.00
$end

@@@

$rem
   scf_convergence   5
   thresh            14
   exchange          omegab97x-v
   basis             6-311G*
   solve_peq         true
   peq_switch        3
   mem_total         75000
   mem_static        8000
$end

$molecule
   0 1
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$peqs
   SoluteCavity        Spherical
   SphereRadius        1.884503499
   InterpolLength      1.0
   NonequilJob         true
   NonequilState       Ionized
   NonequilPartition   Marcus
   BatchSize           1193859
$end

$peqs_grid
   DimX 153 -10.00 10.00
   DimY 153 -10.00 10.00
   DimZ 153 -10.00 10.00
$end