<jupyter_start><jupyter_text>Copyright **<NAME>**, January 2022 Numerical solution of equilibrium problemsLet's consider a generic chemical reaction \begin{equation}\nu_aA + \nu_bB + \dots \leftrightharpoons \nu_xX + \nu_yY + \dots \tag{1}\end{equation}if the concentrations of the species are not the equilibrium ones, the reaction will progress (either forward or backwards) to reach the equilibrium conditions, and the free energy of the system will decrease.An infinitesimal change of the system's free energy can then be written in terms of the chemical potential, $\mu$, of all the species involved in the reaction.\begin{equation}\mathrm{d}G = \sum_i \mu_i \mathrm{d}n_i \tag{2}\end{equation}Because the concentrations of reactants and products are coupled through the reaction equation, they cannot change independently, _i.e._ if any products are formed, some reactants must be consumed.We can therefore define a unique parameter, $\xi$, as the _progress of the reaction_ \begin{equation}\mathrm{d}\xi = \frac{\mathrm{d}n_i}{ \nu_i } \tag{3}\end{equation}where $\nu_i$ is the stoichiometric coefficient of the species $i$ taken as a positive number for products (formed) and as a negative number for reactants (consumed), and $\mathrm{d}n_i$ is an infinitesimal change in the concentration of species $i$.Thus, the free energy change due to the progress of the reaction (in any direction) becomes\begin{equation}\mathrm{d}G = \sum_i \nu_i \mu_i \mathrm{d}\xi \tag{4}\end{equation} At equilibrium, the concentrations of all species are constant, $\mathrm{d}n_i=0$, and the reaction is not _progressing_ anymore, $\mathrm{d}\xi=0$, hence the free energy is not changing anymore.On the other hand, if the system is out of equilibrium, there is a driving force that pulls the system towards equilibrium. A force can always be obtained as the negative of the derivative of the energy with respect to some coordinate.In this case, the driving force that pushes the reaction towards equilibrium can be immediately obtained from equation (4).\begin{equation}\frac{\mathrm{d}G}{\mathrm{d}\xi} = \sum_i \nu_i \mu_i = -force \tag{5}\end{equation}If we then substitute the definition of the chemical potential, $\mu=\mu^0+RT\ln x$, where $x$ is the molar fraction of the species, we obtain\begin{eqnarray}\frac{\mathrm{d} G}{\mathrm{d}\xi} &=& \sum_i\nu_i\mu_i\\ &=& \sum_i\nu_i \big[\mu_i^0 +RT\ln x_i\big]\\&=& \sum_i\nu_i\mu_i^0 +RT\sum_i\nu_i\ln x_i\\&=& \Delta_r G^0 +RT\sum_i\ln x_i^{\nu_i}\\&=& \Delta_r G^0 +RT\ln \prod_i x_i^{\nu_i}\\&=& \Delta_r G^0 +RT\ln Q\end{eqnarray}here $Q$ is the reaction quotient. By remembering the definition of equilibrium constant, $-RT\ln K_{eq}=\Delta_r G^0$ we get\begin{equation}force = -\frac{\mathrm{d} G}{\mathrm{d}\xi} = RT\ln K_{eq} - RT\ln Q = RT\ln\ \frac{K_{eq}}{Q} \tag{6}\end{equation} This equation provides us with a conceptual framework to solve *numerically* any equilibrium problem.Because the above equations are differential, they are valid only for infinitesimal changes in the conditions.Hence each a small change in the concentrations, $\nu_i\delta c$ the driving force may change significantly.Therefore we need to set up an iterative procedure where we compute the driving force, change the concentrations, compute the driving force\dots until the equilibrium is reached.Given the concentration of all species, we can compute the driving force using the equation above and alter the concentrations of all species using\begin{equation}[C]_i^1 = [C]_i^0 + \mathrm{d}n_i \tag{7}\end{equation}where the superscript indicates a the passage of an infinitesimal amount of time, and $\delta n_i$ is\begin{equation}\mathrm{d}n_i = \nu_i \times force \times \delta c  \tag{7}\end{equation}where $\nu_i$ are the stoichiometric coefficients of the species and $\mathrm{d}c$ is an adjustable parameter.$\mathrm{d}c$ must be small enough to maintain the validity of the approximation abover but not too small to allow for a quick convergence of the procedure.We cen then compute the new reaction quotient, the new driving force and update the concentrations again. This procedure can then be repeated a number of times, until the force is negligible, which indicates that the calculation has converged.  Example 1: Dimerisation reactionIn order to elucidate how that procedure works we can take a model dimerisation reaction\begin{equation}2A \leftrightharpoons B\end{equation}whose equilibrium constant can be written as\begin{equation}K_{eq} = \frac{[B]}{[A]^2}\end{equation}We now want to calculate the equilibrium concentrations of $[A]_{eq}$ and $[B]_{eq}$ given their initial concentrations $[A]_{0}$ and $[B]_{0}$.Although this is a simple problem that can be solved analytically, in this workshop we will learn how we can use an iterative method to numerically solve it.We will use a relatively simple minimisation procedure that can be applied to a large number of problems, for which it is not possible or it is too complicated to get an analytic solution.Imagine mixing the reagents and then to be able to monitor the concentration of all the species in the system at very short discrete time intervals (*timesteps*). What you will see is that the concentrations will change and tend to the equilibrium value. As you have learnt in first year, the reaction quotient, $Q$, can be used to decided which way to reaction will proceed, and that at equilibrium the reaction quotient is equal to the equilibrium constant. Hence, as we have discussed in class, the reaction quotient and the equilibrium constant can be use to define a *driving force* that pulls the system towards equilibrium.This *driving force* can then be used in conjunction with an *ICE* table to numerically compute the equilibrium concentration of reactant and products. The working principle of the minimisation procedure that we will employ is1. compute the reaction quotient\begin{equation}Q = \frac{\mathrm{[B]}_0}{\mathrm{[A]}^2_0}\end{equation}2. compute the driving force\begin{equation}F = RT \ln\bigg[\frac{K_{eq}}{Q}\bigg]\end{equation}3. evolve the concentrations using the ICE table|      | [A]        | [B]| :---: | :--------: |:---------:| *I*  | [A]$_0$    | [B]$_0$| *C*  | $-2F\delta c$        | $+F\delta c$| *E*  | [A]$_0$-2$F\delta c$ | [B]$_0$+$F\delta c$4. Repeat until the solution does not change  Let's try to implement this  Working notebooks<jupyter_code>import ipywidgets as ipw
import os

from IPython.display import Javascript
import glob as glob
import nbformat as nbf

label_layout = ipw.Layout(width='300px')

##########
pfiles = ['.protectedFiles.txt' , '../.protectedFiles.txt']
for fff in pfiles:
    if os.path.isfile(fff):
        with open(fff) as f:
            protectedFiles = f.read().splitlines()
##########
def launchNotebook(filename):
    text = "   var name_of_the_notebook = '" + filename + "'"
    vv="""
    var url = window.location.href.split('/')
    var newurl = url[0] + '//'
    for (var i = 1; i < url.length - 1; i++) {
        console.log(url[i], newurl)
        newurl += url[i] + '/'
    }
    newurl += name_of_the_notebook
    window.open(newurl)
    """
    text = text + vv
    display(Javascript(text))

def openNewNotebook(btn):
    if os.path.exists(notebookeName.value):
        print("Filename exists - Please select a different name")
        return
    
    nb = nbf.v4.new_notebook()
    text = """# Click 'Edit App' to start coding"""

    code = """\
# python packages
import pandas as pd # Dataframes and reading CSV files
import numpy as np # Numerical libraries
import matplotlib.pyplot as plt # Plotting library
from lmfit import Model # Least squares fitting library
from scipy.optimize import curve_fit # Alternative curve fittting library"""

    nb['cells'] = [nbf.v4.new_markdown_cell(text),
                   nbf.v4.new_code_cell(code)]

    if notebookeName.value in protectedFiles or notebookeName.value in listOfFiles:
        print("File already exists, select a different filename")
    else:
        with open(notebookeName.value, 'w') as f:
            nbf.write(nb, f)
        launchNotebook(notebookeName.value)

##########
listOfFiles = []
files = glob.glob1("./","*.ipynb")
for f in files:
    if f in protectedFiles:
        continue
    listOfFiles.append(f)

def dropdown_filesHandler(change):
    for i in range(0,len(listOfFiles)):
        if listOfFiles[i] == change.new:
            oldNotebookeName[0] = listOfFiles[i]

def createMenuFiles(data):
    option_list = ["Choose one"]
    option_list.extend(data)

    dropdown = ipw.Dropdown(description="", options=option_list, layout=ipw.Layout(width="300px"))
    dropdown.observe(dropdown_filesHandler, names='value')

    return dropdown

##########
oldNotebookeName = ["None"]
def openOldNotebook(btn):
    if oldNotebookeName[0] == "None":
        print("Please select a filename")
    elif oldNotebookeName[0] in protectedFiles:
        print("Please select a different filename")
    else:
        launchNotebook(oldNotebookeName[0])

##########
actions0 = []

notebookeName = ipw.Text("Empty.ipynb")

btn_new = ipw.Button(description="Create a new notebook", layout=label_layout)
btn_new.on_click(openNewNotebook)

btn_old = ipw.Button(description="Open an old notebook", layout=label_layout)
btn_old.on_click(openOldNotebook)

actions0.append(ipw.HBox([btn_new,notebookeName]))
actions0.append(ipw.HBox([btn_old,createMenuFiles(listOfFiles)]))
ipw.VBox(actions0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11053686320008438, 'fraction_numerical': 0.006750342790844847, 'mean_word_length': 4.384440658716638, 'pattern_counts': {'":': 1, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '173104', 'license': 'MIT', 'n_tokens_mistral': 2839, 'n_tokens_neox': 2714, 'n_words': 1067, 'path': 'week_09_chemicalEquilibrium/equilibrium.ipynb', 'repo_name': 'praiteri/TeachingNotebook', 'size': '9488'}