# 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.

• +10. Beautiful post! I plan to take a closer look when I have more time, but a high-quality post sometimes needs to be acknowledged! Jan 21 at 20:34
• Thanks! FYI it's probably easier to attempt the derivation using the tight form of Buckingham potential and then transform back at the end, but I haven't been able to complete the derivation either way. Jan 22 at 0:18
• Maybe there is some approximation involved in the derivations? Jan 22 at 10:59

I will focus on deriving equations (3) and (4) from equations (6) and (7). In brief: we use equation (7) to derive $$r_j$$ as a function of $$r_i$$, and then differentiate (6) in $$r_i$$ to derive the condition connecting the $$B$$ coefficients.

Then, we use (6) itself at $$r_i$$ to derive the last condition connecting the $$A$$ coefficients. This second part is very gory and full of tricks -- the sort that, if this were a Math Olympiad question, would make me think there has to be a much neater solution to the whole thing, if one can only find the right substitution!

We have

$$U^{\mathrm{rep}}_{ii}(\mathrm{r}) \equiv A_{ii}e^{-B_{ii}\mathrm{r}}$$

(using roman typeface for variable $$\mathrm{r}$$) and thus

$$\partial_{\mathrm{r}} U^{\mathrm{rep}}_{ii}(\mathrm{r}) = -B_{ii}A_{ii}e^{-B_{ii}\mathrm{r}} = -B_{ii} U^{\mathrm{rep}}_{ii}(\mathrm{r}).$$

Substituting into equation (7) gives an expression for $$r_j$$ in terms of $$r_i$$:

$$B_{jj} A_{jj} e^{-2 B_{jj}r_j} = B_{ii} A_{ii} e^{-2 B_{ii}r_i} \Rightarrow r_j = b r_i - \frac{1}{2 B_{jj}} \ln ab,$$

where $$a \equiv A_{ii}/A_{jj}$$ and $$b \equiv B_{ii}/B_{jj}$$.

Now take total derivatives of equation (6) in $$r_i$$:

\begin{align} -(1+b)B_{ij}U^{\mathrm{rep}}_{ij}(r_i + r_j) &\stackrel{\mathrm{d}/\mathrm{d}r_{i}}= -[B_{ii}U_{ii}^{\mathrm{rep}}(2 r_i) + b B_{jj}U_{jj}^{\mathrm{rep}}(2r_j)] \\ &\stackrel{\mathrm{eq}(6)}= -2 B_{ii} U_{ij}^{\mathrm{rep}}(r_i + r_j) \\ \Rightarrow \frac{1}{B_{ij}} = \frac{1}{2} \frac{1+b}{B_{ii}} &= \frac{1}{2} \left(\frac{1}{B_{ii}} + \frac{1}{B_{jj}} \right) \end{align}

which is equation (4).

It remains to tackle equation (6) in full gory detail, substituting our expression for $$r_j$$. The left-hand side expression is:

$$A_{ij} \exp \left(- B_{ij} \left[r_i + b r_i - \frac{1}{2 B_{jj}} \ln ab \right] \right).$$

From earlier, recall that $$(1+b)B_{ij} = 2 B_{ii}$$, and with some work we also have $$B_{ij}/(2B_{jj}) = b/(b+1)$$, so that the above expression becomes

$$(ab)^{\frac{b}{b+1}}A_{ij} e^{-2B_{ii} r_i}.$$

Now the right-hand side expression is:

$$\frac{1}{2} \left[ A_{ii}e^{-2 B_{ii} r_i} + A_{jj} e^{-2 B_{jj} r_j} \right] = \frac{1}{2} \left[ A_{ii}e^{-2 B_{ii} r_i} + A_{jj} ab \,e^{-2 B_{ii} r_i} \right].$$

where we used the definition of $$r_j$$ to simplify. Let's non-dimensionalize our expressions by defining $$\alpha_i \equiv A_{ii}/A_{ij}$$ and $$\alpha_j \equiv A_{jj}/A_{ij}$$, and cancel out factors of $$e^{-2 B_{ii} r_i}$$, to give our simplified equation (6):

$$(ab)^{\frac{b}{b+1}} = \frac{1}{2} \left( \alpha_i + \alpha_j ab \right)$$

Here come two tricks: first, multiply through by $$(ab)^{-1/2}$$ (motivated by the need to symmetrize between $$i$$ and $$j$$, after unceremoniously doing everything in $$i$$ up to now). Then note that

$$\frac{\alpha_i}{\sqrt{a}} = \alpha_j \sqrt{a} = \sqrt{\frac{A_{ii}A_{jj}}{A_{ij}A_{ij}}}$$

and

$$\sqrt{\frac{B_{ii}B_{jj}}{B_{ij}B_{ij}}} = \sqrt{B_{ii}B_{jj}} \times \frac{1}{2} \left(\frac{1}{B_{ii}} + \frac{1}{B_{jj}} \right) = \frac{1}{2} \left(\sqrt{b} + \frac{1}{\sqrt{b}} \right)$$

to get

$$(ab)^{\frac{b}{b+1} - \frac{1}{2}} = \left(\frac{A_{ii}A_{jj}}{A^2_{ij}} \frac{B_{ii}B_{jj}}{B^2_{ij}}\right)^{1/2}.$$

Remembering that

$$ab = (A_{ii}B_{ii})^1 \times (A_{jj}B_{jj})^{-1}$$

grouping powers and rewriting finally gives us

$$A_{ij} B_{ij} = (A_{ii}B_{ii})^{\frac{B_{ij}}{2B_{ii}}} (A_{jj}B_{jj})^{\frac{B_{ij}}{2B_{jj}}}$$

which is equation (3).

• Thank you, really impressive that you figured it out! Jan 24 at 5:01
• I have identified where I failed in my attempts. I did not properly differentiate the function proportional to $r_j$. I forgot rj can be expressed as a function of $r_i$! Jan 24 at 7:17