Magnitude of Lennard-Jones force evaluated with molecular dynamics simulation

Question

The intent is to evaluate the error source of a molecular dynamics (MD) simulation which is only based on the Lennard-Jones interaction between noble gases (helium in this case). What is the magnitude of the force?

For reference: $\varepsilon = \pu{1.49E-22 J}$ and $\sigma = \pu{2.967E-10 m}.$ The equation used to calculate the force is:

$$F = 4\varepsilon\left(\frac{6\sigma^6}{r_\mathrm{abs}^7} - \frac{12\sigma^{12}}{r_\mathrm{abs}^{13}}\right) \sim \pu{E-13 N}$$

Combined with the mass of helium this leads to an acceleration which seems to be too large. Would the above make sense?

Our step size is roughly $\pu{5E-3 s}$ and the velocity has a magnitude of $\pu{E13 m s^-1}.$ This is our problem at the moment because those values are obviously too large. We are using the Verlet algorithm.

Answer 1

The following plots the force for reasonable distances

Intuitively, if you are implementing the verlet algorithm then the step size will be $\delta x \propto a\delta t^2$ (ignoring velocities). If you want a step size of order ~$f \times \sigma$ (where f is a scaling factor) then you'd like to satisfy a relation such as $\delta t < \sqrt{\frac{f\sigma}{a}}$ (if ignoring velocities). The timestep cutoff $\delta t$ computed according to this equation with a force of $\pu{10^-12 N}$ and f=0.1 is 0.3 psec. But 0.3 psec is still a relatively long timestep. Clearly the problem is not with acceleration (provided the timestep is shorter than the cutoff just calculated), but maybe there is one with the velocities. You can proceed in the same way but solving the complete quadratic verlet equation for the timestep in terms of the step size, and inserting a value for the rms speed based on the desired T and step size, to arrive at a timestep cutoff. Using this reasoning I came up with a cutoff $\delta t < \pu{ 15 fsec}$

One problem with a large step size (and discretization of the trajectory in general) is that the energy is not conserved and can easily blow up (along with the temperature): if you take a step of order $\sigma$ from the potential well minimum toward the maximum in the potential then the atoms will suddenly find themselves in the strongly repulsive high energy region.

You may want to check Ref. 1 which discusses this to some extent, and provides this advice:

One wants to find the largest time step that will maintain the conservation of energy.

Reference

1. Sangrak Kim. Issues on the Choice of a Proper Time Step in Molecular Dynamics. Physics Procedia 53 (2014) 60–62. doi: 10.1016/j.phpro.2014.06.027