## Equilibrium simulator

I’m taking a chemistry class right now on acid-base equilibria, and one problem that keeps getting assigned is one where you’re given some initial concentrations, and they want you to find the equilibrium concentrations.

For example, here’s one: Consider 1.00 L of a 0.082 M solution of aqueous formic acid (HCO2H), where $K_a= 1.78 \times 10^{-4}$. What are the equilibrium concentrations of HCO2H, HCO, H3O+, and OH?

There are two relevant reactions (we know their equilibrium constants):

We can set up an ICE table with the initial, change in, and equilibrium concentrations ([X] denotes concentration of chemical X):

[HCO2H] [HCO] [H3O+] [OH]
initial $$0.082$$ $$10^{-7}$$ $$0$$ $$10^{-7}$$
change $$-x$$ $$x + y$$ $$x$$ $$y$$
equilibrium $$0.082-x$$ $$10^{-7}+x + y$$ $$x$$ $$10^{-7}+y$$

And then we can plug these final concentrations into our two equilibrium constant formulas:

We can solve this system by either making some approximations or shoving it all into something like Mathematica. If we do that, we get [HCO2H] = 0.0783 M, [HCO] = 0.00373 M, [H3O+] = 0.00373 M, and [OH] = 2.68 × 10-12 M.

But that’s annoying. Being a programmer, I thought it would be a good idea to try and write a program that numerically solves these problems. It was a fun distraction from my actual homework! It works by pushing each reaction in the right direction until its reaction quotient becomes its equilibrium constant.The below is the result (Chrome is best — Firefox doesn’t support the HTML5 <input type="range"> element!):

Equilibrium simulator Reactions
Concentrations

Here is the code:

HTML

<fieldset id="equilibria">
<legend>Equilibrium simulator <button id="eq-start">Start</button></legend>

<b>Reactions</b>
<ul id="eq-equations"></ul>

<b>Concentrations</b>
<div id="eq-concentrations"></div>
</fieldset>


JavaScript

// Initial concentrations of each chemical (in mol/L)
var concentrations = {
'HCO2H': 0.082,
'HCO2-': 0,
'H3O+': 1e-7,
'OH-': 1e-7
};

// The equations describing the equilibria and their equilibrium constants
var equilibria = [
{
'left': ['H2O'],
'right': ['H3O+', 'OH-'],
'K': 1e-14
},
{
'left': ['HCO2H', 'H2O'],
'right': ['HCO2-', 'H3O+'],
'K': 1.78e-4
},
];

// Whether the simulation is running
var going = false;

// Some helper functions
function p(x) { return -Math.log(x) / Math.LN10; }
function e(x) { return Math.pow(10, -x); }

function product(a, b) { return a * b; }
function min(a, b) { return a < b ? a : b; }
function max(a, b) { return a > b ? a : b; }

// Returns the concentration of a chemical
function c(chem) { return concentrations[chem]; }

// Returns true or false depending on whether the reaction quotient
// depends on this chemical (that is, does its activity = 1?)
function acts(chem) { return chem in concentrations; }

// Takes an equation and pushes it towards equilibrium
function equilibriate(eq) {
var left = eq.left.filter(acts);
var right = eq.right.filter(acts);

// Calculate Q
var quotient = right.map(c).reduce(product, 1) / left.map(c).reduce(product, 1);
var reactants = quotient < eq.K ? left : right;
var products = quotient < eq.K ? right : left;
eq.updateQ(quotient);

// Calculate how close we are to equilibrium and adjust the step size.
// This formula is kind of hacked together, and it breaks on certain reactions.
// Any suggestions for how to implement how big the step size should be?
var closeness = Math.min(Math.abs(eq.K / quotient - 1), 1);
var diff = (reactants.length ?
reactants.map(c).reduce(min) :
Math.pow(products.map(c).reduce(product), 1 / products.length))
/ 2 * closeness * closeness;

// Convert some reactants into products
reactants.forEach(function (chem) { concentrations[chem] -= diff; });
products.forEach(function (chem) { concentrations[chem] += diff; });
}

// Display all the outputs... this is some ugly code; don't judge!

// Write out each equilibrium reaction
equilibria.forEach(function (eq, i) {
var li = document.createElement('li');
var quotient = document.createElement('span');

li.appendChild(document.createTextNode(eq.left.join(' + ') + ' \u21C4 ' + eq.right.join(' + ') + '; K = ' + eq.K + ', Q = '));
li.appendChild(quotient);
document.getElementById('eq-equations').appendChild(li);

eq.updateQ = function (q) { quotient.innerHTML = q; }
});

// Make all the concentration outputs
var updateConcentration = {};
Object.getOwnPropertyNames(concentrations).forEach(function (chem) {
var display = document.createElement('span');
var slider = document.createElement('input');
slider.type = 'range';
slider.min = 0;
slider.max = 14;
slider.step = 0.1;
slider.onchange = function () {
concentrations[chem] = e(this.value);
};

var div = document.createElement('div');
div.appendChild(document.createTextNode('[' + chem + '] = '));
div.appendChild(display);
div.appendChild(document.createTextNode(' mol/L'));
div.appendChild(document.createElement('br'));
div.appendChild(slider);

updateConcentration[chem] = function () {
display.innerHTML = c(chem);
slider.value = p(c(chem));
};

document.getElementById('eq-concentrations').appendChild(div);
});

document.getElementById('eq-start').onclick = function () {
this.innerHTML = going ? 'Start' : 'Stop';
going = !going;
};

// Start the simulation!
setInterval(function () {
if (going)
equilibria.forEach(equilibriate);

for (var chem in concentrations)
updateConcentration[chem]();
}, 50);


I think it’s pretty cool that it gets the right answers! The thing I ran into the most trouble with was how much reactant to convert to product at each step — it’s hard to pick one because the concentrations span many orders of magnitude. If anyone has any ideas, I’d love to know in the comments below!

February 13, 2013, 8:49pm by Casey