Statistical Ensembles

1
Statistical Ensembles

• Classical phase space is 6N variables (pi, qi)
  and a Hamiltonian function H(q,p,t).
• We may know a few constants of motion such as
  energy, number of particles, volume...
• Ergodic hypothesis each state consistent with
  our knowledge is equally likely the
  microcanonical ensemble.
• Implies the average value does not depend on
  initial conditions.
• A system in contact with a heat bath at
  temperature T will be distributed according to
  the canonical ensemble
• exp(-H(q,p)/kBT )/Z
• The momentum integrals can be performed.
• Are systems in nature really ergodic? Not always!

2
Ergodicity

• Fermi- Pasta- Ulam experiment (1954)
• 1-D anharmonic chain V ?(q i1-q i)2? (q
  i1-q i)3
• The system was started out with energy with the
  lowest energy mode. Equipartition would imply
  that the energy would flow into the other modes.
• Systems at low temperatures never come into
  equilibrium. The energy sloshes back and forth
  between various modes forever.
• At higher temperature many-dimensional systems
  become ergodic.
• Area of non-linear dynamics devoted to these
  questions.

3

• Let us say here that the results of our
  computations were, from the beginning, surprising
  us. Instead of a continuous flow of energy from
  the first mode to the higher modes, all of the
  problems show an entirely different behavior.
  Instead of a gradual increase of all the higher
  modes, the energy is exchanged, essentially,
  among only a certain few. It is, therefore, very
  hard to observe the rate of thermalization or
  mixing in our problem, and this wa s the initial
  purpose of the calculation.
• Fermi, Pasta, Ulam (1954)

4

• Equivalent to exponential divergence of
  trajectories, or sensitivity to initial
  conditions. (This is a blessing for numerical
  work. Why?)
• What we mean by ergodic is that after some
  interval of time the system is any state of the
  system is possible.
• Example shuffle a card deck 10 times. Any of the
  52! arrangements could occur with equal
  frequency.
• Aside from these mathematical questions, there is
  always a practical question of convergence. How
  do you judge if your results converged? There is
  no sure way. Why? Only experimental tests.
• Occasionally do very long runs.
• Use different starting conditions.
• Shake up the system.
• Compare to experiment.

5
Statistical ensembles

• (E, V, N) microcanonical, constant volume
• (T, V, N) canonical, constant volume
• (T, P N) constant pressure
• (T, V , ?) grand canonical
• Which is best? It depends on
• the question you are asking
• the simulation method MC or MD (MC better for
  phase transitions)
• your code.
• Lots of work in last 2 decades on various
  ensembles.

6
Definition of Simulation

• What is a simulation?
• An internal state S
• A rule for changing the state Sn1 T (Sn)
• We repeat the iteration many time.
• Simulations can be
• Deterministic (e.g. Newtons equationsMD)
• Stochastic (Monte Carlo, Brownian motion,)
• Typically they are ergodic there is a
  correlation time T. for times much longer than
  that, all non-conserved properties are close to
  their average value. Used for
• Warm up period
• To get independent samples for computing errors.

7
Problems with estimating errors

• Any good simulation quotes systematic and
  statistical errors for anything important.
• Central limit theorem after enough averaging,
  any statistical quantity approaches a normal
  distribution.
• One standard deviation means 2/3 of the time the
  correct answer is within ? of the estimate.
• Problem in simulations is that data is correlated
  in time. It takes a correlation time to be
  ergodic
• We must throw away the initial transient and
  block successive parts to estimate the mean value
  and error.
• The error and mean are simultaneously determined
  from the data. We need at least 20 independent
  data points.

8
Estimating Errors

• Trace of A(t)
• Equilibration time.
• Histogram of values of A ( P(A) ).
• Mean of A (a).
• Variance of A ( v ).
• estimate of the mean ?A(t)/N
• estimate of the variance,
• Autocorrelation of A (C(t)).
• Correlation time (k ).
• The (estimated) error of the (estimated) mean (s
  ).
• Efficiency 1/(CPU time error 2)

9
Statistical thinking is slippery

• Shouldnt the energy settle down to a constant
• NO. It fluctuates forever. It is the overall
  mean which converges.
• My procedure is too complicated to compute
  errors
• NO. Run your whole code 10 times and compute the
  mean and variance from the different runs
• The cumulative energy has converged.
• BEWARE. Even pathological cases have smooth
  cumulative energy curves.
• Data set A differs from B by 2 error bars.
  Therefore it must be different.
• This is normal in 1 out of 10 cases.

10
Characteristics of simulations.

• Potentials are highly non-linear with
  discontinuous higher derivatives either at the
  origin or at the box edge.
• Small changes in accuracy lead to totally
  different trajectories. (the mixing or ergodic
  property)
• We need low accuracy because the potentials are
  not very realistic. Universality saves us a
  badly integrated system is probably similar to
  our original system. This is not true in the
  field of non-linear dynamics or, in studying the
  solar system
• CPU time is totally dominated by the calculation
  of forces.
• Memory limits are not too important.
• Energy conservation is important roughly
  equivalent to time-reversal invariance. allow
  0.01kT fluctuation in the total energy.

11
The Verlet Algorithm

• The nearly universal choice for an integrator is
  the Verlet (leapfrog) algorithm
• r(th) r(t) v(t) h 1/2 a(t) h2 b(t) h3
  O(h4) Taylor expand
• r(t-h) r(t) - v(t) h 1/2 a(t) h2 - b(t) h3
  O(h4) Reverse time
• r(th) 2 r(t) - r(t-h) a(t) h2 O(h4) Add
• v(t) (r(th) - r(t-h))/(2h) O(h2) estimate
  velocities
• Time reversal invariance is built in ? the
  energy does not drift.

8
2
3
4
5
1
6
7
9
10
11
12
13
12
How to set the time step

• Adjust to get energy conservation to 1 of
  kinetic energy.
• Even if errors are large, you are close to the
  exact trajectory of a nearby physical system with
  a different potential.
• Since we dont really know the potential surface
  that accurately, this is satisfactory.
• Leapfrog algorithm has a problem with round-off
  error.
• Use the equivalent velocity Verlet instead
• r(th) r(t) h v(t) (h/2) a(t)
• v(th/2) v(t)(h/2) a(t)
• v(th)v(th/2) (h/2) a(th)

13
Spatial Boundary Conditions

• Important because spatial scales are limited.
  What can we choose?
• No boundaries e.g. droplet, protein in vacuum.
  If droplet has 1 million atoms and surface layer
  is 5 atoms thick? 25 of atoms are on the
  surface.
• Periodic Boundaries most popular choice because
  there are no surfaces (see next slide) but there
  can still be problems.
• Simulations on a sphere
• External potentials
• Mixed boundaries (e.g. infinite in z, periodic in
  x and y)

14
Periodic distances

• Minimum Image Convention take the closest
  distancerM min ( rnL)
• Potential is cutoff so that V(r)0 for rgtL/2
  since force needs to be continuous. How about
  the derivative?
• Image potential
• VI ? v(ri-rjnL)
• For long range potential this leads to the Ewald
  image potential. You need a back ground and
  convergence method (more later)

-L -L/2 0 L/2 L
x
15
Complexity of Force Calculations

• Complexity is defined as how a computer algorithm
  scales with the number of degrees of freedom
  (particles)
• Number of terms in pair potential is N(N-1)/2 ?
  O(N2)
• For short range potential you can use neighbor
  tables to reduce it to O(N)
• (Verlet) neighbor list for systems that move
  slowly
• bin sort list (map system onto a mesh and find
  neighbors from the mesh table)
• Long range potentials with Ewald sums are O(N3/2)
  but Fast Multipole Algorithms are O(N) for very
  large N.

16
Constant Temperature MD

• Problem in MD is how to control the temperature.
  (BC in time.)
• How to start the system? (sample velocities from
  a Gaussian distribution) If we start from a
  perfect lattice as the system becomes disordered
  it will suck up the kinetic energy and cool down.
  (v.v for starting from a gas)
• QUENCH method. Run for a while, compute kinetic
  energy, then rescale the momentum to correct
  temperature, repeat as needed.
• Nose-Hoover Thermostat controls the temperature
  automatically by coupling the microcanonical
  system to a heat bath
• Methods have non-physical dynamics since they do
  not respect locality of interactions. Such
  effects are O(1/N)

17
Quench method

• Run for a while, compute kinetic energy, then
  rescale the momentum to correct temperature,
  repeat as needed.
• Control is at best O(1/N)

18
Nose-Hoover thermostat

• MD in canonical distribution (TVN)
• Introduce a friction force ?(t)

SYSTEM
T Reservoir
Dynamics of friction coefficient to get canonical
ensemble.
Feedback restores makes kinetic energytemperature
Q heat bath mass. Large Q is weak coupling
19
Effect of thermostat

• System temperature fluctuates but how quickly?
• Q1
• Q100

DIMENSION 3 TYPE argon 256 48. POTENTIAL argon
argon 1 1. 1. 2.5 DENSITY 1.05 TEMPERATURE
1.15 TABLE_LENGTH 10000 LATTICE 4 4 4 4 SEED
10 WRITE_SCALARS 25 NOSE 100. RUN MD 2200 .05
20

• Thermostats are needed in non-equilibrium
  situations where there might be a flux of energy
  in or out of the system.
• It is time reversable, deterministic and goes to
  the canonical distribution but
• How natural is the thermostat?
• Interactions are non-local. They propagate
  instantaneously
• Interaction with a single heat bath
  variable-dynamics can be strange. Be careful to
  adjust the mass
• REFERENCES
• S. Nose, J. Chem. Phys. 81, 511 (1984) Mol.
  Phys. 52, 255 (1984).
• W. Hoover, Phys. Rev. A31, 1695 (1985).

21
Constant Pressure

• To generalize MD, follow similar procedure as for
  the thermostat for constant pressure. The size
  of the box is coupled to the internal pressure

• Volume is coupled to virial pressure
• Unit cell shape can also change.

22
Parrinello-Rahman simulation

• 500 KCl ions at 300K
• First P0
• Then P44kB
• System spontaneously changes from rocksalt to
  CsCl structure

23

• Can automatically find new crystal structures
• Nice feature is that the boundaries are flexible
• But one is not guaranteed to get out of local
  minimum
• One can get the wrong answer. Careful free
  energy calculations are needed to establish
  stable structure.
• All such methods have non-physical dynamics since
  they do not respect locality of interactions.
• Non-physical effects are O(1/N)
• REFERENCES
• H. C. Andersen, J. Chem. Phys. 72, 2384 (1980).
• M. Parrinello and A. Rahman, J. Appl. Phys. 52,
  7158 (1981).

24
Brownian dynamics

• Put a system in contact with a heat bath
• Will discuss how to do this later.
• Leads to discontinuous velocities.
• Not necessarily a bad thing, but requires some
  physical insight into how the bath interacts with
  the system in question.
• For example, this is appropriate for a large
  molecule (protein or colloid) in contact with a
  solvent. Other heat baths in nature are given by
  phonons and photons,

25
Monitoring the simulation

• Static properties pressure, specific heat etc.
• Density
• Pair correlation in real space and fourier space.
• Order parameters How to tell a liquid from a
  solid

26
Thermodynamic properties

• Internal energykinetic energy potential energy
• Kinetic energy is kT/2 per momentum
• Specific heat mean squared fluctuation in
  energy
• pressure can be computed from the virial theorem.
• compressibility, bulk modulus, sound speed
• But we have problems for the basic quantities of
  entropy and free energy since they are not ratios
  with respect to the Boltzmann distribution. We
  will discuss this later.

27
(No Transcript)
28
Microscopic Density

• ?(r) lt ?i ?(r-r i) gt
• Or you can take its Fourier Transform
• ? k lt ?i exp(ikri) gt
• (This is a good way to smooth the density.)
• A solid has broken symmetry (order parameter).
  The density is not constant.
• At a liquid-gas transition the density is also
  inhomgeneous.
• In periodic boundary conditions the k-vectors are
  on a grid k2?/L (nx,ny,nz) Long wave length
  modes are absent.
• In a solid Lindemanns ratio gives a rough idea
  of melting
• u2 lt(ri-zi)2gt/d2

29
Order parameters

• A system has certain symmetries translation
  invariance.
• At high temperatures one expect the system to
  have those same symmetries at the microscopic
  scale. (e.g. a gas)
• BUT as the system cools those symmetries are
  broken. (a gas condenses).
• At a liquid gas-transition the density is no
  longer fixed droplets form. The density is the
  order parameter.
• At a liquid-solid transition, both rotational
  symmetry and translational symmetry are broken
• The best way to monitor the transition is to look
  for the dynamics of the order parameter.

30
(No Transcript)
31
Electron Density during exchange2d Wigner
crystal (quantum)
32
Snapshots of densities

• Liquid or crystal or glass? Blue spots are
  defects

33
Density distribution within a helium droplet

• During addition of molecule, it travels from the
  surface to the interior.

Red is high density, blue low density of helium
34
Pair Correlation Function, g(r)

• Primary quantity in a liquid is the probability
  distribution of pairs of particles. Given a
  particle at the origin what is the density of
  surrounding particles
• g(r) lt ?iltj ? (ri-rj-r)gt (2 ?/N2)
• Density-density correlation
• Related to thermodynamic properties.

35
g(r) in liquid and solid helium

• First peak is at inter-particle spacing. (shell
  around the particle)
• goes out to rltL/2 in periodic boundary
  conditions.

36
(The static) Structure Factor S(K)

• The Fourier transform of the pair correlation
  function is the structure factor
• S(k) lt?k2gt/N (1)
• S(k) 1 ? ?dr exp(ikr) (g(r)-1) (2)
• problem with (2) is to extend g(r) to infinity
• This is what is measured in neutron and X-Ray
  scattering experiments.
• Can provide a direct test of the assumed
  potential.
• Used to see the state of a system
• liquid, solid, glass, gas? (much better than g(r)
  )
• Order parameter in solid is ?G

37

• In a perfect lattice S(k) will be non-zero only
  on the reciprocal lattice vectors G S(G) N
• At non-zero temperature (or for a quantum system)
  this structure factor is reduced by the
  Debye-Waller factor
• S(G) 1 (N-1)exp(-G2u2/3)
• To tell a liquid from a crystal we see how S(G)
  scales as the system is enlarged. In a solid,
  S(k) will have peaks that scale with the number
  of atoms.
• The compressibility is given by
• We can use this is detect the liquid-gas
  transition since the compressibility should
  diverge as k approaches 0. (order parameter is
  density)

38
Crystal liquid
39
Here is a snapshot of a binary mixture. What
correlation function would be important?
40

• In a perfect lattice S(k) will be non-zero only
  on the reciprocal lattice vectors S(G) N
• At non-zero temperature (or for a quantum system)
  this structure factor is reduced by the
  Debye-Waller factor
• S(G) 1 (N-1)exp(-G2u2/3)
• To tell a liquid from a crystal we see how S(G)
  scales as the system is enlarged. In a solid,
  S(k) will have peaks that scale with the number
  of atoms.
• The compressibility is given by
• We can use this is detect the liquid gas
  transition since the compressibility should
  diverge. (order parameter is density)

41
Dynamical Properties

• Fluctuation-Dissipation theorem
• Here A e-iwt is a perturbation and ? (w) e-iwt
  is the response of B. We calculate the average on
  the lhs in equilibrium (no external
  perturbation).
• Fluctuations we see in equilibrium are
  equivalent to how a non-equilibrium system
  approaches equilibrium. (Onsager regression
  hypothesis)
• Density-Density response function is S(k,w) can
  be measured by scattering and is sensitive to
  collective motions.

42
Diffusion Coefficient

• Diffusion constant is defined by Ficks law and
  controls how systems mix
• Microscopically we calculate
• Alder-Wainwright discovered long-time tails on
  the velocity autocorrelation function so that the
  diffusion constant does not exist in 2D

43
Mixture simulation with CLAMPS
Initial condition Later
44
Transport Coefficients

• Diffusion Particle flux
• Viscosity Stress tensor
• Heat transport energy current
• Electrical Conductivity electrical current
• These can also be evaluated with non-equilibrium
  simulations use thermostats to control.
• Impose a shear flow
• Impose a heat flow