Distance between two bodies remains constant during numerical integration process

There are some details with your code, the main reason that it isn't running is because of this line return np.sum(xAcc), np.sum(yAcc), you calculated a matrix of accelerations from the interactions between all the particles where you added the acceleration of one particle and the inverse of that acceleration to the other, since the mass of all particles is the same then the accelerations are the same, then you sum all the elements of the matrix, so all the terms cancel out, and instead of returning the acceleration of each particle you return the sum of the accelerations of ALL the particles, and given that they all have the same mass, is just 0, even if it was not 0 it would have been wrong because you where mixing them all, so all the particles would have the same acceleration and move in the same direction.

• 1

In that same function you have also something that can be improved, you have something like

if condition:
    A_statement
    B_statement
else:
    B_statement

that is the same as

if condition:
    A_statement
B_statement

since B_statement is always going to get executed.

• 2

You don't have to use global on the variables you used it, where you use them they are still in scope, basically use it when you are going to use that variable in other functions, basically everything that has the same number of tabulations or more knows about that variable, and that ends in the first line that is less tabulated than that one, that leads to the next point.

• 3

def UpdatexVel(v_x, a_x, a1_x, dt):
    return v_x + 0.5*(a_x + a1_x)*dt

def UpdateyVel(v_y, a_y, a1_y, dt):
    return v_y + 0.5*(a_y + a1_y)*dt

There is no need to make this two identical functions, the only difference is the name, they make the same thing, the name of the parameters of the function go out of scope when you end the function body, so you could reuse the a single function for both x and y.

Some examples

def fun(a):
   return a
x = a # NameError: name 'a' is not defined

a is not in scope when we use it for x, it only existed in the body of fun, so you could reuse that name for other variable or parameter.

a = 1
def fun(a):
  return 2 + a
print( fun(7) ) # = 9

That happens because a is shadowed by the function parameter of the same name, while

a = 1
def fun():
  return a
print( fun() ) # = 1

because when fun() looks for a it doesn't find it in its scope, so it looks at a higher scope, and finds a variable a that holds the value 1

• 4

    return xPositions, yPositions

sim_xpos = RunMD(0.1, 10, x, y)[0]
sim_ypos = RunMD(0.1, 10, x, y)[1]

This is smaller, but python has tuples so this code that repeats the calculations of RunMD (which you don't want), can be simplified to

    return xPositions, yPositions

sim_xpos, sim_ypos = RunMD(0.1, 10, x, y)

where the pair that is returned by RunMD is assined to the pair (a tuple of two elements for python) sim_xpos and sim_ypos.

If I where to rewrite parts of that code, I would remove numpy to test the algorithm, and then bring numpy to vectorize the operations and make the code much more efficient, it would look something like this

import math

def cart_to_polar(x, y):
    rho = math.sqrt(x**2 + y**2)
    phi = math.atan2(y, x)
    return rho, phi

def polar_to_cart(rho, theta):
    x = rho * math.cos(theta)
    y = rho * math.sin(theta)
    return x, y

def GetLJForce(r, epsilon, sigma):
    return 48 * r * epsilon * ( ( ( sigma / r ) ** 14 ) - 0.5 * ( ( sigma / r ) ** 8 ) ) / ( sigma ** 2 )

def GetAcc(x_positions, y_positions, m):
    xAcc = [0]*len(x_positions)
    yAcc = [0]*len(x_positions)

    for i in range(0, len(x_positions)-1):
        for j in range(i+1, len(x_positions)):

            # calculations are made from the point of view of i, and then flipped for j
            delta_x = x_positions[j] - x_positions[i]
            delta_y = y_positions[j] - y_positions[i]

            radius, theta = cart_to_polar(delta_x, delta_y)
            # in case two particles are at the same place
            # this is some times called a cutoff radius
            if radius==0: radius = 1e-10

            force_mag = GetLJForce(radius, 0.84, 2.56)

            # since the polar coordinates are centered in the i particle the result is readilly usable
            # particle j interaction with particle i
            force_x, force_y = polar_to_cart(force_mag, theta)
            xAcc[i] += force_x / m[i]
            yAcc[i] += force_y / m[i]

            # flip the sing of the force to represent the 
            # particle i interaction with particle j
            force_x, force_y = polar_to_cart(-force_mag, theta)
            xAcc[j] += force_x / m[j]
            yAcc[j] += force_y / m[j]

    return xAcc, yAcc

def update_pos(x, v, a, dt):
    return x + v*dt + 0.5 * a * dt ** 2

def update_vel(v, a, dt):
    return v + a * dt

def runMD(dt, x, y, v_x, v_y, m):

    num_particles = len(x)

    a_x, a_y = GetAcc(x, y, m)

    for i in range(num_particles):
        v_x[i] = update_vel(v_x[i], a_x[i], dt)
        v_y[i] = update_vel(v_y[i], a_y[i], dt)

    for i in range(num_particles):
        x[i] = update_pos(x[i], v_x[i], a_x[i], dt)
        y[i] = update_pos(y[i], v_y[i], a_y[i], dt)

    return x, y

# number of particles in the system
num_particles = 3
# mass of the particles
m = [1] * num_particles
# starting positions
x = [1, 9, 15]
y = [16, 22, 26]
# starting velocities
v_x = [0] * num_particles
v_y = [0] * num_particles
# number of steps of the simulation
number_of_steps = 10

# the simulation
for i in range(number_of_steps):
    x, y = runMD(0.1, x, y, v_x, v_y, m)
    print(x, y)

I don't know if the physics is precisely how it should be for molecular dynamics, when I was in physics I only made simulations of dynamical systems, maybe my approach is too precise or classical for your system, either way you can use it to compare.