我们的模型是由奇数个硬球串联而成的自驱动杆，在二维周期性边界介质中以固定速度运动，具备x/y平动扩散和转动扩散三个独立扩散常数，初始取向随机。当不同杆的球体进入相互作用半径时，会产生Lennard-Jones势衍生的作用力，该力通过力矩驱动杆缓慢旋转，最终使所有杆取向一致。

现有代码运行后，杆的旋转效果不符合物理预期，以下是具体修正点：

1. 核心问题：错误的角度更新逻辑 #

当前update_angles直接将杆的角度设置为杆中心合力的方向，完全违背了“扭矩驱动旋转”的物理机制——杆的旋转由每个球体受到的力对杆中心的力矩总和决定，而非直接朝向合力方向。

2. 球体位置生成错误 #

当前get_sphere_positions会重复添加杆的中心球体（i=0时，正负方向的位置均为杆中心），导致杆的球体数量变为偶数（比如设置9个球的杆实际生成10个位置）。

修正代码： #

def get_sphere_positions(start_pos, angle, num_spheres):
    sphere_positions = [tuple(start_pos)]  # 单独添加中心球
    half_length = (num_spheres - 1) // 2
    for i in range(1, half_length + 1):
        # 向杆的负方向生成球体
        x = start_pos[0] - 2 * i * SPHERE_RADIUS * np.cos(angle)
        y = start_pos[1] - 2 * i * SPHERE_RADIUS * np.sin(angle)
        sphere_positions.append((x, y))
        # 向杆的正方向生成球体
        x1 = start_pos[0] + 2 * i * SPHERE_RADIUS * np.cos(angle)
        y1 = start_pos[1] + 2 * i * SPHERE_RADIUS * np.sin(angle)
        sphere_positions.append((x1, y1))
    return sphere_positions

3. 缺失扩散项（平动+转动） #

模型明确要求具备x/y平动扩散和转动扩散，但现有代码完全未实现，导致运动过于规则，不符合介质中自驱动粒子的布朗运动特性。

修正代码： #

# 添加到参数定义区
D_TRANS = 0.1  # 平动扩散常数
D_ROT = 0.05   # 转动扩散常数
DT = 0.1       # 时间步长

def update_positions(positions, angles, velocity, dt):
    theta = angles
    dx = velocity * np.cos(theta) * dt
    dy = velocity * np.sin(theta) * dt
    # 添加平动扩散高斯噪声
    dx += np.sqrt(2 * D_TRANS * dt) * np.random.randn(len(positions))
    dy += np.sqrt(2 * D_TRANS * dt) * np.random.randn(len(positions))
    positions[:, 0] += dx
    positions[:, 1] += dy

    # 周期性边界处理
    positions[:, 0] = positions[:, 0] % WIDTH
    positions[:, 1] = positions[:, 1] % HEIGHT
    return positions

4. 力矩计算与角度更新的正确实现 #

重写力和扭矩计算逻辑，基于力矩驱动角度变化：

def compute_forces_and_torques(positions, angles):
    num_rods = len(positions)
    forces = np.zeros_like(positions)  # 杆中心的合力
    torques = np.zeros(num_rods)       # 杆受到的总扭矩

    for i in range(num_rods):
        sphere_pos_i = np.array(get_sphere_positions(positions[i], angles[i], NUM_SPHERES_PER_ROD))
        for j in range(i + 1, num_rods):
            sphere_pos_j = np.array(get_sphere_positions(positions[j], angles[j], NUM_SPHERES_PER_ROD))
            for idx_i, pos_i in enumerate(sphere_pos_i):
                for idx_j, pos_j in enumerate(sphere_pos_j):
                    # 周期性边界下的最短距离向量
                    dist_vec = pos_j - pos_i
                    dist_vec[0] = np.mod(dist_vec[0] + WIDTH/2, WIDTH) - WIDTH/2
                    dist_vec[1] = np.mod(dist_vec[1] + HEIGHT/2, HEIGHT) - HEIGHT/2
                    dist = np.linalg.norm(dist_vec)
                    if dist < NEIGHBOR_DISTANCE_THRESHOLD and dist > 1e-6:
                        # 计算Lennard-Jones力
                        r6 = (SIGMA / dist) ** 6
                        r12 = r6 ** 2
                        force_mag = 24 * EPSILON * (2 * r12 - r6) / dist
                        force_vec = force_mag * (dist_vec / dist)
                        # 作用与反作用力
                        forces[i] += force_vec
                        forces[j] -= force_vec
                        # 计算对杆中心的力矩（叉乘）
                        r_i = pos_i - positions[i]
                        torques[i] += np.cross(r_i, force_vec)
                        r_j = pos_j - positions[j]
                        torques[j] += np.cross(r_j, -force_vec)

    return forces, torques

def update_angles(angles, torques, dt):
    # 简化假设：转动惯量为1，角加速度=扭矩/转动惯量
    angular_acc = torques
    angles += angular_acc * dt
    # 添加转动扩散噪声
    angles += np.sqrt(2 * D_ROT * dt) * np.random.randn(len(angles))
    # 角度归一化到[0, 2π)
    angles = np.mod(angles, 2 * np.pi)
    return angles

5. 动画循环同步调整 #

修改animate函数适配新逻辑：

def animate(frame):
    global positions, angles
    positions = update_positions(positions, angles, VELOCITY, DT)
    forces, torques = compute_forces_and_torques(positions, angles)
    angles = update_angles(angles, torques, DT)

    ax.clear()
    for pos, angle in zip(positions, angles):
        draw_rod(ax, pos, angle, NUM_SPHERES_PER_ROD)

    ax.set_xlim(0, WIDTH)
    ax.set_ylim(0, HEIGHT)
    ax.set_aspect('equal')

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