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: enter image description here


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.


2 Answers 2


Maybe your numerical solver is instable.

  • $\begingroup$ Thank you for your suggestion However, it does not truly answer my question. I am looking for someone to correct my units calculation or physical understanding. Do you see any errors in my units? Indeed, I am using a position-verlet integrator to reduce long-term error. However, the most important part of my question is really the order of magnitude of a. If it was between 10^(12) and 10^(15), I wouldn't really be worried. But 10^(26) is clearly too large. The smallest time step I found in literature was 10^(-15) seconds, which will at best balance out a values smaller than 10^(16). $\endgroup$ Commented Oct 21, 2021 at 14:38
  • $\begingroup$ My guess is I am missing an angstrom conversion in the force calculation step, but I reviewed the work several times and could not find anything. Essentially, I was expecting the acceleration to be no more than 10^6 m/s^2 = 10^16 angstrom/s^2. $\endgroup$ Commented Oct 21, 2021 at 14:46
  • $\begingroup$ Your pairwise force term is wrong $F = -grad\,U$ $\endgroup$
    – Samson
    Commented Oct 21, 2021 at 15:22
  • $\begingroup$ $F_{ij}=48*\epsilon/\sigma^2*[(\sigma/r_{ij})^{14}-1/2*(\sigma/r_{ij})^8]*r{ij}$ $\endgroup$
    – Samson
    Commented Oct 21, 2021 at 15:30
  • $\begingroup$ @Samson Thank you for the help! Is $r_{ij}$ in the last term a unit vector? $\endgroup$ Commented Oct 21, 2021 at 16:06

I am not familiar with the exact molecular dynamics model you are using, but I have done some simulations with Newton acceleration on mass or Coulomb acceleration on charges. It can be striking what is embedded in those simple accelerations when time-evolved. If the model just simulates Newton's a = G/r2 attraction, we can get Kepler's laws from time evolving its contents. From Coulomb's a = K*q1q2/r2 we get Gauss's law and some others. So, one way is you can perhaps change it to a simpler equation at first and see the results.

  • $\begingroup$ That is a good idea. I will try to implement it as a fallback after trying Samson's answer. You get an upvote from me! $\endgroup$ Commented Oct 21, 2021 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.