physical chemistry - How to obtain product ratio from energy differences via Boltzmann statistics

I am currently calculating some reactions and I would like to express the ratio of the products in per cents. I am using the Boltzmann statistic in the following way:

1. 
   

   The probability $p_j$ of a particle $j$ to be in a state with the energy $\varepsilon_j$ is given as $\text{(1a,a')}$. The statement $\text{(1b)}$ must therefore hold, whereas $\text{(1c)}$ also follows. \begin{align} \text{(a)} && p_j &= \frac{N_j}{N}\\ \text{(a')} && \%P_j &= \frac{N_j}{N}\cdot100\%\\ \text{(b)} && 1 &= \sum_j p_j\\ \text{(c)} && N &= \sum_j N_j\\\tag1 \end{align}

   
   
   


2. 
   

   For any system with $k$ particles, that are in thermal equilibrium, the probability to find particle $j$ can be expressed via the Boltzmann statistics $(2)$, where $Z$ is the canonical partition function, the inverse temperature is $\beta = (\mathcal{k}_\mathrm{B}\cdot T)^{-1}$ and the degeneracy coefficient is $g_j$. $$p_j=g_j\cdot \frac1Z\cdot \exp\left\{-\beta\varepsilon_j\right\}\tag2$$

   
   


3. 
   

   I can assume that a state with the energy $\varepsilon_0$ is the ground state, therefore no other state can have lower energy. I can therefore assign the ground state as the baseline and switch to relative energies, as demonstrated in in $(3)$. \begin{align} \text{def.}&& \varepsilon_0 &= 0.0\\ && \varepsilon_j &= \varepsilon_0 + \Delta\varepsilon_j' \\ \therefore && \varepsilon_j &= \Delta\varepsilon_j'\\\tag3 \end{align}

   
   


4. 
   

   As a consequence I can express the canonical partition function in terms of the ground state probability. \begin{align} && p_0 &= g_0\cdot \frac1Z\cdot \exp\left\{-\beta\varepsilon_0\right\}\\ \therefore && p_0 &= g_0\cdot \frac1Z\\ \equiv && Z &= \frac{g_0}{p_0}\\\tag4 \end{align} The canonical partition function is also defined in $(5)$, but I am quite sure this one would not lead me to my desired goal. $$Z = \sum_l^\infty g_l\cdot\exp\left\{-\beta\varepsilon_l\right\}\tag5$$

   
   



5. 
   

   I can express any probability in terms of the ground state probability following from $(2-4)$. \begin{align} && p_j &= \frac{g_j}{g_0}\cdot p_0\cdot \exp\left\{-\beta\Delta\varepsilon_j'\right\}\\\tag6 \end{align}

   
   


6. 
   

   As a last step I will switch to macroscopic numbers, as my quantum chemical calculations give me $\Delta{}G$ values. \begin{align} \text{def.} && \Delta{}G_j &= \mathcal{N}_\mathrm{A}\cdot \Delta\varepsilon_j'\\ && p_j &= \frac{g_j}{g_0}\cdot p_0\cdot \exp\left\{-\frac{\Delta{}G_j}{\mathcal{R}\cdot{}T}\right\}\\\tag7 \end{align}

   
   

• How do I obtain $p_0$ in a way, that I can use it in spreadsheet software like libreOffice calc, so that I satisfy $\text{(1b)}$? If I insert this equation, I would have to solve it iteratively, as I have the variable to solve for on both sides. I could then go ahead and use $Z$ directly. \begin{align} p_0 &= 1 - p_j -\sum_{i\neq0,j} p_i\\ \color{\red}{p_j} &= \frac{g_j}{g_0}\cdot \left(1 - \color{\red}{p_j} -\sum_{i\neq0,j} p_i\right)\cdot \exp\left\{-\frac{\Delta{}G_j}{\mathcal{R}\cdot{}T}\right\}\\ \end{align}



• Am I making correct assumptions along the way?


• Could the following work? Assign an arbitrary probability to $p_0$, to obtain an arbitrary value for $[p_j]^\ddagger$. What would be reasonable? I picked $0.25$ as an example. Can I then normalize to get the total percentages? \begin{align} \text{def.} && p_0 &= [p_0]^\ddagger\ =\ 0.25\\ && [p_j]^\ddagger &= \frac{g_j}{g_0}\cdot \color{\green}{0.25}\cdot \exp\left\{-\frac{\Delta{}G_j}{\mathcal{R}\cdot{}T}\right\}\\ \text{Norm.} && \%P_j &= \frac{[p_j]^\ddagger}{\sum_i [p_i]^\ddagger}\cdot 100\% \end{align} This would maybe break down the iterative procedure to a two step process. Is this procedure really equivalent to $\text{(1a')}$?

---

Closely related: 4 Compounds in Equilibrium with One Another - Determine Their Equilibrium Concentrations