admin
d930fb558c
对比 WorkBuddy/Claude/Codex 三款 AI 工具对同一代码库的 优化建议，以表格形式评价各自优劣（Bug 发现/代码质量/战略 思维/代码示例），最终整合为 6 阶段实施计划（11 人天）。
16 KiB
Dynamics 项目优化建议
分析日期：2026-06-12
分析范围：D:\Share\Data\aliyun-gitea\dynamics完整代码库
当前算例：examples/case06（120 原子一维链横波，NT=100,000 步，C 引擎）
另见：optimization/workbuddy.md（架构拆分、全局变量封装、测试等宏观建议）
本文档聚焦以下四类问题，不与 workbuddy.md 重复：
- 已确认的 Bug（可能导致运行崩溃或结果错误）
- Python 引擎性能（向量化等具体代码改动）
- 引擎一致性问题（C/C++/Fortran 行为差异）
- 小的代码质量问题（错误文本、重复代码等）
一、已确认的 Bug
B1. run_simulation() 内引用了不在作用域的 config 变量 ⚠️ 严重
# 问题代码（run_simulation 函数内部，config 未传入）
"camera_distance": str(config.get("camera_distance", 40.0)),
"camera_elevation": str(config.get("camera_elevation", 0)),
"camera_azimuth": str(config.get("camera_azimuth", 0)),
run_simulation() 是独立函数，不接收 config 参数，但函数内部直接使用了 run_from_config() 的局部变量 config。只要使用 Python 引擎（engine: python）就会触发 NameError。
修复方案：这三个字段的值已在 run_from_config 里读取并存入全局变量（camera_keyframes_raw 已有），
对 camera_distance/elevation/azimuth 同样声明为全局变量并在 run_from_config 中赋值：
# 在模块顶部全局变量区添加（compute.py ~65 行附近）
camera_distance = 40.0
camera_elevation = 0
camera_azimuth = 0
# 在 run_from_config 中赋值（~756 行附近，已有 use_marker 赋值的位置）
camera_distance = float(config.get("camera_distance", 40.0))
camera_elevation = float(config.get("camera_elevation", 0))
camera_azimuth = float(config.get("camera_azimuth", 0))
# run_simulation 内改为读全局变量
"camera_distance": str(camera_distance),
"camera_elevation": str(camera_elevation),
"camera_azimuth": str(camera_azimuth),
B2. dynamics.py 绘图块引用未定义的变量 ⚠️ 严重
# step_plot=1 时触发，但这些变量从未赋值
time_arr = np.arange(NT) * DT # NT、DT 在此作用域未定义
n_atoms = all_x.shape[1] # all_x 不存在
atom_ids_list = data.get("atom_ids", ...) # data 不存在
这段绘图代码写于旧版（run_simulation 返回完整轨迹时），现在 run_from_config 已不再在此作用域填充这些变量，整块代码是失效的死代码。当前所有案例恰好设置 step_plot: 0 所以不触发，但一旦用户设为 1 就崩溃。
修复方案：从 run_from_config 的返回值读取数据重构绘图代码，或暂时添加保护：
# 临时保护（快速修复）
if not no_plot and config.get("step_plot", 1):
print("[run] 警告: step_plot 功能暂未实现，已跳过")
B3. plot_wave.py 使用旧 JSON 格式字段读取新格式 display.txt ⚠️ 严重
文件：plot_wave.py:22-27 和 plot_wave.py:96-111
# plot_wave.py 的 load_disp_data()
def load_disp_data(output_dir):
return compute.load_text_data(disp_path) # 调用了 JSON 格式的加载函数！
# 然后访问旧 JSON 字段名
n_frames = int(data["n_frames"]) # 新格式没有这个键
x = np.array(data["disp_all_x"]) # 新格式没有这个键
masses = np.array(data["atom_masses"]) # 新格式没有这个键
但 display.txt 现在已是新文本格式（由 save_display_txt 写出），正确的加载函数是 compute.load_display_txt()，且字段名完全不同（frames_x 而非 disp_all_x）。
影响：step_plot_wave: 1 时必然崩溃，case05/case06 的波形动画功能完全失效。
修复方案：将 load_disp_data 改为调用 load_display_txt，并适配字段名：
def load_disp_data(output_dir):
disp_path = os.path.join(output_dir, "display.txt")
if not os.path.exists(disp_path):
raise FileNotFoundError(f"找不到 {disp_path}")
return compute.load_display_txt(disp_path) # 改这里
# 后续访问字段名也需同步修改：
# data["disp_all_x"] → disp_data["frames_x"]
# data["n_frames"] → disp_data["frames_x"].shape[0]
# data["atom_masses"] → 从 header_fields 读取或从 bond/coord 文件读
注：
plot_wave.py依赖atom_masses、bond_pairs等物理量，这些在新的 display.txt 格式中 通过header_fields字符串存储，需要额外解析。较彻底的修复需要在 display.txt header 中保留atom_masses序列（目前atom_radii已以逗号分隔字符串保存，可仿照此方式）。
B4. run_dynamics.py（case06）描述文字错误
文件：examples/case06/run_dynamics.py:34
# 当前（错误）
parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01")
# 应为
parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case06")
B5. draw.py 裸 except 吞掉所有异常
文件：draw.py:97
try:
alpha_list = [float(x) for x in raw_alpha.split(",")]
if len(alpha_list) != 6:
alpha_list = alpha_list * 6
except: # 裸 except：连 KeyboardInterrupt 都会捕获
alpha_list = [float(raw_alpha)] * 6
修复：改为 except (ValueError, AttributeError):，并加 print 提示。
二、Python 引擎性能优化
背景：case06 使用 C 引擎，以下优化主要影响
engine: python路径。 当前 Python vs C 速度比约为 6-8 倍（100,000 步，120 原子）。
P1. 弹簧力计算向量化（最大收益）
# 当前：Python for 循环，119 键 × 100,000 步 ≈ 1200 万次迭代
for bond_idx, (idx_1, idx_2) in enumerate(BOND_PAIRS):
dx = x[idx_2] - x[idx_1]
dy = y[idx_2] - y[idx_1]
dz = z[idx_2] - z[idx_1]
dist = np.sqrt(dx * dx + dy * dy + dz * dz)
if dist <= 1e-12:
continue
stretch = dist - BOND_REST_LENGTHS[bond_idx]
force_mag = BOND_STIFFNESS[bond_idx] * stretch
...
向量化版本：
if ELASTIC_FORCE and BOND_PAIRS is not None and len(BOND_PAIRS) > 0:
i_arr = BOND_PAIRS[:, 0] # (n_bonds,)
j_arr = BOND_PAIRS[:, 1]
dx = x[j_arr] - x[i_arr]
dy = y[j_arr] - y[i_arr]
dz = z[j_arr] - z[i_arr]
dist = np.sqrt(dx * dx + dy * dy + dz * dz)
valid = dist > 1e-12
inv_dist = np.where(valid, 1.0 / np.maximum(dist, 1e-12), 0.0)
stretch = dist - BOND_REST_LENGTHS
force_mag = BOND_STIFFNESS * stretch * inv_dist # 归一化后的分量力幅
fx_bond = force_mag * dx
fy_bond = force_mag * dy
fz_bond = force_mag * dz
np.add.at(fx, i_arr, fx_bond)
np.add.at(fx, j_arr, -fx_bond)
np.add.at(fy, i_arr, fy_bond)
np.add.at(fy, j_arr, -fy_bond)
np.add.at(fz, i_arr, fz_bond)
np.add.at(fz, j_arr, -fz_bond)
预期加速：Python 引擎整体 5-15 倍（对于链状键合体系，弹簧力是瓶颈）。
np.add.at支持重复索引的原子累加，对有分叉键的复杂体系也正确。 若键无重复端点（链状），可改为fx[i_arr] += fx_bond; fx[j_arr] -= fx_bond更快。
P2. apply_fixed_constraints() 减少临时数组分配
# 当前：每步创建 2 个 (n_atoms, 3) 的临时数组
positions = np.column_stack((x, y, z))
velocities = np.column_stack((vx, vy, vz))
positions = np.where(fixed, ATOM_POSITIONS, positions)
velocities = np.where(fixed, 0.0, velocities)
return positions[:, 0], positions[:, 1], positions[:, 2], ...
对于 case06（x/y 全部固定，z 自由），每步都在做 120 原子 × 3 列的 column_stack，但实际上只有 z 方向的原子可以运动。
优化版本（预先计算掩码，直接原地修改）：
# 在模块初始化时（run_from_config 中）预计算固定掩码
_FIXED_MASK_X = ATOM_FIXED[:, 0] != 0 # bool 数组，提前缓存
_FIXED_MASK_Y = ATOM_FIXED[:, 1] != 0
_FIXED_MASK_Z = ATOM_FIXED[:, 2] != 0
def apply_fixed_constraints(x, y, z, vx, vy, vz):
if np.any(_FIXED_MASK_X):
x[_FIXED_MASK_X] = ATOM_POSITIONS[_FIXED_MASK_X, 0]
vx[_FIXED_MASK_X] = 0.0
if np.any(_FIXED_MASK_Y):
y[_FIXED_MASK_Y] = ATOM_POSITIONS[_FIXED_MASK_Y, 1]
vy[_FIXED_MASK_Y] = 0.0
if np.any(_FIXED_MASK_Z):
z[_FIXED_MASK_Z] = ATOM_POSITIONS[_FIXED_MASK_Z, 2]
vz[_FIXED_MASK_Z] = 0.0
return x, y, z, vx, vy, vz
P3. run_simulation() 中冗余的 frame_indices 列表
frame_indices = []
...
frame_indices.append(step) # 每 NSTEP 步追加一次，仅用于最后计算 n_frames_actual
# 实际上等价于：
n_frames_actual = record_steps // NSTEP # 直接计算，O(1)
100,000 步中约有 1,000 次 append，开销小但无必要。删除 frame_indices 列表，用计数器替代。
P4. apply_driving_force() 内部创建临时 t_vec
# 每次调用都分配一个 3 元素 array
t_vec = np.array([t, t, t], dtype=np.float64)
pos_drive = d["amp"] * np.cos(2.0 * np.pi * d["freq"] * t_vec + d["phi"])
t 是标量，d["freq"] 和 d["phi"] 是 (3,) 数组，直接用标量广播即可：
pos_drive = d["amp"] * np.cos(2.0 * np.pi * d["freq"] * t + d["phi"])
vel_drive = -d["amp"] * 2.0 * np.pi * d["freq"] * np.sin(2.0 * np.pi * d["freq"] * t + d["phi"])
P5. GRAVITY_INTERACTION 的 O(N²) 双重循环
# 当前：纯 Python 双重循环
for i in range(n):
for j in range(i + 1, n):
...
向量化方案（适用于 N 不太大时）：
if GRAVITY_INTERACTION:
# 构造所有原子对的差向量
xi = x[:, np.newaxis]; xj = x[np.newaxis, :]
yi = y[:, np.newaxis]; yj = y[np.newaxis, :]
zi = z[:, np.newaxis]; zj = z[np.newaxis, :]
r2 = (xj - xi)**2 + (yj - yi)**2 + (zj - zi)**2
np.fill_diagonal(r2, np.inf) # 排除自身
r = np.sqrt(r2)
mi = m[:, np.newaxis]; mj = m[np.newaxis, :]
# 力幅（标量，上三角）
f_mag = GRAVITY_STRENGTH * mi * mj / r2
fx_ij = f_mag * (xj - xi) / r # (n, n) 矩阵
fy_ij = f_mag * (yj - yi) / r
fz_ij = f_mag * (zj - zi) / r
fx += np.sum(fx_ij - fx_ij.T, axis=1) * 0.5 # 利用反对称性
fy += np.sum(fy_ij - fy_ij.T, axis=1) * 0.5
fz += np.sum(fz_ij - fz_ij.T, axis=1) * 0.5
此向量化在 N ≤ ~500 时比 Python 循环快约 20 倍；N > 1000 时内存占用高（N² 矩阵）， 需要分块或使用 Barnes-Hut 树算法。
三、引擎一致性问题
E1. 弹性力计算：Python 与 C 引擎的差异
文件：compute.py:1273-1277 vs engines/c/main.c
Python 引擎将整根键的弹性势能 ½k(d-d₀)² 记给原子 i（不均分），而 C/C++ 引擎均分给两端。
这不影响力的计算（力相同），但如果未来基于势能做分析，两个引擎的势能分配会不同。
建议注释中明确说明约定。
E2. C++ 引擎 save_trajectory 默认值为 1
// engines/c/param.json: "save_trajectory": 0
// engines/cpp/param.json: "save_trajectory": 1 ← 与 C 不一致
用户切换引擎时行为不同，可能意外生成大文件。建议统一为 0。
E3. Fortran 引擎不支持 save_trajectory=0
Fortran 引擎总是输出 JSON 格式的 trajectory.txt，不支持直接写 display.txt。
workbuddy.md 中已将此列为 P0，此处确认影响范围：
- case 使用
engine: fortran时必然崩溃（Python 侧期待display.txt） - 应参照 C 引擎的
write_display_txt函数实现
E4. 外部引擎每次运行都重新校准
每次运行前跑 min(1000, NT//10) 步测速。对于 NT=100,000 的 case06，校准跑 10,000 步，
占总运算量的 10%。
快速修复（跳过短运行的校准）：
# 当总步数 < 5000 时跳过校准，直接用默认估计
if total_steps < 5000:
_step_time = 1e-5 # 合理的保守估计
_est_total = _step_time * total_steps
else:
# 现有校准逻辑...
彻底修复：缓存校准结果到 engines/{engine}/_calib_cache.json，key 为 (n_atoms, method)，
命中时跳过校准（workbuddy.md §2.1 有详细代码示例）。
四、代码质量问题
Q1. 相机路径解析函数重复
dynamics.py:35-69 的 _load_camera_kf() 与
compute.py:90-136 的 _load_camera_motion() 是功能完全相同的函数，
仅变量名略有差异。
建议：保留 compute.py 中的版本，dynamics.py 改为从 compute 导入：
# dynamics.py
from compute import _load_camera_motion as _load_camera_kf
Q2. load_parameters() 函数已废弃但仍保留
这是老版通过 key = value 文本格式读参数的函数，现在所有案例已改用 YAML。
只有 compute.py 的 main() 还调用它，而 compute.main() 本身也已被 dynamics.py 取代。
整个函数（77 行）加上 main()（25 行）都可以安全删除，或移入 tools/ 目录存档。
Q3. 边界条件逻辑不一致
文件：compute.py:1380-1387 vs compute.py:1420-1428
def Limit_in_box(a, amin, amax, va):
# 反弹：碰壁时速度反向（弹性碰撞）
va = np.where(over | under, -va, va)
def wrap_position(x, y, z):
# 周期性边界：从一侧出去从另一侧进入
x = np.where(x > X_MAX, X_MIN, x)
在 run_simulation 中两者都被调用：先通过 apply_motion_update() 执行 Limit_in_box（反弹），
再执行 wrap_position()（周期）。顺序执行时，Limit_in_box 已把粒子限制在盒内，
wrap_position 实际上永远不会触发——但这一逻辑关系无文档说明，容易误导读者。
建议：在 run_simulation 中加一行注释说明两者的分工，或删除其中一个。
Q4. 案例说明文字未更新（case06）
文件：examples/case06/run_dynamics.py:34
见 B4，在此重申。
五、总结与优先级
| 编号 | 问题 | 优先级 | 预估工作量 |
|---|---|---|---|
| B1 | run_simulation 内 config 未定义 |
🔴 立即 | 15 min |
| B2 | dynamics.py 绘图块死代码（NameError） |
🔴 立即 | 30 min |
| B3 | plot_wave.py 使用旧格式加载新 display.txt |
🔴 立即 | 1-2 h |
| B4 | case06 描述文字写着 case01 | 🟡 低 | 1 min |
| B5 | draw.py 裸 except |
🟡 低 | 5 min |
| P1 | 弹簧力向量化（Python 引擎） | 🟢 中 | 1 h |
| P2 | apply_fixed_constraints 减少临时数组 |
🟢 低 | 30 min |
| P3 | frame_indices 列表改计数器 |
🟢 低 | 5 min |
| P4 | apply_driving_force 去掉 t_vec |
🟢 低 | 5 min |
| P5 | GRAVITY_INTERACTION O(N²) 向量化 |
🟢 中 | 1 h |
| E1 | 势能归属约定加注释 | 🟡 低 | 10 min |
| E2 | C++ save_trajectory 默认值统一 |
🟡 低 | 5 min |
| E3 | Fortran 引擎支持 display.txt 输出 | 🔴 高 | 2-3 h |
| E4 | 外部引擎校准缓存 | 🟡 中 | 1-2 h |
| Q1 | 相机解析函数去重 | 🟡 低 | 15 min |
| Q2 | 废弃的 load_parameters + main() 删除 |
🟡 低 | 15 min |
| Q3 | 边界条件逻辑加注释 | 🟡 低 | 10 min |
| Q4 | case06 描述文字 | 🟢 低 | 1 min |
最高回报的三件事（按投入产出比排序）：
- B3（
plot_wave.py格式不匹配）：修复后波形动画功能恢复，影响 case05/case06 - B1（
config未定义）：修复后 Python 引擎可完整运行 case06 - P1（弹簧力向量化）：修复后 Python 引擎速度提升 5-15 倍，适合大规模测试和教学演示
本文档由 Claude Sonnet 4.6 自动生成，基于代码静态分析，未实际运行测试验证。建议在实施前先运行相关案例确认 Bug 现象。