分子气体模拟中NaN异常及结果错误问题求助

遇到的问题 #

• 使用np.sqrt()计算分子间距时，得到全是nan的列表；
• 替换为()**0.5消除NaN后，结果仍错误：判断dist12 == 0返回0, 0，导致除前4个粒子外，其余粒子位置列表大多为0。

预期代码流程 #

1. 从Czastka类初始化16个粒子及对应初始条件；
2. simulate()函数迭代1000步，通过lennard_jones_forces()计算每个粒子受其他粒子的合力，粒子间距超过2.5*sigma时忽略相互作用；
3. 通过leapforg_single_step()计算新的位置和速度；
4. 调用p1.update_position()更新位置并应用周期性边界条件（通过取模运算self.x = new_x % box_size实现）。

完整代码 #

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 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]

particle_lst = []
for i in range(0, 4):
    for j in range(0, 4):
        particle_lst.append(Czastka(radius, initial_x[j], initial_y[j], initial_vx[j], initial_vy[j]))
#print(len(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
    print(r12[0], ' ', r12[1])
    if dist12 == 0:
        return 0, 0

    # 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):
        print(Force21_x, ' ', Force21_y)
        print("Nan detected in force calculation")

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


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

def animation():
    pass

问题根源与修复方案 #

1. 粒子初始位置重复导致的dist12 == 0问题 #

你的初始化代码中，嵌套循环生成了16个粒子，但所有粒子的初始位置都是[1,3,5,7]的4种组合重复4次，大量粒子位置完全重合。重合粒子的间距为0，触发dist12 == 0返回(0,0)，同时后续力计算会出现除零异常，导致粒子运动停滞。

修复：生成均匀分布的不重复初始位置：

# 替换原有粒子初始化代码
initial_x = np.linspace(1, 7, 4)
initial_y = np.linspace(1, 7, 4)
particle_lst = []
for x in initial_x:
    for y in initial_y:
        # 给每个粒子分配随机初始速度（符合温度设定）
        vx = np.random.normal(np.sqrt(k_B*temp/m), 0.1)
        vy = np.random.normal(np.sqrt(k_B*temp/m), 0.1)
        particle_lst.append(Czastka(radius, x, y, vx, vy))

2. Leapfrog算法实现错误 #

标准Leapfrog算法需要先基于当前力更新半步速度，再更新位置，最后用新位置的力更新完整速度。你的实现中两次使用同一组力（旧位置的力），导致运动计算完全偏离预期。

修复：调整simulate()函数实现正确的Leapfrog流程：

def simulate():
    global box_size
    num_steps = 1000

    # 初始化记录初始状态
    for p in particle_lst:
        p.append_postions(p.x, p.y)
        p.update_velocities(p.vx, p.vy)

    for step in range(num_steps):
        # 第一步：计算所有粒子当前位置的力并存储
        current_forces = []
        for i, p1 in enumerate(particle_lst):
            fx, fy = 0, 0
            for j, p2 in enumerate(particle_lst):
                if i != j:
                    f_x, f_y = lennard_jones_forces(p1.x, p1.y, p2.x, p2.y)
                    fx += f_x
                    fy += f_y
            current_forces.append((fx, fy))
        
        # 第二步：更新每个粒子的位置和速度
        for i, p1 in enumerate(particle_lst):
            fx, fy = current_forces[i]
            # 计算半步速度
            vx_half = p1.vx + 0.5 * (fx/m) * dt
            vy_half = p1.vy + 0.5 * (fy/m) * dt
            # 更新位置并应用周期性边界
            x_next = (p1.x + vx_half * dt) % box_size
            y_next = (p1.y + vy_half * dt) % box_size
            p1.x, p1.y = x_next, y_next
            p1.append_postions(x_next, y_next)
            
            # 计算新位置下的合力
            new_fx, new_fy = 0, 0
            for j, p2 in enumerate(particle_lst):
                if i != j:
                    f_x, f_y = lennard_jones_forces(p1.x, p1.y, p2.x, p2.y)
                    new_fx += f_x
                    new_fy += f_y
            
            # 更新完整速度
            vx_next = vx_half + 0.5 * (new_fx/m) * dt
            vy_next = vy_half + 0.5 * (new_fy/m) * dt
            p1.vx, p1.vy = vx_next, vy_next
            p1.update_velocities(vx_next, vy_next)

3. 周期性边界条件未应用于间距计算 #

LJ势模拟中，周期性边界需要用最小镜像约定计算粒子间距，即取粒子到最近镜像的距离，而非直接坐标差。这也是之前np.sqrt()出现NaN的潜在原因（极端情况下坐标差平方可能因浮点误差出现负值）。

修复：修改lennard_jones_forces()函数的间距计算逻辑：

def lennard_jones_forces(x1, y1, x2, y2):
    global sigma, eps, box_size

    # 应用最小镜像约定
    dx = x2 - x1
    dx = dx - box_size * round(dx / box_size)
    dy = y2 - y1
    dy = dy - box_size * round(dy / box_size)
    
    dist12_sq = dx**2 + dy**2
    dist12 = np.sqrt(dist12_sq)
    
    # 用极小值判断避免浮点精度问题，不用直接等于0
    if dist12 < 1e-10:
        return 0.0, 0.0

    # 提前判断截断条件，减少无效计算
    if dist12 > 2.5 * sigma:
        return 0.0, 0.0

    # 正确计算LJ力分量
    sigma_over_dist = sigma / dist12
    term1 = sigma_over_dist ** 14
    term2 = 0.5 * sigma_over_dist ** 8
    force_magnitude = -(48 * eps / sigma**2) * (term1 - term2) * dist12_sq
    fx = force_magnitude * dx / dist12
    fy = force_magnitude * dy / dist12

    return fx, fy

4. NaN问题彻底解决 #

通过最小镜像约定保证间距平方非负，同时用极小值判断避免除零，此时np.sqrt()不会再返回NaN。

内容的提问来源于stack exchange，提问作者Mark