I am trying to write some codes for my research project. I have 3 different arrays which are `x`, `y`, and `z`. These together represent the positions of the atoms that belong to my system. When I start the simulation, the code prints the positions for every atom. I was able to do this with this code for a 6-atom system via this:

```
for i in range(number_of_steps):
    x, z = RunMD(0.1, x, z, v_x, v_z, m)

    print(num_particles)
    print "Ni"
    print "Ni", x[0], y[0], z[0]
    print "Ni", x[1], y[1], z[1]
    print "Ni", x[2], y[2], z[2]
    print "Ni", x[3], y[3], z[3]
    print "Ni", x[4], y[4], z[4]
    print "Ni", x[5], y[5], z[5]

```

However, now I need to study with much more number of atoms. How can I handle this printing problem for 100-atom system? All the code is below. Thanks

```
import numpy as np
import math

'''
units: eV / A / ps
'''

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

    return rho, phi

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

    return x, z

def GetLJForce(r, epsilon, sigma):
    return 48 * epsilon * np.power(
        sigma, 12) / np.power(
        r, 13) - 24 * epsilon * np.power(
        sigma, 6) / np.power(r, 7)

def GetSpringForce(k, radius):
    return -1 * k * radius

def GetAcc(x_positions, z_positions, m):
    xAcc = [0]*len(x_positions)
    zAcc = [0]*len(z_positions)

    for i in range(0, len(x_positions)-1):
        for j in range(i+1, len(x_positions)):
            delta_x = x_positions[j] - x_positions[i]
            delta_z = z_positions[j] - z_positions[i]

            radius, theta = cart_to_polar(delta_x, delta_z)

            if(radius == 0):
                radius = 1

            force_mag = GetLJForce(radius, 0.84, 2.56) + GetSpringForce(0, radius)

            force_x, force_z = polar_to_cart(force_mag, theta)
            xAcc[i] += force_x / m[i]
            zAcc[i] += force_z / m[i]

            force_x, force_z = polar_to_cart(-force_mag, theta)
            xAcc[j] = force_x / m[j]
            zAcc[j] = force_z / m[j]

    return xAcc, zAcc

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

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

def RunMD(dt, x, z, v_x, v_z, m):
    num_particles = len(x)

    a_x, a_z = GetAcc(x, z, m)

    for i in range(num_particles):
        v_x[i] = UpdateVel(v_x[i], a_x[i], dt)
        v_z[i] = UpdateVel(v_z[i], a_z[i], dt)

    for i in range(num_particles):
        x[i] = UpdatePos(x[i], v_x[i], a_x[i], dt)
        z[i] = UpdatePos(z[i], v_z[i], a_z[i], dt)

        if(z[1], z[2], z[3], z[4], z[5] != 0):
            z[1] = 0
            z[2] = 0
            z[3] = 0
            z[4] = 0
            z[5] = 0

    return x, z

num_particles = 6

m = [1] * num_particles

x = [4, 5, 10, 15, 20, 25]
y = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
z = [5, 0.0, 0.0, 0.0, 0.0, 0.0]

v_x = [0] * num_particles
v_z = [0] * num_particles

v_x[0] = 5

number_of_steps = 100

for i in range(number_of_steps):
    x, z = RunMD(0.1, x, z, v_x, v_z, m)

    print(num_particles)
    print "Ni"
    print "Ni", x[0], y[0], z[0]
    print "Ni", x[1], y[1], z[1]
    print "Ni", x[2], y[2], z[2]
    print "Ni", x[3], y[3], z[3]
    print "Ni", x[4], y[4], z[4]
    print "Ni", x[5], y[5], z[5]

```