I'm trying to implement a simple molecular dynamic simulation in Python for N atoms with periodic boundary conditions while allowing coordinate wrapping to visualize the atoms interacting through boundaries. Here is my code:

import numpy as np
import matplotlib.pyplot as plt

def LJ_VF(r):
    if r == 0:
        return 0, 0  
    V = 4 * epsilon * ( (sigma/r)**(12) - (sigma/r)**6 )
    F = 24 * epsilon * ( 2 * ((sigma**12)/(r**(13))) - ( (sigma**6)/(r**7) )) 
    return V , F

def position_verlet(x, v, f_old, dt):
    return x + v * dt + 0.5 * f_old * dt**2

def velocity_verlet(v, f_old, f_new, dt):
    return v + 0.5 * (f_old + f_new) * dt

def PBC(pos, uc_size):
    return pos - uc_size * np.floor(pos / uc_size)

def MID(r1, r2, uc_size):
    d = r2 - r1
    return d - ( uc_size * np.round(d / uc_size) )

def LJ_N(positions):
    N = len(positions)
    total_F = np.zeros(N)  
    total_V = 0.0  
    for i in range(N):
        for j in range(i + 1, N): 
            r = MID(positions[i], positions[j], L)   
            V, F = LJ_VF(r)  
            total_V += V  
            total_F[i] += -F
            total_F[j] += F 
    return total_V, total_F


def simulate_N_atoms(N_atoms, dt):

    positions = np.linspace(0, L, N_atoms, endpoint=False)
    velocities = np.zeros(N_atoms)

    V_tot, forces = LJ_N(positions)

    pos_list, vel_list = [positions.copy()], [velocities.copy()]
    V_list, K_list, H_list = [V_tot], [0.5 * np.sum(velocities**2)], [V_tot + 0.5 * np.sum(velocities**2)]

    for _ in range(N - 1):
        new_positions = np.zeros(N_atoms)
        new_velocities = np.zeros(N_atoms)

        new_positions = positions + velocities * dt + 0.5 * forces * dt**2

        new_V_tot, new_forces = LJ_N(new_positions)

        new_velocities = velocities + 0.5 * (forces + new_forces) * dt

        positions, velocities, forces = new_positions, new_velocities, new_forces
        pos_list.append(positions.copy())
        vel_list.append(velocities.copy())
        V_list.append(new_V_tot)
        K_list.append(0.5 * np.sum(velocities**2))
        H_list.append(new_V_tot + K_list[-1])

    return np.array(pos_list), np.array(vel_list), np.array(V_list), np.array(K_list), np.array(H_list)

#Constants

epsilon = 0.0103       
sigma = 3.4  
m = 1.0
t0 = 0.0
v0 = 0.0
dt = 0.01 
N = 10
L = 10.0
N_atoms = 5


#Plots

positions, velocities, V_list, K_list, H_list = simulate_N_atoms(N_atoms, dt)
time = np.arange(N) * dt

for i in range(N_atoms):
    plt.plot(time, positions[:, i], label=f"Atom {i+1}")
plt.axhline(y=L, color='black', linestyle="--", label="Box Boundary")
plt.axhline(y=0, color='black', linestyle="--")
plt.xlabel("Time (t)")
plt.ylabel("x Position (Å)")
plt.title("Particle Positions Over Time with PBC")
plt.legend()
plt.grid(True)
plt.show()

I'm trying to do a Bouncing test by plotting the atoms positions as a function of time to make sure the PBC is working, but I'm somehow getting an incorrect plot which shows that they are not bouncing and hence the energy is not conserved!

I have went through the code multiple times to see where I went wrong but it was a failed mission!:( Any advice?