Harmonic oscillators with periodic boundary conditions

Question

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[1]):
            print(f'{data.shape[0]}\n', file=f)
            for i in range(data.shape[0]):
                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[0]) * F1
        v2[i, :] = v2[i, :] + (dt / m[1]) * 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[0] > ((max_x - min_x)/2):
            d_vec[0] = (max_x - min_x) - d_vec[0]                 
        elif d_vec[1] > ((max_y - min_y)/2):
            d_vec[1] = (max_y - min_y) - d_vec[1]

        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!

Answer 1

You seem to be applying periodic boundary corrections wrongly to your dvec, as well as to your vector x2. The x2 problem is just that you've copied part of the formula for x1, without changing 1 to 2. Really, checking your code is not what this site is for, but I couldn't just let this slip by.

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.