I am doing a 2D MD simulations of charge carriers in graphene using the Leapfrog algorithm. I am trying to prove that, in some specific cases (when distance between particles is small), the method is unstable.

The Hamiltonian is given as $H = |\vec{p_{1}}| + |\vec{p_2}| - \frac{\alpha}{r_{12}}$. With $p$ the momentums, $\alpha$ some numerical constant and $r_{12}$ some constants.

The leapfrog iteration would be:

$$ \vec{p_i} \to \vec{p_i} - 0.5 \Delta t \frac{\alpha}{r_{ij}^2} \hat{r_{ij}}\\ \vec{x_i} \to \vec{x_i} + \Delta t \hat{p_i} \\ \vec{p_i} \to \vec{p_i} - 0.5 \Delta t \frac{\alpha}{r_{ij}^2} \hat{r_{ij}} $$

where the vectors with hat's are the direction vectors. I expect the stability to break down when $r_{ij}$ gets small. Yet, I tried using analysis from this (page 7 onward) source and it didn't really help.

Also, a different way to derive leapfrog is to define coordinate shift operators $L_p$ and $L_r$ that move $p \to p +\Delta p$ and $x \to x + \Delta x$. Then the differential equation can be written as: $$ f(p(t), r(t)) = \exp\left(t (L_p + L_r) \right) f(p(0), r(0)). $$ Using Trotter's identity we can write: $$ \exp\left( \Delta t (A + B) \right) \approx \exp\left( 0.5\Delta t A \right) \exp\left( \Delta t B \right) \exp\left(0.5 \Delta t A \right) $$ and this is exactly the Leapfrog algorithm. I tried figuring out if things breakdown in 1D case of two particles going straight at each other. To do this I worked out the commutator term $[L_p,L_r]$ that comes up in Zassenhaus formula. This also didn't help.

Anyone has a fresh look or has some derivations that go beyond basic 1D cases where things go nicely? Any help is appreciated.

  • $\begingroup$ Are you sure about your Hamiltonian? There are squares and a factor 1/2 missing in the kinetic part. $\endgroup$ Jan 14, 2019 at 17:52
  • $\begingroup$ Yes. If velocities are relativistic ($v \to c$) then kinetic energy no longer equals $\frac{p^2}{2m}$. This is the case of $m = 0$. $\endgroup$ Jan 14, 2019 at 18:41
  • 1
    $\begingroup$ @PiotrBenedysiuk This question should be asked on the physics-site. $\endgroup$
    – Peter
    Jan 21, 2019 at 8:35
  • 1
    $\begingroup$ Good point. I mirrored it to over there $\endgroup$ Jan 21, 2019 at 10:02
  • 1
    $\begingroup$ The exact claim, as can be found in the works of Hairer et al., is that the Verlet method in any of its implementation preserves to $O(Δt^4)$ an $O(Δt^2)$ modification of the Hamiltonian. This modification includes derivatives of the Hamiltonian, so that close to singular points of the Hamiltonian the $O(Δt^2)$ modification need not be small at all. // Of course this does not help if the Hamiltonian used is not even a first integral of the system. $\endgroup$ Jan 23, 2019 at 13:58


You must log in to answer this question.

Browse other questions tagged .