Tag Archives: physics

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("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!

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!”

Coercing More Effectively

Coercion’s blade needs a lot of work – I think the video of Ripto bending it 30 or so degrees proved that. The question is how to change it. I know Coercion is a good bot, but I have noticed some continual flaws. First of all, the blade is too heavy. It spins up too slowly and it takes way too much work to get up to speed (thanks to losses all over the place), so the weapon motor runs hot and pulls something in the 30A neighborhood. I thought I had more or less fixed the problem before Bot Blast, but running my robot on 7.4V made it far more sluggish than I anticipated and I still wasn’t getting ideal blade performance. The blade is also so heavy that it tends to flip over its nose, which provided a loss at Boy Blast. I have decided to apply some physics to make Coercion more hard-hitting and hopefully more efficient than ever.

Erotational = Iw^2/2

Where omega (henceforth referred to as w) is angular velocity, E is kinetic energy, and I is moment of inertia.

Awesome, so based on that, we can do one of two things to increase the kinetic energy of the weapon:

  1. Increase the angular velocity (how fast the weapon is spinning). This will be the more effective of the two methods because it is a squared term, so if we double the angular velocity, we quadruple the kinetic energy.
  2. Increase the moment of inertia. MOI is the sum of the masses of the particles that make up an object multiplied by the radii from the center of rotation squared. In other words, the heavier something is, the greater its MOI, and the greater something’s diameter is, the greater its MOI.

So there we have it. To make Coercion’s blade work better, we make it heavier and spin it faster. Great physics lesson, Einstein. The problem here is that to make weight for better wheel guards, a better weapon motor, new batteries and the other improvements I have planned for Coercion, I need to make the blade lighter. Remember that MOI is also based on the distance the particles of an object are from its center of rotation. Thus, if I can concentrate the mass as far away from the center as possible, I can “cheat” (not really) a little bit: the blade will hit harder while not weighing in any more.

My target weight of the blade was 11ish ounces, which is about a big reduction over the old blade. I had to make the inside as skinny as possible while making the outside part as wide as possible. I kept the same blade dimensions – 11.75″x1.50″x0.19″ – and I plan to have this blade waterjetted out of S7 and then hardened (I learned my lesson at Bot Blast!).

New Coercion Blade

10.35oz, far lighter than the old, nearly-one-pound blade.

The plan is to spin this blade at 10,000rpm with a 3S LiPoly and a new weapon motor in the 1500-1900kv range (and a new, 2:1 weapon gear ratio). If I recall correctly, Coercion’s blade was spinning at about 7000rpm at Bot Blast, and according to the MOI differences I saw in Solidworks (the new one has about a 22% lower MOI than the old one), the new blade will actually store about 50% more energy. Not bad for saving three or four ounces, huh?