Some statistical mechanics

1
Some statistical mechanics

• How can we analyze the data generated from an ab
  initio simulation?
• Wed like to calculate thermodynamic, measurable
  quantities free energies, enthalpies, entropies,
  etc.
• From these we can predict, for example, binding
  FEs for drug design, FE pathways for ion transit
  through channels, etc
• FEs are difficult to calculate due to entropic
  contributions

2
Solvation Stat Mech Tutorial

• This lecture will rely heavily on our book The
  Potential Distribution Theorem and Models of
  Molecular Solutions, by Beck, Paulaitis, and
  Pratt.
• Well label the Potl Distn Thm --gt PDT
• The focus is on the chemical potential, which is
  the partial molar Gibbs FE

3
A Theory for Mu

• Originally developed by Widom in 1963
• Well use the Helmholtz FE and give the old
  fashioned derivation from the canonical ensemble.
  You can also derive this from the grand
  canonical ensemble, in our book.

(Well consider a one component system to keep
things simple, easy to generalize to multiple
components)
4
PDT contd

• TD/SM calc of mu (monatomic fluid here)

5
PDT contd

• We can rearrange
• Then

where
The average in the second term means the solute
is not coupled with the rest of the solution
during the averaging. Also we can pick any point
in the solution to do the averaging, or average
over all points.
6
PDT contd

• The first term is the ideal contribution to the
  chemical potential which is a function of the
  number density and the thermal de Broglie
  wavelength
• Notice that in this derivation we utilized
  Boltzmann statistics but there is still the n!
  there which is due to the indistinguishability of
  particles

7
PDT contd

• At equilibrium the chemical potential for a given
  component is constant everywhere in space. But
  the number density and the last term can vary so
  as to add up to that constant (in an
  inhomogeneous system)

8
PDT contd

• The last term is the excess chemical potential
  and has to do with the interactions of the solute
  with its surroundings.
• What happens if we have a molecular solute?

Exercise derive all of the above formulas step
by step
9
PDT contd

• Above, in the first term, is the internal
  partition function for the molecule in the gas
  phase (or vacuum really) at some temperature T
• That can be calculated by quantum chemistry for
  example
• In the last term the double brackets signify
  statistical sampling of the solute in vacuum and
  the solvent in the condensed phase, but
  uncoupled. Then those uncoupled configurations
  are overlaid and the interaction energy is
  computed

10
PDT contd

• Lets go back to an atomic solute, like argon as
  an example. Rearranging the PDT, we see
• Thus, at equilibrium, the density and the
  insertion probability are proportional
• If there are mainly unfavorable (positive)
  interaction energies, then the density will be
  low, or if there are many favorable interactions,
  then the density will be high

Exercise derive the above formula for the density
11
PDT contd

• We could imagine an inhomogeneous system of a
  biological bilayer membrane, with charged or
  polar head groups interacting with water on each
  side and a nonpolar domain in the middle
• Then the free energy profile for an argon
  solute will be unfavorable in water and near the
  head groups, but favorable in the nonpolar region

12
PDT contd

• Formula for inhomogeneous systems
• Here the potential is an external potential
  imposed say by distributions of charges in the
  system

Exercise derive the Nernst formula for the
distribution of charges between two phases kept
at different potentials
13
PDT contd

• Partition function perspective

Exercise derive the inverse formula below
14
PDT contd

• Formula for averages

Exercise prove this formula
15
Solute partitioning

• Imagine a water sample with vapor above it
  containing some Ar atoms, what is the solubility
  of Ar in water?

16
Chemical equilibrium

• Basic equation
• Then plugging in our formula for the chemical
  potentials
• If theres no interactions with the solvent (that
  is were in the gas phase)

Exercise derive these chemical equilibrium
formulas from the PDT
17
Enthalpies and entropies

• How do we get (partial molar) enthalpies and
  entropies from the chemical potential?
  Temperature derivatives

Exercise prove these formulas and show that the
same formulas hold for the excess quantities
18
Quasi-chemical theory

• This gives a way to partition the excess chemical
  potential into various manageable parts
• It partitions the FE into inner-shell and
  outer shell parts. The nice thing is then that
  we might use different approximations for those
  two regions, say quantum chemistry for the inner
  shell and a classical treatment for the outer
  shell

19
QCT contd

• Basic quasi-chemical idea

Molecule
OS
IS
20
QCT contd

• The outer-shell part can be further partitioned
  into a packing part and a long-ranged interaction
  part.
• Important to note, this partitioning is exact, no
  approximations yet
• Next well derive the QCT

21
QCT contd

• We use an often helpful trick in statistical
  mechanics of multiplying and dividing by the same
  thing, then rearranging

Here UN is the interaction energy for the N
solvent molecules. Now multiply by
The interaction energy labeled with HS is for a
hard sphere molecule of the size of the inner
shell radius.
22
QCT contd

• Rearrange to get

Then
Exercise go through all these steps in the
derivation
23
QCT contd

• What are the 3 terms in the excess mu?
• First is the inner shell (IS) term. x0 says
  what is probability that IS region is not
  occupied with any solvent while the solute is in
  the sampling. Well see later that this is like
  a chemical binding contribution.
• p0 is the probability that the whole IS region is
  unoccupied with solvent. This yields the outer
  shell (OS) packing contribution.
• The last term is due to long ranged interactions
  of the solute with the solvent. As written it
  involves sampling with the HS (hard sphere)
  solute involved, that is all solvent molecules
  are pushed away from the solute.

24
QCT contd

• Lets look at the IS term. It is minus the work
  to push the solvent molecules out of the IS
  region away from the solute. We can express it in
  chemical equilibrium language.
• Here S is the solute, W is water, and SWn is the
  complex in solution

25
QCT contd

• Then
• But note that the sum in the denom truncates
  sharply -- we can only pack so many solvent
  molecules in the IS.
• We can rewrite as

26
QCT contd

• We define an equilibrium constant
• So
• And we get
• Thus

Exercise go through the steps of this derivation
27
QCT contd

• Notice we can combine the last two terms

That combination is the total OS excess mu.
So
Exercise check this formula
28
QCT contd

• Now we focus on an ion in water, and suggest that
  likely the solvation structure around that ion is
  dominated by one or two structures. Maybe it is
  a K ion and there are 4-6 waters around it.
• There is a nice trick to view the problem in a
  different way, using our chemical approach
• We look at the IS term and say that one term in
  the sum dominates, and is much greater than 1.

29
QCT contd

• Now lets look at our equilibrium constant

The equilibrium constant in the last line is in
the gas phase
30
QCT contd

• Now we notice that when we assemble all this to
  get the full excess mu, the OS terms cancel
  exactly, and we get
• The only approximation here was to say one term
  in the chemical equilibrium dominates!
• We can try this for several ns and see which is
  most stable in terms of FE.
• What do we need to calculate to implement this
  expression?

Exercise work through all these steps
31
QCT contd

• We need to calculate the gas phase equilibrium
  constant for the complex. That could be obtained
  from a quantum chemistry code such as Gaussian,
  NWCHEM, etc.
• We know the excess chemical potential of water in
  water, about -6 kcal/mol
• BUT we have created another problem, what is the
  excess mu for the ion/water complex?

Asthagiri and Pratt, Chem. Phys. Lett. 371, 613
(2003). This paper examined the Be2 ion in water.
32
QCT contd

• The assumption of Asthagiri and Pratt was that
  the IS chemistry should be handled accurately,
  but the solvation FE of the whole complex in
  water could be treated with a continuum solvation
  model.
• Alternatively that FE could be more accurately
  estimated by MD simulations of the complex in
  water.
• At any rate, Pratt et al have computed very
  accurate solvation FEs of (mainly positive) ions
  with this approach.

33
Be2 data
Appears n4 is most stable complex
34
QCT contd

• We can also calculate x0 and p0 directly from
  simulation, see PRE 68, 041505 (2003) for AIMD
  water mu models.
• What about the long-ranged part of the OS excess
  mu?
• We can develop FE bounds for that term, which
  turn out to be helpful

35
Long ranged term

• Ways to express the long ranged term

So
Exercise derive the above two expression for the
lr part.
36
Cumulant expansions

• Cumulant expansions express approximations to the
  log of the average of an exponential
• Now we expand the log

Here
37
Cumulant expn contd

• We get
• Similarly
• From these two expressions we can derive upper
  and lower bounds for the long ranged part of the
  OS term

Exercise confirm the above cumulant expansions
38
Long ranged bounds

• Approximate expressions for the long ranged part

The fluctuation term is a variance which must
be positive. Thus
Exercise discuss why these averages give bounds.
39
Long ranged bounds

• If the sampling is largely gaussian, then we can
  average these two approximations to get a good
  estimate, since the fluctuation terms cancel

This average of mean-field terms is easy to
calculate and converges quickly. The larger is
the HS radius for the IS, the more gaussian the
sampling becomes for this OS term.
40
Data for CH4 and CF4
D. M. Rogers and T. L. Beck (in preparation)
41
Data for water
D. M. Rogers and T. L. Beck (in preparation)
42
Data for Na and Cl-
D. M. Rogers and T. L. Beck (in preparation)