Attempting to do verlet integration

Question

I have tried adapting euler method code, and am using euler to calculate the first value so there are two for me to use in verlet but when i plot the graph i just get two perpendicular straight lines.

this is the code:

for t in t_array:

    if t == 0:
        x0 = x
        x1 = x0
        a = -k * x1 / m
        x2 = x1 + dt * v
        v = v + dt * a
        x_list.append(x2)
        v_list.append(v)
    else:
        x2 = x1
        x1 = x0
        x2 = 2*x1 - x0 + dt**2*a
        v = (1/dt)*(x2 - x1)
        x_list.append(x2)
        v_list.append(v)

and then i plot a graph using matplotlib.

Answer 1

I have implemented a simple Verlet algorithm implementation. You can find the examples on my own GitHub repository.

Here is implementation for Lennard Jones potential:

import numpy as np

def lennard_jones_potential(r, sigma: float, eps: float) -> np.ndarray:
    """
    Lennard-Jones potential
    """
    return 4 * eps * np.power(sigma, 12) / np.power(r, 12) - 4 * eps * np.power(sigma, 6) / np.power(r, 6)


def lennard_jones_force(r, sigma: float, eps: float) -> np.ndarray:
    """
    Lennard-Jones force (Derivative of LJ potential)
    """
    return 24 * eps * np.power(sigma, 6) / np.power(r, 7) - 48 * eps * np.power(sigma, 12) / np.power(r, 13)


def get_kinetic_energy(m, v):
    """
    Get kinetic energy of particle
    """
    return m * (v ** 2) / 2


def get_verlet_next_x(x, v, f, m, dt):
    """
    Get next position of particle
    """
    a = f / m
    return x + v * dt + a * (dt ** 2) / 2


def get_verlet_next_f(r, rx, sigma: float, eps: float):
    """
    Get next LJ interaction force of particle
    """
    return np.sum(lennard_jones_force(r, sigma, eps) * rx / r)


def get_verlet_next_u(rs, sigma: float, eps: float):
    """
    Get next LJ interaction potential of particle
    """
    return np.sum(lennard_jones_potential(rs, sigma, eps)) / 2


def get_verlet_next_v(v, f, f_next, m, dt):
    """
    Get next velocity of particle using Verlet algorithm
    """
    a, a_next = f / m, f_next / m
    return v + ((a + a_next) * dt) / 2


def get_verlet_next(p: np.ndarray, rx: np.ndarray, ry: np.ndarray, dt: float, sigma: float, eps: float):
    """
    Get next state values of particle using Verlet algorithm
    (force, velocity, potential, kinetic energy)
    """
    r = np.sqrt(np.power(rx, 2) + np.power(ry, 2))

    fx_next = get_verlet_next_f(r, rx, sigma, eps)
    fy_next = get_verlet_next_f(r, ry, sigma, eps)

    vx_next = get_verlet_next_v(p[P_VX], p[P_FX], fx_next, p[P_M], dt)
    vy_next = get_verlet_next_v(p[P_VY], p[P_FY], fy_next, p[P_M], dt)
    v_next = np.sqrt(np.power(vx_next, 2) + np.power(vy_next, 2))

    u_next = get_verlet_next_u(r, sigma, eps)
    ek_next = get_kinetic_energy(p[P_M], v_next)

    return fx_next, fy_next, vx_next, vy_next, v_next, u_next, ek_next

# Indexes of particle properties
P_ID, P_M, P_X, P_Y, P_VX, P_VY, P_FX, P_FY = 0, 1, 2, 3, 4, 5, 6, 7
P_X_NEXT, P_Y_NEXT, P_VX_NEXT, P_VY_NEXT, P_FX_NEXT, P_FY_NEXT = 8, 9, 10, 11, 12, 13

# Indexes of result properties
R_ID, R_M, R_X, R_Y, R_VX, R_VY, R_FX, R_FY, R_U, R_EK, R_ET, R_T = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11

def verlet_algo(particles: np.ndarray, simulation_time: float, dt: float, sigma: float, eps: float) -> np.ndarray:
    """
    Verlet algorithm for N particles in 2D space
    """
    time = np.arange(0, simulation_time + dt, dt)

    result = np.zeros(shape=(particles.shape[0], time.shape[0], 12))

    for i in range(particles.shape[0]):
        result[i, :, R_ID] = np.repeat(particles[i, P_ID], result.shape[1])
        result[i, :, R_M] = np.repeat(particles[i, P_M], result.shape[1])

    result[:, 0, R_X] = particles[:, P_X]
    result[:, 0, R_Y] = particles[:, P_Y]
    result[:, 0, R_VX] = particles[:, P_VX]
    result[:, 0, R_VY] = particles[:, P_VY]
    result[:, 0, R_FX] = particles[:, P_FX]
    result[:, 0, R_FY] = particles[:, P_FY]
    result[:, 0, R_U] = 0
    result[:, 0, R_EK] = 0
    result[:, 0, R_ET] = 0
    result[:, :, R_T] = time

    # Calculate the initial forces and energies
    for i in range(particles.shape[0]):
        p = particles[i]
        p_id = p[P_ID].astype(np.int64)
        neighbours = particles[particles[:, P_ID] != p[P_ID]]

        rx = neighbours[:, P_X] - p[P_X]
        ry = neighbours[:, P_Y] - p[P_Y]

        fx, fy, vx, vy, v, u, ek = get_verlet_next(p, rx, ry, dt, sigma, eps)

        particles[i, P_FX] = fx
        particles[i, P_FY] = fy

        result[p_id, 0, R_FX] = fx
        result[p_id, 0, R_FY] = fy
        result[p_id, 0, R_U] = u
        result[p_id, 0, R_EK] = ek
        result[p_id, 0, R_ET] = ek + u

    for curr in range(len(time) - 1):
        next = curr + 1

        # Get following coordinates of all particles:
        particles[:, P_X_NEXT] = get_verlet_next_x(particles[:, P_X], particles[:, P_VX], particles[:, P_FX],
                                                   particles[:, P_M], dt)
        particles[:, P_Y_NEXT] = get_verlet_next_x(particles[:, P_Y], particles[:, P_VY], particles[:, P_FY],
                                                   particles[:, P_M], dt)

        for i in range(particles.shape[0]):
            p = particles[i]
            p_id = p[P_ID].astype(np.int64)
            neighbours = particles[particles[:, P_ID] != p[P_ID]]

            # Get distances between current particle and all other particles:
            rx = neighbours[:, P_X_NEXT] - p[P_X_NEXT]
            ry = neighbours[:, P_Y_NEXT] - p[P_Y_NEXT]

            # Get next force, velocity and energies for current particle using verlet algorithm:
            fx_next, fy_next, vx_next, vy_next, v_next, u_next, ek_next = get_verlet_next(p, rx, ry, dt, sigma, eps)

            # Save next values of current particle:
            p[P_VX_NEXT] = vx_next
            p[P_VY_NEXT] = vy_next
            p[P_FX_NEXT] = fx_next
            p[P_FY_NEXT] = fy_next

            # Save next values of current particle to result array:
            result[p_id, next, R_X] = p[P_X_NEXT]
            result[p_id, next, R_Y] = p[P_Y_NEXT]
            result[p_id, next, R_VX] = vx_next
            result[p_id, next, R_VY] = vy_next
            result[p_id, next, R_FX] = fx_next
            result[p_id, next, R_FY] = fy_next
            result[p_id, next, R_U] = u_next
            result[p_id, next, R_EK] = ek_next
            result[p_id, next, R_ET] = ek_next + u_next

        # Apply next values to all particles:
        particles[:, P_X] = particles[:, P_X_NEXT]
        particles[:, P_Y] = particles[:, P_Y_NEXT]
        particles[:, P_VX] = particles[:, P_VX_NEXT]
        particles[:, P_VY] = particles[:, P_VY_NEXT]
        particles[:, P_FX] = particles[:, P_FX_NEXT]
        particles[:, P_FY] = particles[:, P_FY_NEXT]

    return result

particles = np.array([
    #id,m, x,  y,  vx,vy,fx,fy,ẋ, ẏ, vẋ,vẏ,fẋ,fẏ (ẋ == x_next, ...)
    [0, 2, 11, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 1, 10, 13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [2, 3, 12, 13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
], dtype=np.float64)

result = verlet_algo(particles, simulation_time=50, dt=0.002, sigma=1, eps=0.7, animate=True)

Energies plot: