Actually, no, it isn’t. It’s pretty darn interesting, though.

Since AP Chemistry ran out of things to do (since we took the AP test), my teacher decided to assign our class a project: come up with a way to relate chemistry to your future career field. *Easy*, I thought. *MIT has Course 6-7 (computational biology), so some sort of chemistry equivalent (computational chemistry) must exist somewhere on the Internet*. As it turned out, I was spot on. Computational chemistry is big.

I started looking at computational chemistry software. All of it carried a steep learning curve and a steeper price tag, and I couldn’t really think of much high-level research, since I just finished AP Chemistry. Instead, I decided to take matters into my own hands. I thought it might be interesting to do some “exploration” around the Ideal Gas Law. At first, the goal of my project was simply to show that PV did indeed equal nRT. After I finished the first step, I started looking into why (hint: gas molecules aren’t very big). After that, it dawned on me that I could calculate an approximate value of R, the Universal Gas Constant.

What’s more, I wrote my program in Python:

#!/usr/bin/env python3 # Computational chemistry: does PV = nRT? # Copyright (c) 2013 John Parsons. All rights reserved. import sys import random import math import copy class Particle: """Gas molecule.""" def __init__(self, vel, pos): """Creates a new object instance. Vel: (xVelocity, yVelocity, zVelocity). Pos = (x, y, z) coordinates.""" # Set up variables self.velocity = vel self.position = pos def getListOfDistances(self, allParticlesInTank): """Gets a sorted list of distances between the current particle and all other particles. allParticlesInTank: list of all other particles.""" particles = [x for x in allParticlesInTank if x.velocity != self.velocity and x.position != self.position] # Get a shallow copy of the particle list and remove the current particle particleDistances = [] # List of tuples that stores distances for particle in particles: # Loop through particles distance = 0 # Quick way of implementing 3D distance formula for x in range(3): distance += (self.position[x] - particle.position[x])** 2 distance = distance ** 0.5 # Take the square root particleDistances.append((particle, distance)) # Add the particle - distance tuple to the list return sorted(particleDistances, key = lambda x:x[1]) # Return particles sorted by their distance from the current particle # Shameless self-promo print("compchem.py -- Copyright (c) 2013 John Parsons") print("Compares simulated gas pressure with the Ideal Gas Law.") # Print usage if user doesn't know what he's doing if len(sys.argv) < 7: print("USAGE: python compchem.py TEMPERATURE_IN_KELVIN MOLAR_MASS_IN_GRAMS NUMBER_OF_PARTICLES ITERATIONS ATOMIC_RADIUS_PICOMETERS INTERMOLECULAR_COLLISIONS_TRUE-FALSE") print("EX: python compchem.py 300 16.00 1000 1000 200 True") exit() # Constants / calculations to make the later code more verbose AVOGADRO_NUMBER = 6.022 * (10 ** 23) CUBE_SIDE_WIDTH = 1000000 # Bounds of simulation, in nanometers SURFACE_AREA_M = (6 * (CUBE_SIDE_WIDTH ** 2)) * 10**-18 NUMBER_OF_PARTICLES = int(sys.argv[3]) # User-defined number of particles to play with GAS_CONSTANT = 8.3144621 # R GAS_TEMPERATURE = float(sys.argv[1]) # T GAS_MOLAR_MASS = float(sys.argv[2]) # M (grams -- note that we must convert to KG for momentum!) GAS_PARTICLE_MASS = GAS_MOLAR_MASS / AVOGADRO_NUMBER # Convenient mass of single particle GAS_SPEED = (3 * GAS_CONSTANT * GAS_TEMPERATURE / (GAS_MOLAR_MASS / 1000)) ** 0.5 # (3rt/m)^.5, root-mean-square velocity of gas TANK_VOLUME = CUBE_SIDE_WIDTH ** 3 # Volume of container in nm^3, used later. SIMULATOR_ITERATIONS = int(sys.argv[4]) ATOMIC_RADIUS = float(sys.argv[5]) # In picometers = 10**-12m INTERMOLECULAR_COLLISIONS = False if sys.argv[6].lower() == "true": # Whether or not to check for collisions between molecules INTERMOLECULAR_COLLISIONS = True # Print all the simulation parameters (calculated and hard-coded) in a nicely formatted manner print("Simulation paramaters:\r\n\tGAS_SPEED = " + str(GAS_SPEED) + " m/s\r\n\tGAS_MOLES = " + str(NUMBER_OF_PARTICLES / AVOGADRO_NUMBER) + " moles\r\n\tTANK_VOLUME = " + str(TANK_VOLUME) + " nm^3\r\n\tATOMIC_RADIUS = " + str(ATOMIC_RADIUS) + " pm\r\n\tCHECK_INTERMOLECULAR_COLLISIONS = " + str(INTERMOLECULAR_COLLISIONS) + "\r\n\tSIMULATION_TIME = " + str(SIMULATOR_ITERATIONS * 10**-8) + " s") # Apply ideal gas law: num particles/Avogadro = n; GAS_CONSTANT = R, GAS_TEMPERATURE = T, tank_volume * 10**-24 = tank volume in m^3 IDEAL_PRESSURE = ((NUMBER_OF_PARTICLES / AVOGADRO_NUMBER) * GAS_CONSTANT * GAS_TEMPERATURE) / (TANK_VOLUME * (10 ** -27)) print("\r\nIdeal Gas Law calculated pressure: " + str(IDEAL_PRESSURE) + " Pa\r\n") def velocityVectorWithSpeed(speed): """Creates a velocity vector with random direction and specified speed.""" x = random.randint(-500, 500) # Create random vector components y = random.randint(-500, 500) z = random.randint(-500, 500) magnitude = (x**2 + y**2 + z**2)**0.5 # Magnitude of vector unitVector = (x/magnitude, y/magnitude, z/magnitude) return ([i * speed for i in unitVector]) # Return vector such that its magnitude is the specified speed def randomPosition(): """Generates a random tuple (x, y, z) that stores a position.""" return ([random.randint(0, CUBE_SIDE_WIDTH) for x in range(3)]) particles = [] # list of particles # Create all our cute little particles for i in range(0, NUMBER_OF_PARTICLES): particles.append(Particle(velocityVectorWithSpeed(GAS_SPEED), randomPosition())) sumOfMomenta = 0 # Used for storing total impulse imparted to walls of container print("Simulating: 0/" + str(SIMULATOR_ITERATIONS) + "\r", end = "") sys.stdout.flush() # Force printing # Now run the simulation -- each iteration will represent 10^-8 seconds for ii in range(0, SIMULATOR_ITERATIONS): # Do some number of 10**-8 second iterations (as defined by user) if ii % 100 == 0: # Make it pretty as hell print("Simulating: " + str(ii) + "/" + str(SIMULATOR_ITERATIONS) + "\r", end="") sys.stdout.flush() # Make it actually print for particle in particles: # Loop through all the particles repositionVector = ([x * 10 for x in particle.velocity]) # x m/s * 10**9 nm/m * 10**-8 s/iteration = nm/iteration # Loop through and change positions for i in range(3): particle.position[i] = particle.position[i] + repositionVector[i] # Adjust particle position if particle.position[i] >= CUBE_SIDE_WIDTH or particle.position[i] <= 0: # Check to see if we need to change direction particle.velocity[i] = particle.velocity[i] * -1 # We do! Now we need to figure out how much impulse is imparted if particle.position[i] >= CUBE_SIDE_WIDTH: # Prevent the particle from being turned around again by putting the particle where it should really be. particle.position[i] = CUBE_SIDE_WIDTH - (particle.position[i] - CUBE_SIDE_WIDTH) else: particle.position[i] = -particle.position[i] # p = mv (mass in kilograms). J = delta p = dp dp = (GAS_PARTICLE_MASS * 2 * math.fabs(particle.velocity[i])) / 1000 sumOfMomenta += dp # Now we do collision detection if INTERMOLECULAR_COLLISIONS: distances = particle.getListOfDistances(particles) for distance in distances: if distance[1] <= ATOMIC_RADIUS * 10**-3: # If they're too close together, we have a collision # We had a collision! print("Collision between molecules.") sys.stdout.flush() else: break # There will be no distances greater than this one print("Simulating: " + str(ii + 1) + "/" + str(SIMULATOR_ITERATIONS), end="") print("\n\r\nSum of momentum changes: " + str(sumOfMomenta) +" m*N") # We ran for (10**-8 * SIMULATOR ITERATIONS) seconds; P=F/A; mv = Ft so P = mv/tA ACTUAL_PRESSURE = sumOfMomenta / ((10**-8 * SIMULATOR_ITERATIONS) * SURFACE_AREA_M) PERCENT_ERROR = (math.fabs(IDEAL_PRESSURE - ACTUAL_PRESSURE) / IDEAL_PRESSURE) * 100 # |actual-experimental|/actual # Let the user know how we did print("Average pressure: " + str(ACTUAL_PRESSURE) + " Pa") print("Pressure error: " + str(PERCENT_ERROR) + "%") if PERCENT_ERROR > 5: print("Percent error was too great; Ideal Gas Law did not hold.") else: print("Percent error was reasonable; Ideal Gas Law held as expected.") print() # Formatting # Now calculate the simulated value of the gas constant # PV/nT = R ACTUAL_GAS_CONSTANT = (ACTUAL_PRESSURE * (TANK_VOLUME * 10**-27)) / ((NUMBER_OF_PARTICLES / AVOGADRO_NUMBER) * GAS_TEMPERATURE) GAS_CONSTANT_ERROR = (math.fabs(GAS_CONSTANT - ACTUAL_GAS_CONSTANT) / GAS_CONSTANT) * 100 # Same error formula as above print("Calculated gas constant: " + str(ACTUAL_GAS_CONSTANT) + " m^3*Pa*K^-1*mol^-1") print("Gas constant error: " + str(GAS_CONSTANT_ERROR) + "%") input("Press enter to exit.") |

I won’t go into any sort of massive description of my code. I will go through a few points, though.

- You will notice from my (likely very improvable) code that adding molecular collisions will add significantly to the running time of the program.
- You will get a significant amount of error when molecules are travelling very quickly. This is expected. The error, however, is on the wrong side of the Ideal Gas Law. It happens because particles literally escape the “test chamber” and thus are unable to impart pressure on the inside of its walls.
- Low temperatures also create significant error.
- Newtonian Mechanics are exclusively used. You won’t see Schrodinger’s Equation anywhere in this program. There are a few reasons for this. First, it’s written in Python, so it’s slow to begin with. Second, I have no idea what a partial differential equation is. I might come back to this issue next year!
- Everything in the program is based on the idea that impulse, or change in momentum, is equal to force times change in time, which is equal to mass times the change in velocity. Thus, mass times change in velocity divided by total change in time is equal to total average force. Divide that by area and you get a pressure! Woo-hoo! Physics!
- My code probably isn’t very “pythonic.” That’s because I don’t know Python very well. Yet. I’ll get better!
- Estimation quality and computing time, as they should be, are inversely related (as opposed to the direct relation to time spent doing computational chemistry and fun). The more molecules you simulate for more time, the better your results will be. Fewer molecules give crappier results. Numbers in the thousands for both iterations and particles seem to do okay.

Here’s the program in action. Note the relatively accurate estimation! Fun-fun-fun!

Oh. By the way. This trial eventually finished. FYI, the accepted value of the Universal Gas Constant is 8.3144621 m^3 Pa mol^-1 K^-1. We were off by about 0.0001–we had four sig figs of perfection!