Background
In the world of atomistic modeling with classical force fields, one is often given a force field defined by like interactions (e.g. argon-argon interactions). If one is working with a system involving multiple species, (e.g. argon + krypton), mixing rules are usually applied to derive the parameters of the cross-interactions (e.g. argon-krypton interactions).
A common non-bonded pairwise force field class is the exp-6 or Buckingham potential, used to model close-range repulsive and dispersive-attractive van der Waals forces, which can be written in its tight form as $$ U(r) = \underbrace{A e^{-Br}}_{U^{\mathrm{rep}}(r)} \underbrace{-\frac{C}{r^6}}_{U^{\mathrm{att}}(r)},\tag{1}\label{BH_tight} $$ where $A$, $B$, and $C$ are empirical parameters and $r$ is the internuclear separation. The potential is often transformed into the equivalent loose form, $$ U(r)=\frac{\varepsilon}{1-\frac{6}{\gamma}}\left[\frac{6}{\gamma} e^{\gamma\left(1-\frac{r}{\sigma}\right)}-\left(\frac{\sigma}{r}\right)^6\right],\tag{2}\label{BH_loose} $$ where $A=\frac{6 \varepsilon e^\gamma}{\gamma-6}$, $B=\frac{\gamma}{\sigma}$ and $C=\frac{\gamma \varepsilon}{\gamma-6} \sigma^6$ define the transform between the two parameter spaces.
While the loose form is perhaps a little uglier, it has the advantage of the length scale ($\sigma$) and energy scale ($\epsilon$) parameters being separate, defining the interaction well position and depth, respectively (assuming $\gamma >6$). The third unitless shape parameter $\gamma$ defines the repulsive wall steepness and well width.
My question is about a commonly cited set of combining rules for this force field, usually referred to as the Kong-Chakrabarty rules, which take on the rather unwieldy form of \begin{align} \left[\frac{\epsilon_{ij} \gamma_{ij} e^{\gamma_{ij}}}{\left(\gamma_{ij}-6\right) \sigma_{ij}}\right]^{2 \sigma_{ij} / \gamma_{ij}} &= \left[\frac{\epsilon_{ii} \gamma_{ii} e^{\gamma_{ii}}}{\left(\gamma_{ii}-6\right) \sigma_{ii}}\right]^{\sigma_{ii} / \gamma_{ii}}\left[\frac{\epsilon_{jj} \gamma_{jj} e^{\gamma_{jj}}}{\left(\gamma_{jj}-6\right) \sigma_{jj}}\right]^{\sigma_{jj} / \gamma_{jj}}\tag{3}\label{rule_1},\\ \frac{\sigma_{ij}}{\gamma_{ij}} &= \frac{1}{2}\left(\frac{\sigma_{ii}}{\gamma_{ii}}+\frac{\sigma_{jj}}{\gamma_{jj}}\right)\tag{4}\label{rule_2},\\ \frac{\epsilon_{ij} \gamma_{ij} \sigma_{ij}^6}{\left(\gamma_{ij}-6\right)} &= \left[\frac{\epsilon_{ii} \gamma_{ii} \sigma_{ii}^6}{\left(\gamma_{ii}-6\right)} \frac{\epsilon_{jj} \gamma_{jj} \sigma_{jj}^6}{\left(\gamma_{jj}-6\right)}\right]^{\frac{1}{2}}\tag{5}\label{rule_3}. \end{align}
According to the original papers (here, and here), these combining rules are derived from the following ansatz:
\begin{align} U_{ij}^{\mathrm{rep}}(r_{i} + r_{j}) &= \frac{1}{2}\left[U_{ii}^{\mathrm{rep}}\left(2 r_i\right)+U_{jj}^{\text {rep }}\left(2 r_j\right)\right]\tag{6}\label{ansatz_1},\\ \left[\frac{d U_{ii}^{\mathrm{rep}}(R) }{ d R} \right]_{R=2 r_i} &= \left[\frac{d U_{jj}^{\mathrm{rep}}(R) }{ d R}\right]_{R=2 r_j}\tag{7}\label{ansatz_2},\\ U_{ij}^{\mathrm{att}}(R) &= \left[U_{ii}^{\mathrm{att}}(R) U_{jj}^{\mathrm{att}}(R)\right]^{1/2},\tag{8}\label{ansatz_3} \end{align} where $U_{ij}^{\mathrm{rep}}$ and $U_{ij}^{\mathrm{att}}$ refers to the repulsive and attractive components of the interaction between species $i$ and $j$, respectively. In general, the subscripts $ii$ and $jj$ refer to like-interactions with empirically derived parameters while the parameters labeled with $ij$ are cross parameters defined by the mixing rules.
My Question
How does one derive combining rules (\ref{rule_1}) and (\ref{rule_2}) from the ansatz (\ref{ansatz_1}), (\ref{ansatz_2}), and possibly (\ref{ansatz_3})? I have attempted the algebra myself to no avail, even going so far as to try with a computer algebra system, but I feel like I'm missing something here. The original papers offer little guidance in deriving these mixing rules, they are simply stated outright. Note that mixing rule (\ref{rule_3}) follows immediately from ansatz (\ref{ansatz_3}).
You might ask why I care so much as to write up this question. Well, I would like to apply the same ansatz (\ref{ansatz_1}), (\ref{ansatz_2}), and (\ref{ansatz_3}) to a related but different force field (the Wang-Buckingham potential), but first would like to understand how it's done to the original Buckingham form. The Wang-Buckingham potential is essentially a damped version of the Buckingham potential, which eliminates the singularity at $r=0$, and takes the form of $$ U(r)=\frac{2 \epsilon}{1-\frac{3}{\gamma+3}}\left(\frac{\sigma^6}{\sigma^6+r^6}\right)\left[\frac{3}{\gamma+3} \mathrm{e}^{\gamma\left(1-\frac{r}{\sigma}\right)}-1\right].\tag{9} $$ Note the three parameters of the Wang-Buckingham potential ($\sigma$, $\epsilon$, and $\gamma$) are not equivalent to the Buckingham potential parameters.
My Attempts on Mathematica
I have attempted to solve this using Wolfram Mathematica as well as pen and paper. You can see my attempts using Mathematica. Download my Wolfram Mathematica notebook here.
