I wrote these notes partially to try out the new plotting, math, and citation support in my blog. These notes are for using coordination number as a collective variable for biasing a molecule dynamics simulation.

Using the pair-pair correlation function, coordination number in a statistical ensemble is:

\begin{equation} \left<N(r_0)\right> = \rho \int_0^R\,dr\left[ 1 - u(r - r_0) \right] 4\pi r^2 g(r) \end{equation}

where is the number density and is the heaviside function

In a molecular dynamics simulation, it’s calculated as

\[ \left<N(r_0)\right> = \frac{1}{N_pT} \sum_{t=0}^T \sum_{i < j} 1 - u(r_{ij} - r_0) \]

where is the number of particle pairs and is the number of frames used in the calculation.

To create a soft-derivative, the following mollified heaviside is used

\[ 1 - u(r - r_0) \approx h(s,n,m) = \left\{ \begin{array}{lr} \frac{1 - s^n}{1 - s^m} & s > 1\\ 1 & s < 1\\ \end{array}\right. \]

where

\[ s = \frac{r - r_0}{w} \]

Moments of Coordination Number

The idea of studying moments of coordination number was probably thought of before, but I used it in one of my recent papers and didn’t see it previously in the literature. We can use higher coordination number moments with this equation using the pair-pair correlation function

\[ \left<r^iN(r_0)\right> = \rho \int_0^R\,dr \left[ 1 - u(r - r_0) \right] 4\pi r^{2+i} g(r) \]

In simulation it’s calculated as

\[ \left<r^iN(r_0)\right> = \frac{1}{N_pT} \sum_{t=0}^T \sum_{i < j} r^i(1 - u(r_{ij} - r_0)) \]

With the mollification, there is simply an extra term added.

Choosing Parameters

There are a few considerations for choosing parameters. Assuming you’re biasing a simulation with coordination number, it’s important to balance good properties of the derivative with accuracy in the function. Read the next paragraph for details, but essentially start with , as half the distance from the start of the rise of the first peak in to the first minimum and as the start of the rise of the first peak. Next, increase as much as possible while maintaining accuracy and then increase as much as possible while still keeping good overlap between and .

The overlap between the derivative, and determines how low the forces will be. More overlap, means low forces distributed across many particles. Less overlap means very high forces on a small number of particles. Controlling is done by increasing , which scales the entire line, or by decreasing , which elongates the tails. The second consideration is accuracy. Sharpening the tails is the single most important consideration because long tails can drastically increase . This is from the term seen in Equation 1. Thus, making large increases accuracy.

A plot of along with its derivatives for is shown below:

n
m
i
Update

Implementation Notes

I first saw this mollification from the PLUMED free energy calculation package (Bonomi et al., 2009), but I don’t know where it originally came from.

There is a removable discontinuity of this function:

\[ h(s=1, n, m) = \frac{n}{m} \]

with a derivative of

\[ h(s=1, n, m) = \frac{n(n - m)}{2m} \]

If is truncated at and , the error is approximately

\[ h(r) - h(R) \lt \epsilon \approx \frac{s^n}{s^m} \]

which means may be chosen for a given tolerance as

\[ R = w\epsilon^{1 / (n - m)} + r_0 \]

or

\[ S = \epsilon^{1 / (n - m)} \]

The derivative may be calculated as

\[ \frac{\partial}{\partial x_i} h(s) = \frac{\partial h(s)}{\partial s_i} \frac{\partial s_i}{\partial x_i} = \frac{ns^n(1 - s^m) - ms^m(1 - s^n)}{s(1 - s^m)^2} \frac{\partial s}{\partial x_i} \] \[ \frac{\partial s}{\partial x_i} = \frac{x_i}{wr} \]

where is the th moment of the pairwise distance vector

Maximizing re-use of pieces of the calculation, the following form can be used:

\[ s = (r - r_0) / w, \, a = s^n, \, b = s^m \] \[ h(s) = (1 - a) / (1 - b) \] \[ \frac{\partial h(s) }{\partial x_i} = \frac{na(1 - b) - mb(1 - a)}{wrs(1 - b)^2} x_i \]

The derivative with moments is slightly different:

\[ \frac{\partial}{\partial x_i} r^ih(s) = \left[ir^{i-1} h(s) + r^i\frac{\partial h(s)}{w\partial s}\right] \frac{x_i}{r} \]

Colvars, another free energy calculation package, enforces only even numerator and denominator exponents () (Fiorin et al., 2013). This allows them to avoid taking square roots in the calculations. That is a nice trick, but it would also require only even moments for the coordination numbers with moments. I’ve found that in practice to be too restrictive.

Cited References

  1. Bonomi, M., Branduardi, D., Bussi, G., Camilloni, C., Provasi, D., Raiteri, P., Donadio, D., Marinelli, F., Pietrucci, F., Broglia, R. A., & Parrinello, M. (2009). PLUMED: A portable plugin for free-energy calculations with molecular dynamics. Computer Physics Communications, 180(10), 1961–1972. https://doi.org/10.1016/j.cpc.2009.05.011
  2. Fiorin, G., Klein, M. L., & Hénin, J. (2013). Using collective variables to drive molecular dynamics simulations. Molecular Physics, 111(22-23), 3345–3362. https://doi.org/10.1080/00268976.2013.813594