How does pairwise energy (not force) affect the simulation?

I’m attempting to implement the following pair style by Hsu and Kremer (DOI: 10.1063/1.5089417):

U(r) = 4\varepsilon \Big[(\frac\sigma r)^{12} - (\frac\sigma r)^{6} + \frac 1 4\Big] for r<2^{1/6}\sigma (WCA potential)
U(r) = \alpha \cos\Big[\pi \frac{r^{2}}{2^{1/3}\sigma^2}\Big] for 2^{1/6}\sigma \leq r < 2^{2/3}\sigma

It’s written in such a way that the energy is discontinuous, but the force is exactly zero at the crossover point. I know that the pressure and the Verlet algorithm use the per-atom forces and not the energy, so they shouldn’t be affected by a constant shift in the potential, but the computed potential energy (e.g., for thermodynamic output or tallying) would be. What would other consequences be in using such a discontinuous potential?

In MD simulations, the energy is just a diagnostic. In minimization it is different.

Thanks for the confirmation!

Also, you can easily make any piece-wise defined potential function continuous by shifting the pieces, i.e you can remove the 1/4 term from the WCA potential or add 1.0 to the cosine function. That will also minimize problems with minimization.

FWIW, this would be a perfect example for using pair style python (see examples/python) and then you can use “pair_write” to create a tabulation and then you even have GPU acceleration, if you want. :wink:
So it is very easy to protoype this.

1 Like