计算科学 - 具有周期性边界条件的谐波振荡器 - 吾爱随笔录

我正在尝试在周期性边界条件下模拟多个谐波振荡器（随后在 VMD 中可视化该过程）。我已经使用 Leapfrog 算法成功地模拟了多个 HO，但是当我尝试对模拟施加周期性边界条件时，开始发生无法解释的事情。下面是我的代码——我是 python 新手，所以我知道它绝对不是最有效的，但是它的编写方式让我很直观，我看不出有什么问题。

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")

如果有人知道出了什么问题，我将不胜感激。谢谢您的帮助！