11.2 Chemical Solvent Models

11.2.8 SM8, SM12, and SMD Models

The SMx models developed by Cramer, Truhlar, and coworkers at the University of Minnesota are a class of implicit solvation models that are designed to be “universal” in the sense that they can be applied to any solvent for which a small of descriptors is known.188 (Note that the x in SMx 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. 635. 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,91 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.635

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. 587 and Bondi’s values for the van der Waals radii,91 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 (i) assuming fixed and model-dependent values of the intrinsic Coulomb radii, and (ii) 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.468 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. 468. Other basis sets should not be used in SM8 calculations.

The SM8 solvation free energy is output at T=298 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 P=1 atm for the gase. To convert 1 M-to1 M solvation free energies at 298 K to a standard state consisting of P=1 atm for the gas and a 1 M concentration in solution, add +1.89 kcal/mol to the computed solvation free energy.

Solution-phase geometry optimizations can be carried out, but basis sets that use spherical harmonic d functions, or angular momentum higher than d (f, g, etc.) are not supported. Since, by definition, the 6-31G*, 6-31+G*, and 6-31+G** basis sets have Cartesian d 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. 1118.

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 SMx 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 SMx solvent
INPUT SECTION: $smx
TYPE:
       STRING
DEFAULT:
       water
OPTIONS:
       Any name from the list of solvents given below.
RECOMMENDATION:
       NONE

1,1,1-trichloroethane bromoethane m-ethylbenzoate
1,1,2-trichloroethane bromooctane m-ethylethanoate
1,1-dichloroethane butanal m-ethylmethanoate
1,2,4-trimethylbenzene butanoicacid m-ethylphenylketone
1,4-dioxane butanone m-ethylpropanoate
1-bromo-2-methylpropane butanonitrile m-ethylbutanoate
1-bromopentane butylethanoate m-ethylcyclohexane
1-bromopropane butylamine m-ethylformamide
1-butanol butylbenzene m-xylene
1-chloropentane carbon disulfide heptane
1-chloropropane carbon tetrachloride hexadecane
1-decanol chlorobenzene hexane
1-fluorooctane chlorotoluene nitrobenzene
1-heptanol cis-1,2-dimethylcyclohexane nitroethane
1-hexanol decalin nitromethane
1-hexene cyclohexane methylaniline
1-hexyne cyclohexanone nonane
1-iodobutane cyclopentane octane
1-iodopentene cyclopentanol pentane
1-iodopropane cyclopentanone o-chlorotoluene
1-nitropropane decane o-cresol
1-nonanol dibromomethane o-dichlorobenzene
1-octanol dibutyl ether o-nitrotoluene
1-pentanol dichloromethane o-xylene
1-pentene diethyl ether pentadecane
1-pentyne diethylsulfide pentanal
1-propanol diethylamine pentanoic acid
2,2,20trifluoroethanol diiodomethane pentylethanoate
2,2,4-trimethylpentane dimethyldisulfide pentylamine
2,4-dimethylpentane dimethylacetamide perfluorobenzene
2,4-dimethylpyridine dimethylformamide phenyl ether
2,6-dimethylpyridine dimethylpyridine propanal
2-bromopropane dimethyl sulfoxide propanoic acid
2-chlorobutane dipropylamine propanonitrile
2-heptanone dodecane propylethanoate
2-hexanone E-1,2-dichloroethene propylamine
2-methylpentane E-2-pentene p-xylene
2-methylpyridine ethanethiol pyridine
2-nitropropane ethanol pyrrolidine
2-octanone ethylethanoate sec-butanol
2-pentanone ethylmethanoate t-butanol
2-propanol ethylphenyl ether t-butylbenzene
2-propen-1-ol ethylbenzene tetrachloroethene
3-methylpyridine ethylene glycol tetrahydrofuran
3-pentanone fluorobenzene tetrahyrothiophenedioxide
4-heptanone formamide tetralin
4-methyl-2-pentanone formic acid thiophene
4-methylpyridine hexadecyliodide thiophenol
5-nonanone hexanoic acid toluene
acetic acid iodobenzene trans-decalin
acetone iodoethane tribromomethane
acetonitrile iodomethane tributylphosphate
aniline isobutanol trichloroethene
anisole isopropyl ether trichloromethane
benzaldehyde isopropylbenzene triethylamine
benzene isopropyltoluene undecane
benzonitrile m-cresol water
benzyl alcohol mesitylene Z-1,2-dichloroethene
bromobenzene methanol other

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:

Dielec dielectric constant, ε, of the solvent
SolN index of refraction at optical frequencies at 293 K, n20D
SolA Abraham’s hydrogen bond acidity, α2H
SolB Abraham’s hydrogen bond basicity, β2H
SolG γ=γm/γ0 (default is 0.0), where γm is the macroscopic surface tension at air/solvent
interface at 298 K, and γ0 is 1 cal mol-1 Å-2 (1 dyne/cm = 1.43932 cal mol-1 Å-2)
SolC aromaticity, ϕ : the fraction of non-hydrogenic solvent atoms that are aromatic
carbon atoms
SolH electronegative “halogenicity”, ψ : the fraction of non-hydrogenic solvent atoms that are
F, Cl or Br

For a desired solvent, these values can be derived from experiment or from interpolation or extrapolation of data available for other solvents. Solvent parameters for common organic solvents are tabulated in the Minnesota Solvent Descriptor Database. The latest version of this database is available at:

http://comp.chem.umn.edu/solvation/mnsddb.pdf

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 CH3OH2+ in water at the B3LYP/6-31G* level;

  • finite-difference frequency (with analytical gradient) calculation for protonated methanol CH3OH2+ 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. 960. 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. 187.

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 nonequilibrium algorithm required for vertical excitation energies578 is not yet available in Q-Chem. (Nonequilibrium versions of PCMs are available instead; see Section 11.2.2.3.)

11.2.8.2 The SM12 Model

The SM12 model633 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,634 which are based on Hirshfeld population analysis, or else charges derived from the electrostatic potential,881, 103 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 charges634 MK Merz-Singh-Kollman charges881 CHELPG ChElPG charges103
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.12  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 model632 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.

In version 5.2 and after, the default surface discretization method is changed from VTN to SWIG in order to ensure the smoothness of potential energy surface. Also, the gas phase SCF calculation that takes place before the SMD/SM8/SM12 calculation is turned off by default. If one wants to obtain the solvation free energy, then the gas phase calculation is still required and it can be turned on by setting SMX_GAS_PHASE = TRUE. Setting this $rem variable to TRUE might also be helpful if directly converging SCF with the SMx models is difficult.

SMX_GAS_PHASE
       Converge the gas-phase SCF first before doing calculations with SMx models
TYPE:
       LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       FALSE Run SMx calculations directly TRUE Run gas-phase calculation first
RECOMMENDATION:
       Use the default unless solvation free energy is needed. Set it to TRUE if the SCF calculation fails to converge otherwise.

Example 11.13  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