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.

  • 1
    $\begingroup$ +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! $\endgroup$ Jan 21 at 20:34
  • $\begingroup$ 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. $\endgroup$
    – Hayden S
    Jan 22 at 0:18
  • $\begingroup$ Maybe there is some approximation involved in the derivations? $\endgroup$
    – wzkchem5
    Jan 22 at 10:59

1 Answer 1


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}}} $$


$$ \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).

  • 1
    $\begingroup$ Thank you, really impressive that you figured it out! $\endgroup$
    – Hayden S
    Jan 24 at 5:01
  • $\begingroup$ 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$! $\endgroup$
    – Hayden S
    Jan 24 at 7:17

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .