Solvation Models

As we have noted several times, both the molecular mechanics calculations performed by PCModel and the MO calculations of SPARTAN, GAMESS, and GAUSSIAN treat the molecules as alone in a vacuum: no neighboring molecules of the same kind, and no solvent.

PCModel does allow us to change the value of the dielectric constant that is used in calculations of dipolar repulsions.

• This will change the relative energies of isomeric structures that have distinctly different dipole moments, but has little effect on the calculated structure itself.

• For example, with D = 1.5, the energy difference between anti- and gauche conformers of 1,2-dichloroethane is calculated to be 1.7 kcal/mol. With D = 80, the difference is 0.8 kcal/mol.

• Calculated structural parameters change not at all: the C-Cl bond length, for example, goes from 1.787 A to 1.788 A, and the Cl-C-C angle changes by 0.3 degrees.

We would expect that when a real molecule passes from the gas phase into solution

• The structure will relax to permit greater charge separation (electronic polarization)

• The solute also will polarize the solvent (electric polarization).

• When the favorable intermolecular consequences of further polarization are overcome by the intramolecular cost of that distortion, relaxation will end.

The energy terms associated with these polarizations are generally referred to as ENP (Electronic, Nuclear, and Polarization), and are primarily a function of the bulk dielectric.

• Other effects exist that are more specifically associated with the first solvation shell
  • The cavitation energy (making a "hole" in the solvent for the solute)

  • Attractive dispersion forces between the solute and solvent molecules

  • Local structural changes in the solvent, such as changes in the extent of hydrogen bonding.

• These effects are referred to by the acronym CDS (for Cavitation, Dispersion, and solvent Structural).

The current methods for solvation modeling generally follow one of two approaches.

• The first involves the explicit construction of a solvent shell, consisting of anywhere from several dozen to several hundred solvent molecules, about the solute.¹
  • This supermolecular system serves as a basis for simulations from which thermodynamic data related to solvation can be extracted.

  • The simulations typically involve either a Monte-Carlo pseudo-random sampling² or solution of Newton's equations of motion to generate molecular dynamics trajectories.³

  • In either case, the size of the system usually prohibits analysis at the quantum mechanical level, and a classical mechanical force field is employed instead.

  • Polarizable solvent models are not general, and polarization of the solute also is not usually considered; thus modeling of the ENP terms is incomplete.

  • By contrast, the explicit representation of the solvent leads automatically to inclusion of the CDS terms in the force field.

• The main alternative to these simulation procedures replaces the explicit solvent molecules with a continuum having the appropriate bulk dielectric constant.
  • Within this much simpler system, quantum mechanical approaches may be employed for the ENP effects.

  • The simplest approach, known as the Kirkwood-Onsager model⁴, uses a Taylor series to represent the classical multipolar expansion of the electronic structure.
    • The expansion is truncated at the dipolar term, thus including only the interactions of atomic charges and dipoles with the medium.

    • This simplification, along with the assumption that the molecule resides in a spherical cavity in the solvent, leads to some inaccuracies, but is simple to implement.

• More recently, continuum approaches have been developed that use a generalized Born⁵ formalism for the interaction of atomic partial charges with a surrounding dielectric.

• In principle, the atom centered monopoles generate all of the multipoles required to represent the electronic distribution.

Neither of these continuum approaches as described attempts to take account of the energetic effects in the CDS terms, thus focusing on a different set of contributions to solvation than the molelcular dynamics methods.

Cramer, Truhlar, and their coworkers, have developed the SMx set of solvation models in an attempt to address these limitations.⁶ Their models are incorporated into SPARTAN, along with refinements developed by Hehre and coworkers at Wavefunction. The C-T models were specifically developed to be associated with the AM1 and PM3 semi-empirical Hamiltonians. (However, this group has very recently published an extension of their methods to some ab initio basis sets.)

These methods use an extended Born treatment of the ENP solvation terms.

• They include atomic parameters to account for the CDS effects by assigning unique surface tensions to the solvent accessible surfaces of various functional groups within the solute molecule.

• The relationship between solvent accessible surface and specific solvation phenomena, including hydrogen bonding, is well established.⁷

• The particular virtue of including atom-specific terms is that it allows fitting to experimental free energies of solution. This in turn permits calculation of absolute free energies of solution (for a standard state taken as 1 M).

As implemented in SPARTAN, however, these methods result in calculation of a solvation enthalpy as the difference in enthalpy of formation between the molecule in the gas phase and the molecule in solution.

• This saves computation time, and represents no great loss of information, since few experimental solvation free energies are available for organic molecules, anyhow.

• Inclusion of solvation does result in a substantial percentage increase in computation time, because of the additional integrals that must be evaluated
  • In general, the increase is proportional to the square of the number of heavy atoms, if full geometry optimization is done.

  • That is, for six heavy atoms, an optimization in water will take 36 times as long as one in the gas phase.

• For this reason, one may wish to do single point calculations in solvent studies, using the optimized gas phase molecular structure.

• This, of course, eliminates consideration of the relaxation effects decribed above, and introduces some error.

These molecules are small enough that you can do full optimizations in both gas phase and water (use either AM1 or PM3 with the C-T model; it gives more negative, and therefore better, energies.) Interestingly, the energy of the pyridinol is strongly dependent upon the orientation of the O-H bond; syn- to N is lower than anti-. (Why?)

In principle, solvation energies calculated by the Cramer-Truhlar method are transferable from semi-empirical to ab initio structures. In this case, one should optimize the geometry with the appropriate ab initio basis set, and then do single point AM1 calculations using the ab initio geometry to evaluate the solvation energy. The result can then be applied to the ab initio energy. (Remember the difference in units!)

GAUSSIAN implements a variety of SCRF (Self-Consistent Reaction Field) solvation models.

• These all model the solvent as a continuum of dielectric constant e

• The solute is placed in a cavity within the solvent. The available methods differ in how they define the cavity and the reaction field
  • The simplest is the Onsager model⁹, obtained by using the keyword SCRF=dipole, which is simply a dipole-induced dipole model. The radius of the cavity and the dielectric constant of the solvent must be supplied as input.

  • The Polarized Continuum Model¹⁰ (SCRF=PCM) defines the cavity as the union of a series of interlocking spheres centered on the atoms, and uses a numerical representation of the polarization of the solvent; integration is numerical rather than analytical.

  • The isodensity PCM model¹¹ (SCRF=IPCM) uses an isodensity surface of the molecule as the cavity, using a converged SCF, and iterating until the cavity shape no longer changes.

  • The Self-Consistent IPCM model¹² (SCRF=SCIPM) includes the effect of solvation of the solution of the SCF problem

• GAMESS includes the Onsager model, by use of the $SCRF group, and the PCM model with the $PCM group. Consult the manual for proper input format if you wish to use these methods.

Clearly, solvation modeling is an extremely complex area, and requires considerable sophistication for proper usage.