I have lots of difficulties trying to make a phase plot for the motion of the particle trapped in Lennard-Jones potential:


where $q$ is particle's position and $q_\mathrm{min}$ is the point where $V$ has minimum value which is $-\epsilon$.

The Hamiltonian equations for the particle are:



For simplicity, I set all my quantities $\epsilon$, $q_\mathrm{min}$, and $m$ based on an Argon particle and do all my calculations in SI units.

import numpy as np
import matplotlib.pyplot as plt

# adjustable constants
epsilon = 1.65532127e-21
sigma = 3.4e-10
q_min = sigma * 2**(1/6)
m = 6.63352088e-26 # mass of argon in kg
N = 10 # number of iterations
h = 0.1 # time step in seconds

def T(p):
    return p**2 / (2*m)

def V(q):
    return epsilon * (pow(q_min / q, 12) - 2 * pow(q_min / q, 6))

def dq_dt(p):
    return p / m

def dp_dt(q):
    return 12 * epsilon * (pow(q_min / q, 12) - pow(q_min / q, 6)) / q

def plot_phase_space(q_0, p_0):
    t = np.arange(0, N, 1)
    q = np.zeros(N)
    p = np.zeros(N)
    E = np.zeros(N)

    q[0] = q_0
    p[0] = p_0
    E[0] = T(p_0) + V(q_0)
    print(f'T[0]={T(p_0)}, V[0]={V(q_0)}')

    for n in range(0, N - 1):
        k1q, k1p = dq_dt(p[n]), dp_dt(q[n])
        k2q, k2p = dq_dt(p[n] + 0.5 * k1p * h), dp_dt(q[n] + 0.5 * k1q * h)
        k3q, k3p = dq_dt(p[n] + 0.5 * k2p * h), dp_dt(q[n] + 0.5 * k2q * h)
        k4q, k4p = dq_dt(p[n] + k3p * h), dp_dt(q[n] + k3q * h)

        print(f'k1q={k1q}, k1p={k1p}')
        print(f'k2q={k2q}, k2p={k2p}')
        print(f'k3q={k3q}, k3p={k3p}')
        print(f'k4q={k4q}, k4p={k4p}')

        q[n+1] = q[n] + (k1q + 2 * k2q + 2 * k3q + k4q) * h / 6
        p[n+1] = p[n] + (k1p + 2 * k2p + 2 * k3p + k4p) * h / 6

        print(f'T[{n+1}]={T(p[n+1])}, V[{n+1}]={V(q[n+1])}')
        E[n+1] = T(p[n+1]) + V(q[n+1])

        print(f'q[{n+1}]={q[n+1]}; p[{n+1}]={p[n+1]}; E[{n+1}]={E[n+1]}\n')

    plt.plot(q, p)
    plt.title('p vs. q')
    plt.xlabel('q (m)')
    plt.ylabel('p (kg*m/s)')
    plt.plot(t, E)
    plt.title('E vs. t')
    plt.xlabel('t (s)')
    plt.ylabel('E (J)')

plot_phase_space(4e-10, 0)
# plot_phase_space(3e-10, 0)

I want to make a graph that looks like the following

Instead, I got

and here is energy vs. time graph

I am not sure where I went wrong. I am pretty sure that I implemented my Runge-Kutta algorithm correctly. Below is the console output for the first 3 steps:

k1q=0.0, k1p=-9.203520182693993e-12
k2q=-6937130634835.653, k2p=-9.203520182693993e-12
k3q=-6937130634835.653, k3p=1.0160565946844496e-157
k4q=1.5317003037524915e-133, k4p=7.937942145972262e-160
T[1]=1.5961505451923683, V[1]=-1.0454038764160691e-147
q[1]=-462475375655.71027; p[1]=-4.601760091346997e-13; E[1]=1.5961505451923683

k1q=-6937130634835.653, k1p=1.3562718338469806e-158
k2q=-6937130634835.653, k2p=2.698238917184523e-160
k3q=-6937130634835.653, k3p=2.698238917184523e-160
k4q=-6937130634835.653, k4p=2.2221157725748937e-161
T[2]=1.5961505451923683, V[2]=-4.2819742778002206e-150
q[2]=-1156188439139.2756; p[2]=-4.601760091346997e-13; E[2]=1.5961505451923683

k1q=-6937130634835.653, k1p=2.2221157725748937e-161
k2q=-6937130634835.653, k2p=3.5413040479903167e-162
k3q=-6937130634835.653, k3p=3.5413040479903167e-162
k4q=-6937130634835.653, k4p=8.27802632963245e-163
T[3]=1.5961505451923683, V[3]=-2.5522555576564188e-151
q[3]=-1849901502622.841; p[3]=-4.601760091346997e-13; E[3]=1.5961505451923683

You can see a massive decrease in the particle's position $q$. Should I suspect that there is an issue with the units or something else?

  • 1
    $\begingroup$ This is of course only the motion without angular dynamic, for two particles moving directly center-to-center in the bounded region of the potential. One would need 3 particle interaction to go from free to bounded pairs. But only with angular dynamic. $\endgroup$ Mar 14, 2023 at 10:47

1 Answer 1


Your time step is far too large. Just looking at the first stage, you have $\dot{p}\approx10^{-12}$, so for the next stage, you will have something like $\dot{q} \approx h\dot{p}/m \approx 10^{-1}(10^{-12})/(10^{-26})\approx 10^{13}$. To keep the solution from blowing up, this suggests that you will need a time step of $\approx 10^{-13}$ to keep the iteration stable. You can also think of it this way: If you want to keep everything in physical units, then you should also simulate your system on a physical time scale. This system operates a much faster frequency than $h=0.1$ seconds.

I ran your code for 1000 time steps with $h = 10^{-14}$ and got reasonable results. Notice that the energy is not preserved. Integratio schemes known as symplectic integrators are designed to preserve this energy better. Runge-Kutta-Nystrom methods are symplectic integrators that are based on Runge-Kutta methods like the one you have implemented here, but you can also just use smaller time steps to reduce this error in the energy with your existing implementation.

enter image description here

enter image description here

  • 3
    $\begingroup$ RK4 tends to slightly dissipate energy if it is in it's stability range I have noticed. $\endgroup$ Mar 12, 2023 at 21:10

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.