Motion of the particle trapped in potential

Question

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

$$V(q)=\epsilon\left[\left(\frac{q_\mathrm{min}}{q}\right)^{12}-2\left(\frac{q_\mathrm{min}}{q}\right)^{6}\right]$$

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:

$$\dot{q}=\frac{p}{m}$$

$$\dot{p}=\frac{12\epsilon}{q}\left[\left(\frac{q_\mathrm{min}}{q}\right)^{12}-\left(\frac{q_\mathrm{min}}{q}\right)^{6}\right]$$

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.show()
    plt.plot(t, E)
    plt.title('E vs. t')
    plt.xlabel('t (s)')
    plt.ylabel('E (J)')
    plt.show()

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?

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.