An Ideal Gas Monte Carlo

1. An Ideal Gas Monte Carlo Vyassa Baratham, Stony Brook University April 20, 2013, 2:20-3:20pm cSplash 2013

2. What is a montecarlo? • Q: If I gave you an unfair coin, how would you determine its bias? • A: Flip it several times, find the ratio of heads to tails • What you’ve done is a Monte Carlo • Not exact, but can get as close as you want (Law of Large Numbers, stay tuned)

3. What is a Monte Carlo? An Example • In the spirit of the name “Monte Carlo,” how would you find the probability of getting a full house with triple aces in Texas Hold’em with 3 other players? • 1. Calculate it from first principles of probability • 2. Deal lots of random games and count the number that gave trip aces and a pair

4. Our Monte Carlo – ideal Gasses • How does pressure vary as a function of volume, number of molecules (N, not n), and temperature in an ideal gas? • 1. Calculate it from first principles of mechanics (statistical mechanics) • 2. Simulate lots of particles moving and interacting according to the assumptions of kinetic theory and the laws of physics

5. Assumptions of Kinetic Theory • Gasses consist of pointlike molecules with negligible volume but significant mass • All molecules have the same mass • The number of particles is large (we can use statistics; the Monte Carlo method is justified) • Molecules interact elastically (with each other, and with the walls of the container) • Elastic = kinetic energy is conserved • Molecules are in random motion, the x, y, and z components of their velocities follow normal distributions http://quantumfreak.com/

6. Normal Distributions (Bell Curve) • There are more particles with velocity near the mean (μ, average) • Histogram the particle velocities, and they will look like the graph below • Spread is determined by standard deviation (σ) • Probability density function (pdf): • Cumulative distribution function (cdf): • When we (randomly) choose the velocity of • our simulated particles, they must follow this • distribution (histogram must look like bell curve) http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc51.htm

7. Law of Large numbers • Sample mean vs population mean • As n increases, the sample mean approaches the population mean • Aside: Central Limit Theorem: • Regardless of population distribution, the histogram of observed sample means follows a normal distribution

8. Initialization • Define physical parameters • chamber size • number of particles • particle mass • temp • Define simulation parameters • length of simulation • size of timestepΔt • Randomly choose particle positions and velocities • Position: uniform distribution for each component • Velocity: normal distribution (bell curve) for each component: , • Note: All the information about the state of the system is entirely contained in the positions and velocities of the particles

9. Initialization - Code chamber_size = .05 #meters a_tot = 24 * chamber_size ** 2 #m^2 (6 sides, each squares of side 2*chamber_size) num_particles = 100000 particle_mass = 6.6335e-26 #kg, mass of argon temp = 300 #kelvin k = 1.3806503e-23 #boltzmann constant, in m^2 kg s^-2 particle_pos = [[uniform(-chamber_size, chamber_size) for i inrange(3)] for i inrange(num_particles)] #meters particle_vel = [[normalvariate(0, sqrt(k*temp/particle_mass)) for i inrange(3)] for i inrange(num_particles)] #m/s particle_vel: [[vx1, vy1, vz1], [vx2, vy2, vz2], ... [vxn, vyn, vzn]] particle_pos: [[x1, y1, z1], [x2, y2, z2], ... [xn, yn, zn]] How do I access the x coordinate of the 84th particle? How do I access the z component of the 4th particle’s velocity?

10. Simulation: What Happens in Δt Seconds? • Each particle moves a distance of (vx Δt) in the x direction, (vyΔt) in the y direction, (vzΔt) in the z direction • Code: particle_pos[i][coord] += particle_vel[i][coord] * delta_t xf = xi + (vx, vy, vz)∙Δt xi

11. Simulation: What Happens in Δt Seconds? • Particles collide with walls, imparting momentum (impulse) to them • Q: What would happen in our simulation? • A: The particle would go through the wall! Fix this: (x’f, y’f) (xf, yf) x’f= xwall – (xf – xwall) = 2*xwall – xf y'f = yf (xi, yi) xwall

12. Simulation: What Happens in Δt Seconds? • Particles collide with walls, imparting momentum (impulse) to them • Q: What is the momentum imparted to the wall by one collision? • A: It is equal to the change in momentum of the particle, ie 2 * m * vx v’xf= -vxf= -vxi v'yf = vyf= vyi Δvx= -2vx Δvy = 0 (v’xf, v’yf) (vxf, vyf) (vxi, vyi) xwall

13. Simulation: What Happens in Δt Seconds? • ifparticle_pos[i][coord] > chamber_size: • #put it back in the chamber • particle_pos[i][coord] = 2 * chamber_size - particle_pos[i][coord] • #reverse its velocity in this direction • particle_vel[i][coord] = -particle_vel[i][coord] • #add the momentum imparted to the wall of the chamber • p_tot+= 2 * particle_mass * abs(particle_vel[i][coord]) Particles collide with walls, imparting momentum (impulse) to them

14. Simulation: What Happens in Δt Seconds? • Particles collide with walls, imparting momentum (impulse) to them • Furthermore, the sum of momentum imparted to the walls by all the particles in one timestep is equal to the total impulse on the walls in that timestep • This total impulse gives us the force on the walls, which gives us the pressure (P=F/A, and we know the area because we chose the size of the container!) #calculate pressure from momentum f_tot= p_tot / delta_t pressure = f_tot / a_tot

15. Kinetic Energy

16. Results P = (4.13E-16)(1/V) With num_particles = 100000, temp = 300K, particle_mass = 3E-30kg Compare with P = NkT (1/V) NkT = 4.14E-16

17. Results P = (1.38E-15)T With num_particles = 100000, volume = .001m3, particle_mass = 3E-30kg Compare with P = (Nk/V)T Nk/V = 1.38E-15

18. Results With num_particles = 100000, volume = .001m3, temp = 500K Compare with P = (Nk/V)T Nk/V = 1.38E-15

19. Thank you! • Thanks for your attention • Please provide feedback, comments, questions, etc: vyassa.baratham@stonybrook.edu