11.2.8 SM8, SM12, and SMD Models

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,⁹¹ 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.⁶³⁵

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,⁹¹ 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.⁴⁶⁸ 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 SM$x$ calculations should be listed in a $smx input section; for SM8, the only additional variable that is required is the name of the solvent:

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:

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 ${\mathrm{CH}}_3{\mathrm{OH}}_2^{+}$ in water at the B3LYP/6-31G* level;

• •

  finite-difference frequency (with analytical gradient) calculation for protonated methanol ${\mathrm{CH}}_3{\mathrm{OH}}_2^{+}$ 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 energies⁵⁷⁸ is not yet available in Q-Chem. (Nonequilibrium versions of PCMs are available instead; see Section 11.2.2.3.)

The SM12 model⁶³³ 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,⁶³⁴ 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.

The SMD model⁶³² 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.