# Harmonic oscillators with periodic boundary conditions

I am trying to simulate multiple harmonic oscillators in periodic boundary conditions (subsequently visualizing the process in VMD). I have successfully simulated multiple HOs by using the Leapfrog algorithm, but when I try to impose the periodic boundary conditions on the simulations, unexplainable things start to happen. Below is my code - I am new to python, so I know it is definitely not the most efficient, but the way it is written makes to me a lot of intuitive sense and I cannot see, what is wrong.

import numpy as np
import numpy.random as npr

# distance vector between two particles
def distance_vec(pos_vec1, pos_vec2):
return pos_vec1 - pos_vec2

# norm of the above distance vector
def distance(dist_vec):
return np.linalg.norm(dist_vec)

# script for writing the .xyz file for visualization in VMD
def write_xyz(data, filename) -> None:
with open(filename, 'w') as f:
for j in range(data.shape):
print(f'{data.shape}\n', file=f)
for i in range(data.shape):
print(f'p{i}   {data[i, j, 0]}   {data[i, j, 1]}   0', file=f)

############################ Simulation parameters #############################

N = 2                                        # Number of particles in the HO
M = 5                                        # Number of HOs
num_step = 10000                             # Number of timesteps
dt = 0.001                                   # Timestep size (accuracy)
dim = 2                                      # Number of dimensions

######################### Periodic boundary conditions #########################

min_x = -10
max_x = 10
min_y = -10
max_y = 10

############################# Spring properties ################################

bond_length = 5                              # Bond length (equilibrial position)
k = 30                                       # Spring constant

########################## Particle properties #################################

m = np.array([1,1])                          # Particle masses

######################## Randomized initial conditions ##########################

x1 = npr.uniform(-10,10,(M, dim))              # Particle 1, initial position
x2 = npr.uniform(-10,10,(M, dim))              # Particle 2, initial position

v1 = npr.uniform(-10,10,(M, dim))              # Particle 1, initial velocity
v2 = npr.uniform(-10,10,(M, dim))              # Particle 2, initial velocity

################################ Simulation #####################################

particles = np.zeros((N*M, num_step, dim))

for i in range(M):

# initial distance calculation
d_vec = distance_vec(x1[i, :], x2[i, :])        # not normalized
d = distance(d_vec)

for j in range(num_step):

# Leapfrog steps 1 and 2
F1 = -k * d_vec * (1 - (bond_length / d))
F2 = k * d_vec * (1 - (bond_length / d))

# Leapfrog step 3
v1[i, :] = v1[i, :] + (dt / m) * F1
v2[i, :] = v2[i, :] + (dt / m) * F2

# Leapfrog step 4
x1[i, :] = x1[i, :] + v1[i, :] * dt
x2[i, :] = x2[i, :] + v2[i, :] * dt

#Periodic boundaries in x-axis
if x1[i, 0] < min_x:
x1[i, 0] = max_x - (min_x - x1[i, 0])
elif x1[i, 0] > min_x:
x1[i, 0] = min_x + (x1[i, 0] - max_x)

if x2[i, 0] < min_x:
x2[i, 0] = max_x - (min_x - x1[i, 0])
elif x2[i, 0] > min_x:
x2[i, 0] = min_x + (x1[i, 0] - max_x)

#Periodic boundaries in y-axis
if x1[i, 1] < min_y:
x1[i, 1] = max_y - (min_y - x1[i, 1])
elif x1[i, 1] > min_y:
x1[i, 1] = min_y + (x1[i, 1] - max_y)

if x2[i, 1] < min_y:
x2[i, 1] = max_y - (min_y - x1[i, 1])
elif x2[i, 1] > min_y:
x2[i, 1] = min_y + (x1[i, 1] - max_y)

# Saving positions/timestep
particles[2 * i, j, :] = x1[i, :]
particles[2 * i + 1, j, :] = x2[i, :]

# Distance calculation in periodic boundaries
d_vec = distance_vec(x1[i, :], x2[i, :])

if d_vec > ((max_x - min_x)/2):
d_vec = (max_x - min_x) - d_vec
elif d_vec > ((max_y - min_y)/2):
d_vec = (max_y - min_y) - d_vec

d = distance(d_vec)

# Visualizing trajectory
write_xyz(particles, str(M) + "xHO_2D.xyz")


If anyone has an idea about what is going wrong, I would really appreciate it. Thank you for your help!

The formula for dvec should be very similar to the one you use to correct absolute positions. In words, you want to add or subtract $$L$$, (or in general, multiples of $$L$$, but this is not essential if you are always correcting the positions), where $$L$$ is the box length, so as to bring the components of the vector into the range $$\pm \frac{1}{2}L$$. That's not what you are doing: there are mistakes in the tests that you are doing (the if ... elif tests, which mix up what you are doing to the two independent components of dvec), and then the actual algebra performed on the components, which does not match the formula you are applying to x1 and x2. This is the main problem with your approach.