as Chris Cramer has shifted the thread from a special Gaussian to a
general continuum solvation context, and since he suggests his post
"as a useful archival reference for CCL users present and future" I
need to make a few remarks on Chris' comprehensive entry:
- Chris is starting his entry with a reference to a paper by Ribero,
Marenich, Cramer and Truhlar (RMCT), which he considers as a
dispositive paper. Unfortunately that paper already in the abstract
explicitly contradicts to our recent paper (Ho, Klamt, Coote,
Comment on the correct use of continuum solvent models,
, 13442). This clearly
shows that apparently some aspects are under controversial
discussion between the experts, while Chris is describing things as
if their view would be generally accepted. Without any question, I
could agree with 70% of his statements, but I disagree with at least
20%.
- It should be especially noted that in the theoretical derivation
of the _expression_ for the effective vibrational potential of the
solute given in the RMCT paper there is an important either
mathematical or conceptional error: They are starting from a
correct, classical formulation of the partition function, incl.
integrals over the solute coordinates R and their momenta P, and
over the external (solvent) coordinates r and momenta p. Then they
separate the integral over R and define the remainder to be an
entity W(R), and in the next sentence they identify W(R) with the
potential of mean force (PMF(R)). But the remainder of the partition
function includes the integration over the momenta P of the solute,
while the PMF(R) at a given set of internal coordinates R is just
the integral over all states of the outer system at these fixed
coordinates R, i.e. PMF(R) does not include an integration over the
kinetic states P of the solute. Hence W(R) and PMF(R) cannot be the
same. By assuming that they are the same RMCT are tacitly omitting
the P integral, but it is just this integral which is the
problematic part in the partition function of a solute in a solvent.
Hence all the further statement in the RMCT paper are based on a
wrong derivation.
- Also the statements of Chris regarding the liberation free energy
are only partly right: Yes, the volume integral gives the same
(almost trivial) RTln(V/V0) result as in the gas phase, but Ben Naim
has no _expression_ for the integral over the kinetic states. Hence
nobody, neither Ben Naim, nor Chris Cramer, nor me has a formally
correct theory for the complete partition function of a solute in a
solvent with respect to the translational, stiff-rotational, and
vibration/rotation contributions. The solvation free energy
contribution resulting from the differences of partition function
compared to the ideal gas partition function are all effectively
parameterized into the non-electrostatic contributions of our
solvation models, e.g. in to the cavitation free energy, etc.
Lukily, these contributions seem to be sufficiently systematic (most
likely to first order surface proportional) in order to be subsumed
in the functional forms of the non-electrostatic contributions.
Otherwise all of our continuum solvation models would not work.
- But as a consequence of this - if we do not know how to formulate
the partition function in solution, and since the differences are
already empirically included in the solvation free energy of the
continuum solvation mode - it is worthless, yes even
counterproductive, to try to calculate the vibrational part of this
partition function from the frequencies in the continuum solvent, as
RMCT are suggesting in their paper.
I am absolute convinced, and I have stated this several times before
in the CCL, and we worked it out in more detail in the Ho, Klamt,
Coote paper) that the most consistent way to calculate the free
energy in solution is
- to first calculate the total free energy in the gas phase as
accurately as we can
- and then add the effective free energy of solvation of the solute
as calculated with a good continuum solvation model (electrostatic +
non-electrostatic, or the dG_solv in COSMO-RS)
Stay out of the calculation of the frequencies in solution for th
purpose of a free energy calculation, unless you first write down a
consistent theory for the partition function in solution, which
would be highly appreciated, but seems to be out of reach.
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
(website
includes information about the
textbook "Essentials
of
Computational Chemistry: Theories
and Models, 2nd Edition")