12.2 Chemical Solvent Models

12.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 Chipman742, 743, 744, 745 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 12.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,961 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

ΔGCMIRS(ρs)=ΔGSS(V)PE(ρs)+ΔGDEFESR(ρs), (12.10)

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

ΔGDEFESR=ΔGdisp+ΔGexch+ΔGFESR. (12.11)

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 Zhi-Qiang You and John Herbert.1036 In the course of this work, a serious error was discovered in the the original implementation of ΔGdisp 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. (12.11) changes significantly.1036 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 ΔG compare very favorably to those of the SMx models that are described in Section 12.2.8, e.g., mean unsigned errors <0.7 kcal/mol in benzene and cyclohexane and <1.5 kcal/mol in water. The latter statistic includes challenging ionic solutes; errors for charge-neutral aqueous solutes are smaller still.1036

Solvent SolvRho A B C D Gamma
Benzene 0.0421 -0.00522 0.01294
Cyclohexane 0.0396 -0.00938 0.03184
DMSO 0.05279 -0.00951 0.044791 -162.07 4.1
CH3CN 0.03764 -0.008178 0.045278 -0.33914 1.3
Water 0.05 -0.006736 0.032698 -1249.6 -21.405 3.7
Table 12.4: Optimized CMIRS parameters (from Ref. 1036) using RHOISO = 0.0010 a.u. and, for ions, the proton reference ΔG value from Ref. 463.
Solvent SolvRho A B C D Gamma
Benzene 0.0421 -0.00572 0.01116
Cyclohexane 0.0396 -0.00721 0.05618
DMSO 0.05279 -0.002523 0.011757 -817.93 4.3
CH3CN 0.03764 -0.003805 0.03223 -0.44492 1.2
Water 0.05 -0.006496 0.050833 -566.7 -30.503 3.2
Table 12.5: Optimized CMIRS parameters (from Ref. 1036) using RHOISO = 0.0005 a.u. and, for ions, the proton reference ΔG value from Ref. 463.

The current implementation of CMIRS in Q-Chem computes the electrostatic energy using Chipman’s isodensity SS(V)PE module. The resulting isodensity 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 = 1 in an isodensity SS(V)PE calculation (see Section 12.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 δ=7 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 40 grid points for efficient integration with accuracy. Optimized parameters for the supported solvents are listed in Tables 12.4 and 12.5 for two different values of ρs.1036

Example 12.11  CMIRS calculation for methane in the water solvent.

$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
   MEM_TOTAL       = 4000
   MEM_STATIC      = 400
   SYM_IGNORE      = true
   XC_GRID         = 000096000974
$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         = 000096000974
$end

$molecule
   read
$end

$svp
   RHOISO=0.001, DIELST=78.36, NPTLEB=974,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