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^{26} \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^{-14}$ 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.