Solving Chemical Equilibria on the Computer

Any chemistry student will tell you: solving chemical equilibria can be a real pain, especially when you have complex sets of equations. As a chemical engineer, I always wondered if I could tell the computer to solve them for me... and that's what I will be talking about today. A few weeks ago, I was at the university library consulting a few books on rocket propulsion and I ran into an old book (almost in decomposition state...) describing a method to solve chemical equilibria on a computer. The method was described on one or two pages and gave me a strong starting point but I have quickly realized that the method described was a bit incomplete regarding to several computational issues; I will cover them here. Also, I have lost the book title and author... so if you find it please send me an e-mail so I can update my post with the reference.

For those of you who are not familiar with chemical equilibria, you can see them as reactions which always tend to respect some proportions between the products and reactants of a solution. These proportions are fixed by thermodynamical properties that I will not cover here. That is, we will consider a given proportion law (available in any chemistry handbook) and see how we can find the actual products and reactants concentrations at the end of the reaction based on that law.

For example, you may be familiar with the combustion of hydrogen which is the way we send rockets into space:

which tell us that hydrogen (H₂) reacts with oxygen (O₂) to form water (H₂O). So far, we don't care about the state of the chemicals (either solid, liquid or gaseous) but only on how atoms change their configuration. In the case of hydrogen combustion, hydrogen and oxygen atoms bindings are broken and new binding are formed which gives a water molecule. The nature of the chemical change during reaction but the atom stays the same. The same thing happens for humans: we change chemicals into our bodies but the atoms stay the same. In fact, the atoms that compose your body have been unchanged for billions of years!

If you look at the previous reaction, you will see that the reactants (left chemicals of the equation) are appended with numerals (called stoichiometric constants or stoichiometric numbers) which tells how much oxygen will react with how much hydrogen to produce (right chemicals of the equation) how much water. You may multiply these numbers by any positive real value because all what matters is the relative proportions of chemicals and not their overall amount.

The previous reaction also has an implicit meaning that it is complete: all the hydrogen and oxygen will react to form water and no reactant will be left if the proportion were right. Obviously, if you have put an excess of hydrogen (for example), there will be some hydrogen remaining. This is a valid assumption up to some point but it is not the best description of what is happening. In fact, the reaction acts as an equilibrium which will keep defined proportions of reactants and products:

The reaction is reversible and can go in one direction or the other. When some conditions are met, water will spontaneously decompose into hydrogen and oxygen gas. The thermodynamical condition which tells in which direction the reaction will be driven is the free energy of the system. Free energy not only tell which direction the reaction will takes but also how strong the reaction wants to go in that direction. For instance, at room temperature, hydrogen combustion is mainly driven to complete reaction to the right of the equation and almost no hydrogen nor oxygen remains. This is why we make the former assumption that the reaction undergoes completion such that we wrote at the beginning of the post.

By doing a bit of math, free energy gives us the proportion law we were talking about. The product of all the products concentration of the reaction raised to the power of their respective stoichiometric number divided by the product of the reactants concentration raised to the power of their respective stoichiometric number should always be equal to a well defined value called the equilibrium constant that you can find into chemistry handbooks or compute through the free energy, but this is off topic.

So the proportions of hydrogen, oxygen and water of the previous reaction obey to the equilibrium law:

where the concentration of a species X is noted \[X\] and is defined as the number of moles (1 mole = 6.02e²³ molecules) per volume. So one mole of oxygen contained in a 1 litre bottle has a concentration of 1 mole/l. The relation between moles and mass is given through the molar mass of the chemical that you also find in handbooks. For example, 1 mole of oxygen weights 32 grams and 1 mole of hydrogen weights 2 grams.

Using the proportion law we may compute the final values of a reaction knowing K and the initial concentrations of hydrogen, oxygen and water. This is a classical high-school exam question for chemistry classes and it is not that hard but still takes a few lines of maths. The problem is when you start adding a second equilibrium reaction, let's say we not only have oxygen and hydrogen but we also have nitrogen (N₂) in the mixture. Nitrogen reacts with oxygen to product nitrous oxide (NO):

For which we can write a second proportion law:

So a part of oxygen will react with nitrogen and a part will react with hydrogen. The two proportion laws allow finding the final concentration of oxygen, hydrogen, nitrogen, water and nitrogen oxide but we cannot solve these equilibria independently: all the equilibria should be satisfied at the same time. To do this, we have to solve a system of five equations which is relatively troublesome. But wait, it's not over yet. Nitrous oxide also reacts with oxygen to form nitric oxide (NO₂):

which also have its proportion law:

All of this just to solve the combustion of hydrogen in air! So imagine when you add other chemicals: carbon, sulphur, chlorine... It soon becomes extremely complex. And when you think you have got the picture, things become even worse: when the chemicals react, they produce heat which increase the temperature depending on the chemicals presents and all the proportion laws constant depends on the temperature! Clearly, we cannot solve this by hand.

"Ok, ok" you will tell me. "All of this just to say that we should solve the equation sets by a computer. But that's easy!" Well... no, it's not.

There are indeed tools to solve systems of linear equations but our system here mixes linear and non-linear equations. Let us consider our former example with the combustion of hydrogen. We have three unknowns (the amount of oxygen, hydrogen and water) so we should be able to write a system of three equations to solve the problem. We have only one equation with the proportion law, the two other are mass conservation equations. The technique I'm using here might look strange to chemical students at first sight but it is easier to expand on more complex system. All we have to do is to state that the numbers of oxygen and hydrogen atoms at the beginning and at the end of the reaction are the same:

with n₀(X) the number of moles of X at the initial conditions and n_f(X) the number of moles of X at the end of the reaction.

We may linearize the proportion law by computing its logarithm but we cannot mix log of moles with moles to write our system. The technique explained in the book consists of approaching the solution step by step by computing the amount of change from one step to the other.

Rather than considering n(X) we will consider the number of moles at step i+1 as a small increase of the number of moles at step i:

The advantage is that taking the logarithm of that writing can be approximated by a linear function:

Because of the Taylor series truncated at the first term:

Using this, we may rewrite our previous system

and so,

Where I have defined

and used a short-hand notation for the total number of atoms:

The system can be re-arranged as a matrix:

What we are looking for in that system is the values of the various ξ_i(X) to compute the next estimation of n_i+1(X). If we are lucky, the values of n_i+1(X) should be slightly closer to the solution of the problem and all we have to do is to find the next ξ_i(X) until we get close enough to return a solution to the problem.

If we write the system as AX=F, we can find the value of X by inverting the matrix A: X=A^-1F. F can be seen as a driving force which gets smaller as we get closer to the solution of the equilibria set. So we may create a criterion of convergence as:

with ε some fixed tolerance. In the case of our system, this would translate to:

Now that the theory is in place, I will cover various issues I had with this framework.

Always clamp the n_i(X)

Let us consider a simple partition equilibria between a chemical into two different solvents. I will call A the chemical in the first solvent and B the chemical in the second solvent. The equilibrium is such that:

with K=B/A.

We may write the system:

And using our methodology:

or as a matrix:

which we may inverse:

You may see that something bad may happen if n_i(A)+n_i(B)=0. Also, if either n_i(A) or n_i(B) are equal to zero, we get into trouble with the logQ term. It might look like an odd case but, in my experience, it happens extremely fast if you don't pay attention to clamp the values of n_i(A) and n_i(B).

Clamping is a simple operation; all we have to do is to check that the values of n_i(X) are within a given boundary. So when updating these value we just have to add:

with:

Always saturate the ξ_i(X)

If you remember how we get our approximation to linearize the log function, we have made some implicit hypothesis that the ξ_i(X) were sufficiently small such that:

If we don't pay attention to the values of ξ_i(X) returned when solving the system, we may end up well beyond the accuracy of our linearization. We should always saturate the values of ξ_i(X) before updating the n_i(X):

with:

I recommend keeping δ well beyond 1.

Use damping

If you run the hydrogen combustion example, you will notice that the system may oscillate without ever reaching a solution. This is because the values of the ξ_i(X) are not sufficiently bounded. A solution might be to reduce δ but this is not the only way to solve the issue. A better solution, in my experience, is to use damping when updating the n_i(X) through the usage of a damping coefficient λ:

where

The damping coefficient reduces artificially the value of the ξ_i(X) by a factor that changes over time. Using the law given above, it will go from λ₀ at the first iteration to λ_∞ after some time. The larger τ the more aggressive the damping will be. Slower damping may be achieved with smaller τ values.

Playing with solubility

When we put salt (NaCl) into water, it will split into two ions Na⁺ and Cl^-. This is true for a small addition of salt only because at some point the solution will become saturated and no more ions will be able to dissolve. On the contrary, if we mix two highly concentrated solutions of Na⁺ and Cl^-, NaCl salt will immediately precipitates such that the equilibrium constant of solubility is met:

with

This proportion law is a bit different from the former and is due to the fact that we should not add pure solid nor pure liquids terms because their activities are unity. Don't look mad at me, it's a special case of the free energy thing I promised not to talk about.

Anyway, we can use our methodology to process this just like we did previously:

with:

But wait... this is true only if there is an excess of NaCl, not if we are below its saturation value.

If we are below its saturation value (when Q_i<K), the reaction acts as a complete decomposition of the salt:

for which the condition becomes:

because we know that no reactants should be left. We may then write the new equation set as:

which translates to the following matrix form:

So we have to switch from one equation system to the other depending on the condition Q_i<K.

You may wonder why I have chosen to write the reaction n_i(NaCl)=0 as n_i(NaCl)*ξ_i(NaCl)=-n_i(NaCl) and not simply as ξ_i(NaCl)=-1 which would have had the same driving force effect for the system. Actually, I did in a first place but changed to the final version to have the convergence criterion working. If we put ξ_i(NaCl)=-1, then the convergence criterion never get smaller than 1 because:

but with the new system, we get:

which may go down to zero and let us appreciate the convergence of the system.

Towards a systematic approach

So far we have derived the systems by hand but a more rigorous approach is required to have something fully automated.

We will consider a set of n chemicals X_j with r reactions between them. These reactions may be either equilibria, complete reactions or solubility equilibria. Each reaction k have a stoichiometric constant for chemical j written ν_k,j which may be null if they do not play any role in the reaction. These stoichiometric constants are defined positive for products and negative for reactants. For instance, in the combustion of hydrogen we have 3 chemicals (0=H₂, 1=O₂, 2=H₂O) with one equilibrium reaction such that ν_0,0=-1, ν_0,1=-½ and ν_0,2=+1.

We also consider a set of m atom types A_p such that the chemical X_j contains μ_j,p atoms of type p. In our hydrogen combustion example, we have 2 atom types (0=H, 1=O) such that μ_0,0=2, μ_1,0=0, μ_2,0=2, μ_0,1=0, μ_1,1=2 and μ_2,1=1.

The first thing to check is that r+m=n. If the previous equality is not met then we probably have too few reactions to be able to solve the system.

The goal is now to write a matrix Mξ=F such that we can solve for ξ using Gauss triangularization theorem (or any other method you like but it is much faster than inverting the matrix). We know that the size of M must be nxn because we have n chemicals involved. ξ will be a column matrix with ξ_j corresponding to the advance of chemical X_j.

We may first compute the coefficient from the mass balances. For each atom type p, we know that the coefficient M_j,p of the matrix is given by the number of atoms type A_p in the chemicals X_j multiplied by the number of moles of X_j:

The driving force term F_j is simply the difference between the current condition and the initial condition:

we can now append the system with the r reactions, depending on the nature of the reaction.

For equilibrium reactions, we simply put the value of M_j,m+k to ν_k,j and the term F_m+k to:

where:

and K_k is the thermodynamic constant for the reaction k which might as well be recomputed at every iteration if it is susceptible to change (effect of temperature or activity).

For complete reactions, we have to put the value of M_j,m+k and F_m+k to n_i(Xj) for the elements where ν_k,j<0 (reactants):

Finally, solubility reactions are a mix of the two previous cases:

Then, we compute the updated value n_i+1 of every chemical X_j using:

until the convergence criterion gets above the fixed threshold: