Tag Archives: Computational

Computational Chemistry is the Spice of Life

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("EX: python compchem.py 300 16.00 1000 1000 200 True")
# 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
if sys.argv[6].lower() == "true": # Whether or not to check for collisions between molecules
# 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
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) 
                    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
            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.")
                    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
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) + "%")
    print("Percent error was too great; Ideal Gas Law did not hold.")
    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
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!

Mmmmmm. Semiconductor chemistry simulating physics. I wonder if the transistors in my CPU feel like actors.

Mmmmmm. Semiconductor chemistry simulating more chemistry. I wonder if the transistors in my CPU feel like actors.

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!

Not bad for a rough approximation that uses "10^-20 moles of gas" and runs for "0.0001 seconds!"

Not bad for a rough approximation that uses “10^-20 moles of gas” and runs for “0.0001 seconds!”