Subsections

• Potential cut-off in MD

Lennard-Jones potential revisited

Before we proceed developing program for the MD simulations of particles interacting via the Lennard-Jones (or any other) potential we have to derive the expression for the force, and for the characteristic time scale, to determine the integration time step $\delta t$. 

$$U_{LJ}(r_{ij})=-4\epsilon \left(\left( \frac{\sigma}{r_{ij}} \right)^6 - \left(\frac{\sigma}{r_{ij}}\right)^{12} \right)$$ (9)

 

$$\mathbf{F}_{LJ}(\mathbf{r}_{ij})_{\alpha}=\frac{24\epsilon}... ...{\sigma}{r_{ij}}\right)^{12} \right) {\mathbf{r}_{ij}}_{\alpha}$$ (10)

The time scale can be determined by calculating the characteristic oscillations frequency for two particles interacting via the Lennard-Jones potential. Consider two particles of the same mass that are bound together by the Lennard-Jones potential with parameters $\sigma$ and $\epsilon$. Distance between them is close to the equilibrium distance $r_{ij}^0=2^{1/6}\sigma$.As the force is central we can consider a one-dimensional problem

If we assume, that they oscillate around the r_ij⁰ with a very small amplitude, we can expand the potential energy around r_ij⁰ and find for the force:

$$\begin{aligned} F_i=&-\frac{\partial U(r_i-r_j)}{\partial r_... ...2}\right\rvert_{r^0} =\frac{72 \epsilon}{{r^0}^2} \end{aligned}$$ (11)

Equations of motion and characteristic frequency then are:

$$\begin{aligned} m \ddot{x_1} =& - k (x_1 - x_2 -r^0) \ m \d... ...2 = & \frac{k}{m} = \frac{72 \epsilon}{m{r^0}^2} \end{aligned}$$

$\omega$ is the characteristic frequency of the oscillations, which sets a time scale to choose the time step. Usually for the Verlet algorithm $\delta t \sim 0.01 \omega^{-1}$.

Potential cut-off in MD

One technical remark should be made before we proceed to the program code. Same as in Monte-Carlo, the potential cut-off is used in the MD simulations. Two important consequences:

• energy has a discontinuity at the cut-off distance r_c, and
• force is then undefined at r_c

This introduces the error in the numerical integration and energy conservation. The cure is to use spline, that would make the potential and force continuous. This, however precludes the cut-off correction or makes it difficult. Simple way to cope with it is just to ignore the error in the force if the r_c is large enough, and to use the shifted potential energy U = U - U(r_c) to control the energy conservation.