# Computational Thermo Equations¶

## $$pK_a$$ Equation¶

\begin{align} pK_a\text{ from }K_a:\\ \ce{HA &<=> A- + H+}\\ \ce{K_a}&=\frac{\ce{[A-][H+]}}{\ce{[HA]}}\\ \ce{pK_a}&=\ce{-log(K_a)}\\ pK_a\text{ from }pH:\\ \ce{pK_a}&=\text{pH}\ce{-log}\bigg(\frac{\ce{[A-]}}{\ce{[HA]}}\bigg) \end{align}

## $$K_{eq}$$ From $$pK_a$$ of Reactant Species¶

$pK_{eq}=pK_a^{\text{left acid}}-pK_a^{\text{right acid}}\\ K_{eq}=10^{-pK}$
e.g. Will ammonia and acetonitrile react?
$\ce{NH3 + H2C#N- -> NH2- + H3C#N}\\ 36 \hskip{1.7cm} \ce{->} \hskip{1.5cm} 25$

\begin{align} pK_{eq}&=36-25\\ &=11\\ k_{eq}&=1\times10^{-11}\\ k_{eq}&<<<1 \end{align}

$$k_{eq}$$ is much less than 1, so we have much more reactant than product. No reaction, but ammonium will:

$\ce{NH4+ + H2C#N- -> NH3 + H3C#N}\\ 9.3 \hskip{1.7cm} \ce{->} \hskip{1.5cm}25$

\begin{align} pK_{eq}&=9.3-25\\ &=11\\ k_{eq}&=5\times10^{15}\\ k_{eq}&>>>1 \end{align}

As a shortcut, we can also say that kinetcs favours the acid with the higher $$pK_a$$, and as there’s a big difference between the two, the reactants will be favoured over the reagents

## Calculating $$pK_a$$ From Thermodynamics Calculations¶

Using ‘the proton exchange method’ from A universal approach for continuum solvent $$pK_a$$ calculations: are we there yet? (Ho and Coote 2010), we can use a reference species with a known $$pK_a$$to act as our proton acceptor, this will give us the reference $$pK_a(\text{HRef})$$.

\begin{align} \text{For reaction:}\\ \ce{HA_{(g)} + Ref_{(g)} &-> A-_{(g)} + HRef_{(g)}}\\ \text{We can use the equation:}\\ \Delta G_g&= (\ce{A-_{(g)} + HRef_{(g)}})-(\ce{HA_{(g)} + Ref_{(g)}})\\ \text{And susequent solvation equtions}\\ \Delta G_{solv}(\ce{X})&=\Delta G(\ce{X_{(aq)}})-\Delta G(\ce{X_{(g)}}) \end{align}

From this point, we need to calculate $$\Delta G_{sol}$$ (not to be mixed with $$\Delta G_{solv}$$

$\Delta G_{soln}^*=\Delta G_g^*+\Delta G_{solv}^*(\ce{A-})+\Delta G_{solv}^*(\ce{HRef})-\Delta G_{solv}^*(\ce{HA})-\Delta G_{solv}^*(\ce{Ref})$

Warning

Watch your units. You’ll need your $$\Delta G_{soln}^*$$ in $$J\cdot mol^{-1}$$

And now we can use the reference value to correct our calculated $$pK_a$$

$pK_a=\frac{\Delta G^*_{soln}}{RT\ln(10)}+pK_a(\ce{HRef})$

Another method if using water as a proton acceptor is to use this equation, which accounts for the bulk concentration of water. Change it to $$+log[\ce{H2O}]$$ if using $$\ce{OH-}$$ as your proton acceptor.

$pK_a=\frac{\Delta G^*_{soln}}{RT\ln(10)}-\log[H2O]$

Warning

If using the $$\ce{H2O/H3O+}$$ cycle as your HRef, calculations seem to struggle with $$\Delta G_{solv}(\ce{H3O+})$$, so the experimental value of -0.1756148 can be used instead, along with $$pK_a(\ce{HRef})=-1.7$$. You could also use $$\ce{OH-/H2O}\:\:[pK_a(\ce{HRef})=14]$$ instead,

## ZPVE and Thermodynamics¶

The ZPVE can be calculated from the frequencies as such:

$E_{ZPVE}=\frac{1}{2}\sum\limits_ih\nu_i$

Which can then be used to calculate $$G^\circ$$ from $$E^\circ$$

Where:

• $$q=$$ Partition function
• $$Q_c=$$ Reaction quotient $$\bigg(Q_c=\frac{[C]^c[D]^d}{[A]^a[B]^b} \text{ where: } \ce{aA + bB -> cC + dC}\bigg)$$
\begin{align} \Delta G^\circ&=\Delta H^\circ-T\Delta S^\circ\\ \Delta G^\circ&=E^\circ+E_{ZPVE}-RT\ln(q)\\ \Delta U&=\Delta H - RT\Delta n\\ \\ \Delta U&=E^\circ+E_{ZPVE}\\ \Delta H^\circ&=\Delta U+P\Delta V\\ \\ \Delta_r G&=\Delta_r G^\circ+RT\ln Q \end{align}

## Reaction Constants from Thermodynamics¶

### Transition State Theory¶

• $$k =$$ Rate constant
• $$T =$$ Temperature (K)
• $$\Delta G =\Delta G_f-\Delta G_i$$
• $$R =$$ Gas constant
• $$h =$$ Planck constant
• $$\kappa=$$ Transmission coefficient (see To Account for Tunnelling)
\begin{align} \text{Between two minima:}\\ k&=exp\bigg(\frac{-\Delta G}{RT}\bigg)\\ \text{Between minima and TS (Eyring–Polanyi):}\\ k&=\bigg(\frac{\kappa k_b T}{h}\bigg)exp\bigg({\frac{-\Delta G ^{\circ\ddagger}}{RT}}\bigg)\\ \end{align}

Note

If not accounting for tunnelling, $$\kappa=1$$

### Eyring $$\big(\ln\big(\frac{k}{T}\big)$$ over $$\frac{1}{T}\big)$$¶

$\ln\bigg(\frac{k}{T}\bigg)=\frac{-\Delta H^{\circ\ddagger}}{RT}+\frac{\Delta S^{\circ\ddagger}}{R}+\ln{\bigg(\frac{\kappa k_B}{h}\bigg)}\\$

### Arrhenius $$\big(\ln(k)$$ over $$\frac{1}{T}\big)$$¶

\begin{align} \ln(k)&=\frac{-E_a}{RT}+\ln(A)\\ &\hskip{1cm}or\\ k&=A\exp\bigg(\frac{-E_a}{RT}\bigg)\\\\ Where:\\ E_a&=\Delta H^{\circ\ddagger}+RT\\ A&=\frac{k_BT}{h}\exp{\bigg(\frac{1+\Delta S^{\circ\ddagger}}{R}\bigg)} \end{align}

## To Account for Tunnelling¶

Note

I haven’t managed to make this work in plotting, I think I’m messing up something with how I’m handling $$\text{Im}(\nu^\ddagger)$$.
Edit #1: I have figured this out… it turns out I just need to retake a basic spectroscopy unit
Edit #2: Okay, I’m not sure anymore, since the barrier height seems to not do anything
Edit #3: I made it better! I switched my constants to be in units of $$Eh$$ and made the input to be in $$\kjmol$$ which I convert to $$Eh$$. Since the function is outputting $$\kappa(T)$$, the internal units are irrelevant.
Edit #4: While experimenting with some other things in this plot, I noticed that the discontinuation between the two pieces of the piecewise functions was because I had incorrectly transcribed the function (in one piece I wrote $$\frac{1}{\alpha-\beta}$$ where I should have put a $$\frac{1}{\beta-\alpha}$$, so they didn’t line up properly)

### Skodje-Truhlar¶

Where:

• $$V^\ddagger=$$ Activation energy (this is formally $$V_0=E^\ddagger-E_i$$)
• $$V=\max(V_1,0)$$ (formally $$V_2$$); where $$V_1=\max\big([E_f-E_i],0\big)$$
• $$\text{Im}(\nu^\ddagger)=$$ Magnitude of the imaginary frequency of the transition state
• Note that this needs to be the frequency in $$Hz$$ or $$s^{-1}$$, so you’ll need to convert it from $$\tilde\nu(cm^{-1})$$ as $$\nu(s^{-1})=\frac{c\cdot100}{\tilde\nu(cm^{-1})^{-1}}$$
• Also note that the frequency is temperature dependent, so this can’t be extrapolated beyond the temperature of the calculation.
\begin{align} \kappa(T)&=\begin{matrix} \frac{\beta}{\beta-\alpha}\bigg[\exp\bigg((\beta-\alpha)(\Delta V^\ddagger-V)\bigg)-1\bigg], & \alpha\leq\beta \\ \frac{\beta\pi/\alpha}{\sin(\beta\pi/\alpha)}-\frac{\beta}{\beta-\alpha}\exp\bigg((\beta-\alpha)(\Delta V^\ddagger-V)\bigg), & \alpha\geq\beta\\ \end{matrix}\\ \\ \alpha&=\frac{2\pi}{h\text{Im}(\nu^\ddagger)} \\ \beta&=\frac{1}{k_BT} \end{align}

This gives your transmission coefficient ($$\kappa(T)$$) which can be inserted back into the Eyring–Polanyi equation

Example

A working?? example of this can be played with here on Desmos

## Boltzmann Distributions (Fractional Population Distribution)¶

Note

If you’re dealing with absolute energies, the energy could be relative to each other. e.g. $$\max(w_i-w_n)$$. This will solve issues with numbers overflowing if they get too big.

\begin{align} \text{Calculate the weighting of each species you want to compare:}\\ w_n&=exp\bigg({\frac{−\Delta G_n}{k_B T}}\bigg)\\ \text{Calculate the sum of their weights (partition function):}\\ q &= \sum\limits_i w_i\\ \text{Find the proportion of each species that would exist:}\\ PD&=\frac{w_n}{q}\\ \end{align}