# Confusion on molecular dynamics (MD) simulation units leading to absurd acceleration values

I am trying to code up a very simple MD simulation in order to learn more about it. I am using the Leonard-Jones potential, expressed as $$V=4\epsilon (\frac{\sigma}{r}^{12}-\frac{\sigma}{r}^{6})$$ The pairwise force is therefore $$\frac{1}{r}\frac{dV}{dr}=\frac{24\epsilon}{r}((\frac{\sigma}{r})^6-2(\frac{\sigma}{r})^{12})$$

I want to use parameter values as listed in Li's publication: http://li.mit.edu/A/Papers/05/Li05-2.8.pdf

However, when I work through the math, I end up with unreasonable values as shown in the following calculation:

Suppose two argon atoms are separated by 4 angstroms. The ambient temperature is 300K. In Joules, $$\epsilon=119.8k_B T=119.8\cdot 1.3806\cdot 10^{-23} \frac{J}{K} \cdot 300K=4.9619\cdot 10^{-19} J$$.

The pairwise force is then $$4.9619\cdot 10^{-19}J * \frac{24}{4\cdot 10^{-10} m} ((\frac{3.405Å}{4.0Å})^6-2(\frac{3.405Å}{4.0Å})^{12}) = 2.26\frac{J}{m}=2.71\cdot 10^{-9}N$$.

The mass of one a.m.u. in kg is $$1.6605\cdot 10^{-27} kg$$.

$$f=\frac{m}{a}\implies a=\frac{f}{m}=\frac{2.71\cdot 10^{-9} N}{39.948*1.6605\cdot 10^{-27} kg}=4.0854\cdot 10^{16}\frac{m}{s^2}$$

My simulation has units of angstroms, so I usually convert the acceleration to angstroms.

$$a=4.0854 \cdot 10^{26} \frac{Å}{s^2}$$.

This number is absurd. Even with a timestep of 1ps, the particles will fly out of the bounding box in just a few simulation steps.

Where did I go wrong?

EDIT: Thank you to @Samson for providing the correct pairwise force equation, $$\vec{F_{IJ}=\frac{48\epsilon}{\sigma^2}[(\frac{\sigma}{r_{IJ}})^{14}-0.5\cdot(\frac{\sigma}{r_{IJ}})^{8}]}\vec{r_{IJ}}$$. Here, $$\sigma$$ has units of angstroms, $$\epsilon$$ units of Joules, $$r_{IJ}$$ units of angstroms. $$\vec{r_{IJ}}$$ is the vector difference of the pair's positions, i.e. $$\vec{r_I}-\vec{r_J}$$. If you use that equation, and handle the units correctly, you get acceleration on the order of $$10^{25} \frac{Å}{s^2}$$. Regardless of what integrator you use, the total displacement after one timestep will be proportional to the timestep squared. If the timestep is $$10^{-13}$$ seconds or so, the order of magnitude of displacement in each step will be less than one angstrom, which is an agreeable number to work with. In my implementation, I expressed the time-step as multiples of picoseconds (in some cases, a fraction of a picosecond) and did the unit conversions on paper. This avoids using extreme floating-point numbers and leads to better numerical stability.

• Your pairwise force term is wrong $F = -grad\,U$ Oct 21, 2021 at 15:22
• $F_{ij}=48*\epsilon/\sigma^2*[(\sigma/r_{ij})^{14}-1/2*(\sigma/r_{ij})^8]*r{ij}$ Oct 21, 2021 at 15:30
• @Samson Thank you for the help! Is $r_{ij}$ in the last term a unit vector? Oct 21, 2021 at 16:06