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
multipolar expansion method, also known as a
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.
(This should not be confused with the Onsager model, in which a dipole approximation is used,
and to avoid confusion the multipolar expansion method is perhaps preferable.)
More sophisticated models, which use a molecule-shaped cavity and
the full molecular electrostatic potential, include the
conductor-like PCM (C-PCM),
J. Phys. Chem. A
(1998), 102, pp. 1995. , 226 J. Comput. Chem.
(2003), 24, pp. 669. the conductor-like screening model (COSMO), the “surface and simulation of volume polarization for electrostatics” [SS(V)PE] model, 199 J. Chem. Phys.
(2000), 112, pp. 5558. and the closely-related “integral equation formalism” (IEF-PCM). 152 J. Chem. Phys.
(1997), 107, pp. 3032. , 153 J. Chem. Phys.
(2001), 114, pp. 4744. For an overview of all of these methods and their interconnections, see the review by Herbert.
The C-PCM, IEF-PCM, and SS(V)PE are examples of what are called “apparent
surface charge” SCRF models, and the term polarizable continuum
models (PCMs), as popularized by Tomasi and coworkers,
(2005), 106, pp. 2999. 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,606, 607, 609, 446, 604 which 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. 606 J. Phys. Chem. Lett.
(2010), 1, pp. 556. , 607 J. Chem. Phys.
(2010), 133, pp. 244111. 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”
J. Phys. Chem. B
(1997), 101, pp. 5583. , 317 J. Phys. Chem. B
(1999), 103, pp. 10282. as well as versions 8 and 12 of the SM models, and the SMD model, developed at the University of Minnesota. 733 J. Chem. Theory Comput.
(2007), 3, pp. 2011. , 731 J. Chem. Theory Comput.
(2013), 9, pp. 609. , 730 J. Phys. Chem. B
(2009), 113, pp. 6378. SM8 and SM12 are based upon the generalized Born method for electrostatics, augmented with atomic surface tensions intended to capture nonelectrostatic effects (cavitation, dispersion, exchange repulsion, and changes in solvent structure). Empirical corrections of this sort are also available for the PCMs mentioned above, although these have not always been particularly successful. In contrast, the nonelectrostatic interactions in SM8 and SM12 are carefully parameterized to reproduce experimental free energies of solvation. The SMD model (in which the “D” is for “density”) combines IEF-PCM with similar nonelectrostatic corrections. Statistically speaking, SMD is not any more or less accurate than other SM models, but has the advantage that its electrostatic interactions are density-based and this model can be used in arbitrary basis sets. The SM8 model uses Mulliken atomic charges and thus can only be used in certain small basis sets for which it was parameterized, while the SM12 model replaces these with Hirshfeld-based atomic charges that are more robust with respect to the choice of basis set.
|Langevin Dipoles||atomic spheres||dipoles in||no||all|
|C-PCM||atomic spheres||point charges or||user-||all|
|SS(V)PE/||atomic spheres||point charges or||user-||all|
|Isodensity SS(V)PE||isodensity contour||point charges||none||all|
|CMIRS||isodensity contour||point charges||automatic||all|
|Generalized Born electrostatic model; does not require cavity construction.|
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, which the caveat that only certain basis sets are supported for SM8. 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. This represents a “zeroth-order” inclusion of the solvent effects at the correlated level of theory, but is perfectly adequate for many applications as higher-order corrections are usually small. Table 11.2 also summarizes the availability of analytical energy gradients for implicit solvent models. (Finite-difference gradients and Hessians are requested automatically for calculations where the requisite analytic derivatives are not available.)
Note: The 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.
|Energy Derivatives||C-PCM||SS(V)PE &||CMIRS||COSMO||SM8||SM12||SMD|
|SCF energy gradient||yes||yes||no||yes||yes||no||yes|
|SCF energy Hessian||yes||no||no||yes||no||no||no|
|CIS/TDDFT energy gradient||yes||no||— unsupported —|
|CIS/TDDFT energy Hessian||yes||no||— unsupported —|
|MP2 & double-hybrid energy||— unsupported —|
|Coupled cluster methods||— unsupported —|
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, SS(V)PE, and IEF-PCM. The latter two models are formally
equivalent at the level of integral equations,
J. Chem. Phys.
(2001), 114, pp. 4744. but exhibit some some differences in their numerical implementation. 609 Chem. Phys. Lett.
(2011), 509, pp. 77. 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”), 554 J. Chem. Phys.
(1996), 105, pp. 9972. although later work cast doubt on the theoretical justification for this and other ad hoc charge renormalization procedures. 198 J. Chem. Phys.
(1999), 110, pp. 8012. [The IEF-PCM and SS(V)PE methods already contain an implicit correction for outlying charge, as does the C-PCM method that is derived as an approximation to these models, 198 J. Chem. Phys.
(1999), 110, pp. 8012. although this was not recognized at the time that COSMO was formulated. See Ref. Herbert:2021b for a historical discussion of these developments.] In any case, COSMO is described 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 in
conjunction with appropriate job-control variables in the $pcm input section)
uses a solute cavity constructed from atom-centered spheres, in keeping with 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.
Theor. Chem. Acc.
(2002), 107, pp. 90. , 192 J. Chem. Phys.
(2003), 119, pp. 10289. , 201 J. Chem. Phys.
(2006), 124, pp. 224111. This is an appealing, one-parameter cavity construction that avoid many of the problems with cusps in “van der Waals” cavity surfaces that are constructed from atom-centered spheres, and the isodensity implementation of SS(V)PE forms the basis of a physics-based continuum solvation model called CMIRS that is described in Section 11.2.6. 879 J. Chem. Theory Comput.
(2011), 7, pp. 3952. , 880 J. Phys. Chem. A
(2013), 117, pp. 5812. , 881 J. Chem. Theory Comput.
(2014), 10, pp. 211. , 882 J. Phys. Chem. A
(2015), 119, pp. 5173. , 1224 J. Chem. Theory Comput.
(2016), 12, pp. 4338. CMIRS is competitive in accuracy (for solvation free energies) with the best-available SM solvation models, despite using far fewer fitting parameters. However, analytic energy gradients are not available for the isodensity cavity construction and therefore not available for CMIRS, whereas these gradients are available when the cavity surface is constructed from atom-centered spheres.
Regarding the accuracy of these models for solvation free energies (),
the SM models generally achieve
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.
Acc. Chem. Res.
(2008), 41, pp. 760. (Note that the SM12 and SMD models generally do not improve too much upon SM8 in any statistical sense, 730 J. Phys. Chem. B
(2009), 113, pp. 6378. , 731 J. Chem. Theory Comput.
(2013), 9, pp. 609. but do extend these models to arbitrary basis sets whereas SM8 is limited to a few small basis sets.) To achieve accuracy comparable to the SM models within the PCM class of solvent models, nonelectrostatic terms must be added. Among the various PCMs described above, the only one that constitutes a “black box” model for solvation energies is CMIRS, which slightly outperforms the SM models in a statistical sense, 1224 J. Chem. Theory Comput.
(2016), 12, pp. 4338. although it is only available for a few solvents.
The following sections provide more details regarding theory and job control for the various implicit solvent models that are available in Q-Chem. Ref. Herbert:2021b contains both formal comparisons amongst these models as well as a side-by-side comparison of the accuracy of solvation free energies.