# Motion of the particle trapped in potential

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

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

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?

• 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. Mar 14, 2023 at 10:47

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.