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! D. Ceperley Simulation Methods

2. Ergodicity • Fermi- Pasta- Ulam “experiment” (1954) • 1-D anharmonic chain: V= [(q i+1-q i)2+ (q i+1-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. D. Ceperley Simulation Methods

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) D. Ceperley Simulation Methods

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. D. Ceperley Simulation Methods

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. D. Ceperley Simulation Methods

6. Definition of Simulation • What is a simulation? An internal state “S” A rule for changing the state Sn+1 = T (Sn) We repeat the iteration many time. • Simulations can be • Deterministic (e.g. Newton’s equations=MD) • 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. D. Ceperley Simulation Methods

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. D. Ceperley Simulation Methods

8. Estimating Errors • Trace of A(t): • Equilibration time. • Histogramof values of A ( P(A) ). • Meanof 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)] D. Ceperley Simulation Methods

9. Statistical thinking is slippery • “Shouldn’t 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. D. Ceperley Simulation Methods

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. D. Ceperley Simulation Methods

11. The Verlet Algorithm The nearly universal choice for an integrator is the Verlet (leapfrog) algorithm r(t+h) = 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(t+h) = 2 r(t) - r(t-h) + a(t) h2 + O(h4) Add v(t) = (r(t+h) - r(t-h))/(2h) + O(h2) estimate velocities Time reversal invariance is built in the energy does not drift. 8 1 2 3 4 5 6 7 9 10 11 12 13 D. Ceperley Simulation Methods

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 don’t 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(t+h) = r(t) +h [ v(t) +(h/2) a(t)] v(t+h/2) = v(t)+(h/2) a(t) v(t+h)=v(t+h/2) + (h/2) a(t+h) D. Ceperley Simulation Methods

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) D. Ceperley Simulation Methods

14. Periodic distances Minimum Image Convention: take the closest distance:|r|M = min ( r+nL) • Potential is cutoff so that V(r)=0 for r>L/2 since force needs to be continuous. How about the derivative? • Image potential VI =  v(ri-rj+nL) • 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 D. Ceperley Simulation Methods

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. D. Ceperley Simulation Methods

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) D. Ceperley Simulation Methods

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) D. Ceperley Simulation Methods

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 energy=temperature Q= “heat bath mass”. Large Q is weak coupling D. Ceperley Simulation Methods

19. Effect of thermostat System temperature fluctuates but how quickly? Q=1 Q=100 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 D. Ceperley Simulation Methods

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). D. Ceperley Simulation Methods

21. TPN 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. D. Ceperley Simulation Methods

22. Parrinello-Rahman simulation • 500 KCl ions at 300K • First P=0 • Then P=44kB • System spontaneously changes from rocksalt to CsCl structure D. Ceperley Simulation Methods

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). D. Ceperley Simulation Methods

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,… D. Ceperley Simulation Methods

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 D. Ceperley Simulation Methods

26. Thermodynamic properties • Internal energy=kinetic 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. D. Ceperley Simulation Methods

27. D. Ceperley Simulation Methods

28. Microscopic Density (r) = < i (r-r i) > Or you can take its Fourier Transform: k = < i exp(ikri) > (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: k=2/L (nx,ny,nz)Long wave length modes are absent. • In a solid Lindemann’s ratio gives a rough idea of melting: u2= <(ri-zi)2>/d2 D. Ceperley Simulation Methods

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. D. Ceperley Simulation Methods

30. D. Ceperley Simulation Methods

31. Electron Density during exchange2d Wigner crystal (quantum) D. Ceperley Simulation Methods

32. Snapshots of densities Liquid or crystal or glass? Blue spots are defects D. Ceperley Simulation Methods

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 D. Ceperley Simulation Methods

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) = < i<j  (ri-rj-r)> (2 /N2) Density-density correlation Related to thermodynamic properties. D. Ceperley Simulation Methods

35. g(r) in liquid and solid helium • First peak is at inter-particle spacing. (shell around the particle) • goes out to r<L/2 in periodic boundary conditions. D. Ceperley Simulation Methods

36. (The static) Structure Factor S(K) • The Fourier transform of the pair correlation function is the structure factor S(k) = <|k|2>/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 D. Ceperley Simulation Methods

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) D. Ceperley Simulation Methods

38. Crystal liquid D. Ceperley Simulation Methods

39. Here is a snapshot of a binary mixture. What correlation function would be important? D. Ceperley Simulation Methods

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) D. Ceperley Simulation Methods

41. Dynamical Properties • Fluctuation-Dissipation theorem: HereA 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. D. Ceperley Simulation Methods

42. Diffusion Coefficient • Diffusion constant is defined by Fick’s 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 D. Ceperley Simulation Methods

43. Mixture simulation with CLAMPS Initial condition Later D. Ceperley Simulation Methods

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 D. Ceperley Simulation Methods