CCL:G: Gaussian Delta G of solution calculation?



This question, in various forms, does just keep coming up...

While I'm going to take the opportunity to refer to a detailed, formal, and I hope dispositive paper just published by the Minnesota group (doi 10.1021/jp205508z), I'll also make some effort here to be more cookbooky in an explanation here. Pedagogy in front! Not just for Gaussian, but more general.

The free energy of solvation is, clearly, defined as the DIFFERENCE in the free energy of a species in the gas phase and in solution. Free energy is an ensemble property, not a molecular property, so we are immediately faced with the need to make some approximations in order to render the modeling tractable.

In the gas phase, those approximations are by now fairly standard and nearly universal. We make the ideal-gas approximation (so that the partition function of a mole of molecules is simply the product of Avogadro's number of molecular partition functions), and we make the rigid rotator and harmonic oscillator approximations to simplify the rovibrational partition functions, and, voila, we have a means to compute a standard state free energy (Gibbs free energy -- free enthalpy for my colleagues in (far more logical) German-speaking countries). Pretty much every electronic structure program on earth will do this for you when you run a vibrational frequency calculation in the gas phase. Poof -- you get E, H, S (third law), and G (and Cv too, if you care).

Note that you should be careful to recognize that (i) unless you overrode it, the program chose some default standard state concentration (1 bar? 1 atm? you should know...) and temperature (nearly always 298 K) and scale factor for vibrational frequencies (typically 1, but that might not be the best value for a given level of theory) and (ii) the harmonic oscillator approximation is catastrophically bad for super-low frequency vibrations (below, say, 50 wavenumbers, to pick an arbitrary value). There are fixes for the latter problem, but I'll let someone else post about that.

Now, what about G of a solute in solution? Well, to begin, since we're not dealing with a pure component anymore (at least, not if the solvent is different than the solute), we need to assume that we have an ideal solution, and we should recognize that we're talking about a partial molar quantity (more often referred to in the rigorous literature as a chemical potential than a free energy, but that's a matter of tradition). In any case, we need to assemble a free energy as a sum of

electronic energy (i.e., potential and kinetic energy of electrons for a fixed set of nuclear positions)

coupling to medium (which includes electrostatic and non-electrostatic components, although, it being chemistry, EVERYTHING is electrostatic... -- in practice, however, with continuum models, electrostatic means what you get by assuming the molecule is a charge distribution in a cavity embedded in a classical dielectric medium in which case one can apply the Poisson equation -- non-electrostatic is everything else -- dispersion, cavitation, covalent components of hydrogen bonding, hydrophobic effects, you name it)

temperature dependent translational (?), rotational (?), vibrational, and conformational contributions -- the question marks indicate conceptual issues

So, a few points to bear in mind. The optimal geometry in solution is unlikely to be the same as that in the gas phase -- but it might be close. You just have to decide for yourself if you want to reoptimize or not.

Same for the vibrational and conformational contributions to free energy in solution -- they might be very, very close to those in the gas phase -- or they might not. If you assume that they ARE the same, you avoid having to do a frequency calculation in solution (and you avoid wondering what it means to do vibrational frequencies in a continuum, which in principle means a surrounding that is in equilibrium with the solute -- but how can a medium composed of molecules be fully in equilibrium with a molecular solute on the timescale of the solute's vibrations, since the solvent vibrations are on the same timescale?)

Note that if I assume no change in the various parts of the free energy EXCEPT for the electronic energy and solute-solvent coupling (more on that momentarily), life is pretty easy. The free energy of solvation is the difference in the self-consistent reaction-field (SCRF) energy INCLUDING non-electrostatic effects, and the electronic energy in the gas phase. That is, I look at the expectation value of <H+(1/2)V> for the solvated wave function, where V is the reaction field operator, add non-electrostatic effects (typically NOT dependent on the quantum wave function, so added post facto, although there are a few exceptions in the literature), subtract the expectation value of <H> for the gas-phase wave function (note that you might have done the two expectation values at different geometries, or you might have used the gas-phase geometry for both -- your choice -- the former is more "physical", certainly, but the latter is a useful approximation in many instances), and you are done. You've got the free energy of solvation FOR IDENTICAL STANDARD-STATE CONCENTRATIONS. That is, the number you have in hand assumes no change in standard-state concentration. However, many experimental solvation free energies are tabulated for, say, 1 atm gaseous standard states and 1 M solution standard states. To compare the computed value to the tabulated value, one needs to correct for the standard-state concentration difference.

In the interest of the cookbook, let me be more practical. Thus, let's say that I compute a gas-phase G value, including all contributions, electronic and otherwise of

-3.000 00 a.u.

and, let's say that the electronic energy alone in the gas phase is

-3.020 00 a.u.  (so ZPVE and thermal contributions to G are +0.020 00 a.u.)

and, finally, let's say that my SCRF calculation provides an electronic energy INCLUDING non-electrostatic effects of

-3.030 00 a.u.

In that case, my free energy of solvation is -0.010 00 a.u. (difference of -3.030 00 and -3.020 00 a.u.) And, if I want to think about my free energy in solution, I can make the assumption that there is no change in the ZPVE and thermal contributions, in which case I would have G in solution equals -3.010 00 a.u. (which is gas-phase G of -3.000 00 a.u. plus free energy of solvation -0.010 00 a.u.)

But, just to be clear, if my gas phase G referred to a 1 atm standard state concentration, for an ideal gas at 298 K and 1 atm, that implies a molarity of 1/24.5 M. If I want my G in solution to be for a 1 M standard state, I need to pay the entropy penalty to compress my concentration from 1/24.5 M to 1 M, which is about 1.9 kcal/mol (the proof is left to the reader...) So, my 1 M free energy in solution is not -3.010 00 but rather about -3.006 99 a.u.

The above is an example of how almost all free energies of solvation and free energies in solution are computed in the literature using continuum solvation models (at least if they're done properly!)

Lots of important details glossed over a bit above (in the interests of clarity, I demur). But, to be more thorough, let's note:

1)  Why was it (1/2)V in the SCRF calculation? -- the 1/2 comes from linear response theory and assumes that you spend precisely half of the favorable coupling energy organizing the medium so that it provides a favorable reaction field.

2)  How can there be a translational partition function for a solute in solution? There isn't one -- but there is something called a liberational free energy associated with accessible volume, and Ben-Naim showed some time ago that the value is identical to that for the a particle-in-a-box having the same standard-state concentration -- i.e., there is no change on going from gas-phase translational partition function to liberational partition function for the same standard-state concentration for an ideal solution. When there are issues with non-accessible volume, however, account must be taked (cf. Flory-Huggins theory).

3)  How can there be a rotational partition function for a solute in solution? There isn't one -- solute rotations become librations that are almost certainly intimately coupled with first-solvation-shell motions. In essence, assuming no change in "rotational partition function" implies assuming no free energetic consequence associated with moving from rotations to librations. This remains a poorly resolved question, but, in practice, since most continuum solvation models are semiempirical in nature (having been parameterized against experimental data) any actual changes in free energy have been absorbed in the parameterization as best as possible. If you find that unsatisfying, hey, feel free not to use continuum solvent models -- it's certainly ok by me...

4)  Where did those non-electrostatic effects come from? Every model is different in that regard, and I won't attempt to summarize a review's-worth of material in an email. Lots of nice Chem. Rev. articles over the years on continuum models if you want to catch up.

Finally, what is described above is a popular approach for computing solvation free energies and free energies in solution, but by no means the only approach out there. A non-exhaustive list to compute either or both solvation free energies or free energies in solution includes free-energy perturbation from explicit simulations, RISM-based models, fragment-based models derived > from a statistical mechanical approach (including COSMO-RS and variations on that theme), fragment-based models from expert learning, and models relying on alternative physicochemical approaches to computing interaction energies (e.g., SPARC). These alternative models can be quantal, classical, or SMILESal (which is to say, more in the realm of chemoinformatics than physical chemistry). Let a thousand flowers bloom.

I hope that this post serves as a useful archival reference for CCL users present and future. Best wishes to all for a peaceful winter solstice (or summer, for my antipodeal colleagues).

Chris



On Nov 30, 2011, at 1:46 PM, Close, David M. CLOSED#,#mail.etsu.edu wrote:

  Does anyone know how Gaussian calculates deltaG of solvation?  This was in G98 and was automatic.  In G03 one had to add SCFVAC in the scrf(CPCM,read) read list to get the delta G result.  I have several questions about the methods used.  I presume that during the SCRF calculation the program has to have a separate SCF calculation on the input molecule  in the gas phase, along with a frequency calculation to know the free energy of the gas phase structure.  Is this correct?  If so, are the ZPE and free energy correction terms multiplied by a scale factor (0.92 for example?).  Also to do solvation energy calculations one need to convert the gas phase reference state from 1 atmos. to 1 M.  Does the program report this corrected delta G value? 
  My reason for asking is that I have done these actual calculations separately with energy/frequency calculations on the neutral and anion, and I get slightly difference answers from those answers using SCFVAC, and I need to know why this is?
  Regards, Dave Close.

--


Christopher J. Cramer

Elmore H. Northey Professor

University of Minnesota

Department of Chemistry

207 Pleasant St. SE

Minneapolis, MN 55455-0431

--------------------------

Phone:  (612) 624-0859 || FAX:  (612) 626-7541

Mobile: (952) 297-2575

email:  cramer|a|umn.edu

jabber:  cramer|a|jabber.umn.edu

http://pollux.chem.umn.edu

(website includes information about the textbook "Essentials

    of Computational Chemistry:  Theories and Models, 2nd Edition")