Lennard-Jones力3原子系统数值积分中r_x恒为0问题求助

我帮你仔细排查了代码里的问题，发现几个关键错误导致r_x始终为0且原子位置没有更新，下面逐个分析并给出修正方案：

1. 加速度计算传入了错误的位置数组 #

在RunMD函数中，你初始化了xPositions = np.zeros((number_of_steps, 3))后，直接用这个全0数组调用GetAcc计算初始加速度——这完全偏离了你的初始原子位置x和y。后续循环里你依然传入未更新的xPositions，导致计算时用的全是0值位置，自然r_x始终为0。

2. 原子对循环范围错误，漏掉大部分相互作用 #

你在GetAcc里的循环写了：

for i in range(0, xPositions.shape[0]-1):
    for j in range(i+1, xPositions.shape[0]-1):

对于3原子系统，xPositions.shape[0]是3，这会导致i只能取0、1，j只能取1（当i=0时），完全漏掉了原子0-2、1-2这两对的相互作用，加速度计算完全不完整。

3. 全局变量滥用+距离判断逻辑错误 #

• r_x、rmag这类临时变量没必要设为全局，局部变量足够用，全局变量反而容易引发意外覆盖。
• rmag是单个原子对的距离标量，但你用了rmag[0]==0这种数组索引判断，会直接报错，正确的做法是判断距离是否接近0（避免除以0）。

4. 初始加速度和速度设置错误 #

初始加速度基于全0位置计算，结果为0；初始速度也设为0，导致位置更新公式x = x + v_x*dt + 0.5*a_x*dt²完全没有变化，原子位置一直停留在初始值。

import numpy as np
np.seterr(invalid="ignore")

m = 1
epsilon = 0.84
sigma = 2.56

def GetLJForce(r):
    return 48 * epsilon * (sigma**12) / (r**13) - 24 * epsilon * (sigma**6) / (r**7)

def GetAcc(x_current, y_current):
    num_atoms = x_current.size
    xAcc = np.zeros((num_atoms, num_atoms), dtype=np.float64)
    yAcc = np.zeros((num_atoms, num_atoms), dtype=np.float64)
    
    # 遍历所有i<j的原子对，避免重复计算
    for i in range(num_atoms):
        for j in range(i+1, num_atoms):
            r_x = x_current[j] - x_current[i]
            r_y = y_current[j] - y_current[i]
            rmag = np.sqrt(r_x**2 + r_y**2)
            
            # 避免除以0，设置极小值替代
            if rmag < 1e-10:
                rmag = 1e-10
            
            force_scalar = GetLJForce(rmag)
            force_x = force_scalar * r_x / rmag
            force_y = force_scalar * r_y / rmag
            
            # 牛顿第三定律，相互作用力大小相等方向相反
            xAcc[i,j] = force_x / m
            xAcc[j,i] = -force_x / m
            yAcc[i,j] = force_y / m
            yAcc[j,i] = -force_y / m
    
    # 对每个原子的所有相互作用加速度求和
    return np.sum(xAcc, axis=1), np.sum(yAcc, axis=1)

def UpdatexPos(x, v_x, a_x, dt):
    return x + v_x*dt + 0.5*a_x*dt**2

def UpdateyPos(y, v_y, a_y, dt):
    return y + v_y*dt + 0.5*a_y*dt**2

def UpdatexVel(v_x, a_x, a1_x, dt):
    return v_x + 0.5*(a_x + a1_x)*dt

def UpdateyVel(v_y, a_y, a1_y, dt):
    return v_y + 0.5*(a_y + a1_y)*dt

def RunMD(dt, number_of_steps, x_init, y_init):
    xPositions = np.zeros((number_of_steps, 3), dtype=np.float64)
    yPositions = np.zeros((number_of_steps, 3), dtype=np.float64)
    
    # 初始化当前位置和速度
    x_current = x_init.copy()
    y_current = y_init.copy()
    v_x = np.zeros_like(x_current)
    v_y = np.zeros_like(y_current)
    
    # 用初始位置计算初始加速度
    a_x, a_y = GetAcc(x_current, y_current)
    
    for i in range(number_of_steps):
        # 更新位置
        x_current = UpdatexPos(x_current, v_x, a_x, dt)
        y_current = UpdateyPos(y_current, v_y, a_y, dt)
        
        # 计算新的加速度
        a1_x, a1_y = GetAcc(x_current, y_current)
        
        # 更新速度
        v_x = UpdatexVel(v_x, a_x, a1_x, dt)
        v_y = UpdateyVel(v_y, a_y, a1_y, dt)
        
        # 更新加速度为最新值
        a_x, a_y = a1_x, a1_y
        
        # 保存当前位置到结果数组
        xPositions[i, :] = x_current
        yPositions[i, :] = y_current
    
    return xPositions, yPositions

# 初始位置
x = np.array([1, 9, 15], dtype=np.float64)
y = np.array([16, 22, 26], dtype=np.float64)

# 运行模拟
sim_xpos, sim_ypos = RunMD(0.1, 10, x, y)
print(sim_xpos)

1. 传入正确的位置计算加速度：初始加速度用真实的初始位置计算，后续循环用当前步的原子位置调用GetAcc，确保相互作用计算正确。
2. 修复原子对循环范围：遍历所有i<j的原子对，覆盖3原子系统的全部3对相互作用。
3. 移除冗余全局变量：临时变量改为局部，避免变量污染。
4. 修复除以0问题：用极小值1e-10替代接近0的距离，既避免报错，又不会严重影响力的计算。
5. 优化计算逻辑：每次循环只调用一次GetAcc，避免冗余计算；用np.sum(xAcc, axis=1)正确求和每个原子的总加速度。

运行修正后的代码，你会看到r_x不再为0，原子位置也会随着模拟步骤正常更新。

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