Python 分子气体模拟，产生恒定值

Question

我正在编写一个脚本，该脚本应该使用伦纳德-琼斯势来模拟气体分子的相互作用。 问题是，尽管代码有效，但它返回的结果表明粒子几乎没有移动。但实际情况并非如此！

所以代码流程可能有问题。但我无法找到问题所在。因此，让我首先简要描述一下预期的代码流程。

1. 粒子初始化，使用函数

   class Czastka

   从

   initialization()

   初始化16个粒子。

2. simulate()

   函数会迭代

   num_step = 1000

   遍历粒子列表，并计算由于与其他粒子相互作用而作用在每个粒子上的力的总和。为此，使用

   lennard_jones_forces(x1, y1, x2, y2)

   函数。所应用的条件是，如果粒子相距超过

   2.5*sigma

   ，则它们不会感受到任何吸引力。

3. 使用函数

   leapforg_single_step(particle, force_x_comp, force_y_comp):

   计算新的位置和速度。

4. 然后使用

   p1.update_position(p1_x_next, p1_y_next, box_size)

   将其附加到粒子位置列表中。 这里还检查周期性边界条件。其背后的想法是，粒子是否会离开模拟/包含气体的框架/盒子。从另一面看来，这是通过模除来实现的

   self.x = new_x % box_size

   。

正如我所说，我不明白为什么代码会产生几乎恒定的值。也许这与周期性边界条件有关。

import numpy as np
import random


#Inicjalizacja paramatrów
particle_num, box_size, eps, sigma = 16, 8.0, 1.0, 1.0
dt, temp, k_B, m, M = 0.0001, 2.5, 1, 1, 1
radius = sigma/2



class Czastka:
    #Inicjalizacja danych cząstki
    def __init__(self, radius, x_pos, y_pos, x_vel, y_vel):
        #Składowe wartości wektorowych
        self.x, self.y = x_pos, y_pos
        self.vx, self.vy = x_vel, y_vel

        """
        NOTE: The goal of the simulation in the end is to get only these 4 lists
        """
        #Positions lists
        self.x_positions_lst = []
        self.y_positions_lst = []

        #Velocity lists
        self.x_velocity_lst = []
        self.y_velocity_lst = []

    def append_postions(self, x_pos, y_pos):
        self.x_positions_lst.append(x_pos)
        self.y_positions_lst.append(y_pos)

    def update_velocities(self, x_vel, y_vel):
        self.x_velocity_lst.append(x_vel)
        self.y_velocity_lst.append(y_vel)

    # Stosujemy periodyczne warunki brzegowe poprzez liczenie modułowe
    def update_position(self, new_x, new_y, box_size):
        self.x = new_x % box_size
        self.y = new_y % box_size
        self.append_postions(self.x, self.y)

# Inicjalizacja cząstek
initial_x, initial_y = [1, 3, 5, 7], [1, 3, 5, 7]
initial_vx, initial_vy = [1, 1, 2, 0.5], [1, 0.5, 0.25, 4]


#Inicjalizacja cząstek wraz z przynaleznymi wartościami początkowymi
particle_lst = []

def initialization():
    n, m = 0, 0
    global particle_lst

    # Inicjalizacja cząstek
    initial_x, initial_y = [1, 3, 5, 7], [1, 3, 5, 7]
    initial_vx, initial_vy = [1, 1, 2, 0.5], [1, 0.5, 0.25, 4]

    # Generuje losowe prędkości
    initial_vx_random = [random.uniform(-2, 2) for _ in range(16)]
    initial_vy_random = [random.uniform(-2, 2) for _ in range(16)]

    # Inicjalizuje cząstki
    for i in range(0, 4):
        for j in range(0, 4):
            x_pos = initial_x[i]
            y_pos = initial_y[j]
            vx = initial_vx_random[n]
            vy = initial_vy_random[n]

            #print('x-pos: {}, y-pos: {}, vx: {}, vy: {}'.format(x_pos, y_pos, vx, vy))

            Czastka(radius, x_pos, y_pos, vx, vy)
            particle_lst.append(Czastka(radius, x_pos, y_pos, vx, vy))

            n += 1

    # Dodaję początkowe położenia i prędkości cząstek do ich własnych list
    for i, p1 in enumerate(particle_lst):
        m = i % 4
        for k in range(0, 4):
            continue
        p1.append_postions(initial_x[m], initial_y[k])
        p1.update_velocities(initial_vx_random[i], initial_vy_random[i])

initialization()


def simulation_results(target_lst):
    for i, particle in enumerate(target_lst):
        print('Particle number: {}. List lenght: {} and x-positions list: {} and vx-velocity list: {}'.format(i, len(particle.x_positions_lst), particle.x_positions_lst, particle.x_velocity_lst))

#simulation_results(particle_lst)



#Siła jaka działa na skutek odziaływań w potencjale Lennard-Jonesa
def lennard_jones_forces(x1, y1, x2, y2):
    global sigma, eps

    #Obliczenia pomocnicze
    r12 = [x2-x1, y2-y1]
    dist12 = (r12[0]**2 + r12[1]**2)**0.5
    if dist12 == 0:
        return 0, 0


    # Siła jaką cząstka 2 wywiera na cząstekę 1
    Force21 = -(48*eps)/(sigma**2) * ((sigma/dist12)**14 - 0.5*(sigma/dist12)**8)
    Force21_x = Force21*r12[0]
    Force21_y = Force21*r12[1]


    #W tym momecie nasza funkcja jest gotowa ALE musimy sprawdzić czy zachodzi warunek odcięcia
    #Żeby zwiększyć wydajnośc obliczniową powinno się dokonać tego sprwdzenia, przed obliczniem sił
    if dist12 > (2.5*sigma):
        #print('Cut off')
        Force21_x, Force21_y = 0, 0
        return Force21_x, Force21_y
    else:
        #print('Normal operation')
        return Force21_x, Force21_y


# Obliczanie poprzez wykorzystanie algorytmu żabki nowych współrzędnych
def leapforg_single_step(particle, force_x_comp, force_y_comp):
    global dt, m

    #Obliczanie pół-krokowych prędkości
    vx_half = particle.vx + 0.5 * (force_x_comp / m) * dt
    vy_half = particle.vy + 0.5 * (force_y_comp / m) * dt

    #Obliczanie pół-krokowych położeń
    x_next = particle.x + vx_half * dt
    y_next = particle.y + vy_half * dt

    #Obliczanie nowych prędkości
    vx_next = vx_half + 0.5 * (force_x_comp / m) * dt
    vy_next = vy_half + 0.5 * (force_y_comp / m) * dt

    return x_next, y_next, vx_next, vy_next



def simulate():
    global box_size
    num_steps = 1000

    #Cała symulacja składa się z num_steps kroków
    for step in range(num_steps):
        #Obliczmy sumę sił działającą na każdą cząstkę
        for i, p1 in enumerate(particle_lst):
            force_x, force_y = 0, 0
            for j, p2 in enumerate(particle_lst):
                if i != j:
                    # Obliczam sume sił działająca na cząsteki p1 (bez interakcji z samym sobą)
                    f_x, f_y = lennard_jones_forces(p1.x, p1.y, p2.x, p2.y)
                    force_x += f_x
                    force_y += f_y

            p1_x_next, p1_y_next, p1_vx_next, p1_vy_next = leapforg_single_step(p1, force_x, force_y)
            p1.update_position(p1_x_next, p1_y_next, box_size)
            p1.update_velocities(p1_vx_next, p1_vy_next)

simulate()
simulation_results(particle_lst)


def animation():
    pass

任何帮助将不胜感激！干杯！

Answer 1

零距离是非物理的。您有一个可能不正确的初始迭代，仅从 4 个初始位置值生成 16 个粒子。当我解决这个问题时，即使我通过引发异常来检测相同的噪音，也不会出现 NaN。这里还有其他需要改进的地方（你根本没有正确使用 Numpy）。

import numpy as np

# Inicjalizacja paramatrów
particle_num, box_size, eps, sigma = 16, 8.0, 1.0, 1.0
dt, temp, k_B, m, M = 0.0001, 2.5, 1, 1, 1
radius = sigma/2


class Particle:
    # Inicjalizacja danych cząstki
    def __init__(self, radius: float, x_pos: float, y_pos: float, x_vel: float, y_vel: float) -> None:
        # Składowe wartości wektorowych
        self.x, self.y = x_pos, y_pos
        self.vx, self.vy = x_vel, y_vel

        """
        NOTE: The goal of the simulation in the end is to get only these 4 lists
        """
        # Positions lists
        self.x_positions_lst: list[float] = []
        self.y_positions_lst: list[float] = []

        # Velocity lists
        self.x_velocity_lst: list[float] = []
        self.y_velocity_lst: list[float] = []

    def append_postions(self, x_pos: float, y_pos: float) -> None:
        self.x_positions_lst.append(x_pos)
        self.y_positions_lst.append(y_pos)

    def update_velocities(self, x_vel: float, y_vel: float) -> None:
        self.x_velocity_lst.append(x_vel)
        self.y_velocity_lst.append(y_vel)

    # Stosujemy periodyczne warunki brzegowe poprzez liczenie modułowe
    def update_position(self, new_x: float, new_y: float, box_size: float) -> None:
        self.x = new_x % box_size
        self.y = new_y % box_size
        self.append_postions(self.x, self.y)

    # Obliczanie poprzez wykorzystanie algorytmu żabki nowych współrzędnych
    def leapfrog_single_step(self, force_x_comp: float, force_y_comp: float) -> tuple[
        float,  # x next
        float,  # y next
        float,  # vx next
        float,  # vy next
    ]:
        # Obliczanie pół-krokowych prędkości
        vx_half = self.vx + 0.5 * (force_x_comp / m) * dt
        vy_half = self.vy + 0.5 * (force_y_comp / m) * dt

        # Obliczanie pół-krokowych położeń
        x_next = self.x + vx_half * dt
        y_next = self.y + vy_half * dt

        # Obliczanie nowych prędkości
        vx_next = vx_half + 0.5 * (force_x_comp / m) * dt
        vy_next = vy_half + 0.5 * (force_y_comp / m) * dt

        return x_next, y_next, vx_next, vy_next


# Siła jaka działa na skutek odziaływań w potencjale Lennard-Jonesa
def lennard_jones_forces(
    x1: float, y1: float, x2: float, y2: float,
) -> tuple[float, float]:
    # Obliczenia pomocnicze
    r12 = [x2-x1, y2-y1]
    dist12 = (r12[0]**2 + r12[1]**2)**0.5
    # print(r12[0], ' ', r12[1])
    if dist12 == 0:
        raise ValueError('A zero distance is non-physical')

    # Calculate Lennard-Jones force
    sigma_over_dist12 = sigma / dist12
    sigma_over_dist12_14 = sigma_over_dist12 ** 14
    sigma_over_dist12_8 = sigma_over_dist12 ** 8

    Force21 = -(48 * eps / (sigma ** 2)) * (sigma_over_dist12_14 - 0.5 * sigma_over_dist12_8)
    Force21_x = Force21 * r12[0]
    Force21_y = Force21 * r12[1]

    # W tym momecie nasza funkcja jest gotowa ALE musimy sprawdzić czy zachodzi warunek odcięcia
    # Żeby zwiększyć wydajnośc obliczniową powinno się dokonać tego sprwdzenia, przed obliczniem sił
    if np.isnan(Force21_x) or np.isnan(Force21_y):
        raise ValueError(f"Nan detected in force calculation: {Force21_x} {Force21_y}")

    if dist12 > (2.5*sigma):
        # print('Cut off')
        Force21_x, Force21_y = 0, 0
        return Force21_x, Force21_y
    else:
        # print('Normal operation')
        return Force21_x, Force21_y


def simulate(particle_lst: list[Particle]) -> None:
    num_steps = 1000

    # Cała symulacja składa się z num_steps kroków
    for step in range(num_steps):
        # Obliczmy sumę sił działającą na każdą cząstkę
        for i, p1 in enumerate(particle_lst):
            force_x, force_y = 0, 0
            for j, p2 in enumerate(particle_lst):
                if i != j:
                    # Obliczam sume sił działająca na cząsteki p1 (bez interakcji z samym sobą)
                    f_x, f_y = lennard_jones_forces(p1.x, p1.y, p2.x, p2.y)
                    force_x += f_x
                    force_y += f_y

            p1_x_next, p1_y_next, p1_vx_next, p1_vy_next = p1.leapfrog_single_step(force_x, force_y)
            p1.update_position(p1_x_next, p1_y_next, box_size)
            p1.update_velocities(p1_vx_next, p1_vy_next)


def simulation_results(particle_lst: list[Particle]) -> None:
    for i, particle in enumerate(particle_lst):
        print(
            f'Particle number: {i}. '
            f'List length: {len(particle.x_positions_lst)} and '
            f'x-positions list: {particle.x_positions_lst}'
        )


def main() -> None:
    # Inicjalizacja cząstek
    initial_x, initial_y = [1, 3, 5, 7], [1, 3, 5, 7]
    initial_vx, initial_vy = [1, 1, 2, 0.5], [1, 0.5, 0.25, 4]

    particle_lst = []
    for x, y, vx, vy in zip(initial_x, initial_y, initial_vx, initial_vy):
        particle_lst.append(Particle(radius, x, y, vx, vy))
    # print(len(particle_lst))

    simulate(particle_lst)
    simulation_results(particle_lst)


if __name__ == '__main__':
    main()