feat: 真蛙跳法重构（Python/C/C++/Fortran 四引擎统一）

- 新增 compute_accel_conservative / accel_conservative：
  保守力加速度（弹簧+重力+原子间引力），不含阻尼，供蛙跳专用
- 重写 leapfrog_step / leapfrog_full：
  - 无阻尼：纯辛积分器，每步 1 次力计算（原 Velocity-Verlet 需 2 次）
  - 有阻尼：半隐式处理 v(t+dt/2)=[v(t-dt/2)*(1-α)+a_c*dt]/(1+α)，无条件稳定
- 主循环加初始化反向半步 v(-dt/2)=v(0)-0.5*a_c(0)*dt
- 修复 C/C++ number of frames 字段写采样帧数而非总积分步数的 bug
- Python 引擎：新增 display.npz 二进制格式，draw.py/plot_wave.py 优先读取
- 编译参数统一为 -O3 -march=native -ffast-math

.gitignore

@@ -99,3 +99,7 @@ desktop.ini
 *.zip
 *.tar.gz
 *.tar.bz2
+output_test/
+
+optimization
+optimization/*

compute.py

@@ -5,11 +5,12 @@ compute.py
 
 功能：
   1. 运行 NT 步物理模拟（ kinematics / dynamics 等运动模式）
-  2. 将每一步的 (x, y, z, vx, vy, vz) 保存到 output/trajectory.txt
-  3. 同时保存所有模拟参数元数据
+  2. 按 NSTEP 抽帧，输出 output/display.txt（新文本格式）
+  3. 可选（save_trajectory=1）保留完整轨迹 output/trajectory.txt（JSON）
 """
 
 import json
+import re
 import numpy as np
 import os
 import platform
@@ -63,6 +64,14 @@ ball_color_b = None
 box_color_r = None
 box_color_g = None
 box_color_b = None
+use_marker = 0
+camera_keyframes_raw = ""
+camera_distance = 40.0
+camera_elevation = 0.0
+camera_azimuth = 0.0
+FIXED_MASK_X = None
+FIXED_MASK_Y = None
+FIXED_MASK_Z = None
 
 # 力开关
 GRAVITY_FIELD = 1         # 均匀重力场
@@ -84,6 +93,55 @@ Z_MIN = None
 Z_MAX = None
 
 
+def _load_camera_motion(path):
+    """读取 move_camera.txt（速度段格式），返回 JSON 字符串。
+    
+    格式：每行是一个运动段
+      start-end  vx=f1 vy=f2 vz=f3  rx=d1 ry=d2 rz=d3
+    示例：
+      1-60   vx=1.0  rx=10
+      30-90  vy=2.0  ry=20  rz=10
+    
+    返回 JSON: [{"start":N,"end":N,"v":[x,y,z],"r":[x,y,z]},...]
+    """
+    import re
+    if not os.path.exists(path):
+        print(f"[compute] 警告: 未找到 {path}，跳过运动相机")
+        return ""
+    segments = []
+    with open(path, "r", encoding="utf-8") as f:
+        for line in f:
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            # 解析帧范围：支持 "all"（全程）或 "N-M"（区间）
+            if line.lower().startswith("all") or re.match(r'^\s*all\s', line, re.IGNORECASE):
+                start, end = 0, 10**9
+            else:
+                m = re.match(r'(\d+)\s*-\s*(\d+)', line)
+                if not m:
+                    continue
+                start, end = int(m.group(1)), int(m.group(2))
+            v = [0.0, 0.0, 0.0]
+            r = [0.0, 0.0, 0.0]
+            # 解析 vx=, vy=, vz=
+            for i, axis in enumerate(['x', 'y', 'z']):
+                m2 = re.search(r'v' + axis + r'\s*=\s*([-\d.]+)', line)
+                if m2:
+                    v[i] = float(m2.group(1))
+            # 解析 rx=, ry=, rz=
+            for i, axis in enumerate(['x', 'y', 'z']):
+                m2 = re.search(r'r' + axis + r'\s*=\s*([-\d.]+)', line)
+                if m2:
+                    r[i] = float(m2.group(1))
+            if any(v) or any(r):
+                segments.append({"start": start, "end": end, "v": v, "r": r})
+    if not segments:
+        return ""
+    import json
+    return json.dumps(segments)
+
+
 def _to_text_value(value):
     """Convert numpy-heavy objects into JSON-friendly plain Python values."""
     if isinstance(value, np.ndarray):
@@ -117,7 +175,7 @@ def _from_text_value(value):
 
 
 def save_text_data(path, data):
-    """Save structured simulation data as formatted JSON text."""
+    """Save structured simulation data as formatted JSON text (旧格式，保留兼容)."""
     os.makedirs(os.path.dirname(path), exist_ok=True)
     with open(path, "w", encoding="utf-8") as f:
         json.dump(_to_text_value(data), f, ensure_ascii=False, indent=2)
@@ -125,12 +183,208 @@ def save_text_data(path, data):
 
 
 def load_text_data(path):
-    """Load structured simulation data from JSON text."""
+    """Load structured simulation data from JSON text (旧格式)."""
     with open(path, "r", encoding="utf-8") as f:
         data = json.load(f)
     return _from_text_value(data)
 
 
+# ========================================================================
+# 新 display.txt 格式：纯文本，按帧分块
+#   第1行: number of frames: N
+#   第2行: number of particles: M
+#   第3行: frame:  1
+#   第4行: n   x            y            z           vx         vy         vz
+#   第5+行: 数据行（每个原子一行）
+#   重复第3-5行直到所有帧
+# ========================================================================
+
+def save_display_txt(path, frames_x, frames_y, frames_z,
+                     frames_vx, frames_vy, frames_vz,
+                     atom_ids, n_total_frames, n_total_particles,
+                     header_fields=None):
+    """Write display.txt in new text format.
+    
+    Args:
+        path: 输出文件路径
+        frames_x/y/z/vx/vy/vz: (n_frames, n_atoms) 数组
+        atom_ids: (n_atoms,) 原子编号数组
+        n_total_frames: 总帧数（含未采样）
+        n_total_particles: 总粒子数
+        header_fields: 可选的额外元数据字典（写入文件头之后）
+    """
+    os.makedirs(os.path.dirname(path), exist_ok=True)
+    n_frames = frames_x.shape[0]
+    n_atoms = frames_x.shape[1]
+    
+    # 格式化辅助：固定宽度，6位小数
+    def fmt(v): return f"{v:13.6f}"
+    
+    with open(path, "w", encoding="utf-8") as f:
+        f.write(f"number of frames: {n_total_frames}\n")
+        f.write(f"number of particles: {n_total_particles}\n")
+        # 写入额外元数据
+        if header_fields:
+            for k, v in header_fields.items():
+                f.write(f"{k}: {v}\n")
+        
+        for fr in range(n_frames):
+            f.write(f"\nframe:  {fr + 1}\n")
+            f.write(f"n   x            y            z           vx         vy         vz\n")
+            for a in range(n_atoms):
+                f.write(f"{atom_ids[a]:d}"
+                        f"{fmt(frames_x[fr, a])}"
+                        f"{fmt(frames_y[fr, a])}"
+                        f"{fmt(frames_z[fr, a])}"
+                        f"{fmt(frames_vx[fr, a])}"
+                        f"{fmt(frames_vy[fr, a])}"
+                        f"{fmt(frames_vz[fr, a])}\n")
+    return path
+
+
+def load_display_txt(path):
+    """Read display.txt new text format into numpy arrays（快速版）.
+    
+    Returns dict with keys: frames_x/y/z/vx/vy/vz, atom_ids,
+           n_total_frames, n_total_particles, header_fields
+    """
+    import re
+    
+    with open(path, "r", encoding="utf-8") as f:
+        raw = f.read()
+    
+    # 解析 header 行
+    header_fields = {}
+    n_total_frames = 0
+    n_total_particles = 0
+    
+    lines = raw.splitlines()
+    data_start = 0
+    for i, line in enumerate(lines):
+        line_stripped = line.strip()
+        if line_stripped.startswith("number of frames:"):
+            n_total_frames = int(line_stripped.split(":")[1].strip())
+        elif line_stripped.startswith("number of particles:"):
+            n_total_particles = int(line_stripped.split(":")[1].strip())
+        elif line_stripped.startswith("frame:"):
+            data_start = i
+            break
+        else:
+            if ":" in line_stripped:
+                k, v = line_stripped.split(":", 1)
+                header_fields[k.strip()] = v.strip()
+    
+    # 快速定位所有数据行：跳过 frame header 和 column header
+    # 数据行格式：每行 7 个字段（n x y z vx vy vz），固定宽度列
+    data_text = []
+    i = data_start
+    n_frames = 0
+    while i < len(lines):
+        line = lines[i].strip()
+        if line.startswith("frame:"):
+            n_frames += 1
+            i += 2  # 跳过 "frame: N" 和列头行
+            continue
+        if line:
+            data_text.append(line)
+        i += 1
+    
+    if n_frames == 0 or not data_text:
+        raise ValueError(f"{path} 中没有有效帧数据")
+    
+    # 用 numpy 批量解析所有数据行（远比逐行 split+float 快）
+    data_array = np.genfromtxt(data_text, dtype=np.float64)
+    # data_array shape: (n_frames * n_atoms, 7) — 列: n, x, y, z, vx, vy, vz
+    
+    n_atoms = n_total_particles
+    atoms_per_frame = len(data_text) // n_frames
+    
+    # 提取原子ID（第一帧即可）
+    atom_ids = data_array[0:n_atoms, 0].astype(np.int64)
+    
+    # 重塑为 (n_frames, n_atoms, 6) — 去掉第0列(原子ID)
+    all_data = data_array[:, 1:].reshape(n_frames, n_atoms, 6)
+    
+    return {
+        "frames_x":  all_data[:, :, 0],
+        "frames_y":  all_data[:, :, 1],
+        "frames_z":  all_data[:, :, 2],
+        "frames_vx": all_data[:, :, 3],
+        "frames_vy": all_data[:, :, 4],
+        "frames_vz": all_data[:, :, 5],
+        "atom_ids": atom_ids,
+        "n_total_frames": n_total_frames,
+        "n_total_particles": n_total_particles,
+        "header_fields": header_fields,
+    }
+
+
+def save_display_npz(path, frames_x, frames_y, frames_z,
+                     frames_vx, frames_vy, frames_vz,
+                     atom_ids, header_fields=None):
+    """Write display.npz binary format alongside display.txt.
+
+    Stores all frame arrays in compressed NumPy format.  Metadata is
+    serialised as a JSON byte-string stored in the 'meta' array so the
+    file stays self-contained.
+
+    Args:
+        path: output path ending in '.npz'
+        frames_x/y/z/vx/vy/vz: (n_frames, n_atoms) float64 arrays
+        atom_ids: (n_atoms,) int array
+        header_fields: dict of metadata key-value strings (same as for
+                       save_display_txt)
+    """
+    meta_bytes = np.bytes_(json.dumps(header_fields or {}, ensure_ascii=False))
+    stem = path[:-4] if path.endswith(".npz") else path
+    np.savez_compressed(
+        stem,
+        frames_x=frames_x,
+        frames_y=frames_y,
+        frames_z=frames_z,
+        frames_vx=frames_vx,
+        frames_vy=frames_vy,
+        frames_vz=frames_vz,
+        atom_ids=atom_ids,
+        meta=np.array(meta_bytes),
+    )
+    return path
+
+
+def load_display_npz(path):
+    """Read display.npz and return the same dict as load_display_txt."""
+    data = np.load(path, allow_pickle=False)
+    header_fields = json.loads(data["meta"].item().decode("utf-8"))
+    frames_x = data["frames_x"]
+    n_total_frames = int(header_fields.get("number_of_frames",
+                                           frames_x.shape[0]))
+    n_total_particles = frames_x.shape[1]
+    return {
+        "frames_x":  frames_x,
+        "frames_y":  data["frames_y"],
+        "frames_z":  data["frames_z"],
+        "frames_vx": data["frames_vx"],
+        "frames_vy": data["frames_vy"],
+        "frames_vz": data["frames_vz"],
+        "atom_ids":  data["atom_ids"],
+        "n_total_frames": n_total_frames,
+        "n_total_particles": n_total_particles,
+        "header_fields": header_fields,
+    }
+
+
+def load_display(output_dir):
+    """Load display data, preferring binary .npz over text .txt."""
+    npz_path = os.path.join(output_dir, "display.npz")
+    txt_path = os.path.join(output_dir, "display.txt")
+    if os.path.exists(npz_path):
+        return load_display_npz(npz_path)
+    if os.path.exists(txt_path):
+        return load_display_txt(txt_path)
+    raise FileNotFoundError(
+        f"找不到 display.npz 或 display.txt in {output_dir}")
+
+
 def get_output_dir(base_dir=None):
     """Return the output directory used for generated artifacts."""
     override = os.environ.get("DYNAMICS_OUTPUT_DIR")
@@ -445,9 +699,8 @@ def apply_driving_force(x, y, z, vx, vy, vz, t, step, drivers, dt):
             period_steps = None  # 全程驱动
 
         # 当前驱动力下的位置 / 速度
-        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"])
-        vel_drive = -d["amp"] * 2.0 * np.pi * d["freq"] * np.sin(2.0 * np.pi * d["freq"] * t_vec + d["phi"])
+        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"])
 
         x[idx] = pos_drive[0]
         y[idx] = pos_drive[1]
@@ -510,6 +763,8 @@ def run_from_config(config, out_dir=None):
     global X_MIN, X_MAX, Y_MIN, Y_MAX, Z_MIN, Z_MAX
     global ball_radius, ball_color_r, ball_color_g, ball_color_b
     global box_color_r, box_color_g, box_color_b
+    global use_marker, camera_keyframes_raw, camera_distance, camera_elevation, camera_azimuth
+    global FIXED_MASK_X, FIXED_MASK_Y, FIXED_MASK_Z
     global warmup_steps, sample_start, sample_end
     global GRAVITY_FIELD, GRAVITY_INTERACTION, ELASTIC_FORCE, DAMPING_FORCE, GRAVITY_STRENGTH
     global DRIVING_FORCE, DRIVER_DATA
@@ -526,6 +781,9 @@ def run_from_config(config, out_dir=None):
         coord_path = os.path.join(out_dir, coord_path)
     (ATOM_IDS, ATOM_MASSES, ATOM_RADII, ATOM_POSITIONS,
      ATOM_VELOCITIES, ATOM_FIXED) = load_coord_file(coord_path)
+    FIXED_MASK_X = ATOM_FIXED[:, 0] != 0
+    FIXED_MASK_Y = ATOM_FIXED[:, 1] != 0
+    FIXED_MASK_Z = ATOM_FIXED[:, 2] != 0
     BOND_CONNECTION_FILE = str(config.get("connection_file", os.path.join("input", "connection.txt")))
     BOND_PARAMETER_FILE = str(config.get("bond_file", os.path.join("input", "bond.txt")))
     connection_path = BOND_CONNECTION_FILE
@@ -570,6 +828,10 @@ def run_from_config(config, out_dir=None):
     box_color_r = float(config.get("box_color_r", 0.8))
     box_color_g = float(config.get("box_color_g", 0.8))
     box_color_b = float(config.get("box_color_b", 0.85))
+    use_marker = int(config.get("use_marker", 0))
+    camera_distance = float(config.get("camera_distance", 40.0))
+    camera_elevation = float(config.get("camera_elevation", 0.0))
+    camera_azimuth = float(config.get("camera_azimuth", 0.0))
 
     # 力开关
     global GRAVITY_FIELD, GRAVITY_INTERACTION, ELASTIC_FORCE, DAMPING_FORCE, GRAVITY_STRENGTH
@@ -590,13 +852,25 @@ def run_from_config(config, out_dir=None):
             driver_path = os.path.join(out_dir, driver_rel)
         DRIVER_DATA = load_driver_file(driver_path, ATOM_IDS)
 
+    # 加载运动相机关键帧
+    camera_keyframes_raw = ""
+    move_camera = int(config.get("move_camera", 0))
+    camera_keyframes_raw = ""
+    if move_camera:
+        cam_rel = str(config.get("move_camera_file", os.path.join("input", "move_camera.txt")))
+        cam_path = cam_rel
+        if out_dir is not None and not os.path.isabs(cam_rel):
+            cam_path = os.path.join(out_dir, cam_rel)
+        camera_keyframes_raw = _load_camera_motion(cam_path)
+
     print(f"[compute] 使用算法: {METHOD}")
     print(f"[compute] 已加载成键信息: {len(BOND_PAIRS)} 条键")
     if config.get("_skip_run", False):
         return None, None, None, None, None, None
     t_start = time.time()
     t_start_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
-    traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = run_simulation()
+    traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = run_simulation(
+        save_trajectory=int(config.get("save_trajectory", 0)))
     t_end = time.time()
     t_end_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
     elapsed = t_end - t_start
@@ -698,6 +972,20 @@ def run_engine(engine, input_dir, output_dir, config):
         "damping_force": int(config.get("damping_force", 0)),
         "gravity_strength": float(config.get("gravity_strength", 1.0)),
         "driving_force": int(config.get("driving_force", 0)),
+        "save_trajectory": int(config.get("save_trajectory", 0)),
+        # 渲染参数（用于 display.txt header）
+        "alpha": config.get("alpha", 0.2),
+        "ball_radius": float(config.get("ball_radius", 0.5)),
+        "ball_color": [float(config.get("ball_color_r", 0.9)),
+                       float(config.get("ball_color_g", 0.2)),
+                       float(config.get("ball_color_b", 0.2))],
+        "box_color": [float(config.get("box_color_r", 0.8)),
+                      float(config.get("box_color_g", 0.8)),
+                      float(config.get("box_color_b", 0.85))],
+        "use_marker": int(config.get("use_marker", 0)),
+        "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)),
     }
     param_path = os.path.join(script_dir, "engines", engine, "param.json")
     os.makedirs(os.path.dirname(param_path), exist_ok=True)
@@ -711,31 +999,55 @@ def run_engine(engine, input_dir, output_dir, config):
     print(f"[compute]   input:  {input_dir}")
     print(f"[compute]   output: {output_dir}")
 
-    # ── 预校准：跑 NT=1000 步测速 ────────────────────────────────
+    # ── 预校准：跑少量步测速（结果缓存，避免每次重跑）────────────
     total_steps = int(config["NT"]) - int(config.get("warmup_steps", 0))
     _calib_nt = min(1000, max(100, total_steps // 10))
-    _calib_param = dict(param_json)
-    _calib_param["NT"] = _calib_nt
-    _calib_path = os.path.join(script_dir, "engines", engine, "_calib.json")
-    with open(_calib_path, "w", encoding="utf-8") as _cf:
-        json.dump(_calib_param, _cf, indent=2)
-    _calib_outdir = os.path.join(script_dir, "engines", engine, "_calib_out")
-    os.makedirs(_calib_outdir, exist_ok=True)
-    _ct0 = time.time()
-    subprocess.run(
-        [engine_path, os.path.abspath(input_dir), _calib_outdir, _calib_path],
-        capture_output=True, timeout=60)
-    # 清理校准临时文件
-    for _f in os.listdir(_calib_outdir):
-        try: os.remove(os.path.join(_calib_outdir, _f))
+    n_atoms_calib = len(ATOM_IDS) if ATOM_IDS is not None else 0
+
+    # 尝试读取缓存；当 n_atoms 相同且 NT 在 50% 范围内时视为有效
+    _cache_path = os.path.join(script_dir, "engines", engine, "_calib_cache.json")
+    _step_time = None
+    try:
+        with open(_cache_path, encoding="utf-8") as _cf:
+            _cache = json.load(_cf)
+        if (_cache.get("n_atoms") == n_atoms_calib and
+                abs(_cache.get("nt", 0) - total_steps) / max(total_steps, 1) < 0.5):
+            _step_time = float(_cache["step_time"])
+            print(f"[compute] 使用校准缓存: {_step_time*1e6:.2f} μs/步")
+    except Exception:
+        pass
+
+    if _step_time is None:
+        _calib_param = dict(param_json)
+        _calib_param["NT"] = _calib_nt
+        _calib_path = os.path.join(script_dir, "engines", engine, "_calib.json")
+        with open(_calib_path, "w", encoding="utf-8") as _cf:
+            json.dump(_calib_param, _cf, indent=2)
+        _calib_outdir = os.path.join(script_dir, "engines", engine, "_calib_out")
+        os.makedirs(_calib_outdir, exist_ok=True)
+        _ct0 = time.time()
+        subprocess.run(
+            [engine_path, os.path.abspath(input_dir), _calib_outdir, _calib_path],
+            capture_output=True, timeout=60)
+        # 清理校准临时文件
+        for _f in os.listdir(_calib_outdir):
+            try: os.remove(os.path.join(_calib_outdir, _f))
+            except OSError: pass
+        try: os.rmdir(_calib_outdir)
         except OSError: pass
-    try: os.rmdir(_calib_outdir)
-    except OSError: pass
-    os.remove(_calib_path)
-    _calib_elapsed = max(time.time() - _ct0, 0.001)
-    _overhead = _calib_elapsed * 0.15
-    _step_time = max(_calib_elapsed - _overhead, 0.0001) / _calib_nt
-    _est_total = max(_calib_elapsed, _overhead + _step_time * total_steps)
+        os.remove(_calib_path)
+        _calib_elapsed = max(time.time() - _ct0, 0.001)
+        _overhead = _calib_elapsed * 0.15
+        _step_time = max(_calib_elapsed - _overhead, 0.0001) / _calib_nt
+        # 保存缓存
+        try:
+            with open(_cache_path, "w", encoding="utf-8") as _cf:
+                json.dump({"n_atoms": n_atoms_calib, "nt": total_steps,
+                           "step_time": _step_time}, _cf)
+        except OSError:
+            pass
+
+    _est_total = _step_time * total_steps
 
     t_start = time.time()
     t_start_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
@@ -760,17 +1072,28 @@ def run_engine(engine, input_dir, output_dir, config):
                 _line = _line.strip()
                 if _line:
                     _engine_lines.append(_line)
+                    # 读取外部引擎真实进度：格式 "[xxx-engine] progress: N/total"
+                    _prog_match = re.search(r'progress:\s*(\d+)/(\d+)', _line)
+                    if _pbar is not None and _prog_match:
+                        _prog_done = int(_prog_match.group(1))
+                        _prog_total = int(_prog_match.group(2))
+                        if _prog_total > 0:
+                            _pbar.n = min(_prog_done, total_steps)
+                            _pbar.refresh()
             if _p.poll() is not None:
                 if _p.stdout:
                     for _r in _p.stdout:
                         _r = _r.strip()
                         if _r:
                             _engine_lines.append(_r)
+                            # 读取残留在管道中的进度消息，避免 20%→100% 跳变
+                            _prog_match = re.search(r'progress:\s*(\d+)/(\d+)', _r)
+                            if _pbar is not None and _prog_match:
+                                _p_done = min(int(_prog_match.group(1)), total_steps)
+                                _pbar.n = max(_pbar.n, _p_done)
+                    if _pbar is not None: _pbar.refresh()
                 break
-            if _pbar is not None and _est_total > 0:
-                _pbar.n = int(min((time.time() - t_start) / _est_total, 0.99) * total_steps)
-                _pbar.refresh()
-            time.sleep(0.2)
+            time.sleep(0.05)
     finally:
         if _pbar is not None:
             _pbar.n = total_steps
@@ -801,7 +1124,32 @@ def run_engine(engine, input_dir, output_dir, config):
         except OSError:
             pass
 
-    # 加载输出的 trajectory.txt
+    # 加载输出的 trajectory.txt / display.txt
+    save_traj = int(config.get("save_trajectory", 0))
+    if not save_traj:
+        # save_trajectory=0：引擎只写 display.txt
+        disp_path = os.path.join(os.path.abspath(output_dir), "display.txt")
+        if not os.path.exists(disp_path):
+            raise FileNotFoundError(f"引擎未生成 display.txt: {disp_path}")
+        print(f"[compute] 引擎已生成 {disp_path}")
+        # 将 display.txt 转换为 display.npz（加快后续 draw.py / plot_wave.py 加载）
+        _npz_path = disp_path.replace("display.txt", "display.npz")
+        try:
+            _d = load_display_txt(disp_path)
+            save_display_npz(
+                _npz_path,
+                _d["frames_x"], _d["frames_y"], _d["frames_z"],
+                _d["frames_vx"], _d["frames_vy"], _d["frames_vz"],
+                _d["atom_ids"],
+                header_fields={**_d["header_fields"],
+                                "number_of_frames": str(_d["n_total_frames"]),
+                                "number_of_particles": str(_d["n_total_particles"])})
+            print(f"[compute] display.npz 已生成: {_npz_path}")
+        except Exception as _e:
+            print(f"[compute] 警告: display.npz 生成失败({_e})，将使用 display.txt")
+        return None, None, None, None, None, None
+    
+    # save_trajectory=1：加载完整 trajectory.txt
     traj_path = os.path.join(os.path.abspath(output_dir), "trajectory.txt")
     if not os.path.exists(traj_path):
         raise FileNotFoundError(f"引擎未生成 trajectory.txt: {traj_path}")
@@ -1020,30 +1368,28 @@ def compute_force(x, y, z, vx, vy, vz, m, g, b):
         fz -= b[2] * vz
 
     if ELASTIC_FORCE and BOND_PAIRS is not None and len(BOND_PAIRS) > 0:
-        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
+        idx_1 = BOND_PAIRS[:, 0]
+        idx_2 = BOND_PAIRS[:, 1]
+        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)
+        valid = dist > 1e-12
+        if np.any(valid):
+            force_scale = np.zeros_like(dist)
+            stretch = dist[valid] - BOND_REST_LENGTHS[valid]
+            force_scale[valid] = BOND_STIFFNESS[valid] * stretch / dist[valid]
 
-            stretch = dist - BOND_REST_LENGTHS[bond_idx]
-            force_mag = BOND_STIFFNESS[bond_idx] * stretch
-            ux = dx / dist
-            uy = dy / dist
-            uz = dz / dist
+            fx_bond = force_scale * dx
+            fy_bond = force_scale * dy
+            fz_bond = force_scale * dz
 
-            fx_bond = force_mag * ux
-            fy_bond = force_mag * uy
-            fz_bond = force_mag * uz
-
-            fx[idx_1] += fx_bond
-            fy[idx_1] += fy_bond
-            fz[idx_1] += fz_bond
-            fx[idx_2] -= fx_bond
-            fy[idx_2] -= fy_bond
-            fz[idx_2] -= fz_bond
+            np.add.at(fx, idx_1, fx_bond)
+            np.add.at(fx, idx_2, -fx_bond)
+            np.add.at(fy, idx_1, fy_bond)
+            np.add.at(fy, idx_2, -fy_bond)
+            np.add.at(fz, idx_1, fz_bond)
+            np.add.at(fz, idx_2, -fz_bond)
 
     if GRAVITY_INTERACTION:
         n = len(m)
@@ -1073,6 +1419,50 @@ def compute_acceleration(x, y, z, vx, vy, vz, m, g, b):
     return fx / m, fy / m, fz / m
 
 
+def compute_accel_conservative(x, y, z, m, g):
+    """保守力加速度（弹簧 + 重力 + 原子间引力），不含阻尼。
+    供真蛙跳法使用：阻尼项由 leapfrog_staggered_step 半隐式处理。
+    """
+    _v0 = np.zeros_like(x)
+    _b0 = np.zeros(3)
+    fx, fy, fz = compute_force(x, y, z, _v0, _v0, _v0, m, g, _b0)
+    return fx / m, fy / m, fz / m
+
+
+def leapfrog_staggered_step(x, y, z, vx_h, vy_h, vz_h, dt, m, g, b):
+    """真蛙跳一步：x(t), v(t-dt/2) → x(t+dt), v(t+dt/2)
+
+    - 无阻尼 (DAMPING_FORCE=False 或 b≈0)：
+        纯保守蛙跳，每步仅 1 次力计算，辛积分器。
+        v(t+dt/2) = v(t-dt/2) + a_c(t)·dt
+
+    - 有阻尼 (DAMPING_FORCE=True 且 b≠0)：
+        半隐式处理阻尼，仍只 1 次力计算，对任意阻尼无条件稳定。
+        利用 v(t) ≈ [v(t-dt/2) + v(t+dt/2)] / 2 解析求解：
+        v(t+dt/2) = [v(t-dt/2)·(1 - α) + a_c(t)·dt] / (1 + α)
+        其中 α = b·dt / (2·m)，每个方向独立。
+    """
+    ax_c, ay_c, az_c = compute_accel_conservative(x, y, z, m, g)
+
+    has_damping = DAMPING_FORCE and np.any(b != 0)
+    if has_damping:
+        alpha_x = b[0] * dt / (2.0 * m)
+        alpha_y = b[1] * dt / (2.0 * m)
+        alpha_z = b[2] * dt / (2.0 * m)
+        vx_h_new = (vx_h * (1.0 - alpha_x) + ax_c * dt) / (1.0 + alpha_x)
+        vy_h_new = (vy_h * (1.0 - alpha_y) + ay_c * dt) / (1.0 + alpha_y)
+        vz_h_new = (vz_h * (1.0 - alpha_z) + az_c * dt) / (1.0 + alpha_z)
+    else:
+        vx_h_new = vx_h + ax_c * dt
+        vy_h_new = vy_h + ay_c * dt
+        vz_h_new = vz_h + az_c * dt
+
+    x_new = x + vx_h_new * dt
+    y_new = y + vy_h_new * dt
+    z_new = z + vz_h_new * dt
+    return x_new, y_new, z_new, vx_h_new, vy_h_new, vz_h_new
+
+
 def Explicit_Euler_Method(x, y, z, vx, vy, vz, dt, m, g, b):
     ax, ay, az = compute_acceleration(x, y, z, vx, vy, vz, m, g, b)
     x = x + vx * dt
@@ -1176,15 +1566,16 @@ def apply_motion_update(x, y, z, vx, vy, vz, dt, m, g, b):
 
 def apply_fixed_constraints(x, y, z, vx, vy, vz):
     """Keep fixed degrees of freedom at their initial coordinate with zero speed."""
-    fixed = ATOM_FIXED != 0
-    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],
-        velocities[:, 0], velocities[:, 1], velocities[:, 2],
-    )
+    if FIXED_MASK_X is not None and np.any(FIXED_MASK_X):
+        x[FIXED_MASK_X] = ATOM_POSITIONS[FIXED_MASK_X, 0]
+        vx[FIXED_MASK_X] = 0.0
+    if FIXED_MASK_Y is not None and np.any(FIXED_MASK_Y):
+        y[FIXED_MASK_Y] = ATOM_POSITIONS[FIXED_MASK_Y, 1]
+        vy[FIXED_MASK_Y] = 0.0
+    if FIXED_MASK_Z is not None and 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
 
 def wrap_position(x, y, z):
     """边界回绕（ dynamics 模式）。"""
@@ -1201,12 +1592,13 @@ def wrap_position(x, y, z):
 # 主计算流程
 # ===========================================================================
 
-def run_simulation():
-    """计算 NT 步轨迹，返回位置/速度数组。
+def run_simulation(save_trajectory=0):
+    """计算 NT 步轨迹，直接抽帧并保存 display.txt。
 
     步骤控制：
     - warmup_steps: 预热步数，跳过不记录（用于稳定初始状态）
-    - 实际记录步数 = NT - warmup_steps
+    - 按 NSTEP 抽帧保存到 display.txt（新格式）
+    - save_trajectory=1 时额外保存完整 trajectory.txt（JSON 旧格式）
     """
     # 预热阶段
     x = ATOM_POSITIONS[:, 0].copy()
@@ -1216,15 +1608,27 @@ def run_simulation():
     vy = ATOM_VELOCITIES[:, 1].copy()
     vz = ATOM_VELOCITIES[:, 2].copy()
     x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz)
-    # 初始时刻驱动力（t=0 时原子 1 的位置由驱动力决定而非 coord.txt）
     x, y, z, vx, vy, vz = apply_driving_force(x, y, z, vx, vy, vz, 0.0, 0, DRIVER_DATA, DT)
 
+    # 真蛙跳初始化：从 v(0) 反推 v(-dt/2)，使后续每步只需 1 次力计算。
+    # 其他方法（euler/midpoint）vx/vy/vz 始终存整步速度，不受影响。
+    if METHOD == "leapfrog":
+        _ax0, _ay0, _az0 = compute_accel_conservative(x, y, z, ATOM_MASSES, G)
+        vx = vx - 0.5 * _ax0 * DT
+        vy = vy - 0.5 * _ay0 * DT
+        vz = vz - 0.5 * _az0 * DT
+        # 从此 vx/vy/vz 存 v(t-dt/2)（蛙跳半步速度）
+
     if warmup_steps is not None and warmup_steps > 0:
         print(f"[compute] 预热阶段: 前 {warmup_steps} 步不记录")
         for step in trange(warmup_steps, desc="[compute] 预热"):
             t = (step + 1) * DT
             x, y, z, vx, vy, vz = apply_driving_force(x, y, z, vx, vy, vz, t, step, DRIVER_DATA, DT)
-            x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
+            if METHOD == "leapfrog":
+                x, y, z, vx, vy, vz = leapfrog_staggered_step(
+                    x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
+            else:
+                x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
             x, y, z = wrap_position(x, y, z)
             x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz)
         print(
@@ -1232,33 +1636,153 @@ def run_simulation():
             f"({x[PLOT_ATOM_ROW]:.4f}, {y[PLOT_ATOM_ROW]:.4f}, {z[PLOT_ATOM_ROW]:.4f})"
         )
 
-    # 记录阶段
+    # 记录阶段 - 按 NSTEP 抽帧保存
     record_steps = NT - (warmup_steps or 0)
     n_atoms = len(ATOM_IDS)
-    traj_x  = np.zeros((record_steps, n_atoms), dtype=np.float64)
-    traj_y  = np.zeros((record_steps, n_atoms), dtype=np.float64)
-    traj_z  = np.zeros((record_steps, n_atoms), dtype=np.float64)
-    traj_vx = np.zeros((record_steps, n_atoms), dtype=np.float64)
-    traj_vy = np.zeros((record_steps, n_atoms), dtype=np.float64)
-    traj_vz = np.zeros((record_steps, n_atoms), dtype=np.float64)
+    n_frames = (record_steps + NSTEP - 1) // NSTEP
+    # 按 NSTEP 抽帧的临时缓冲区（远小于全量轨迹）
+    sampled_x  = np.zeros((n_frames, n_atoms), dtype=np.float64)
+    sampled_y  = np.zeros((n_frames, n_atoms), dtype=np.float64)
+    sampled_z  = np.zeros((n_frames, n_atoms), dtype=np.float64)
+    sampled_vx = np.zeros((n_frames, n_atoms), dtype=np.float64)
+    sampled_vy = np.zeros((n_frames, n_atoms), dtype=np.float64)
+    sampled_vz = np.zeros((n_frames, n_atoms), dtype=np.float64)
+
+    # 如果 save_trajectory，准备完整轨迹缓冲区
+    if save_trajectory:
+        traj_x  = np.zeros((record_steps, n_atoms), dtype=np.float64)
+        traj_y  = np.zeros((record_steps, n_atoms), dtype=np.float64)
+        traj_z  = np.zeros((record_steps, n_atoms), dtype=np.float64)
+        traj_vx = np.zeros((record_steps, n_atoms), dtype=np.float64)
+        traj_vy = np.zeros((record_steps, n_atoms), dtype=np.float64)
+        traj_vz = np.zeros((record_steps, n_atoms), dtype=np.float64)
 
     for step in trange(record_steps, desc="[compute] 计算中"):
         t = (step + (warmup_steps or 0)) * DT
-        # 先施加驱动力（受驱原子的位置覆盖初始/固定约束，为弹簧力提供正确参考）
         x, y, z, vx, vy, vz = apply_driving_force(x, y, z, vx, vy, vz, t, step, DRIVER_DATA, DT)
 
-        traj_x[step]  = x
-        traj_y[step]  = y
-        traj_z[step]  = z
-        traj_vx[step] = vx
-        traj_vy[step] = vy
-        traj_vz[step] = vz
+        if save_trajectory:
+            traj_x[step]  = x
+            traj_y[step]  = y
+            traj_z[step]  = z
+            traj_vx[step] = vx
+            traj_vy[step] = vy
+            traj_vz[step] = vz
 
-        x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
+        # 抽帧：NSTEP 间隔保存
+        if step % NSTEP == 0:
+            fi = step // NSTEP
+            sampled_x[fi]  = x
+            sampled_y[fi]  = y
+            sampled_z[fi]  = z
+            sampled_vx[fi] = vx
+            sampled_vy[fi] = vy
+            sampled_vz[fi] = vz
+        if METHOD == "leapfrog":
+            x, y, z, vx, vy, vz = leapfrog_staggered_step(
+                x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
+        else:
+            x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
         x, y, z = wrap_position(x, y, z)
         x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz)
 
-    return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz
+    # 写入 display.txt（新格式）
+    output_dir = get_output_dir()
+    disp_path = os.path.join(output_dir, "display.txt")
+    n_frames_actual = (record_steps + NSTEP - 1) // NSTEP
+    save_display_txt(
+        disp_path,
+        sampled_x[:n_frames_actual], sampled_y[:n_frames_actual], sampled_z[:n_frames_actual],
+        sampled_vx[:n_frames_actual], sampled_vy[:n_frames_actual], sampled_vz[:n_frames_actual],
+        np.array(ATOM_IDS), record_steps, n_atoms,
+        header_fields={"DT": str(DT), "NSTEP": str(NSTEP), "method": str(METHOD),
+                       "warmup_steps": str(warmup_steps or 0),
+                       "dynamic_steps": str(record_steps),
+                       "T_total": str(NT * DT),
+                       "X_MAX": str(X_MAX), "X_MIN": str(X_MIN),
+                       "Y_MAX": str(Y_MAX), "Y_MIN": str(Y_MIN),
+                       "Z_MAX": str(Z_MAX), "Z_MIN": str(Z_MIN),
+                       "ball_radius": str(ball_radius),
+                       "ball_color_r": str(ball_color_r),
+                       "ball_color_g": str(ball_color_g),
+                       "ball_color_b": str(ball_color_b),
+                       "box_color_r": str(box_color_r),
+                       "box_color_g": str(box_color_g),
+                       "box_color_b": str(box_color_b),
+                       "gravity_field": str(GRAVITY_FIELD),
+                       "gravity_interaction": str(GRAVITY_INTERACTION),
+                       "elastic_force": str(ELASTIC_FORCE),
+                       "damping_force": str(DAMPING_FORCE),
+                       "gravity_strength": str(GRAVITY_STRENGTH),
+                       "driving_force": str(DRIVING_FORCE),
+                       "use_marker": str(use_marker),
+                       "alpha": ",".join(str(a) for a in (alpha if isinstance(alpha, list) else [alpha])),
+                       "atom_masses": json.dumps([float(v) for v in ATOM_MASSES]),
+                       "atom_positions": json.dumps(ATOM_POSITIONS.tolist()),
+                       "bond_pairs": json.dumps(BOND_PAIRS.tolist() if BOND_PAIRS is not None else []),
+                       "bond_stiffness": json.dumps(BOND_STIFFNESS.tolist() if BOND_STIFFNESS is not None else []),
+                       "bond_rest_lengths": json.dumps(BOND_REST_LENGTHS.tolist() if BOND_REST_LENGTHS is not None else []),
+                       "G": json.dumps(G.tolist() if G is not None else [0.0, 0.0, 0.0]),
+                       "atom_radii": ",".join(str(r) for r in ATOM_RADII),
+                       "camera_distance": str(camera_distance),
+                       "camera_elevation": str(camera_elevation),
+                       "camera_azimuth": str(camera_azimuth),
+                       "camera_keyframes": str(camera_keyframes_raw)}
+    )
+    print(f"[compute] display.txt 已保存至: {disp_path} ({n_frames_actual} 帧)")
+
+    # 同时保存二进制 display.npz（I/O 速度提升 5-10x）
+    _npz_path = os.path.join(output_dir, "display.npz")
+    _hdr = {"DT": str(DT), "NSTEP": str(NSTEP), "method": str(METHOD),
+            "warmup_steps": str(warmup_steps or 0),
+            "dynamic_steps": str(record_steps),
+            "T_total": str(NT * DT),
+            "X_MAX": str(X_MAX), "X_MIN": str(X_MIN),
+            "Y_MAX": str(Y_MAX), "Y_MIN": str(Y_MIN),
+            "Z_MAX": str(Z_MAX), "Z_MIN": str(Z_MIN),
+            "ball_radius": str(ball_radius),
+            "ball_color_r": str(ball_color_r),
+            "ball_color_g": str(ball_color_g),
+            "ball_color_b": str(ball_color_b),
+            "box_color_r": str(box_color_r),
+            "box_color_g": str(box_color_g),
+            "box_color_b": str(box_color_b),
+            "gravity_field": str(GRAVITY_FIELD),
+            "gravity_interaction": str(GRAVITY_INTERACTION),
+            "elastic_force": str(ELASTIC_FORCE),
+            "damping_force": str(DAMPING_FORCE),
+            "gravity_strength": str(GRAVITY_STRENGTH),
+            "driving_force": str(DRIVING_FORCE),
+            "use_marker": str(use_marker),
+            "alpha": ",".join(str(a) for a in (alpha if isinstance(alpha, list) else [alpha])),
+            "atom_masses": json.dumps([float(v) for v in ATOM_MASSES]),
+            "atom_positions": json.dumps(ATOM_POSITIONS.tolist()),
+            "bond_pairs": json.dumps(BOND_PAIRS.tolist() if BOND_PAIRS is not None else []),
+            "bond_stiffness": json.dumps(BOND_STIFFNESS.tolist() if BOND_STIFFNESS is not None else []),
+            "bond_rest_lengths": json.dumps(BOND_REST_LENGTHS.tolist() if BOND_REST_LENGTHS is not None else []),
+            "G": json.dumps(G.tolist() if G is not None else [0.0, 0.0, 0.0]),
+            "atom_radii": ",".join(str(r) for r in ATOM_RADII),
+            "camera_distance": str(camera_distance),
+            "camera_elevation": str(camera_elevation),
+            "camera_azimuth": str(camera_azimuth),
+            "camera_keyframes": str(camera_keyframes_raw),
+            "number_of_frames": str(record_steps),
+            "number_of_particles": str(n_atoms)}
+    save_display_npz(_npz_path,
+                     sampled_x[:n_frames_actual], sampled_y[:n_frames_actual],
+                     sampled_z[:n_frames_actual],
+                     sampled_vx[:n_frames_actual], sampled_vy[:n_frames_actual],
+                     sampled_vz[:n_frames_actual],
+                     np.array(ATOM_IDS), header_fields=_hdr)
+    print(f"[compute] display.npz 已保存至: {_npz_path}")
+
+    # 可选：保存完整 trajectory.txt
+    if save_trajectory:
+        save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz)
+        print(f"[compute] trajectory.txt 已保存（完整轨迹）")
+
+    return sampled_x[:n_frames_actual], sampled_y[:n_frames_actual], sampled_z[:n_frames_actual], \
+           sampled_vx[:n_frames_actual], sampled_vy[:n_frames_actual], sampled_vz[:n_frames_actual]
 
 
 def save_trajectory_table_txt(txt_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT):
@@ -1301,10 +1825,13 @@ def main():
     output_dir = get_output_dir(script_dir)
 
     print(f"[compute] 开始计算  NT={NT}  DT={DT} ")
-    traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = run_simulation()
+    traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = run_simulation(save_trajectory=0)
     print(f"[compute] 计算完成，共 {NT} 步")
+    print(f"[compute] display.txt 已在 run_simulation 中保存")
 
-    save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, script_dir)
+    # 如果需要完整轨迹，以上传 save_trajectory=1 重新运行
+    # 以下旧函数保留兼容但不再自动调用
+    # save_trajectory_txt(...)
 
     # 同时保存为逐行表格，便于直接查看
     txt_path = os.path.join(output_dir, "trajectory_table.txt")

draw.py

@@ -1,9 +1,8 @@
 """VisPy 演示：加载预计算轨迹数据，驱动小球运动动画。
 
 计算与显示完全分离：
-  1. 先运行 compute.py  → 生成 output/trajectory.txt（全量 NT 步轨迹）
-  2. 再运行 sample.py  → 从 output/trajectory.txt 抽帧生成 output/display.txt
-  3. 本文件加载 output/display.txt，按帧播放动画
+  1. 运行 run_dynamics.py  → 生成 output/display.txt（新格式，直接抽帧）
+  2. 本文件加载 output/display.txt，按帧播放动画
 
 用法：
     python draw.py                                          # 使用 dynamics 根目录下的 output/
@@ -11,6 +10,7 @@
 """
 
 import numpy as np
+import json
 import os
 import sys
 from vispy import app, scene
@@ -30,88 +30,99 @@ else:
     output_dir = compute.get_output_dir(script_dir)
 os.environ["DYNAMICS_OUTPUT_DIR"] = output_dir
 disp_path = os.path.join(output_dir, "display.txt")
+npz_path  = os.path.join(output_dir, "display.npz")
 
-if not os.path.exists(disp_path):
+if not os.path.exists(npz_path) and not os.path.exists(disp_path):
     raise FileNotFoundError(
-        f"找不到 display.txt！\n"
-        f"期望路径: {disp_path}\n"
+        f"找不到 display.npz 或 display.txt！\n"
+        f"期望路径: {output_dir}\n"
         f"请先运行 compute.py 计算轨迹，再运行 sample.py 生成显示数组。\n"
         f"用法: python draw.py [案例输出目录]"
     )
 
-disp_data = compute.load_text_data(disp_path)
+# 优先读二进制 npz（加载速度约快 5-10x）
+if os.path.exists(npz_path):
+    disp_data = compute.load_display_npz(npz_path)
+else:
+    disp_data = compute.load_display_txt(disp_path)
+h = disp_data["header_fields"]
 
-# 单原子数据（plot_atom：用于信息显示）
-DISP_X    = disp_data["disp_x"]
-DISP_Y    = disp_data["disp_y"]
-DISP_Z    = disp_data["disp_z"]
-DISP_VX   = disp_data["disp_vx"]
-DISP_VY   = disp_data["disp_vy"]
-DISP_VZ   = disp_data["disp_vz"]
+# 全原子帧数据
+DISP_ALL_X  = disp_data["frames_x"]       # (n_frames, n_atoms)
+DISP_ALL_Y  = disp_data["frames_y"]
+DISP_ALL_Z  = disp_data["frames_z"]
+DISP_ALL_VX = disp_data["frames_vx"]
+DISP_ALL_VY = disp_data["frames_vy"]
+DISP_ALL_VZ = disp_data["frames_vz"]
 
-# 全原子数据（用于多球绘制）
-DISP_ALL_X  = disp_data["disp_all_x"]     # (n_frames, n_atoms)
-DISP_ALL_Y  = disp_data["disp_all_y"]
-DISP_ALL_Z  = disp_data["disp_all_z"]
-DISP_ALL_VX = disp_data["disp_all_vx"]
-DISP_ALL_VY = disp_data["disp_all_vy"]
-DISP_ALL_VZ = disp_data["disp_all_vz"]
+# 第一个原子的轨迹（用于信息显示）
+DISP_X    = DISP_ALL_X[:, 0]
+DISP_Y    = DISP_ALL_Y[:, 0]
+DISP_Z    = DISP_ALL_Z[:, 0]
+DISP_VX   = DISP_ALL_VX[:, 0]
+DISP_VY   = DISP_ALL_VY[:, 0]
+DISP_VZ   = DISP_ALL_VZ[:, 0]
 
-DISP_T    = disp_data["disp_t"]
-DISP_STEP = disp_data["disp_step"]
-N_FRAMES  = int(disp_data["n_frames"])
-NT        = int(disp_data["NT"])
-DT        = float(disp_data["DT"])
-NSTEP     = int(disp_data["NSTEP"])
+N_FRAMES  = DISP_ALL_X.shape[0]
+NT        = int(disp_data["n_total_frames"])
+N_ATOMS   = int(disp_data["n_total_particles"])
+DT        = float(h.get("DT", 0.001))
+NSTEP     = int(h.get("NSTEP", 1))
+DISP_STEP = np.arange(N_FRAMES) * NSTEP
+DISP_T    = DISP_STEP * DT
 
 # 原子信息
-ATOM_IDS  = disp_data.get("atom_ids", np.array([1]))
-ATOM_RADII = disp_data.get("atom_radii", np.array([float(disp_data["ball_radius"])]))
-N_ATOMS   = len(ATOM_IDS)
-PLOT_ATOM_ROW = int(disp_data.get("plot_atom_row", 0))
-PLOT_ATOM_ID  = int(disp_data.get("plot_atom_id", ATOM_IDS[0]))
-BOND_PAIRS = disp_data.get("bond_pairs", [])
+ATOM_IDS  = disp_data["atom_ids"]
+# 优先使用 per-atom 半径，否则用统一的 ball_radius
+_raw_radii = h.get("atom_radii", "")
+if _raw_radii.strip():
+    ATOM_RADII = np.array([float(x) for x in _raw_radii.split(",")])
+else:
+    ATOM_RADII = np.full(N_ATOMS, float(h.get("ball_radius", 0.5)))
+PLOT_ATOM_ROW = 0
+PLOT_ATOM_ID  = int(ATOM_IDS[0])
+BOND_PAIRS = []  # display 格式不含成键信息，从原始数据加载
 
 # 渲染方式：0=Sphere(网格球体), 1=Marker(GPU点精灵)
-USE_MARKER = int(disp_data.get("use_marker", 0))
+USE_MARKER = int(h.get("use_marker", 0))
 
 if N_FRAMES <= 0:
     raise ValueError(
         "output/display.txt 中没有可播放的帧，请检查 sample_start/sample_end/NSTEP 配置。")
 
 # 保留模拟边界常量（用于场景缩放、相机等），从 output/display.txt 中读取
-X_MIN = float(disp_data["X_MIN"]); X_MAX = float(disp_data["X_MAX"])
-Y_MIN = float(disp_data["Y_MIN"]); Y_MAX = float(disp_data["Y_MAX"])
-Z_MIN = float(disp_data["Z_MIN"]); Z_MAX = float(disp_data["Z_MAX"])
-X0 = float(disp_data["X0"]); Y0 = float(disp_data["Y0"]); Z0 = float(disp_data["Z0"])
-raw_alpha = disp_data["alpha"]
-if isinstance(raw_alpha, (list, tuple, np.ndarray)):
-    alpha_list = [float(a) for a in raw_alpha]
+X_MIN = float(h.get("X_MIN", -10)); X_MAX = float(h.get("X_MAX", 10))
+Y_MIN = float(h.get("Y_MIN", -10)); Y_MAX = float(h.get("Y_MAX", 10))
+Z_MIN = float(h.get("Z_MIN", -10)); Z_MAX = float(h.get("Z_MAX", 10))
+raw_alpha = h.get("alpha", "0.2")
+try:
+    alpha_list = [float(x) for x in raw_alpha.split(",")]
     if len(alpha_list) != 6:
-        raise ValueError(f"alpha 数组长度须为 6，实际为 {len(alpha_list)}")
-else:
+        alpha_list = alpha_list * 6
+except (ValueError, AttributeError):
     alpha_list = [float(raw_alpha)] * 6
+
 # 绘图参数
-ball_radius = float(disp_data["ball_radius"])
-ball_color_r = float(disp_data["ball_color_r"])
-ball_color_g = float(disp_data["ball_color_g"])
-ball_color_b = float(disp_data["ball_color_b"])
-box_color_r = float(disp_data["box_color_r"])
-box_color_g = float(disp_data["box_color_g"])
-box_color_b = float(disp_data["box_color_b"])
+ball_radius = float(h.get("ball_radius", 0.5))
+ball_color_r = float(h.get("ball_color_r", 0.9))
+ball_color_g = float(h.get("ball_color_g", 0.2))
+ball_color_b = float(h.get("ball_color_b", 0.2))
+box_color_r = float(h.get("box_color_r", 0.8))
+box_color_g = float(h.get("box_color_g", 0.8))
+box_color_b = float(h.get("box_color_b", 0.85))
 
 
 # ===========================================================================
 # 图形界面无关的几何参数（不参与物理计算，仅用于场景外观）
 # ===========================================================================
 
-info_margin  = 36
+info_margin  = 8
 axis_length  = 10.0
 
 initial_camera = {
-    "distance": 40.0,
-    "elevation": 0,
-    "azimuth": 0,
+    "distance": float(h.get("camera_distance", 40.0)),
+    "elevation": float(h.get("camera_elevation", 0)),
+    "azimuth": float(h.get("camera_azimuth", 0)),
     "center": (0, 0, 0),
 }
 
@@ -255,7 +266,7 @@ camera_info = scene.visuals.Text(
 
 # 左上角：小球信息
 ball_info = scene.visuals.Text(
-    text="", color=(0.2, 1.0, 0.2, 1.0), font_size=28,
+    text="", color=(0.2, 1.0, 0.2, 1.0), font_size=18,
     pos=(0, 0), anchor_x="left", anchor_y="top",
     face="黑体", bold=True, parent=canvas.scene)
 
@@ -508,8 +519,9 @@ def handle_mouse_press(event):
 def _update_atom_positions(f_idx):
     """更新所有原子到第 f_idx 帧的位置。"""
     if USE_MARKER:
-        for i in range(N_ATOMS):
-            marker_pos[i] = [DISP_ALL_X[f_idx, i], DISP_ALL_Y[f_idx, i], DISP_ALL_Z[f_idx, i]]
+        marker_pos[:, 0] = DISP_ALL_X[f_idx]
+        marker_pos[:, 1] = DISP_ALL_Y[f_idx]
+        marker_pos[:, 2] = DISP_ALL_Z[f_idx]
         balls.set_data(pos=marker_pos)
     else:
         for i in range(N_ATOMS):
@@ -549,9 +561,77 @@ print(f"[draw] 渲染方式: {mode_str}")
 print(f"[draw] 绘图参数: ball_radius={ball_radius}, box_color=({box_color_r:.2f},{box_color_g:.2f},{box_color_b:.2f}), alpha={alpha_list}")
 
 
-# ===========================================================================
-# 每帧回调：仅推进帧索引，从预存数组读取位置，零物理计算
-# ===========================================================================
+# 运动相机（速度段驱动）：优先读取 move_camera.txt，其次用 display.txt header 缓存
+import re as _re
+
+def _load_move_camera_txt():
+    """直接读取 input/move_camera.txt（与 output 同级的 input 目录）。"""
+    input_dir = os.path.join(os.path.dirname(output_dir), "input")
+    cam_path = os.path.join(input_dir, "move_camera.txt")
+    if not os.path.exists(cam_path):
+        return None
+    segs = []
+    with open(cam_path, "r", encoding="utf-8") as f:
+        for line in f:
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            # 解析帧范围：支持 "all"（全程）或 "N-M"（区间）
+            if line.lower().startswith("all") or _re.match(r'^\s*all\s', line, _re.IGNORECASE):
+                start, end = 0, 10**9  # 用极大值表示全程
+            else:
+                m = _re.match(r'(\d+)\s*-\s*(\d+)', line)
+                if not m:
+                    continue
+                start, end = int(m.group(1)), int(m.group(2))
+            v, r = [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]
+            for i, axis in enumerate(['x', 'y', 'z']):
+                m2 = _re.search(r'v' + axis + r'\s*=\s*([-\d.]+)', line)
+                if m2: v[i] = float(m2.group(1))
+                m2 = _re.search(r'r' + axis + r'\s*=\s*([-\d.]+)', line)
+                if m2: r[i] = float(m2.group(1))
+            if any(v) or any(r):
+                segs.append({"start": start, "end": end, "v": v, "r": r})
+    return segs if segs else None
+
+# 先试 move_camera.txt 直读，没有则用 display.txt 缓存
+_CAM_MOTION = _load_move_camera_txt()
+if not _CAM_MOTION:
+    _CAM_MOTION = json.loads(h.get("camera_keyframes", "null")) if h.get("camera_keyframes") else None
+if _CAM_MOTION:
+    _cam_center = [0.0, 0.0, 0.0]
+    _cam_elev = initial_camera["elevation"]
+    _cam_azim = initial_camera["azimuth"]
+    _cam_dist = initial_camera["distance"]
+    src = "move_camera.txt" if _load_move_camera_txt() else "display.txt header"
+    print(f"[draw] 运动相机已启用（数据来源: {src}，{len(_CAM_MOTION)} 段）")
+
+
+def _update_motion_camera(f_idx):
+    """速度段驱动：每帧累加平移/旋转。
+    
+    时间交叠时所有段同时生效，按文件中出现的顺序依次作用。
+    矩阵操作不具有对易性，排在前面的段优先作用于相机位置。
+    """
+    if not _CAM_MOTION:
+        return
+    global _cam_center, _cam_elev, _cam_azim, _cam_dist
+    # 找当前帧所有活动的段（时间交叠=同时作用），按文件顺序依次应用
+    for seg in _CAM_MOTION:
+        if seg["start"] <= f_idx < seg["end"]:
+            _cam_center[0] += seg["v"][0]
+            _cam_center[1] += seg["v"][1]
+            _cam_center[2] += seg["v"][2]
+            _cam_elev += seg["r"][0]
+            _cam_azim += seg["r"][1]
+            # rz 预留
+
+    view.camera.center    = tuple(_cam_center)
+    view.camera.distance  = _cam_dist
+    view.camera.elevation = _cam_elev
+    view.camera.azimuth   = _cam_azim
+
+
 def update(event):
     global frame_idx
     frame_idx = (frame_idx + 1) % N_FRAMES     # 循环播放
@@ -563,6 +643,9 @@ def update(event):
     if bond_lines is not None and len(BOND_PAIRS) > 0:
         _update_bond_positions(frame_idx)
 
+    # 运动相机：速度段驱动
+    _update_motion_camera(frame_idx)
+
     # 信息面板显示 plot_atom 的数据
     x = float(DISP_X[frame_idx])
     y = float(DISP_Y[frame_idx])

dynamics.py

@@ -14,6 +14,7 @@ import sys
 import subprocess
 import time
 import argparse
+import json
 from contextlib import contextmanager
 from pathlib import Path
 
@@ -25,6 +26,57 @@ sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
 import compute
 
 
+def _fmt_alpha(v):
+    """将 alpha 值格式化为逗号分隔字符串，兼容 numpy 数组/list/标量。"""
+    if isinstance(v, (list, tuple, np.ndarray)):
+        return ",".join(str(float(x)) for x in v)
+    return str(float(v))
+
+
+def _json_field(value):
+    """Serialize arrays/lists for display header metadata."""
+    if isinstance(value, np.ndarray):
+        value = value.tolist()
+    return json.dumps(value, ensure_ascii=False)
+
+
+def _load_camera_kf(config, runtime_base):
+    """加载 move_camera.txt（速度段格式）→ JSON 字符串。"""
+    import re, json
+    if not int(config.get("move_camera", 0)):
+        return ""
+    cam_file = str(config.get("move_camera_file",
+                    os.path.join("input", "move_camera.txt")))
+    cam_path = cam_file
+    if not os.path.isabs(cam_file):
+        cam_path = os.path.join(runtime_base, cam_file)
+    if not os.path.exists(cam_path):
+        return ""
+    segments = []
+    with open(cam_path, "r", encoding="utf-8") as f:
+        for line in f:
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            # 解析帧范围：支持 "all"（全程）或 "N-M"（区间）
+            if line.lower().startswith("all") or re.match(r'^\s*all\s', line, re.IGNORECASE):
+                start, end = 0, 10**9
+            else:
+                m = re.match(r'(\d+)\s*-\s*(\d+)', line)
+                if not m:
+                    continue
+                start, end = int(m.group(1)), int(m.group(2))
+            v, r = [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]
+            for i, axis in enumerate(['x', 'y', 'z']):
+                m2 = re.search(r'v' + axis + r'\s*=\s*([-\d.]+)', line)
+                if m2: v[i] = float(m2.group(1))
+                m2 = re.search(r'r' + axis + r'\s*=\s*([-\d.]+)', line)
+                if m2: r[i] = float(m2.group(1))
+            if any(v) or any(r):
+                segments.append({"start": start, "end": end, "v": v, "r": r})
+    return json.dumps(segments) if segments else ""
+
+
 def read_optional_index(data, key, default_value):
     """Read an optional integer index from structured txt metadata."""
     if key not in data:
@@ -94,16 +146,6 @@ def build_sample_indices(total_steps, sample_step, sample_start, sample_end):
     return indices
 
 
-def save_display_txt(data, out_dir=None):
-    """将抽帧数据保存到 output/display.txt（含所有参数元数据）。"""
-    if out_dir is None:
-        out_dir = os.path.dirname(os.path.abspath(__file__))
-    disp_path = os.path.join(compute.get_output_dir(out_dir), "display.txt")
-    compute.save_text_data(disp_path, data)
-    print(f"[sample] 显示数组已保存至: {disp_path}")
-    return disp_path
-
-
 def run_case(config_path, runtime_base, input_dir="input", output_dir="output", no_plot=False):
     """Run one case with explicit program path, input path, and output path."""
     runtime_base = Path(runtime_base).resolve()
@@ -193,273 +235,107 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
             _t0 = _time.time()
 
             if engine == "python":
-                traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = compute.run_from_config(config, str(runtime_base))
-                compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base))
+                compute.run_from_config(config, str(runtime_base))
             else:
-                # 外部引擎：先加载配置到全局变量，再运行引擎，再用 save_trajectory_txt 补全 metadata
+                # 外部引擎：先加载配置到全局变量，再运行引擎
                 config["_skip_run"] = True
                 compute.run_from_config(config, str(runtime_base))
                 config.pop("_skip_run", None)
                 input_dir_abs = str(input_dir_path.resolve())
                 output_dir_abs = str(output_dir_path.resolve())
+                # 外部引擎写完整 trajectory.txt，后续抽帧
                 traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = compute.run_engine(
                     engine, input_dir_abs, output_dir_abs, config)
-                compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base))
+                if int(config.get("save_trajectory", 0)):
+                    compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base))
 
             _elapsed = _time.time() - _t0
             print(f"[run] 引擎: {engine}  计算完成: {record_steps} 步  {_elapsed:.3f} s")
         else:
-            print("[run] 步骤 [模拟] 已跳过，直接加载已有轨迹")
+            print("[run] 步骤 [模拟] 已跳过")
+
+        # 3. 检查/生成 display.txt
+        disp_path_new = os.path.join(output_dir_abs, "display.txt")
+        save_traj = int(config.get("save_trajectory", 0))
+        
+        if os.path.exists(disp_path_new):
+            # Python 引擎或新版外部引擎（save_trajectory=0）已直接写入
+            print(f"[run] 发现已有 display.txt（引擎直接抽帧）")
+        elif engine != "python" and os.path.exists(os.path.join(output_dir_abs, "trajectory.txt")):
+            # 旧版外部引擎：从 trajectory.txt 抽帧
+            traj_path = os.path.join(output_dir_abs, "trajectory.txt")
             if not os.path.exists(traj_path):
-                print(f"[run] 错误: trajectory.txt 不存在，无法跳过模拟")
+                print(f"[run] 错误: 找不到 trajectory.txt 或 display.txt")
                 sys.exit(1)
+            data = compute.load_text_data(traj_path)
+            NT = int(data["NT"]); DT = float(data["DT"]); NSTEP = int(data.get("NSTEP", 1))
+            record_steps = NT - int(data.get("warmup_steps", 0))
+            n_atoms = len(data["atom_ids"])
+            sample_start = 0
+            sample_end = NT
+            indices = np.arange(0, record_steps, NSTEP, dtype=np.int64)
+            if len(indices) == 0:
+                indices = np.array([0])
                 
-        # 3. 抽帧 → output/display.txt
-        traj_path = os.path.join(output_dir_abs, "trajectory.txt")
-        data = compute.load_text_data(traj_path)
+            traj_x = data["traj_x"]; traj_y = data["traj_y"]; traj_z = data["traj_z"]
+            traj_vx = data["traj_vx"]; traj_vy = data["traj_vy"]; traj_vz = data["traj_vz"]
             
-        NT = int(data["NT"]); DT = float(data["DT"]); NSTEP = int(data["NSTEP"])
-        warmup_steps = int(data.get("warmup_steps", 0))
-        plot_atom_row = int(data["plot_atom_row"]) if "plot_atom_row" in data else 0
-        plot_atom_id = int(data["plot_atom_id"]) if "plot_atom_id" in data else int(data["atom_ids"][plot_atom_row])
+            # 构建 header_fields
+            hf = {"DT": str(DT), "NSTEP": str(NSTEP), "method": str(data.get("method", "")),
+                  "warmup_steps": str(data.get("warmup_steps", 0)),
+                  "dynamic_steps": str(record_steps),
+                  "T_total": str(NT * DT),
+                  "X_MAX": str(data.get("X_MAX", 10)), "X_MIN": str(data.get("X_MIN", -10)),
+                  "Y_MAX": str(data.get("Y_MAX", 10)), "Y_MIN": str(data.get("Y_MIN", -10)),
+                  "Z_MAX": str(data.get("Z_MAX", 10)), "Z_MIN": str(data.get("Z_MIN", -10)),
+                  "ball_radius": str(data.get("ball_radius", 0.5)),
+                  "ball_color_r": str(data.get("ball_color_r", 0.9)),
+                  "ball_color_g": str(data.get("ball_color_g", 0.2)),
+                  "ball_color_b": str(data.get("ball_color_b", 0.2)),
+                  "box_color_r": str(data.get("box_color_r", 0.8)),
+                  "box_color_g": str(data.get("box_color_g", 0.8)),
+                  "box_color_b": str(data.get("box_color_b", 0.85)),
+                  "gravity_field": str(data.get("gravity_field", 1)),
+                  "gravity_interaction": str(data.get("gravity_interaction", 0)),
+                  "elastic_force": str(data.get("elastic_force", 1)),
+                  "damping_force": str(data.get("damping_force", 0)),
+                  "gravity_strength": str(data.get("gravity_strength", 1.0)),
+                  "driving_force": str(data.get("driving_force", 0)),
+                  "use_marker": str(config.get("use_marker", 0)),
+                  "alpha": _fmt_alpha(data.get("alpha", 0.2)),
+                  "atom_masses": _json_field(data.get("atom_masses", [])),
+                  "atom_positions": _json_field(data.get("atom_positions", [])),
+                  "bond_pairs": _json_field(data.get("bond_pairs", [])),
+                  "bond_stiffness": _json_field(data.get("bond_stiffness", [])),
+                  "bond_rest_lengths": _json_field(data.get("bond_rest_lengths", [])),
+                  "G": _json_field(data.get("G", [0.0, 0.0, 0.0])),
+                  "atom_radii": _fmt_alpha(data.get("atom_radii", [])),
+                  "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)),
+                  "camera_keyframes": _load_camera_kf(config, str(runtime_base))}
             
-        # 抽帧范围控制
-        sample_start = read_optional_index(data, "sample_start", 0)
-        sample_end = read_optional_index(data, "sample_end", NT)
+            n_frames = len(indices)
+            compute.save_display_txt(
+                disp_path_new,
+                traj_x[indices], traj_y[indices], traj_z[indices],
+                traj_vx[indices], traj_vy[indices], traj_vz[indices],
+                np.array(data["atom_ids"]), n_frames, n_atoms,
+                header_fields=hf)
+            print(f"[run] 从 trajectory.txt 抽帧生成 display.txt ({n_frames} 帧)")
 
-        indices = build_sample_indices(NT, NSTEP, sample_start, sample_end)
-        n_frames = len(indices)
-
-        print(f"[run] 抽帧范围: [{sample_start}, {sample_end}), 共 {n_frames} 帧")
-
-        traj_x = data["traj_x"]
-        traj_y = data["traj_y"]
-        traj_z = data["traj_z"]
-        traj_vx = data["traj_vx"]
-        traj_vy = data["traj_vy"]
-        traj_vz = data["traj_vz"]
-
-        if traj_x.ndim == 1:
-            selected_x = traj_x
-            selected_y = traj_y
-            selected_z = traj_z
-            selected_vx = traj_vx
-            selected_vy = traj_vy
-            selected_vz = traj_vz
-            all_x = traj_x[:, None]
-            all_y = traj_y[:, None]
-            all_z = traj_z[:, None]
-            all_vx = traj_vx[:, None]
-            all_vy = traj_vy[:, None]
-            all_vz = traj_vz[:, None]
-        else:
-            selected_x = traj_x[:, plot_atom_row]
-            selected_y = traj_y[:, plot_atom_row]
-            selected_z = traj_z[:, plot_atom_row]
-            selected_vx = traj_vx[:, plot_atom_row]
-            selected_vy = traj_vy[:, plot_atom_row]
-            selected_vz = traj_vz[:, plot_atom_row]
-            all_x = traj_x
-            all_y = traj_y
-            all_z = traj_z
-            all_vx = traj_vx
-            all_vy = traj_vy
-            all_vz = traj_vz
-
-        if config.get("step_sample", 1):
-            disp_data = {
-                "disp_x": selected_x[indices],
-                "disp_y": selected_y[indices],
-                "disp_z": selected_z[indices],
-                "disp_vx": selected_vx[indices],
-                "disp_vy": selected_vy[indices],
-                "disp_vz": selected_vz[indices],
-                "disp_all_x": all_x[indices],
-                "disp_all_y": all_y[indices],
-                "disp_all_z": all_z[indices],
-                "disp_all_vx": all_vx[indices],
-                "disp_all_vy": all_vy[indices],
-                "disp_all_vz": all_vz[indices],
-                "disp_t": indices * DT,
-                "disp_step": indices,
-                "n_frames": n_frames,
-                "NT": NT, "DT": DT, "NSTEP": NSTEP,
-                "plot_atom_id": plot_atom_id,
-                "plot_atom_row": plot_atom_row,
-                "method": str(data["method"]) if "method" in data else "explicit_euler",
-                "coord_file": str(data["coord_file"]) if "coord_file" in data else os.path.join("input", "coord.txt"),
-                "atom_ids": data["atom_ids"] if "atom_ids" in data else np.array([1]),
-                "atom_masses": data["atom_masses"] if "atom_masses" in data else np.array([float(data["M"])]),
-                "atom_radii": data["atom_radii"] if "atom_radii" in data else np.array([float(data["ball_radius"])]),
-                "atom_positions": data["atom_positions"] if "atom_positions" in data else np.array([[float(data["X0"]), float(data["Y0"]), float(data["Z0"])]]),
-                "atom_velocities": data["atom_velocities"] if "atom_velocities" in data else np.array([[float(data["VX0"]), float(data["VY0"]), float(data["VZ0"])]]),
-                "atom_fixed": data["atom_fixed"] if "atom_fixed" in data else np.array([[0, 0, 0]]),
-                "bond_pairs": data.get("bond_pairs", np.zeros((0, 2), dtype=np.int64)).tolist(),
-                "bond_stiffness": data.get("bond_stiffness", np.zeros(0, dtype=np.float64)).tolist(),
-                "bond_rest_lengths": data.get("bond_rest_lengths", np.zeros(0, dtype=np.float64)).tolist(),
-                "warmup_steps": warmup_steps,
-                "sample_start": sample_start,
-                "sample_end": sample_end,
-                "X_MIN": float(data["X_MIN"]), "X_MAX": float(data["X_MAX"]),
-                "Y_MIN": float(data["Y_MIN"]), "Y_MAX": float(data["Y_MAX"]),
-                "Z_MIN": float(data["Z_MIN"]), "Z_MAX": float(data["Z_MAX"]),
-                "X0": float(data["X0"]), "Y0": float(data["Y0"]), "Z0": float(data["Z0"]),
-                "VX0": float(data["VX0"]), "VY0": float(data["VY0"]), "VZ0": float(data["VZ0"]),
-                "M": float(data["M"]) if "M" in data else 1.0,
-                "alpha": data["alpha"],
-                "ball_radius": float(data["ball_radius"]),
-                "ball_color_r": float(data["ball_color_r"]),
-                "ball_color_g": float(data["ball_color_g"]),
-                "ball_color_b": float(data["ball_color_b"]),
-                "box_color_r": float(data["box_color_r"]),
-                "box_color_g": float(data["box_color_g"]),
-                "box_color_b": float(data["box_color_b"]),
-                "gravity_field": int(data.get("gravity_field", 1)),
-                "gravity_interaction": int(data.get("gravity_interaction", 0)),
-                "elastic_force": int(data.get("elastic_force", 1)),
-                "damping_force": int(data.get("damping_force", 0)),
-                "gravity_strength": float(data.get("gravity_strength", 1.0)),
-                "driving_force": int(config.get("driving_force", 0)),
-                "use_marker": int(config.get("use_marker", 0)),
-            }
-            save_display_txt(disp_data, str(runtime_base))
-            print(f"[run] 抽帧完成: {sample_end - sample_start} 步 -> {n_frames} 帧")
-        else:
-            print("[run] 步骤 [抽帧] 已跳过")
+            # save_trajectory=0 时清理 trajectory.txt
+            if not save_traj:
+                try:
+                    os.remove(traj_path)
+                    print(f"[run] save_trajectory=0，已删除 {traj_path}")
+                except OSError:
+                    pass
 
         # 4. 绘图（可选）
         if not no_plot and config.get("step_plot", 1):
-            try:
-                import matplotlib.pyplot as plt
-
-                # 配置中文字体支持
-                plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'WenQuanYi Micro Hei', 'DejaVu Sans']
-                plt.rcParams['axes.unicode_minus'] = False
-
-                time_arr = np.arange(NT) * DT
-                n_atoms = all_x.shape[1]
-                atom_ids_list = data.get("atom_ids", np.arange(n_atoms) + 1)
-
-                fig, axes = plt.subplots(3, 3, figsize=(15, 13))
-                fig.suptitle("轨迹与能量分析", fontsize=16)
-
-                # ── 位置 / 速度 6 子图（前 2 行，每行 3 列） ──
-                plot_configs = [
-                    (axes[0, 0], all_x, "x - 时间"),
-                    (axes[0, 1], all_y, "y - 时间"),
-                    (axes[0, 2], all_z, "z - 时间"),
-                    (axes[1, 0], all_vx, "vx - 时间"),
-                    (axes[1, 1], all_vy, "vy - 时间"),
-                    (axes[1, 2], all_vz, "vz - 时间"),
-                ]
-
-                colors = plt.cm.tab10(np.linspace(0, 1, n_atoms))
-
-                for ax, data_arr, title in plot_configs:
-                    for i in range(n_atoms):
-                        atom_id = int(atom_ids_list[i])
-                        ax.plot(time_arr, data_arr[:, i], color=colors[i], linewidth=1.5, label=f"原子 {atom_id}")
-                    ax.set_title(title)
-                    ax.set_xlabel("时间 (s)")
-                    ax.grid(True, alpha=0.3)
-                    ax.legend()
-
-                # ── 能量计算 ─────────────────────────────────────
-                masses = np.array(data["atom_masses"])                     # (n_atoms,)
-                G_vec = np.array(data.get("G", [0.0, 0.0, -9.8]))        # [gx, gy, gz]
-                gravity_field_enabled = int(data.get("gravity_field", 1))
-                gravity_interaction_enabled = int(data.get("gravity_interaction", 0))
-                gravity_strength = float(data.get("gravity_strength", 1.0))
-                elastic_force_enabled = int(data.get("elastic_force", 1))
-                damping_force_enabled = int(data.get("damping_force", 0))
-
-                # 动能 Ek = ½ m v²
-                ek = 0.5 * masses[np.newaxis, :] * (all_vx**2 + all_vy**2 + all_vz**2)
-
-                # 均匀重力场势能 Ug = -m G·r
-                ug = np.zeros_like(ek)
-                if gravity_field_enabled:
-                    ug = -masses[np.newaxis, :] * (
-                        G_vec[0] * all_x + G_vec[1] * all_y + G_vec[2] * all_z
-                    )
-
-                # 弹性势能 Us = ½ k (d - d₀)²
-                us = np.zeros_like(ek)
-                bond_pairs = data.get("bond_pairs")
-                bond_stiffness = data.get("bond_stiffness")
-                bond_rest_lengths = data.get("bond_rest_lengths")
-                if elastic_force_enabled and bond_pairs is not None and len(bond_pairs) > 0:
-                    for b_idx in range(len(bond_pairs)):
-                        i, j = bond_pairs[b_idx]
-                        dx = all_x[:, j] - all_x[:, i]
-                        dy = all_y[:, j] - all_y[:, i]
-                        dz = all_z[:, j] - all_z[:, i]
-                        dist = np.sqrt(dx**2 + dy**2 + dz**2)
-                        stretch = dist - bond_rest_lengths[b_idx]
-                        us_each = 0.5 * bond_stiffness[b_idx] * stretch**2
-                        us[:, i] += us_each          # 将整根键的势能记给 i
-
-                # 万有引力势能 Ug_grav = -G_grav * m_i * m_j / r
-                ug_grav = np.zeros_like(ek)
-                if gravity_interaction_enabled:
-                    n_atoms_en = len(masses)
-                    for i in range(n_atoms_en):
-                        for j in range(i + 1, n_atoms_en):
-                            dx = all_x[:, j] - all_x[:, i]
-                            dy = all_y[:, j] - all_y[:, i]
-                            dz = all_z[:, j] - all_z[:, i]
-                            dist = np.sqrt(dx**2 + dy**2 + dz**2)
-                            dist = np.maximum(dist, 1e-12)
-                            pair_pe = -gravity_strength * masses[i] * masses[j] / dist
-                            ug_grav[:, i] += 0.5 * pair_pe
-                            ug_grav[:, j] += 0.5 * pair_pe
-
-                # 各原子总能量
-                e_total = ek + ug + us + ug_grav                           # (NT, n_atoms)
-
-                # 系统能量分量
-                ek_sys = np.sum(ek, axis=1)
-                ug_sys = np.sum(ug, axis=1)
-                us_sys = np.sum(us, axis=1)
-                ug_grav_sys = np.sum(ug_grav, axis=1)
-                e_sys = ek_sys + ug_sys + us_sys + ug_grav_sys
-
-                # ── 第 3 行左：各原子总能量 ──
-                ax_e = axes[2, 0]
-                for i in range(n_atoms):
-                    aid = int(atom_ids_list[i])
-                    ax_e.plot(time_arr, e_total[:, i], color=colors[i], linewidth=1.5, label=f"原子 {aid}")
-                ax_e.set_title("各原子总能量")
-                ax_e.set_xlabel("时间 (s)")
-                ax_e.set_ylabel("能量")
-                ax_e.grid(True, alpha=0.3)
-                ax_e.legend(loc="upper right")
-
-                # ── 第 3 行右：系统总能量 ──
-                ax_sys = axes[2, 1]
-                ax_sys.plot(time_arr, ek_sys, 'b-', linewidth=1.5, label="系统动能")
-                ax_sys.plot(time_arr, ug_sys, 'g-', linewidth=1.5, label="均匀重力势能")
-                if elastic_force_enabled and bond_pairs is not None and len(bond_pairs) > 0:
-                    ax_sys.plot(time_arr, us_sys, color='orange', linewidth=1.5, label="系统弹性势能")
-                if gravity_interaction_enabled:
-                    ax_sys.plot(time_arr, ug_grav_sys, color='purple', linewidth=1.5, label="万有引力势能")
-                ax_sys.plot(time_arr, e_sys, 'r--', linewidth=1.5, label="系统总能量")
-                ax_sys.set_title("系统总能量")
-                ax_sys.set_xlabel("时间 (s)")
-                ax_sys.set_ylabel("能量")
-                ax_sys.grid(True, alpha=0.3)
-                ax_sys.legend(loc="upper right")
-
-                # 隐藏第 3 行第 3 列空子图
-                axes[2, 2].set_visible(False)
-
-                plt.tight_layout(rect=[0, 0.03, 1, 0.95])
-                plot_path = os.path.join(output_dir_abs, "trajectory_plots.png")
-                plt.savefig(plot_path, dpi=300, bbox_inches="tight")
-                print(f"[run] 轨迹图表已保存至: {plot_path}")
-                plt.show()
-            except ImportError:
-                print("[run] 警告: 未安装 matplotlib，跳过绘图")
+            print("[run] 注意: 旧版 step_plot 绘图路径依赖完整轨迹局部变量，当前已暂时跳过。")
+            print("[run] 如需波形/能量动画，请使用 step_plot_wave: 1。")
 
         print(f"[run] 完成！输出目录: {output_dir_abs}")
 

engines/c/main.c

@@ -38,6 +38,13 @@ typedef struct {
     int     damping_force;       /* 阻尼开关 */
     double  gravity_strength;    /* 万有引力强度 */
     int     driving_force;       /* 驱动力开关 */
+    int     save_trajectory;     /* 是否保存完整轨迹文件 */
+    double  alpha[6];            /* 盒子透明度 */
+    double  ball_radius;
+    double  ball_color[3];
+    double  box_color[3];
+    int     use_marker;
+    double  camera_distance, camera_elevation, camera_azimuth;
 } SimParams;
 
 /* ========================================================================
@@ -270,6 +277,20 @@ static void json_read_double3(const char *json, const char *key, double out[3])
     }
 }
 
+static void json_read_double6(const char *json, const char *key, double out[6]) {
+    char search[256];
+    snprintf(search, sizeof(search), "\"%s\"", key);
+    const char *p = strstr(json, search);
+    if (!p) { for (int i=0;i<6;i++) out[i]=0; return; }
+    p = strchr(p, '[');
+    if (!p) { for (int i=0;i<6;i++) out[i]=0; return; }
+    p++;
+    for (int i = 0; i < 6; i++) {
+        while (*p == ' ' || *p == '\t' || *p == '\n' || *p == ',' || *p == ']') p++;
+        out[i] = strtod(p, (char**)&p);
+    }
+}
+
 /* 读取 param.json */
 static int    g_gravity_field = 1;
 static int    g_gravity_interaction = 0;
@@ -304,6 +325,16 @@ static SimParams read_params(const char *path) {
     p.damping_force       = json_read_int(buf, "damping_force");
     p.gravity_strength    = json_read_double(buf, "gravity_strength");
     p.driving_force       = json_read_int(buf, "driving_force");
+    p.save_trajectory     = json_read_int(buf, "save_trajectory");
+    /* 渲染参数 */
+    json_read_double6(buf, "alpha", p.alpha);
+    p.ball_radius         = json_read_double(buf, "ball_radius");
+    json_read_double3(buf, "ball_color", p.ball_color);
+    json_read_double3(buf, "box_color", p.box_color);
+    p.use_marker          = json_read_int(buf, "use_marker");
+    p.camera_distance     = json_read_double(buf, "camera_distance");
+    p.camera_elevation    = json_read_double(buf, "camera_elevation");
+    p.camera_azimuth      = json_read_double(buf, "camera_azimuth");
     g_gravity_field       = p.gravity_field;
     g_gravity_interaction = p.gravity_interaction;
     g_elastic_force       = p.elastic_force;
@@ -389,6 +420,7 @@ static BondData read_bonds(const char *input_dir, const AtomData *atoms) {
     char bond_name[256];
     while (fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name) == 3) n_lines++;
     rewind(f);
+    fgets(line, sizeof(line), f);  // 再次跳过表头
 
     if (n_lines == 0) { fclose(f); return b; }
 
@@ -501,6 +533,20 @@ static void compute_acceleration(
     }
 }
 
+/* 保守力加速度（不含阻尼），供真蛙跳法专用。
+   通过传入零速度调用 compute_acceleration，阻尼项 -B*v/m 自动为零。 */
+static void compute_accel_conservative(
+    int n, const double *x, const double *y, const double *z,
+    const double *m, const double G[3],
+    const BondData *bonds,
+    double *ax, double *ay, double *az)
+{
+    double *v0 = (double*)alloca(n * sizeof(double));
+    for (int i = 0; i < n; i++) v0[i] = 0.0;
+    double Bzero[3] = {0.0, 0.0, 0.0};
+    compute_acceleration(n, x, y, z, v0, v0, v0, m, G, Bzero, bonds, ax, ay, az);
+}
+
 /* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
 static void limit_in_box(double *pos, double *vel, double lo, double hi) {
     if (*pos > hi)  { *pos = hi;  *vel = -*vel; }
@@ -618,6 +664,15 @@ static void midpoint_step(
 }
 
 /* ── 蛙跳法（Velocity-Verlet）── */
+/* 真蛙跳一步：x(t), v(t-dt/2) → x(t+dt), v(t+dt/2)
+ *
+ * 无阻尼：纯保守蛙跳，每步 1 次力计算，辛积分器。
+ *   v(t+dt/2) = v(t-dt/2) + a_c(t)·dt
+ *
+ * 有阻尼：半隐式处理，仍 1 次力计算，对任意阻尼无条件稳定。
+ *   利用 v(t) ≈ [v(t-dt/2) + v(t+dt/2)] / 2 解析求解：
+ *   v(t+dt/2) = [v(t-dt/2)·(1-α) + a_c(t)·dt] / (1+α)，α = B·dt/(2m)
+ */
 static void leapfrog_step(
     int n, double *x, double *y, double *z,
     double *vx, double *vy, double *vz,
@@ -627,47 +682,30 @@ static void leapfrog_step(
     double *ax = (double*)alloca(n * sizeof(double));
     double *ay = (double*)alloca(n * sizeof(double));
     double *az = (double*)alloca(n * sizeof(double));
-    compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az);
 
-    /* 半推速度：v_half = v + 0.5*a*dt */
+    /* 1 次保守力计算（不含阻尼） */
+    compute_accel_conservative(n, x, y, z, m, G, bonds, ax, ay, az);
+
+    int has_damping = g_damping_force && (B[0] != 0.0 || B[1] != 0.0 || B[2] != 0.0);
+
     for (int i = 0; i < n; i++) {
         if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        vx[i] += ax[i] * dt * 0.5;
-        vy[i] += ay[i] * dt * 0.5;
-        vz[i] += az[i] * dt * 0.5;
-    }
-
-    /* 全推位置（不含边界）*/
-    for (int i = 0; i < n; i++) {
-        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        x[i] += vx[i] * dt;   /* vx 此时是 v_half */
+        if (has_damping) {
+            double alphax = B[0] * dt / (2.0 * m[i]);
+            double alphay = B[1] * dt / (2.0 * m[i]);
+            double alphaz = B[2] * dt / (2.0 * m[i]);
+            vx[i] = (vx[i] * (1.0 - alphax) + ax[i] * dt) / (1.0 + alphax);
+            vy[i] = (vy[i] * (1.0 - alphay) + ay[i] * dt) / (1.0 + alphay);
+            vz[i] = (vz[i] * (1.0 - alphaz) + az[i] * dt) / (1.0 + alphaz);
+        } else {
+            vx[i] += ax[i] * dt;
+            vy[i] += ay[i] * dt;
+            vz[i] += az[i] * dt;
+        }
+        x[i] += vx[i] * dt;
         y[i] += vy[i] * dt;
         z[i] += vz[i] * dt;
     }
-
-    /* 显式预测器：v_pred = v_half + 0.5*a_old*dt，用第一次加速度外推半步
-       包含重力+阻尼+弹簧的所有贡献（标准 Velocity-Verlet 预测步）*/
-    double *pred_vx = (double*)alloca(n * sizeof(double));
-    double *pred_vy = (double*)alloca(n * sizeof(double));
-    double *pred_vz = (double*)alloca(n * sizeof(double));
-    for (int i = 0; i < n; i++) {
-        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        pred_vx[i] = vx[i] + 0.5 * ax[i] * dt;
-        pred_vy[i] = vy[i] + 0.5 * ay[i] * dt;
-        pred_vz[i] = vz[i] + 0.5 * az[i] * dt;
-    }
-
-    /* 用新位置 + 预测速度重算加速度 */
-    compute_acceleration(n, x, y, z, pred_vx, pred_vy, pred_vz, m, G, B, bonds, ax, ay, az);
-
-    /* 速度后半步：v = v_half + 0.5*a_next*dt
-       vx 仍为 v_half（未被覆盖）*/
-    for (int i = 0; i < n; i++) {
-        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        vx[i] += ax[i] * dt * 0.5;
-        vy[i] += ay[i] * dt * 0.5;
-        vz[i] += az[i] * dt * 0.5;
-    }
 }
 
 /* ── 驱动力（与 Python apply_driving_force 一致）──────────────── */
@@ -789,12 +827,14 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
     const char *names[] = {"traj_x","traj_y","traj_z","traj_vx","traj_vy","traj_vz"};
     double *arrs[] = {traj->x, traj->y, traj->z, traj->vx, traj->vy, traj->vz};
 
+    printf("[C-engine] 正在写入轨迹数据…\n");
+    fflush(stdout);
     for (int a = 0; a < 6; a++) {
         fprintf(f, "  \"%s\": [\n", names[a]);
         for (int t = 0; t < traj->n_steps; t++) {
             fprintf(f, "    [");
             for (int i = 0; i < traj->n_atoms; i++) {
-                fprintf(f, "%.15g", arrs[a][t * traj->n_atoms + i]);
+                fprintf(f, "%.8g", arrs[a][t * traj->n_atoms + i]);
                 if (i < traj->n_atoms - 1) fputc(',', f);
             }
             fprintf(f, "]");
@@ -808,12 +848,12 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
 
     /* 标量参数 */
     fprintf(f, "  \"NT\": %d,\n", params->NT);
-    fprintf(f, "  \"DT\": %.15g,\n", params->DT);
+    fprintf(f, "  \"DT\": %.8g,\n", params->DT);
     fprintf(f, "  \"NSTEP\": %d,\n", params->NSTEP);
     fprintf(f, "  \"method\": \"%s\",\n", params->method);
     fprintf(f, "  \"warmup_steps\": %d,\n", params->warmup_steps);
-    fprintf(f, "  \"G\": [%.15g, %.15g, %.15g],\n", params->G[0], params->G[1], params->G[2]);
-    fprintf(f, "  \"B\": [%.15g, %.15g, %.15g],\n", params->B[0], params->B[1], params->B[2]);
+    fprintf(f, "  \"G\": [%.8g, %.8g, %.8g],\n", params->G[0], params->G[1], params->G[2]);
+    fprintf(f, "  \"B\": [%.8g, %.8g, %.8g],\n", params->B[0], params->B[1], params->B[2]);
 
     fprintf(f, "  \"atom_ids\": [");
     for (int i = 0; i < atoms->n_atoms; i++) {
@@ -825,7 +865,7 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
     fprintf(f, "  \"atom_masses\": [");
     for (int i = 0; i < atoms->n_atoms; i++) {
         if (i > 0) fputc(',', f);
-        fprintf(f, "%.15g", atoms->masses[i]);
+        fprintf(f, "%.8g", atoms->masses[i]);
     }
     fprintf(f, "],\n");
 
@@ -839,14 +879,14 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
     fprintf(f, "  \"bond_stiffness\": [");
     for (int b = 0; b < bonds->n_bonds; b++) {
         if (b > 0) fputc(',', f);
-        fprintf(f, "%.15g", bonds->stiffness[b]);
+        fprintf(f, "%.8g", bonds->stiffness[b]);
     }
     fprintf(f, "],\n");
 
     fprintf(f, "  \"bond_rest_lengths\": [");
     for (int b = 0; b < bonds->n_bonds; b++) {
         if (b > 0) fputc(',', f);
-        fprintf(f, "%.15g", bonds->rest_lengths[b]);
+        fprintf(f, "%.8g", bonds->rest_lengths[b]);
     }
     fprintf(f, "],\n");
     fprintf(f, "  \"driving_force\": %d\n", params->driving_force);
@@ -855,6 +895,65 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
     fclose(f);
 }
 
+static void write_display_txt(const char *path, const Trajectory *traj,
+                               const SimParams *params, const AtomData *atoms)
+{
+    FILE *f = fopen(path, "w");
+    if (!f) die("无法写入 display.txt");
+
+    int n_frames    = traj->n_steps;   /* 实际采样帧数，用于下面的帧循环 */
+    int n_particles = traj->n_atoms;
+    int dynamic_steps = params->NT - params->warmup_steps;
+    double T_total = dynamic_steps * params->DT;
+
+    /* number of frames 写总积分步数（与 draw.py NT 对应），不是采样帧数 */
+    fprintf(f, "number of frames: %d\n", dynamic_steps);
+    fprintf(f, "number of particles: %d\n", n_particles);
+    fprintf(f, "DT: %.16g\n", params->DT);
+    fprintf(f, "NSTEP: %d\n", params->NSTEP);
+    fprintf(f, "method: %s\n", params->method);
+    fprintf(f, "warmup_steps: %d\n", params->warmup_steps);
+    fprintf(f, "dynamic_steps: %d\n", dynamic_steps);
+    fprintf(f, "T_total: %.16g\n", T_total);
+    fprintf(f, "box_a: %.16g\n", params->box_a);
+    fprintf(f, "alpha: %.16g,%.16g,%.16g,%.16g,%.16g,%.16g\n",
+            params->alpha[0], params->alpha[1], params->alpha[2],
+            params->alpha[3], params->alpha[4], params->alpha[5]);
+    fprintf(f, "ball_radius: %.16g\n", params->ball_radius);
+    fprintf(f, "ball_color_r: %.16g\n", params->ball_color[0]);
+    fprintf(f, "ball_color_g: %.16g\n", params->ball_color[1]);
+    fprintf(f, "ball_color_b: %.16g\n", params->ball_color[2]);
+    fprintf(f, "box_color_r: %.16g\n", params->box_color[0]);
+    fprintf(f, "box_color_g: %.16g\n", params->box_color[1]);
+    fprintf(f, "box_color_b: %.16g\n", params->box_color[2]);
+    fprintf(f, "use_marker: %d\n", params->use_marker);
+    fprintf(f, "camera_distance: %.16g\n", params->camera_distance);
+    fprintf(f, "camera_elevation: %.16g\n", params->camera_elevation);
+    fprintf(f, "camera_azimuth: %.16g\n", params->camera_azimuth);
+    fprintf(f, "\n");
+
+    if (params->driving_force) {
+        fprintf(f, "driving_force: 1\n");
+    } else {
+        fprintf(f, "driving_force: 0\n");
+    }
+    fprintf(f, "\n");
+
+    for (int t = 0; t < n_frames; t++) {
+        fprintf(f, "frame: %3d\n", t + 1);
+        fprintf(f, "n   x            y            z           vx         vy         vz\n");
+        for (int i = 0; i < n_particles; i++) {
+            int idx = t * n_particles + i;
+            fprintf(f, "%4d %12.6f %12.6f %12.6f %10.6f %10.6f %10.6f\n",
+                    atoms->atom_ids[i],
+                    traj->x[idx], traj->y[idx], traj->z[idx],
+                    traj->vx[idx], traj->vy[idx], traj->vz[idx]);
+        }
+    }
+
+    fclose(f);
+}
+
 // ========================================================================
 // 主函数
 // ========================================================================
@@ -899,17 +998,43 @@ int main(int argc, char **argv) {
         vz[i] = atoms.vel_0[i*3+2];
     }
 
-    /* 分配轨迹缓冲区：改用 record_steps */
+    /* 分配轨迹缓冲区 */
     int record_steps = params.NT - params.warmup_steps;
     Trajectory traj;
-    traj.n_steps = record_steps;
     traj.n_atoms = n;
-    traj.x  = (double*)xmalloc(record_steps * n * sizeof(double) * 6);
-    traj.y  = traj.x  + record_steps * n;
-    traj.z  = traj.y  + record_steps * n;
-    traj.vx = traj.z  + record_steps * n;
-    traj.vy = traj.vx + record_steps * n;
-    traj.vz = traj.vy + record_steps * n;
+    if (params.save_trajectory) {
+        traj.n_steps = record_steps;
+        traj.x  = (double*)xmalloc(record_steps * n * sizeof(double) * 6);
+        traj.y  = traj.x  + record_steps * n;
+        traj.z  = traj.y  + record_steps * n;
+        traj.vx = traj.z  + record_steps * n;
+        traj.vy = traj.vx + record_steps * n;
+        traj.vz = traj.vy + record_steps * n;
+    } else {
+        int sampled_steps = (record_steps + params.NSTEP - 1) / params.NSTEP;
+        if (sampled_steps < 1) sampled_steps = 1;
+        traj.n_steps = sampled_steps;
+        traj.x  = (double*)xmalloc(sampled_steps * n * sizeof(double) * 6);
+        traj.y  = traj.x  + sampled_steps * n;
+        traj.z  = traj.y  + sampled_steps * n;
+        traj.vx = traj.z  + sampled_steps * n;
+        traj.vy = traj.vx + sampled_steps * n;
+        traj.vz = traj.vy + sampled_steps * n;
+    }
+
+    /* 真蛙跳初始化：v(0) 反推 v(-dt/2) = v(0) - 0.5*a_c(0)*dt */
+    if (strcmp(params.method, "leapfrog") == 0) {
+        double *ax0 = (double*)alloca(n * sizeof(double));
+        double *ay0 = (double*)alloca(n * sizeof(double));
+        double *az0 = (double*)alloca(n * sizeof(double));
+        compute_accel_conservative(n, x, y, z, atoms.masses, params.G, &bonds, ax0, ay0, az0);
+        for (int i = 0; i < n; i++) {
+            if (atoms.fixed[i*3+0] && atoms.fixed[i*3+1] && atoms.fixed[i*3+2]) continue;
+            vx[i] -= 0.5 * ax0[i] * params.DT;
+            vy[i] -= 0.5 * ay0[i] * params.DT;
+            vz[i] -= 0.5 * az0[i] * params.DT;
+        }
+    }
 
     /* 预热 */
     /* 初始时刻 t=0 驱动力（与 Python run_simulation 一致）*/
@@ -925,16 +1050,28 @@ int main(int argc, char **argv) {
     }
 
     /* 记录 */
+    int _prog_interval = record_steps / 100;
+    if (_prog_interval < 1) _prog_interval = 1;
+    int sample_idx = 0;
     for (int s = 0; s < record_steps; s++) {
+        if (s % _prog_interval == 0 && s > 0) {
+            printf("[C-engine] progress: %d/%d\n", s, record_steps);
+            fflush(stdout);
+        }
         double t = (s + params.warmup_steps) * params.DT;
         if (params.driving_force) apply_driving_force(n, x, y, z, vx, vy, vz, t, s, params.DT, &drivers);
-        for (int i = 0; i < n; i++) {
-            traj.x[ s * n + i]  = x[i];
-            traj.y[ s * n + i]  = y[i];
-            traj.z[ s * n + i]  = z[i];
-            traj.vx[s * n + i]  = vx[i];
-            traj.vy[s * n + i]  = vy[i];
-            traj.vz[s * n + i]  = vz[i];
+        int do_record = params.save_trajectory || (s % params.NSTEP == 0);
+        if (do_record) {
+            int idx = params.save_trajectory ? s : sample_idx;
+            for (int i = 0; i < n; i++) {
+                traj.x[ idx * n + i]  = x[i];
+                traj.y[ idx * n + i]  = y[i];
+                traj.z[ idx * n + i]  = z[i];
+                traj.vx[idx * n + i]  = vx[i];
+                traj.vy[idx * n + i]  = vy[i];
+                traj.vz[idx * n + i]  = vz[i];
+            }
+            if (!params.save_trajectory) sample_idx++;
         }
         apply_step(params.method, n, x, y, z, vx, vy, vz,
                    atoms.masses, params.G, params.B, &bonds, atoms.fixed,
@@ -943,8 +1080,13 @@ int main(int argc, char **argv) {
     }
 
     char out_path[512];
-    snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir);
-    write_trajectory_json(out_path, &traj, &params, &atoms, &bonds);
+    if (params.save_trajectory) {
+        snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir);
+        write_trajectory_json(out_path, &traj, &params, &atoms, &bonds);
+    } else {
+        snprintf(out_path, sizeof(out_path), "%s/display.txt", output_dir);
+        write_display_txt(out_path, &traj, &params, &atoms);
+    }
 
     clock_t t1 = clock();
     double elapsed = (double)(t1 - t0) / CLOCKS_PER_SEC;

engines/cpp/main.cpp

@@ -43,6 +43,15 @@ struct SimParams {
     int damping_force = 0;
     double gravity_strength = 1.0;
     int driving_force = 0;
+    int save_trajectory = 1;
+    double alpha[6] = {0,0,0,0,0,0};
+    double ball_radius = 0.5;
+    double ball_color[3] = {0.9, 0.2, 0.2};
+    double box_color[3] = {0.8, 0.8, 0.85};
+    int use_marker = 0;
+    double camera_distance = 40.0;
+    double camera_elevation = 0;
+    double camera_azimuth = 0;
 };
 
 // ========================================================================
@@ -131,6 +140,11 @@ static std::string json_read_string(const std::string &json, const std::string &
     return result;
 }
 
+/* 检查 JSON 中是否存在某个 key */
+static bool json_has_key(const std::string &json, const std::string &key) {
+    return json.find("\"" + key + "\"") != std::string::npos;
+}
+
 /* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */
 static void json_read_double3(const std::string &json, const std::string &key, double out[3]) {
     auto pos = json.find("\"" + key + "\"");
@@ -146,6 +160,21 @@ static void json_read_double3(const std::string &json, const std::string &key, d
     }
 }
 
+/* 读取 JSON 数组到 double[6] */
+static void json_read_double6(const std::string &json, const std::string &key, double out[6]) {
+    auto pos = json.find("\"" + key + "\"");
+    if (pos == std::string::npos) { for (int i=0;i<6;i++) out[i]=0; return; }
+    pos = json.find('[', pos);
+    if (pos == std::string::npos) { for (int i=0;i<6;i++) out[i]=0; return; }
+    pos++;
+    for (int i = 0; i < 6; i++) {
+        while (pos < json.size() && (json[pos]==' '||json[pos]=='\t'||json[pos]=='\n'||json[pos]==','||json[pos]==']')) pos++;
+        char *end;
+        out[i] = std::strtod(json.c_str() + pos, &end);
+        pos = end - json.c_str();
+    }
+}
+
 /* 解析 param.json */
 static SimParams read_params(const std::string &path) {
     std::string buf = read_file(path);
@@ -165,6 +194,17 @@ static SimParams read_params(const std::string &path) {
     p.damping_force       = json_read_int(buf, "damping_force");
     p.gravity_strength    = json_read_double(buf, "gravity_strength");
     p.driving_force       = json_read_int(buf, "driving_force");
+    // save_trajectory 默认 1（全量），仅在 JSON 中存在该 key 时覆盖
+    if (json_has_key(buf, "save_trajectory"))
+        p.save_trajectory = json_read_int(buf, "save_trajectory");
+    json_read_double6(buf, "alpha", p.alpha);
+    p.ball_radius         = json_read_double(buf, "ball_radius");
+    json_read_double3(buf, "ball_color", p.ball_color);
+    json_read_double3(buf, "box_color", p.box_color);
+    p.use_marker          = json_read_int(buf, "use_marker");
+    p.camera_distance     = json_read_double(buf, "camera_distance");
+    p.camera_elevation    = json_read_double(buf, "camera_elevation");
+    p.camera_azimuth      = json_read_double(buf, "camera_azimuth");
     return p;
 }
 
@@ -380,6 +420,26 @@ static void compute_acceleration(
     }
 }
 
+/* 保守力加速度（不含阻尼），供真蛙跳法专用。
+   传入 damping_force=0 使 compute_acceleration 跳过阻尼项。 */
+static void compute_accel_conservative(
+    int n,
+    const double *x, const double *y, const double *z,
+    const double *m, const double G[3],
+    const BondData &bonds,
+    int gravity_field, int gravity_interaction,
+    int elastic_force, double gravity_strength,
+    double *ax, double *ay, double *az)
+{
+    std::vector<double> v0(n, 0.0);
+    double Bzero[3] = {0.0, 0.0, 0.0};
+    compute_acceleration(n, x, y, z, v0.data(), v0.data(), v0.data(),
+                         m, G, Bzero, bonds,
+                         gravity_field, gravity_interaction,
+                         elastic_force, 0, gravity_strength,
+                         ax, ay, az);
+}
+
 /* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
 static void limit_in_box(double &pos, double &vel, double lo, double hi) {
     if (pos > hi)  { pos = hi;  vel = -vel; }
@@ -503,6 +563,14 @@ static void midpoint_step(
 }
 
 /* ── 蛙跳法（Velocity-Verlet）——与 Python Leapfrog_Method 一致 ── */
+/* 真蛙跳一步：x(t), v(t-dt/2) → x(t+dt), v(t+dt/2)
+ *
+ * 无阻尼：纯保守蛙跳，每步 1 次力计算，辛积分器。
+ *   v(t+dt/2) = v(t-dt/2) + a_c(t)·dt
+ *
+ * 有阻尼：半隐式处理，仍 1 次力计算，对任意阻尼无条件稳定。
+ *   v(t+dt/2) = [v(t-dt/2)·(1-α) + a_c(t)·dt] / (1+α)，α = B·dt/(2m)
+ */
 static void leapfrog_full_step(
     int n, double *x, double *y, double *z,
     double *vx, double *vy, double *vz,
@@ -512,54 +580,33 @@ static void leapfrog_full_step(
     int elastic_force, int damping_force,
     double gravity_strength)
 {
-    // 第一次加速度
     std::vector<double> ax(n), ay(n), az(n);
-    compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds,
-                         gravity_field, gravity_interaction,
-                         elastic_force, damping_force, gravity_strength,
-                         ax.data(), ay.data(), az.data());
 
-    // 半推速度：v_half = v + 0.5*a*dt  (存入 vx, vy, vz)
+    // 1 次保守力计算（不含阻尼）
+    compute_accel_conservative(n, x, y, z, m, G, bonds,
+                               gravity_field, gravity_interaction,
+                               elastic_force, gravity_strength,
+                               ax.data(), ay.data(), az.data());
+
+    bool has_damping = damping_force && (B[0] != 0.0 || B[1] != 0.0 || B[2] != 0.0);
     for (int i = 0; i < n; i++) {
         if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        vx[i] += ax[i] * dt * 0.5;
-        vy[i] += ay[i] * dt * 0.5;
-        vz[i] += az[i] * dt * 0.5;
-    }
-
-    // 全推位置（不含边界，边界在外层统一处理）
-    for (int i = 0; i < n; i++) {
-        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        x[i] += vx[i] * dt;   // vx 此时是 v_half
+        if (has_damping) {
+            double alphax = B[0] * dt / (2.0 * m[i]);
+            double alphay = B[1] * dt / (2.0 * m[i]);
+            double alphaz = B[2] * dt / (2.0 * m[i]);
+            vx[i] = (vx[i] * (1.0 - alphax) + ax[i] * dt) / (1.0 + alphax);
+            vy[i] = (vy[i] * (1.0 - alphay) + ay[i] * dt) / (1.0 + alphay);
+            vz[i] = (vz[i] * (1.0 - alphaz) + az[i] * dt) / (1.0 + alphaz);
+        } else {
+            vx[i] += ax[i] * dt;
+            vy[i] += ay[i] * dt;
+            vz[i] += az[i] * dt;
+        }
+        x[i] += vx[i] * dt;
         y[i] += vy[i] * dt;
         z[i] += vz[i] * dt;
     }
-
-    // 显式预测器：v_pred = v_half + 0.5*a_old*dt，用第一次加速度外推半步
-    // 包含所有力的贡献（标准 Velocity-Verlet 预测步）
-    std::vector<double> pred_vx(n), pred_vy(n), pred_vz(n);
-    for (int i = 0; i < n; i++) {
-        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        pred_vx[i] = vx[i] + 0.5 * ax[i] * dt;
-        pred_vy[i] = vy[i] + 0.5 * ay[i] * dt;
-        pred_vz[i] = vz[i] + 0.5 * az[i] * dt;
-    }
-
-    // 用新位置 + 预测速度重算加速度
-    compute_acceleration(n, x, y, z, pred_vx.data(), pred_vy.data(), pred_vz.data(),
-                         m, G, B, bonds,
-                         gravity_field, gravity_interaction,
-                         elastic_force, damping_force, gravity_strength,
-                         ax.data(), ay.data(), az.data());
-
-    // 速度后半步：v = v_half + 0.5*a_next*dt
-    // vx 仍为 v_half（未被覆盖），直接加上 0.5*a_next*dt
-    for (int i = 0; i < n; i++) {
-        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
-        vx[i] += ax[i] * dt * 0.5;
-        vy[i] += ay[i] * dt * 0.5;
-        vz[i] += az[i] * dt * 0.5;
-    }
 }
 
 /* ── 分发器：调用对应积分方法 + 边界条件（与 Python apply_motion_update 一致）── */
@@ -681,6 +728,67 @@ static void apply_driving_force(
     }
 }
 
+// ========================================================================
+// display.txt 输出（save_trajectory=0 时使用）
+// ========================================================================
+
+static void write_display_txt(
+    const std::string &path,
+    const std::vector<double> &x, const std::vector<double> &y,
+    const std::vector<double> &z, const std::vector<double> &vx,
+    const std::vector<double> &vy, const std::vector<double> &vz,
+    int n_steps, int n_atoms,
+    const SimParams &params, const AtomData &atoms)
+{
+    std::ofstream f(path);
+    if (!f) die("无法写入 " + path);
+    std::cout << "[Cpp-engine] 正在写入显示数据…" << std::endl;
+
+    int dynamic_steps = params.NT - params.warmup_steps;
+    double T_total = dynamic_steps * params.DT;
+
+    /* number of frames 写总积分步数（与 draw.py NT 对应），不是采样帧数 */
+    f << "number of frames: " << dynamic_steps << "\n";
+    f << "number of particles: " << n_atoms << "\n";
+    f << "DT: " << params.DT << "\n";
+    f << "NSTEP: " << params.NSTEP << "\n";
+    f << "method: " << params.method << "\n";
+    f << "warmup_steps: " << params.warmup_steps << "\n";
+    f << "dynamic_steps: " << dynamic_steps << "\n";
+    f << "T_total: " << T_total << "\n";
+    f << "box_a: " << params.box_a << "\n";
+    f << "driving_force: " << params.driving_force << "\n";
+    f << "alpha: " << params.alpha[0] << "," << params.alpha[1] << "," << params.alpha[2] << ","
+                   << params.alpha[3] << "," << params.alpha[4] << "," << params.alpha[5] << "\n";
+    f << "ball_radius: " << params.ball_radius << "\n";
+    f << "ball_color_r: " << params.ball_color[0] << "\n";
+    f << "ball_color_g: " << params.ball_color[1] << "\n";
+    f << "ball_color_b: " << params.ball_color[2] << "\n";
+    f << "box_color_r: " << params.box_color[0] << "\n";
+    f << "box_color_g: " << params.box_color[1] << "\n";
+    f << "box_color_b: " << params.box_color[2] << "\n";
+    f << "use_marker: " << params.use_marker << "\n";
+    f << "camera_distance: " << params.camera_distance << "\n";
+    f << "camera_elevation: " << params.camera_elevation << "\n";
+    f << "camera_azimuth: " << params.camera_azimuth << "\n\n";
+
+    f << std::fixed << std::setprecision(6);
+    for (int t = 0; t < n_steps; t++) {
+        f << "frame:  " << (t + 1) << "\n";
+        f << "n   x            y            z           vx         vy         vz\n";
+        for (int i = 0; i < n_atoms; i++) {
+            int base = t * n_atoms + i;
+            f << std::setw(3) << atoms.ids[i]
+              << std::setw(12) << x[base]
+              << std::setw(12) << y[base]
+              << std::setw(12) << z[base]
+              << std::setw(12) << vx[base]
+              << std::setw(12) << vy[base]
+              << std::setw(12) << vz[base] << "\n";
+        }
+    }
+}
+
 // ========================================================================
 // JSON 输出
 // ========================================================================
@@ -695,7 +803,8 @@ static void write_trajectory_json(
 {
     std::ofstream f(path);
     if (!f) die("无法写入 " + path);
-    f << std::setprecision(15);
+    std::cout << "[Cpp-engine] 正在写入轨迹数据…" << std::endl;
+    f << std::setprecision(8);
 
     f << "{\n";
 
@@ -816,12 +925,29 @@ int main(int argc, char **argv) {
 
     // 分配轨迹缓冲区
     int record_steps = params.NT - params.warmup_steps;
-    std::vector<double> traj_x(record_steps * n);
-    std::vector<double> traj_y(record_steps * n);
-    std::vector<double> traj_z(record_steps * n);
-    std::vector<double> traj_vx(record_steps * n);
-    std::vector<double> traj_vy(record_steps * n);
-    std::vector<double> traj_vz(record_steps * n);
+    int nstep_sampling = (params.save_trajectory == 0) ? params.NSTEP : 1;
+    int buf_steps = (params.save_trajectory == 0) ? (record_steps / nstep_sampling) : record_steps;
+    std::vector<double> traj_x(buf_steps * n);
+    std::vector<double> traj_y(buf_steps * n);
+    std::vector<double> traj_z(buf_steps * n);
+    std::vector<double> traj_vx(buf_steps * n);
+    std::vector<double> traj_vy(buf_steps * n);
+    std::vector<double> traj_vz(buf_steps * n);
+
+    // 真蛙跳初始化：v(0) 反推 v(-dt/2) = v(0) - 0.5*a_c(0)*dt
+    if (params.method == "leapfrog") {
+        std::vector<double> ax0(n), ay0(n), az0(n);
+        compute_accel_conservative(n, x.data(), y.data(), z.data(), atoms.masses.data(), params.G, bonds,
+                                   params.gravity_field, params.gravity_interaction,
+                                   params.elastic_force, params.gravity_strength,
+                                   ax0.data(), ay0.data(), az0.data());
+        for (int i = 0; i < n; i++) {
+            if (atoms.fixed[i*3+0] && atoms.fixed[i*3+1] && atoms.fixed[i*3+2]) continue;
+            vx[i] -= 0.5 * ax0[i] * params.DT;
+            vy[i] -= 0.5 * ay0[i] * params.DT;
+            vz[i] -= 0.5 * az0[i] * params.DT;
+        }
+    }
 
     // 预热
     // 初始时刻 t=0 驱动力（与 Python run_simulation 一致）
@@ -843,18 +969,27 @@ int main(int argc, char **argv) {
     }
 
     // 记录
+    int _prog_int = record_steps / 100;
+    if (_prog_int < 1) _prog_int = 1;
+    int si = 0;  // 采样帧索引
     for (int s = 0; s < record_steps; s++) {
+        if (s % _prog_int == 0 && s > 0) {
+            std::cout << "[Cpp-engine] progress: " << s << "/" << record_steps << std::endl;
+        }
         double t = (s + params.warmup_steps) * params.DT;
         if (params.driving_force)
             apply_driving_force(n, x.data(), y.data(), z.data(), vx.data(), vy.data(), vz.data(), t, s, params.DT, drivers);
-        // 保存当前帧
-        for (int i = 0; i < n; i++) {
-            traj_x[s * n + i]  = x[i];
-            traj_y[s * n + i]  = y[i];
-            traj_z[s * n + i]  = z[i];
-            traj_vx[s * n + i] = vx[i];
-            traj_vy[s * n + i] = vy[i];
-            traj_vz[s * n + i] = vz[i];
+        // 保存当前帧（采样模式仅每 NSTEP 步保存一次）
+        if (s % nstep_sampling == 0) {
+            for (int i = 0; i < n; i++) {
+                traj_x[si * n + i]  = x[i];
+                traj_y[si * n + i]  = y[i];
+                traj_z[si * n + i]  = z[i];
+                traj_vx[si * n + i] = vx[i];
+                traj_vy[si * n + i] = vy[i];
+                traj_vz[si * n + i] = vz[i];
+            }
+            si++;
         }
 
         apply_step(params.method, n, x.data(), y.data(), z.data(),
@@ -868,9 +1003,15 @@ int main(int argc, char **argv) {
     }
 
     // 输出轨迹
-    std::string out_path = output_dir + "/trajectory.txt";
-    write_trajectory_json(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz,
-                          record_steps, n, params, atoms, bonds);
+    if (params.save_trajectory == 0) {
+        std::string out_path = output_dir + "/display.txt";
+        write_display_txt(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz,
+                          si, n, params, atoms);
+    } else {
+        std::string out_path = output_dir + "/trajectory.txt";
+        write_trajectory_json(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz,
+                              record_steps, n, params, atoms, bonds);
+    }
 
     auto t1 = std::chrono::high_resolution_clock::now();
     double elapsed = std::chrono::duration<double>(t1 - t0).count();

engines/fortran/main.f90

@@ -34,7 +34,7 @@ program dynamics_f90
     double precision, allocatable :: bond_stiffness(:), bond_rest_lengths(:)
 
     ! 驱动力数据
-    integer :: n_drivers
+    integer :: n_drivers, prog_step
     integer, allocatable :: drv_atom_idx(:)
     double precision, allocatable :: drv_amp_x(:), drv_amp_y(:), drv_amp_z(:)
     double precision, allocatable :: drv_freq_x(:), drv_freq_y(:), drv_freq_z(:)
@@ -107,6 +107,24 @@ program dynamics_f90
     allocate(traj_x(record_steps, n), traj_y(record_steps, n), traj_z(record_steps, n))
     allocate(traj_vx(record_steps, n), traj_vy(record_steps, n), traj_vz(record_steps, n))
 
+    ! 真蛙跳初始化：v(0) 反推 v(-dt/2) = v(0) - 0.5*a_c(0)*dt
+    if (trim(method) == 'leapfrog') then
+        block
+            double precision :: ax0(n), ay0(n), az0(n)
+            integer :: ii
+            call accel_conservative(n, x, y, z, masses, G, &
+                                    n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths, &
+                                    gravity_field, gravity_interaction, &
+                                    elastic_force, gravity_strength, ax0, ay0, az0)
+            do ii = 1, n
+                if (fixed(ii,1) /= 0 .and. fixed(ii,2) /= 0 .and. fixed(ii,3) /= 0) cycle
+                vx(ii) = vx(ii) - 0.5d0 * ax0(ii) * DT
+                vy(ii) = vy(ii) - 0.5d0 * ay0(ii) * DT
+                vz(ii) = vz(ii) - 0.5d0 * az0(ii) * DT
+            end do
+        end block
+    end if
+
     ! 预热
     ! 初始时刻 t=0 驱动力（与 Python run_simulation 一致）
     if (driving_force /= 0 .and. n_drivers > 0) then
@@ -139,7 +157,12 @@ program dynamics_f90
     end do
 
     ! 记录
+    prog_step = record_steps / 100
+    if (prog_step < 1) prog_step = 1
     do s = 1, record_steps
+        if (mod(s, prog_step) == 0 .and. s > 0) then
+            write(*, '("[Fortran-engine] progress: ", i0, "/", i0)') s, record_steps
+        end if
         if (driving_force /= 0 .and. n_drivers > 0) then
             tw = ((s-1 + warmup_steps) * 1.0d0) * DT
             call apply_driving(n, x, y, z, vx, vy, vz, tw, s-1, DT, &
@@ -161,6 +184,7 @@ program dynamics_f90
     end do
 
     ! 输出轨迹
+    write(*, '("[Fortran-engine] 正在写入轨迹数据…")')
     call write_json(output_dir, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, &
                     record_steps, n_atoms, atom_ids, masses, &
                     NT, DT, NSTEP, warmup_steps, method, G, B, &
@@ -518,6 +542,24 @@ pure subroutine accel(n, x, y, z, vx, vy, vz, m, G, B, &
     end if
 end subroutine accel
 
+! 保守力加速度（不含阻尼），供真蛙跳法专用。
+! 传入零速度、零 B 调用 accel，阻尼项 -B*v/m 自动为零。
+subroutine accel_conservative(n, x, y, z, m, G, nb, bp, bk, br, &
+                              gravity_field, gravity_interaction, &
+                              elastic_force, gravity_strength, ax, ay, az)
+    integer, intent(in) :: n, nb, bp(nb, 2)
+    integer, intent(in) :: gravity_field, gravity_interaction, elastic_force
+    double precision, intent(in) :: x(n), y(n), z(n), m(n), G(3)
+    double precision, intent(in) :: bk(nb), br(nb), gravity_strength
+    double precision, intent(out) :: ax(n), ay(n), az(n)
+    double precision :: v0(n), B0(3)
+    v0 = 0.0d0
+    B0 = 0.0d0
+    call accel(n, x, y, z, v0, v0, v0, m, G, B0, nb, bp, bk, br, &
+               gravity_field, gravity_interaction, &
+               elastic_force, 0, gravity_strength, ax, ay, az)
+end subroutine accel_conservative
+
 ! 边界条件：clamp 位置 + 速度反转
 subroutine limit_in_box(pos, vel, lo, hi)
     double precision, intent(inout) :: pos, vel
@@ -646,7 +688,13 @@ subroutine midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, &
     end do
 end subroutine midpoint_step
 
-! ── 蛙跳法（Velocity-Verlet）──
+! 真蛙跳一步：x(t), v(t-dt/2) → x(t+dt), v(t+dt/2)
+!
+! 无阻尼：纯保守蛙跳，每步 1 次力计算，辛积分器。
+!   v(t+dt/2) = v(t-dt/2) + a_c(t)*dt
+!
+! 有阻尼：半隐式处理，仍 1 次力计算，对任意阻尼无条件稳定。
+!   v(t+dt/2) = [v(t-dt/2)*(1-α) + a_c(t)*dt] / (1+α)，α = B*dt/(2m)
 subroutine leapfrog_full(n, x, y, z, vx, vy, vz, m, G, B, &
                          nb, bp, bk, br, fixed, dt, &
                          gravity_field, gravity_interaction, &
@@ -656,53 +704,36 @@ subroutine leapfrog_full(n, x, y, z, vx, vy, vz, m, G, B, &
     double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
     double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), dt, gravity_strength
     double precision :: ax(n), ay(n), az(n)
-    double precision :: dmp_vx(n), dmp_vy(n), dmp_vz(n)
-    double precision :: gx, gy, gz
+    double precision :: alphax, alphay, alphaz
+    logical :: has_damping
     integer :: i
 
-    call accel(n, x, y, z, vx, vy, vz, m, G, B, nb, bp, bk, br, &
-               gravity_field, gravity_interaction, &
-               elastic_force, damping_force, gravity_strength, ax, ay, az)
+    ! 1 次保守力计算（不含阻尼）
+    call accel_conservative(n, x, y, z, m, G, nb, bp, bk, br, &
+                            gravity_field, gravity_interaction, &
+                            elastic_force, gravity_strength, ax, ay, az)
+
+    has_damping = (damping_force /= 0) .and. &
+                  (B(1) /= 0.0d0 .or. B(2) /= 0.0d0 .or. B(3) /= 0.0d0)
 
-    ! 速度半步推
-    do i = 1, n
-        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
-        vx(i) = vx(i) + ax(i) * dt * 0.5d0
-        vy(i) = vy(i) + ay(i) * dt * 0.5d0
-        vz(i) = vz(i) + az(i) * dt * 0.5d0
-    end do
-
-    ! 全推位置（不含边界）
     do i = 1, n
         if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
+        if (has_damping) then
+            alphax = B(1) * dt / (2.0d0 * m(i))
+            alphay = B(2) * dt / (2.0d0 * m(i))
+            alphaz = B(3) * dt / (2.0d0 * m(i))
+            vx(i) = (vx(i) * (1.0d0 - alphax) + ax(i) * dt) / (1.0d0 + alphax)
+            vy(i) = (vy(i) * (1.0d0 - alphay) + ay(i) * dt) / (1.0d0 + alphay)
+            vz(i) = (vz(i) * (1.0d0 - alphaz) + az(i) * dt) / (1.0d0 + alphaz)
+        else
+            vx(i) = vx(i) + ax(i) * dt
+            vy(i) = vy(i) + ay(i) * dt
+            vz(i) = vz(i) + az(i) * dt
+        end if
         x(i) = x(i) + vx(i) * dt
         y(i) = y(i) + vy(i) * dt
         z(i) = z(i) + vz(i) * dt
     end do
-
-    ! 隐式阻尼处理（用临时数组 dmp_v，不覆盖 vx/vy/vz）
-    do i = 1, n
-        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then
-            dmp_vx(i) = 0; dmp_vy(i) = 0; dmp_vz(i) = 0; cycle
-        end if
-        gx = B(1) / m(i); gy = B(2) / m(i); gz = B(3) / m(i)
-        dmp_vx(i) = (vx(i) + 0.5d0 * G(1) * dt) / (1.0d0 + 0.5d0 * gx * dt)
-        dmp_vy(i) = (vy(i) + 0.5d0 * G(2) * dt) / (1.0d0 + 0.5d0 * gy * dt)
-        dmp_vz(i) = (vz(i) + 0.5d0 * G(3) * dt) / (1.0d0 + 0.5d0 * gz * dt)
-    end do
-
-    ! 用新位置 + 阻尼速度重算加速度
-    call accel(n, x, y, z, dmp_vx, dmp_vy, dmp_vz, m, G, B, nb, bp, bk, br, &
-               gravity_field, gravity_interaction, &
-               elastic_force, damping_force, gravity_strength, ax, ay, az)
-
-    ! 速度后半步：v = v_half + 0.5*a_next*dt（vx 仍为 v_half）
-    do i = 1, n
-        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
-        vx(i) = vx(i) + ax(i) * dt * 0.5d0
-        vy(i) = vy(i) + ay(i) * dt * 0.5d0
-        vz(i) = vz(i) + az(i) * dt * 0.5d0
-    end do
 end subroutine leapfrog_full
 
 ! ── 分发器 + 边界条件 + 自由度约束 ──
@@ -1066,7 +1097,7 @@ subroutine json_arr(u, vals, n, has_next, indent)
     write(u, '(a)', advance='no') indent // '['
     do i = 1, n
         if (i > 1) write(u, '(a)', advance='no') ','
-        write(u, '(g0)', advance='no') vals(i)
+        write(u, '(g0.8)', advance='no') vals(i)
     end do
     if (has_next) then
         write(u, '(a)') '],'
@@ -1103,7 +1134,7 @@ subroutine write_dbl_arr(u, name, arr, n, has_next)
     write(u, '(a)', advance='no') '  "' // trim(name) // '": ['
     do i = 1, n
         if (i > 1) write(u, '(a)', advance='no') ','
-        write(u, '(g0)', advance='no') arr(i)
+        write(u, '(g0.8)', advance='no') arr(i)
     end do
     if (has_next) then
         write(u, '(a)') '],'

examples/case06/input/bond.txt

@@ -1,2 +1,2 @@
 bond_name k rest_length
-k1 1.0 1.0
+k1 50.0 1.0

examples/case06/input/driver.txt

@@ -1,2 +1,2 @@
 n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period
-1     0     0   2.0      0      0    0.1     0     0    90    all
+1     0     0   2.0      0      0   0.05     0     0    90    all

examples/case06/input/input.txt

@@ -5,18 +5,21 @@
 # ── 流程控制 ──────────────────────────────────
 # 每步用 0/1 单独开关，1=执行，0=跳过
 # 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧
-step_simulate:  1   # 运行物理模拟 → output/trajectory.txt
-step_sample:    1   # 抽帧 → output/display.txt
+step_simulate:  1   # 运行物理模拟 → output/display.txt（引擎直接抽帧）
+step_sample:    0   # （旧版）从 trajectory.txt 重新抽帧，默认0=不执行
 step_plot:      0   # 绘制轨迹/能量图 → output/trajectory_plots.png
 step_animation: 1   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
 step_plot_wave: 0   # 绘制波形能量动画
-force_calc: 0       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
+force_calc: 1       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
 plot_wave_save_gif: 1  # 输出波形 GIF（需 step_plot_wave=1）
 plot_wave_save_mp4: 1  # 输出波形 MP4（需 step_plot_wave=1）
 
+# ── 文件保存 ──────────────────────────────────
+save_trajectory: 0  # 0=不保留完整轨迹文件, 1=保留 trajectory.txt（用于后续单独抽帧）
+
 # ── 计算引擎 ──────────────────────────────────
 # 可选: python, c, cpp, fortran, java
-engine: fortran  # 默认使用 python 引擎
+engine: python  # 默认使用 python 引擎
 
 # ── 盒子 ──────────────────────────────────────
 box_a: 80.0         # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
@@ -63,13 +66,13 @@ warmup_steps: 0           # 默认 0（立即开始记录）
 
 # 总模拟时间（秒），程序自动计算 NT = T_total / DT
 # 如果同时指定了 NT，以 NT 为准
-T_total:  10.0
+T_total:  100.0
 
 # 抽帧间隔（每 NSTEP 步取一帧用于动画）
-NSTEP: 100
+NSTEP: 10
 
 # ── 时间步长 ──────────────────────────────────
-DT:    0.001         # 时间步长 (s)
+DT:    0.01         # 时间步长 (s)
 
 # 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧
 sample_start: null        # null 表示从头开始（帧索引从 0 起）
@@ -85,7 +88,7 @@ use_marker: 1
 
 # ── 显示参数 ──────────────────────────────────
 # 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序
-alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5]
+alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 
 # 小球颜色
 # 小球半径从 coord_file 的 radius 列读取
@@ -97,3 +100,9 @@ ball_color_b: 0.90  # B 分量
 box_color_r:  0.80
 box_color_g:  0.80
 box_color_b:  0.85
+
+# ── 摄像机初始位置 ────────────────────────────
+camera_distance:  40.0   # 摄像机到场景中心的距离
+camera_elevation: 0      # 俯仰角（度），负值=俯视
+camera_azimuth:   0      # 方位角（度）
+move_camera:      1      # 0=固定视角, 1=按 move_camera.txt 运动

examples/case06/input/move_camera.txt

@@ -0,0 +1,9 @@
+# move_camera.txt — 摄像机速度段驱动
+# 格式: start-end  vx=f vy=f vz=f  rx=d ry=d rz=d
+#   vx/vy/vz: 平移速度（每帧移动单位）
+#   rx/ry/rz: 旋转速度（每帧度数）
+#     rx → elevation（俯仰）, ry → azimuth（方位）, rz → (预留)
+#
+# 示例：前60帧向右平移+绕x旋转，30-90帧向上平移+绕y绕z旋转
+all   vx=0.02
+# 30-90  vy=0.02  ry=1  rz=1

examples/case06/run_dynamics.py

@@ -31,7 +31,7 @@ def load_dynamics_module(module_path: Path):
 
 
 def main():
-    parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01")
+    parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case06")
     parser.add_argument("--no-plot", action="store_true", help="跳过 matplotlib 绘图")
     args = parser.parse_args()
 

optimization/claude.md

@@ -0,0 +1,450 @@
+# Dynamics 项目优化建议
+
+> 分析日期：2026-06-12  
+> 分析范围：`D:\Share\Data\aliyun-gitea\dynamics` 完整代码库  
+> 当前算例：`examples/case06`（120 原子一维链横波，NT=100,000 步，C 引擎）  
+> 另见：[`optimization/workbuddy.md`](workbuddy.md)（架构拆分、全局变量封装、测试等宏观建议）
+
+本文档聚焦以下四类问题，不与 workbuddy.md 重复：
+1. **已确认的 Bug**（可能导致运行崩溃或结果错误）
+2. **Python 引擎性能**（向量化等具体代码改动）
+3. **引擎一致性问题**（C/C++/Fortran 行为差异）
+4. **小的代码质量问题**（错误文本、重复代码等）
+
+---
+
+## 一、已确认的 Bug
+
+### B1. `run_simulation()` 内引用了不在作用域的 `config` 变量 ⚠️ 严重
+
+**文件**：[`compute.py:1548-1550`](../compute.py)
+
+```python
+# 问题代码（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` 中赋值：
+
+```python
+# 在模块顶部全局变量区添加（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` 绘图块引用未定义的变量 ⚠️ 严重
+
+**文件**：[`dynamics.py:330-451`](../dynamics.py)
+
+```python
+# 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` 的返回值读取数据重构绘图代码，或暂时添加保护：
+
+```python
+# 临时保护（快速修复）
+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) 和 [`plot_wave.py:96-111`](../plot_wave.py)
+
+```python
+# 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`，并适配字段名：
+
+```python
+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`](../examples/case06/run_dynamics.py)
+
+```python
+# 当前（错误）
+parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01")
+
+# 应为
+parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case06")
+```
+
+---
+
+### B5. `draw.py` 裸 `except` 吞掉所有异常
+
+**文件**：[`draw.py:97`](../draw.py)
+
+```python
+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. 弹簧力计算向量化（最大收益）
+
+**文件**：[`compute.py:1253-1277`](../compute.py)
+
+```python
+# 当前：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
+    ...
+```
+
+**向量化版本**：
+
+```python
+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()` 减少临时数组分配
+
+**文件**：[`compute.py:1408-1418`](../compute.py)
+
+```python
+# 当前：每步创建 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 方向的原子可以运动。
+
+**优化版本**（预先计算掩码，直接原地修改）：
+
+```python
+# 在模块初始化时（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` 列表
+
+**文件**：[`compute.py:1470,1510`](../compute.py)
+
+```python
+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`
+
+**文件**：[`compute.py:630`](../compute.py)
+
+```python
+# 每次调用都分配一个 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,) 数组，直接用标量广播即可：
+
+```python
+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²) 双重循环
+
+**文件**：[`compute.py:1279-1297`](../compute.py)
+
+```python
+# 当前：纯 Python 双重循环
+for i in range(n):
+    for j in range(i + 1, n):
+        ...
+```
+
+**向量化方案**（适用于 N 不太大时）：
+
+```python
+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`](../compute.py) vs [`engines/c/main.c`](../engines/c/main.c)
+
+Python 引擎将整根键的弹性势能 `½k(d-d₀)²` 记给原子 `i`（不均分），而 C/C++ 引擎均分给两端。
+这不影响力的计算（力相同），但如果未来基于势能做分析，两个引擎的势能分配会不同。
+建议注释中明确说明约定。
+
+### E2. C++ 引擎 `save_trajectory` 默认值为 1
+
+**文件**：[`engines/cpp/param.json`](../engines/cpp/param.json)
+
+```json
+// engines/c/param.json:    "save_trajectory": 0
+// engines/cpp/param.json:  "save_trajectory": 1   ← 与 C 不一致
+```
+
+用户切换引擎时行为不同，可能意外生成大文件。建议统一为 0。
+
+### E3. Fortran 引擎不支持 `save_trajectory=0`
+
+**文件**：[`engines/fortran/main.f90`](../engines/fortran/main.f90)
+
+Fortran 引擎总是输出 JSON 格式的 `trajectory.txt`，不支持直接写 `display.txt`。
+workbuddy.md 中已将此列为 P0，此处确认影响范围：
+- case 使用 `engine: fortran` 时必然崩溃（Python 侧期待 `display.txt`）
+- 应参照 C 引擎的 `write_display_txt` 函数实现
+
+### E4. 外部引擎每次运行都重新校准
+
+**文件**：[`compute.py:924-948`](../compute.py)
+
+每次运行前跑 `min(1000, NT//10)` 步测速。对于 `NT=100,000` 的 case06，校准跑 10,000 步，
+占总运算量的 10%。
+
+**快速修复**（跳过短运行的校准）：
+
+```python
+# 当总步数 < 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`](../dynamics.py) 的 `_load_camera_kf()` 与
+[`compute.py:90-136`](../compute.py) 的 `_load_camera_motion()` 是功能完全相同的函数，
+仅变量名略有差异。
+
+**建议**：保留 `compute.py` 中的版本，`dynamics.py` 改为从 `compute` 导入：
+
+```python
+# dynamics.py
+from compute import _load_camera_motion as _load_camera_kf
+```
+
+### Q2. `load_parameters()` 函数已废弃但仍保留
+
+**文件**：[`compute.py:1123-1199`](../compute.py)
+
+这是老版通过 `key = value` 文本格式读参数的函数，现在所有案例已改用 YAML。
+只有 `compute.py` 的 `main()` 还调用它，而 `compute.main()` 本身也已被 `dynamics.py` 取代。
+
+整个函数（77 行）加上 `main()`（25 行）都可以安全删除，或移入 `tools/` 目录存档。
+
+### Q3. 边界条件逻辑不一致
+
+**文件**：[`compute.py:1380-1387`](../compute.py) vs [`compute.py:1420-1428`](../compute.py)
+
+```python
+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`](../examples/case06/run_dynamics.py)
+
+见 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 |
+
+**最高回报的三件事（按投入产出比排序）**：
+
+1. **B3**（`plot_wave.py` 格式不匹配）：修复后波形动画功能恢复，影响 case05/case06
+2. **B1**（`config` 未定义）：修复后 Python 引擎可完整运行 case06
+3. **P1**（弹簧力向量化）：修复后 Python 引擎速度提升 5-15 倍，适合大规模测试和教学演示
+
+---
+
+*本文档由 Claude Sonnet 4.6 自动生成，基于代码静态分析，未实际运行测试验证。建议在实施前先运行相关案例确认 Bug 现象。*

optimization/claude_v1.md

@@ -0,0 +1,459 @@
+# Dynamics 项目优化方案综合报告
+
+> 整合日期：2026-06-12  
+> 来源文档：[`claude.md`](claude.md)（Claude Sonnet 4.6 分析）、[`workbuddy.md`](workbuddy.md)（WorkBuddy 分析）  
+> 项目路径：`D:\Share\Data\aliyun-gitea\dynamics`
+
+---
+
+## 第一部分：两份分析报告的评价与对比
+
+### 1.1 WorkBuddy 分析报告评价（`workbuddy.md`）
+
+#### 优势
+
+**① 架构视野宏观、分层清晰**  
+WorkBuddy 从整体软件工程角度出发，明确指出 `compute.py` 的职责混乱问题（1618 行同时承担物理引擎、文件 I/O、参数加载、外部引擎管理五种职责），并给出了完整的模块拆分方案：
+
+```
+compute/
+├── core.py        # 物理引擎
+├── io.py          # 文件 I/O
+├── params.py      # 参数加载
+├── runner.py      # 运行管理
+├── engine_helper.py  # 外部引擎
+└── main.py        # 主入口
+```
+
+这是一个系统性重构方向，有助于长期维护。
+
+**② 关注工程规范与测试**  
+WorkBuddy 专门开辟了"测试与质量"一节，强调了缺少单元测试的风险，并给出了类型标注的示例代码——这两点是 claude.md 未涉及的。对于一个教学/研究用项目，有测试才能安全重构。
+
+**③ 配置管理建议具体**  
+提出了将 `ball_color_r/g/b` 三个键合并为 `ball_color: [r,g,b]` 数组的配置统一方案，以及6个案例 `input.txt` 格式不统一的问题（如 `save_trajectory`、`camera_*` 字段缺失），并给出了统一模板，实操性强。
+
+**④ 引擎一致性梳理全面**  
+清晰列出了 C/C++/Fortran 三引擎在 `save_trajectory`、`box_a` 默认值等参数上的差异表格，对用户切换引擎有直接指导价值。
+
+**⑤ 区分"已完成"与"待做"**  
+文档中标注了"✅ 已修复"、"✅ 已完成"条目，帮助读者快速定位当前状态，避免重复分析。
+
+#### 劣势
+
+**① 缺少具体 Bug 定位**  
+WorkBuddy 的分析主要停留在架构和规范层面，对于已经导致运行崩溃的 Bug（如 `plot_wave.py` 格式不匹配、`run_simulation` 中 `config` 未定义）没有识别和标注。用户按此文档操作时，可能先花时间做重构，但基础功能仍然崩溃。
+
+**② 性能优化建议较抽象**  
+对 Python 引擎性能仅建议"使用 Numba JIT"或"numpy.vectorize"，缺少具体的向量化代码示例。文档注明"Python vs C 慢约 6-8 倍"，但未分析哪个函数是瓶颈（实际上是弹簧力的 Python for 循环）。
+
+**③ 部分建议颗粒度不足**  
+"消除全局变量"建议封装为 `SimulationState` 数据类，但未说明：封装后接口如何调整、外部引擎的 `param.json` 协议是否需要同步更新、draw.py 的全局变量是否属于同一重构范围。
+
+**④ 优先级排列有待商榷**  
+将"全局变量封装"列为 P0（与 Fortran 引擎修复同级），但封装全局变量是一项重构工作（3-4h），而 Fortran 引擎崩溃是已存在的功能阻断性问题，两者紧迫性不同。
+
+---
+
+### 1.2 Claude 分析报告评价（`claude.md`）
+
+#### 优势
+
+**① Bug 定位精准，附有复现条件**  
+每个 Bug 都说明了触发条件（"只要使用 `engine: python`"、"`step_plot_wave: 1` 时"），以及为什么当前没有爆发（"当前所有案例恰好设置 `step_plot: 0`"）。这对于开发者复现和验证问题极有价值。
+
+**② 代码级修复方案完整**  
+每个 Bug 都给出了完整可运行的修复代码，包括需要修改的行号、修改前后的对比。特别是 B3（`plot_wave.py` 格式不匹配）不仅指出问题，还解释了为何修复较复杂（`atom_masses` 等物理量在新格式中缺失，需要同步扩展 header 字段），给出了分层的修复路径（快速修复 vs 彻底修复）。
+
+**③ 弹簧力向量化方案具体可行**  
+给出了完整的 NumPy 向量化代码，包括处理重复索引的 `np.add.at` 用法，以及链状体系可进一步优化为直接索引的提示。预期加速 5-15 倍的估算也有根据（119 键 × 100,000 步的量级分析）。
+
+**④ 引擎差异的量化说明**  
+指出了 `engines/cpp/param.json` 中 `save_trajectory` 默认值为 1（与 C 引擎的 0 不一致），并说明了实际影响（切换引擎时可能意外生成大文件）。
+
+#### 劣势
+
+**① 覆盖面刻意与 workbuddy.md 错开，造成视野盲区**  
+claude.md 开篇声明"不与 workbuddy.md 重复"，因此主动跳过了架构、测试、配置管理等方向。但这导致文档读者如果只看 claude.md，会误以为架构是没问题的。两份文档需要结合阅读，割裂感较强。
+
+**② 部分建议依赖全局变量方案**  
+claude.md 提出的 Bug B1 修复方案（将 `camera_distance` 等加入全局变量）是在全局变量模式下的补丁式修复，与 workbuddy.md 建议的"封装全局变量为类"方向相反。若未来执行 workbuddy.md 的重构，B1 的修复需要再次调整。
+
+**③ 某些结论尚未验证**  
+文档末尾注明"基于代码静态分析，未实际运行测试验证"。例如 B3（`plot_wave.py` 格式不匹配）的结论基于函数调用链分析，如果 `load_text_data` 内部有容错逻辑，实际情况可能不同。
+
+**④ 缺乏对 draw.py 的深入分析**  
+draw.py 是用户交互最频繁的模块（3D 动画播放），但 claude.md 仅指出了一处裸 `except`（B5），未分析动画帧率优化、大原子数渲染性能等用户体验层面的问题。
+
+---
+
+### 1.3 两份报告对比总结
+
+| 维度 | WorkBuddy | Claude |
+|------|-----------|--------|
+| **Bug 识别** | ❌ 未发现运行时 Bug | ✅ 识别 5 个确认 Bug，附触发条件 |
+| **性能优化** | 🟡 方向正确但缺代码 | ✅ 完整向量化代码示例 |
+| **架构设计** | ✅ 完整模块拆分方案 | ❌ 刻意回避，未覆盖 |
+| **测试建议** | ✅ pytest 方案具体 | ❌ 未涉及 |
+| **配置管理** | ✅ 6 案例统一模板 | ❌ 未涉及 |
+| **引擎一致性** | ✅ 参数差异表格 | ✅ 补充了势能归属约定 |
+| **修复代码** | 🟡 仅有架构示意 | ✅ 逐条给出可运行代码 |
+| **优先级合理性** | 🟡 P0 安排有误 | ✅ 按紧迫性分层清晰 |
+| **可独立阅读** | ✅ 自成体系 | ❌ 依赖读者已读 workbuddy.md |
+| **当前状态标注** | ✅ 已完成项有标注 | ❌ 未区分 |
+
+**结论**：WorkBuddy 擅长系统性架构分析和工程规范，Claude 擅长 Bug 精确定位和具体修复代码。两份报告互补，单独使用任何一份都存在明显盲区。
+
+---
+
+## 第二部分：综合优化方案
+
+### 原则
+
+基于以上分析，综合方案遵循以下原则：
+1. **先止血再手术**：运行时 Bug 优先于架构重构，不能让教学演示功能处于崩溃状态
+2. **快速收益优先**：同等重要性时，工作量小的先做
+3. **向量化先于重构**：Python 引擎提速不依赖架构改动，可独立进行
+4. **渐进式重构**：全局变量封装等重构分阶段进行，避免一次性大改引入新 Bug
+
+---
+
+### 阶段一：紧急修复（总工作量约 4 小时）
+
+> 目标：恢复所有已有功能到正常可用状态
+
+#### 1. 修复 `plot_wave.py` 格式不匹配（B3）— 1.5h
+
+**根本原因**：`display.txt` 格式从 JSON 迁移到新文本格式时，`plot_wave.py` 未同步更新。
+
+**修复步骤**：
+
+① 修改加载函数（`plot_wave.py:22-27`）：
+```python
+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)   # 改为新格式加载
+```
+
+② 在 `save_display_txt`（`compute.py:217`）的 header 中补充 `atom_masses` 字段：
+```python
+# 在 header_fields 中添加
+"atom_masses": ",".join(str(m) for m in ATOM_MASSES),
+```
+
+③ 修改 `plot_wave.py` 中的字段访问：
+```python
+# 旧 → 新
+data["n_frames"]       → disp_data["frames_x"].shape[0]
+data["disp_all_x"]     → disp_data["frames_x"]
+data["atom_masses"]    → np.array([float(x) for x in h.get("atom_masses","").split(",") if x])
+data["atom_ids"]       → disp_data["atom_ids"]
+data["bond_pairs"]     → []   # display.txt 不含成键信息，能量计算中弹性势能项跳过
+```
+
+**验证**：运行 case05（`step_plot_wave: 1`），确认生成 `wave_animation.gif`。
+
+---
+
+#### 2. 修复 `run_simulation()` 中的 `config` 未定义（B1）— 15min
+
+在 `compute.py` 模块顶部（约第 68 行，紧随 `camera_keyframes_raw`）添加全局变量：
+
+```python
+camera_keyframes_raw = ""
+camera_distance  = 40.0    # 新增
+camera_elevation = 0       # 新增
+camera_azimuth   = 0       # 新增
+```
+
+在 `run_from_config`（约第 756 行，`use_marker` 赋值附近）添加：
+```python
+use_marker = int(config.get("use_marker", 0))
+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_from_config` 的 `global` 声明行追加这三个变量：
+```python
+global use_marker, camera_keyframes_raw, camera_distance, camera_elevation, camera_azimuth
+```
+
+在 `run_simulation`（第 1548-1550 行）改为读全局变量：
+```python
+"camera_distance": str(camera_distance),
+"camera_elevation": str(camera_elevation),
+"camera_azimuth": str(camera_azimuth),
+```
+
+**验证**：用 `engine: python` 运行 case01，确认生成正确的 `display.txt`。
+
+---
+
+#### 3. 修复 `dynamics.py` 绘图块死代码（B2）— 30min
+
+`step_plot` 功能依赖完整轨迹数据，但当前架构中 `run_from_config` 不再在主流程中保留这些数据。
+快速修复：添加 "功能暂不可用" 保护，避免用户误开后崩溃：
+
+```python
+# dynamics.py 约 322 行
+if not no_plot and config.get("step_plot", 1):
+    print("[run] 注意：step_plot 绘图功能需要 save_trajectory=1 时的完整轨迹，")
+    print("[run]       当前版本暂未支持，已跳过。如需绘图请联系开发者。")
+```
+
+**长期修复**（阶段三再做）：从 `display.txt` 读取抽帧数据重建绘图逻辑。
+
+---
+
+#### 4. 修复 Fortran 引擎输出格式（E3）— 2h
+
+Fortran 引擎总是输出 JSON 格式 `trajectory.txt`，与 Python 侧期待的 `display.txt` 不兼容，
+使用 `engine: fortran` 时必然崩溃。
+
+参照 C 引擎的 `write_display_txt` 函数（`engines/c/main.c`）在 Fortran 中实现：
+1. 读取 `param.json` 中的 `save_trajectory` 参数
+2. 新增 `write_display_txt_f90` 子程序，按 NSTEP 抽帧写文本格式
+3. 根据 `save_trajectory` 决定是否额外写 `trajectory.txt`
+
+**同时**：将 `engines/cpp/param.json` 中的 `"save_trajectory": 1` 改为 `0`，与 C 保持一致。
+
+---
+
+#### 5. 其他 1-5 分钟的小修复
+
+```python
+# B4: examples/case06/run_dynamics.py:34
+description="运行 Dynamics 示例案例 case06"   # 改 case01 → case06
+
+# B5: draw.py:97
+except (ValueError, AttributeError):           # 改裸 except
+
+# Q1: dynamics.py 顶部
+from compute import _load_camera_motion as _load_camera_kf   # 删除重复实现
+
+# E2: engines/cpp/param.json
+"save_trajectory": 0                            # 统一默认值
+```
+
+---
+
+### 阶段二：性能优化（总工作量约 3 小时）
+
+> 目标：Python 引擎提速 5-10 倍，外部引擎减少启动开销
+
+#### 6. 弹簧力向量化（P1）— 1h
+
+将 `compute.py:1253-1277` 的 Python for 循环替换为 NumPy 向量化版本（见 claude.md §P1 完整代码）。
+
+**链状体系（case01-06）的进一步优化**：  
+由于链状体系每根键的两端没有重复（原子 i 和 j 各只出现一次），可以用更快的直接索引替代 `np.add.at`：
+
+```python
+# 比 np.add.at 快约 3 倍（无冲突索引时）
+fx[i_arr] += fx_bond
+fx[j_arr] -= fx_bond
+fy[i_arr] += fy_bond
+fy[j_arr] -= fy_bond
+fz[i_arr] += fz_bond
+fz[j_arr] -= fz_bond
+```
+
+对于有分叉键的复杂体系仍需用 `np.add.at`，可根据 `BOND_PAIRS` 检测是否有重复端点自动选择路径。
+
+---
+
+#### 7. `apply_fixed_constraints` 预计算掩码（P2）— 30min
+
+将 `compute.py:1408-1418` 的 `column_stack` 方案替换为预计算 bool 掩码，见 claude.md §P2 完整代码。  
+对 case06（120 原子，x/y 全固定）每步节省 2 次 `(120, 3)` 数组分配。
+
+---
+
+#### 8. 外部引擎校准缓存（E4）— 1h
+
+在 `compute.py:924` 校准逻辑前添加缓存命中检查：
+
+```python
+import hashlib, json as _json
+
+def _calib_cache_key(engine, n_atoms, method, dt):
+    s = f"{engine}:{n_atoms}:{method}:{dt}"
+    return hashlib.md5(s.encode()).hexdigest()[:8]
+
+_calib_cache_path = os.path.join(script_dir, "engines", engine, "_calib_cache.json")
+_calib_key = _calib_cache_key(engine, len(ATOM_IDS), config.get("method","leapfrog"), float(config["DT"]))
+_cached = {}
+if os.path.exists(_calib_cache_path):
+    with open(_calib_cache_path) as _cf:
+        _cached = _json.load(_cf)
+
+if _cached.get("key") == _calib_key and time.time() - _cached.get("ts", 0) < 86400:
+    # 缓存命中（1天内有效），跳过校准
+    _step_time = _cached["step_time"]
+    _overhead  = _cached["overhead"]
+else:
+    # 执行校准（现有逻辑）...
+    # 校准完成后写缓存
+    with open(_calib_cache_path, "w") as _cf:
+        _json.dump({"key": _calib_key, "ts": time.time(),
+                    "step_time": _step_time, "overhead": _overhead}, _cf)
+```
+
+**预期收益**：case06 第二次运行节省约 10 秒（校准 10,000 步）。
+
+---
+
+#### 9. 其他微优化（P3/P4）— 30min
+
+```python
+# P3: run_simulation 中删除 frame_indices 列表
+# 删除: frame_indices = []  和  frame_indices.append(step)
+# 替换: n_frames_actual = record_steps // NSTEP  （直接计算）
+
+# P4: apply_driving_force 中删除 t_vec 数组创建
+t_vec = np.array([t, t, t], dtype=np.float64)    # 删除此行
+pos_drive = d["amp"] * np.cos(2.0 * np.pi * d["freq"] * t + d["phi"])   # 直接用标量 t
+```
+
+---
+
+### 阶段三：架构重构（总工作量约 15-20 小时）
+
+> 目标：提升长期可维护性，支持并发模拟、单元测试
+
+> **重要说明**：此阶段是 workbuddy.md 建议的核心，以下是对其建议的细化和补充。  
+> 建议在阶段一、二完成并通过全案例验证后再开始。
+
+#### 10. 全局变量封装为 `SimulationState`（workbuddy.md §3.1）— 3-4h
+
+封装后，阶段一中 B1 的补丁修复（新增三个全局变量）可以一并纳入 `SimulationState`，不需要额外保留。
+
+封装时注意：外部引擎路径（`run_engine`）通过 `param.json` 传参，不受全局变量影响，可独立进行。
+
+#### 11. `compute.py` 模块拆分（workbuddy.md §1.1）— 4-6h
+
+拆分时建议**先拆 io.py**（文件读写），因为它与物理算法完全解耦，风险最低。  
+拆分顺序建议：`io.py` → `params.py` → `core.py` → `runner.py`。
+
+**补充**：同步删除已废弃的 `load_parameters()` 函数（102 行）和已被 `dynamics.py` 取代的 `compute.main()`（25 行），见 claude.md §Q2。
+
+#### 12. 添加单元测试（workbuddy.md §5.1）— 3-5h
+
+优先覆盖以下三类：
+1. 物理算法正确性：以简谐振子解析解验证 leapfrog、midpoint 方法（参照 `tools/NumericalMethods.py` 已有实现）
+2. 文件 I/O 往返一致性：`save_display_txt → load_display_txt` 数值不变
+3. 参数验证边界：质量为 0、键长为负、NT 为 0 等异常输入
+
+#### 13. 案例 `input.txt` 格式统一（workbuddy.md §4.1）— 1-2h
+
+将 `case01-05` 的 `input.txt` 补充缺失字段（`save_trajectory`、`camera_*`、`move_camera`），
+统一 `step_sample: 0`（新版引擎已内置抽帧），清理旧格式注释。
+
+#### 14. `dynamics.py` 绘图功能完整重建（B2 长期修复）— 2h
+
+从 `display.txt` 的抽帧数据重建绘图逻辑（替代依赖完整轨迹的旧实现）：
+
+```python
+disp_data = compute.load_display_txt(disp_path)
+frames_x = disp_data["frames_x"]   # (n_frames, n_atoms)
+# 对每帧计算能量，绘制时序图
+```
+
+能量计算需要质量等物理量，依赖阶段一 B3 修复中在 header 补充的 `atom_masses` 字段。
+
+---
+
+### 三个引擎的综合评价
+
+#### C 引擎（`engines/c/main.c`）
+
+| 维度 | 评价 |
+|------|------|
+| **功能完整性** | ✅ 最完整，支持 `display.txt` 直接输出，`save_trajectory` 控制 |
+| **性能** | ✅ 基准引擎，比 Python 快 6-8 倍 |
+| **一致性** | ✅ `save_trajectory` 默认值为 0，与 Python 行为一致 |
+| **代码质量** | 🟡 JSON 解析为自行实现（非 cJSON 库），维护成本高 |
+| **建议** | 作为参考实现，其他引擎向 C 引擎看齐 |
+
+#### C++ 引擎（`engines/cpp/main.cpp`）
+
+| 维度 | 评价 |
+|------|------|
+| **功能完整性** | ✅ 功能与 C 引擎相当 |
+| **性能** | 🟡 编译产物 3.3 MB（远大于 C 的 504 KB），有优化空间 |
+| **一致性** | ⚠️ `save_trajectory` 默认值为 1（与 C 不一致），已在阶段一修复 |
+| **代码质量** | 🟡 物理算法与 C 引擎重复，应提取到公共头文件 |
+| **建议** | 统一 `param.json` 默认值；考虑与 C 引擎共用物理算法头文件 |
+
+#### Fortran 引擎（`engines/fortran/main.f90`）
+
+| 维度 | 评价 |
+|------|------|
+| **功能完整性** | ❌ 不支持 `display.txt` 直接输出，使用时必然崩溃 |
+| **性能** | ✅ 编译性能与 C 相当，适合大规模计算 |
+| **一致性** | ❌ 输出 JSON 格式 trajectory.txt，与其他引擎行为完全不同 |
+| **代码质量** | 🟡 代码结构清晰，但与 C/C++ 引擎物理算法重复 |
+| **建议** | 阶段一 E3 修复为最高优先级；修复后可成为高性能替代引擎 |
+
+---
+
+### 综合优先级总表
+
+| 阶段 | 编号 | 问题 | 来源 | 优先级 | 预估工时 |
+|------|------|------|------|--------|---------|
+| 一 | B3 | `plot_wave.py` 格式不匹配（波形动画失效） | Claude | 🔴 P0 | 1.5h |
+| 一 | E3 | Fortran 引擎不支持 display.txt | 两者 | 🔴 P0 | 2h |
+| 一 | B1 | `run_simulation` 内 `config` 未定义 | Claude | 🔴 P0 | 15min |
+| 一 | B2 | `dynamics.py` 绘图块死代码 | Claude | 🔴 P0 | 30min |
+| 一 | B4/B5 | 小错误修复（case06 文字、裸 except 等） | Claude | 🟡 P1 | 15min |
+| 一 | Q1/E2 | 重复函数去除、C++ 默认值统一 | Claude | 🟡 P1 | 20min |
+| 二 | P1 | 弹簧力向量化（Python 引擎 5-15x） | Claude | 🟢 P1 | 1h |
+| 二 | E4 | 外部引擎校准缓存 | 两者 | 🟢 P1 | 1h |
+| 二 | P2-P4 | apply_fixed_constraints 等微优化 | Claude | 🟢 P2 | 1h |
+| 三 | 架构 | 全局变量封装为 SimulationState | WorkBuddy | 🔵 P2 | 3-4h |
+| 三 | 架构 | compute.py 模块拆分 | WorkBuddy | 🔵 P2 | 4-6h |
+| 三 | 测试 | 添加 pytest 单元测试 | WorkBuddy | 🔵 P2 | 3-5h |
+| 三 | 配置 | 6 案例 input.txt 格式统一 | WorkBuddy | 🔵 P3 | 1-2h |
+| 三 | 绘图 | dynamics.py 绘图功能重建 | Claude | 🔵 P3 | 2h |
+| 三 | 引擎 | C/C++/Fortran 共用物理算法头文件 | WorkBuddy | 🔵 P3 | 2-4h |
+
+---
+
+## 附录：快速执行清单
+
+### 今天可以完成（总计约 30 分钟的小修复）
+
+```
+[ ] examples/case06/run_dynamics.py:34  "case01" → "case06"
+[ ] draw.py:97  裸 except → except (ValueError, AttributeError):
+[ ] dynamics.py  删除 _load_camera_kf()，改为从 compute 导入
+[ ] engines/cpp/param.json  "save_trajectory": 1 → 0
+[ ] compute.py:1548-1550  config.get → 全局变量（B1 修复）
+```
+
+### 本周可以完成（核心功能恢复）
+
+```
+[ ] plot_wave.py 格式适配（B3）
+[ ] dynamics.py 绘图块保护（B2 临时修复）
+[ ] Fortran 引擎 display.txt 支持（E3）
+[ ] 弹簧力向量化（P1）
+[ ] 引擎校准缓存（E4）
+```
+
+### 下阶段规划（架构改善）
+
+```
+[ ] 全局变量 → SimulationState 类
+[ ] compute.py 拆分为子模块
+[ ] pytest 单元测试套件
+[ ] 案例 input.txt 格式统一
+```
+
+---
+
+*本文档综合 Claude Sonnet 4.6（`claude.md`）与 WorkBuddy（`workbuddy.md`）的分析报告生成。  
+所有代码改动建议均基于静态分析，实施前请先在独立分支验证。*

optimization/codex.md

@@ -0,0 +1,351 @@
+# Dynamics 项目优化建议
+
+本文基于对以下文件的静态检查整理：
+
+- `compute.py`
+- `dynamics.py`
+- `draw.py`
+- `engines/c/main.c`
+- `README.md`
+
+目标不是泛泛地“提速”，而是优先找出这个项目当前最可能影响性能、内存占用、可维护性和后续扩展效率的点，并给出按优先级排序的改进方向。
+
+## 一、整体判断
+
+这个项目现在的架构已经有一个很好的基础：`dynamics.py` 做流程编排，`compute.py` 做 Python 参考实现，`engines/c/main.c` 提供高性能引擎，`draw.py` 负责可视化。
+
+从代码结构看，当前主要瓶颈不在“外层调度”，而集中在三类地方：
+
+1. `compute.py` 中逐步积分时的 Python 层循环。
+2. `compute.py` 中弹簧力/粒子间引力的双层或逐键循环。
+3. `display.txt` / `trajectory.txt` 的文本格式读写。
+
+如果后续案例规模继续增大，Python 参考引擎会很快被这三处放大；而且即便切到 C 引擎，文本 I/O 仍然会成为新的主要瓶颈。
+
+## 二、最高优先级建议
+
+## 1. 先把“计算核心”和“输出格式”分开优化
+
+当前项目已经有多语言引擎，这是很正确的方向。建议把优化目标拆成两条线分别推进：
+
+- 计算性能：优先优化 `compute.py` 中的力计算和积分主循环。
+- I/O 性能：优先替换或补充 `display.txt` / `trajectory.txt` 的文本格式。
+
+原因：
+
+- 计算核心的瓶颈主要是 CPU。
+- 轨迹输出的瓶颈主要是字符串格式化、磁盘写入、加载解析。
+- 这两类优化手段完全不同，混在一起做容易看不出收益来源。
+
+建议先增加一个简单的 profiling 开关，例如输出：
+
+- 总模拟时间
+- 力计算时间
+- 抽帧时间
+- 保存 `display.txt` 时间
+- 保存 `trajectory.txt` 时间
+- 加载 `display.txt` 时间
+
+这样后面每做一次优化，都能知道收益来自哪一段。
+
+## 2. 优先向量化 `compute_force()` 里的弹簧键计算
+
+`compute.py:1230` 开始的 `compute_force()` 是 Python 引擎最关键的热点。
+
+其中弹簧键部分目前是：
+
+- 遍历 `BOND_PAIRS`
+- 每根键单独算 `dx/dy/dz`
+- 每根键单独回写两个原子的受力
+
+这在键数增大后会变成明显瓶颈。
+
+建议做法：
+
+- 把 `BOND_PAIRS[:, 0]` 和 `BOND_PAIRS[:, 1]` 拆成两个索引数组 `i_idx`、`j_idx`。
+- 用 NumPy 一次性计算所有键的 `dx/dy/dz/dist/stretch`。
+- 用 `np.add.at()` 或等价的 scatter 累加方式，把键力一次性回写到 `fx/fy/fz`。
+
+预期收益：
+
+- 中等规模体系下，Python 参考引擎会有很可观的提速。
+- 这也是最不改变项目结构、最容易验证正确性的优化。
+
+## 3. 把原子间万有引力视为“可选高成本模块”，不要默认走 Python 双层循环
+
+`compute.py:1281` 开始的 `GRAVITY_INTERACTION` 采用 `for i` / `for j` 双层循环，复杂度是 `O(N^2)`。
+
+这段代码在教学上没问题，但在稍大的粒子数下会非常慢。建议：
+
+- 明确把它标记为“仅适合小规模案例”。
+- 默认推荐用户用 `engine: c` 跑这类 case。
+- 如果未来仍想保留 Python 版本，优先考虑：
+  - 小规模时保持现状。
+  - 中大规模时改成 Numba 或 C 扩展。
+  - 更进一步再考虑 Barnes-Hut、cell list、neighbor list 之类近似/加速算法。
+
+如果只是当前项目阶段，我更建议：
+
+- 不要先在 Python 里硬做复杂天体优化。
+- 先把这类高复杂度场景明确导向 C 引擎。
+
+这是投入产出比更高的路线。
+
+## 4. `display.txt` 文本格式要补一个二进制版本
+
+`compute.py:196` 的 `save_display_txt()` 和 `compute.py:239` 的 `load_display_txt()` 已经比 JSON 好很多，但本质仍然是大文本。
+
+当前问题：
+
+- 保存时每个数都在做字符串格式化。
+- 读取时虽然用了 `np.genfromtxt`，但源数据仍是文本。
+- 帧数和粒子数继续增大后，I/O 会成为明显瓶颈。
+
+建议保留 `display.txt` 作为“人可读格式”，同时新增一个高性能格式，例如：
+
+- `display.npz`
+- 或 `display.npy` + `meta.json`
+
+推荐方案：
+
+- `frames_x/y/z/vx/vy/vz` 存到 `np.savez_compressed`
+- 元信息单独存一个轻量 `meta.json`
+- `draw.py` 优先读二进制，不存在时再回退到 `display.txt`
+
+这样有几个好处：
+
+- 调试和教学仍可保留文本文件。
+- 大规模运行时可以直接避开文本解析开销。
+- Python 和 C 引擎都可以逐步迁移，不需要一次切完。
+
+## 三、中优先级建议
+
+## 5. `run_simulation()` 里完整轨迹缓存要改成“按需流式写出”
+
+`compute.py:1481` 起，如果 `save_trajectory=1`，会一次性申请：
+
+- `traj_x`
+- `traj_y`
+- `traj_z`
+- `traj_vx`
+- `traj_vy`
+- `traj_vz`
+
+这意味着内存复杂度接近 `O(steps * atoms)`，而且还是 6 份 `float64` 数组。
+
+这在教学小案例里没问题，但一旦：
+
+- `NT` 很大
+- 粒子数上来
+- 同时还保留抽帧缓存
+
+内存会膨胀得很快。
+
+建议：
+
+- 如果只是为了最终导出，改成边算边写。
+- 如果还需要后续随机访问，可以改成 `memmap` 或分块写入 `npz/hdf5/zarr`。
+
+推荐优先级：
+
+- 先做“保存完整轨迹时改为 chunked binary writer”。
+- `trajectory.txt` 保留为兼容输出，而不是默认主输出。
+
+## 6. `draw.py` 的 Marker 更新应一次性切片赋值，避免逐原子 Python 循环
+
+`draw.py:517` 附近的 `_update_atom_positions()` 在 `USE_MARKER` 模式下仍然逐原子循环：
+
+- 先 `for i in range(N_ATOMS)`
+- 再逐个写 `marker_pos[i]`
+- 然后 `balls.set_data(pos=marker_pos)`
+
+建议直接改成：
+
+- `marker_pos[:, 0] = DISP_ALL_X[f_idx]`
+- `marker_pos[:, 1] = DISP_ALL_Y[f_idx]`
+- `marker_pos[:, 2] = DISP_ALL_Z[f_idx]`
+
+这样更符合 Marker 模式“批量更新”的初衷。
+
+同理，成键线 `_update_bond_positions()` 也可以进一步尝试批量索引构造，而不是逐键更新。
+
+虽然这部分通常不如计算核心慢，但在大粒子数动画里会直接影响帧率。
+
+## 7. `apply_fixed_constraints()` 每步构造 `column_stack`，可以改成原地掩码写回
+
+`compute.py:1408` 的实现每步都会：
+
+- `np.column_stack((x, y, z))`
+- `np.column_stack((vx, vy, vz))`
+- `np.where(...)`
+
+这会产生额外临时数组。
+
+建议改成：
+
+- 预先缓存 `fixed_x/fixed_y/fixed_z` 三个布尔掩码
+- 直接对 `x/y/z/vx/vy/vz` 原地赋值
+
+例如思路上改成：
+
+- `x[fixed_x] = ATOM_POSITIONS[fixed_x, 0]`
+- `vx[fixed_x] = 0.0`
+
+这类优化单项收益不一定最大，但由于它在每一步都执行，累计下来是有价值的。
+
+## 8. `apply_driving_force()` 有重复小数组分配
+
+`compute.py:599` 的驱动力逻辑中，每个 driver、每一步都构造：
+
+- `t_vec = np.array([t, t, t], dtype=np.float64)`
+
+这是典型的小对象重复分配。
+
+建议：
+
+- 直接分别计算三个轴，不需要生成 `t_vec`
+- 或预先把 driver 参数整理成矩阵，批量更新受驱原子
+
+如果 driver 数量很少，这不是最大瓶颈；但这类微优化很容易做，而且不会增加复杂度。
+
+## 四、架构与可维护性建议
+
+## 9. 减少 `compute.py` 的全局变量依赖，逐步收敛到状态对象
+
+目前 `compute.py` 依赖大量全局变量，例如：
+
+- `ATOM_IDS`
+- `ATOM_MASSES`
+- `BOND_PAIRS`
+- `METHOD`
+- `NT`
+- `DT`
+
+这让代码在以下场景下会越来越难维护：
+
+- 并行运行多个 case
+- 写单元测试
+- 替换不同 force model
+- 将来做 GUI 或服务化封装
+
+建议中期做一个 `SimulationState` / `SimulationConfig` / `SystemData` 分层：
+
+- 配置类：步长、方法、开关、输出选项
+- 系统类：原子、键、边界、驱动参数
+- 状态类：当前 `x/y/z/vx/vy/vz`
+
+不需要一次性重构完，先从最核心的 `compute_force()`、`run_simulation()` 入手即可。
+
+## 10. Python 参考引擎和 C 引擎的“功能等价层”建议更明确
+
+README 里已经强调了多语言对比，这是项目亮点。为了后续更稳，建议把“等价层”标准化：
+
+- 相同输入
+- 相同输出语义
+- 相同采样规则
+- 相同边界行为
+- 相同 driver 行为
+
+然后做一个最小一致性测试集：
+
+- case01-case06 都能跑
+- Python 与 C 输出在容差内一致
+- 不同 method 的结果有基准对照
+
+这样后续做任何优化时，都更容易大胆改，不怕悄悄改坏物理行为。
+
+## 五、需要尽快处理的正确性/维护风险
+
+这些不一定直接是“性能问题”，但会影响后续优化效率，建议优先修一下。
+
+## 11. `dynamics.py` 的绘图分支存在明显变量依赖不完整的风险
+
+在 `dynamics.py:331` 一带，绘图代码直接使用：
+
+- `all_x`
+- `all_y`
+- `all_z`
+- `all_vx`
+- `all_vy`
+- `all_vz`
+- `data`
+
+但从当前文件上下文看，这些变量只在某些分支里才会存在，尤其 Python 引擎路径下很可能未定义。
+
+这会带来两个问题：
+
+- 某些 case 可能直接在绘图阶段报错。
+- 优化时很难判断问题来自性能还是流程分支。
+
+建议：
+
+- 在进入绘图逻辑前统一构造标准数据对象。
+- 不要依赖分支里“顺便留下来的局部变量”。
+
+## 12. `README.md` 描述与当前实现已经有部分不一致
+
+例如 README 仍强调：
+
+- `trajectory.txt (JSON, 统一格式)`
+- `sample.py -> display.txt`
+
+但当前实现里：
+
+- Python 路径已经直接写 `display.txt`
+- `save_trajectory` 变成可选
+- `sample.py` 在主流程中的角色已经下降
+
+文档不一致本身不会拖慢程序，但会拖慢后续协作和排障效率，尤其在你继续演进输出格式时会更明显。
+
+建议在做 I/O 优化时顺手更新 README，避免认知分叉。
+
+## 六、推荐实施顺序
+
+如果按“最少改动获得最大收益”的原则，我建议这样排：
+
+1. 给 `dynamics.py` / `compute.py` 增加基础计时统计。
+2. 向量化 `compute_force()` 的弹簧键计算。
+3. 把 `draw.py` 的 Marker 更新改成切片赋值。
+4. 新增 `display.npz`，`draw.py` 优先读取二进制。
+5. 把 `save_trajectory=1` 改成分块二进制输出，而不是全量内存缓存。
+6. 修复 `dynamics.py` 绘图分支的数据来源问题。
+7. 再考虑 `compute.py` 全局变量收敛和更深层的结构重构。
+
+## 七、如果只做三件事，最值得做什么
+
+如果你现在只想投入一小轮精力，我建议只做这三项：
+
+1. 向量化 `compute_force()` 的弹簧键部分。
+2. 增加 `display.npz` 二进制输出与读取。
+3. 修正 `dynamics.py` 绘图阶段的数据流一致性。
+
+原因：
+
+- 第 1 项直接优化 Python 引擎核心热点。
+- 第 2 项直接优化所有引擎共享的 I/O 瓶颈。
+- 第 3 项能降低后续改动时的维护风险。
+
+这三项一起做，收益通常比零散微优化更明显。
+
+## 八、结论
+
+这个项目最值得肯定的地方，是已经天然分成了：
+
+- 参考实现
+- 高性能引擎
+- 统一可视化管线
+
+这意味着它非常适合做“分层优化”，不需要推倒重来。
+
+从当前代码看，后续最有效的路线不是继续在外围加流程判断，而是：
+
+- 把 Python 热点循环尽量向量化
+- 把文本轨迹格式逐步替换为二进制主格式
+- 把数据流和状态管理收紧
+
+如果按这个方向推进，这个项目会同时得到：
+
+- 更好的性能
+- 更低的内存占用
+- 更稳定的多引擎一致性
+- 更容易继续扩展新的物理项和可视化形式

optimization/codex_v1.md

@@ -0,0 +1,341 @@
+# Dynamics 优化方案综合评估
+
+> 分析日期：2026-06-12  
+> 评估对象：`optimization/claude.md`、`optimization/codex.md`、`optimization/workbuddy.md`  
+> 说明：用户问题中出现两次 `claude.md`，但目录内实际为三份不同文档，因此本文按 **Claude / Codex / Workbuddy** 三个工具的输出进行综合分析。
+
+## 一、结论先行
+
+这三份建议并不是互相冲突，而是分别代表了三种不同但互补的优化视角：
+
+- **Claude**：最像“代码审计员”，擅长找出会崩溃、会出错、会不一致的具体问题。
+- **Codex**：最像“性能工程师”，擅长定位热点路径、区分 CPU 与 I/O 瓶颈、给出高收益优化顺序。
+- **Workbuddy**：最像“架构负责人”，擅长从模块拆分、全局变量治理、配置统一、测试体系等角度规划长期演进。
+
+如果只选一个方案：
+
+- 短期救火，优先参考 **Claude**
+- 中期提速，优先参考 **Codex**
+- 长期重构，优先参考 **Workbuddy**
+
+如果要形成真正适合这个项目的落地方案，最佳做法不是三选一，而是：
+
+1. 先按 **Claude** 修 correctness bug
+2. 再按 **Codex** 做性能主线优化
+3. 最后按 **Workbuddy** 做结构化治理
+
+---
+
+## 二、三款工具的评价
+
+## 1. Claude
+
+### 优势
+
+- **最强的 bug 发现能力**  
+  `claude.md` 明确指出了几个高价值问题，例如：
+  - `compute.py` 里 `run_simulation()` 直接引用未传入的 `config`
+  - `dynamics.py` 绘图分支使用未定义变量
+  - `plot_wave.py` 仍按旧格式读取新 `display.txt`
+
+  这些问题都属于“不是慢，而是可能直接错或崩”的问题，优先级非常高。
+
+- **问题定位具体，能直接改代码**  
+  Claude 的输出不是抽象建议，而是带着明确文件位置、触发条件和修复方向，适合立即进入修复。
+
+- **对多引擎一致性比较敏感**  
+  它注意到了 Python/C/C++/Fortran 之间在 `save_trajectory`、势能分配、输出格式支持等方面的不一致，这对多语言框架很重要。
+
+### 劣势
+
+- **整体视角偏“局部修补”**  
+  它非常擅长发现局部 bug，但对“整体性能主瓶颈怎么排优先级”没有 Codex 那么系统。
+
+- **长期架构建议相对弱一些**  
+  虽然也提到一些结构问题，但没有 Workbuddy 那种模块拆分、配置治理、测试体系的完整路线。
+
+- **偏静态审计，偏保守**  
+  更适合“先修对”，不一定最适合“先做最值的性能投资”。
+
+### 适合的角色
+
+- 第一轮排雷
+- 回归前 bug 清单
+- 多引擎一致性核对
+
+---
+
+## 2. Codex
+
+### 优势
+
+- **性能视角最清晰**  
+  `codex.md` 非常明确地区分了三类瓶颈：
+  - Python 计算循环
+  - 力学计算中的逐键/双层循环
+  - 文本 I/O
+
+  这种分层对项目很有价值，因为它避免把“算法慢”和“文件慢”混成一个问题。
+
+- **优化顺序合理，投入产出比高**  
+  它提出的几个优先项都很务实：
+  - 弹簧力向量化
+  - `display.txt` 增加二进制格式
+  - `draw.py` Marker 批量更新
+
+  这些建议的共同特点是：
+  - 不需要推翻现有架构
+  - 可验证
+  - 回报高
+
+- **能识别共享瓶颈**  
+  它特别强调：即使换成 C 引擎，文本 I/O 仍然会拖后腿。这种判断比只盯 Python 引擎更深入。
+
+### 劣势
+
+- **对当前代码中的直接 bug 关注度不如 Claude**  
+  虽然也指出了 `dynamics.py` 数据流和 README 不一致等风险，但没有 Claude 那么系统地清点“哪些地方一开功能就报错”。
+
+- **工程治理不如 Workbuddy 完整**  
+  它提到全局变量和状态对象，但没有把模块拆分、测试体系、案例配置统一讲得那么全面。
+
+- **更偏策略，不总是落到立即可改的具体 patch**  
+  对项目 owner 来说这很好，但如果是要当天就修 bug，Claude 会更直接。
+
+### 适合的角色
+
+- 性能优化主方案
+- 技术债排序
+- 中期研发路线设计
+
+---
+
+## 3. Workbuddy
+
+### 优势
+
+- **最强的工程化和长期治理视角**  
+  `workbuddy.md` 关注的重点包括：
+  - `compute.py` 拆模块
+  - 消除全局变量
+  - 配置结构统一
+  - 测试体系建设
+  - Fortran/C/C++ 引擎功能对齐
+
+  这类建议对项目走向稳定维护很关键。
+
+- **适合把个人项目变成可持续项目**  
+  它不是只盯某一个性能点，而是在思考这个项目未来怎么更好地维护、测试、扩展。
+
+- **对“项目管理成本”比较敏感**  
+  比如配置格式不统一、废弃文件残留、测试缺失，这些短期不一定影响运行，但长期一定影响开发效率。
+
+### 劣势
+
+- **短期收益不如另外两者明显**  
+  比如模块拆分、类型标注、状态封装，这些很重要，但如果当前主要痛点是“功能坏了”或“运行太慢”，它们不会立刻让用户感受到变化。
+
+- **有些建议偏大改**  
+  像 `compute.py` 模块拆分、全局变量全面状态化，工作量不小，如果现在马上做，容易把项目带入较长重构周期。
+
+- **对现存 bug 的敏锐度不如 Claude，对热点性能的聚焦不如 Codex**  
+  它更像制定制度和框架的人，而不是先冲到一线止血的人。
+
+### 适合的角色
+
+- 中长期架构治理
+- 重构路线设计
+- 代码库标准化与测试建设
+
+---
+
+## 三、三者对比总结
+
+| 维度 | Claude | Codex | Workbuddy |
+|------|--------|-------|-----------|
+| 找 bug | 最强 | 中等 | 中等 |
+| 找性能瓶颈 | 强 | 最强 | 中等 |
+| 架构治理 | 中等 | 强 | 最强 |
+| 可直接落地修复 | 最强 | 强 | 中等 |
+| 中期路线规划 | 中等 | 最强 | 强 |
+| 长期工程化 | 中等 | 强 | 最强 |
+| 风险意识 | 最强 | 强 | 强 |
+| 适合当前项目阶段 | 很适合当前排雷 | 很适合当前提速 | 很适合下一阶段治理 |
+
+可以概括为：
+
+- **Claude 更像 QA/审计**
+- **Codex 更像性能工程师**
+- **Workbuddy 更像架构经理**
+
+---
+
+## 四、综合判断
+
+从这三份文档综合看，`Dynamics` 项目当前不是单一问题，而是三类问题叠加：
+
+1. **已有功能正确性问题**
+   - 某些路径会直接崩溃
+   - 某些脚本已经和新格式脱节
+   - 多引擎行为存在不一致
+
+2. **性能与 I/O 问题**
+   - Python 引擎热点循环还没做足够向量化
+   - `display.txt`/`trajectory.txt` 文本格式开始成为共享瓶颈
+   - 动画路径里还有逐原子更新
+
+3. **工程治理问题**
+   - `compute.py` 过大
+   - 全局变量过多
+   - 配置不统一
+   - 缺少测试
+
+因此，最合理的综合方案必须是 **分阶段方案**，而不是试图一次性全做。
+
+---
+
+## 五、推荐的综合方案
+
+## 阶段 A：先修正确性，避免“带病优化”
+
+这一阶段以 Claude 的结论为主，目标是先把明显会崩溃、会错的地方修掉。
+
+建议优先做：
+
+1. 修复 `compute.py` 中 `run_simulation()` 对 `config` 的非法引用
+2. 修复 `dynamics.py` 绘图分支中未定义变量的死代码问题
+3. 修复 `plot_wave.py` 对新 `display.txt` 格式的不兼容
+4. 统一外部引擎在 `save_trajectory`、输出格式上的行为
+5. 清理少量明显错误与危险写法
+   - case06 文案错误
+   - `draw.py` 裸 `except`
+
+这一阶段的目标不是提速，而是让系统进入“功能可信、路径可跑”的状态。
+
+---
+
+## 阶段 B：再做高收益性能优化
+
+这一阶段以 Codex 的结论为主，重点是先做收益大、改动相对可控的优化。
+
+建议优先做：
+
+1. 给 `compute.py` / `dynamics.py` 增加基础 profiling
+   - 总耗时
+   - 力计算耗时
+   - 抽帧耗时
+   - I/O 耗时
+
+2. 向量化 `compute_force()` 中的弹簧键计算
+
+3. 在 `draw.py` 中将 Marker 更新改成整列切片赋值
+
+4. 新增二进制显示格式
+   - 推荐 `display.npz`
+   - `draw.py` 优先读二进制，回退读文本
+
+5. 将完整轨迹保存改为分块/流式写出
+
+这一步完成后，项目会得到最直观的收益：
+
+- Python 引擎更快
+- C 引擎的后处理更快
+- 大规模 case 的内存压力更低
+
+---
+
+## 阶段 C：最后做工程治理与结构化重构
+
+这一阶段以 Workbuddy 的结论为主，目标是让项目从“能用”变成“好维护”。
+
+建议优先做：
+
+1. 拆分 `compute.py`
+   - `core.py`
+   - `io.py`
+   - `params.py`
+   - `runner.py`
+   - `engine_helper.py`
+
+2. 引入 `SimulationState` / `SimulationConfig` / `SystemData`
+   - 逐步消灭模块级全局变量
+
+3. 统一案例配置格式
+   - 六个案例字段保持一致
+   - 将颜色、相机等零散字段收束成结构化配置
+
+4. 建立最小测试体系
+   - 数值方法单元测试
+   - I/O 一致性测试
+   - Python/C 输出一致性测试
+
+5. 清理废弃脚本和旧逻辑
+   - 尤其是已经不再处于主流程中的脚本
+
+这一步不会像阶段 B 那样马上体现性能收益，但它会显著降低后续维护成本。
+
+---
+
+## 六、最终推荐方案
+
+如果只给一个“综合最优”的执行方案，我建议采用下面这个顺序：
+
+1. **先按 Claude 修 bug**
+   - 因为错误路径不先修，后面性能优化和重构都会建立在不稳定基础上
+
+2. **再按 Codex 做性能主线**
+   - 因为当前性能瓶颈已经很清楚，尤其是 Python 热点和文本 I/O
+
+3. **最后按 Workbuddy 做结构治理**
+   - 因为工程化重构最适合在功能稳定、热点清楚之后进行
+
+换句话说：
+
+- **Claude 决定“先修什么”**
+- **Codex 决定“先优化什么”**
+- **Workbuddy 决定“最后把系统整理成什么样”**
+
+这是三份建议里最不冲突、也最符合实际开发节奏的组合方式。
+
+---
+
+## 七、简版执行清单
+
+如果要把这份综合方案压缩成一份实际待办，我建议如下：
+
+### P0：立即处理
+
+- 修 `run_simulation()` 的 `config` 作用域错误
+- 修 `dynamics.py` 绘图分支未定义变量
+- 修 `plot_wave.py` 与新 `display.txt` 格式不兼容
+- 校正多引擎 `save_trajectory` / 输出行为不一致
+
+### P1：本轮优化主线
+
+- 增加 profiling
+- 向量化弹簧力计算
+- `draw.py` Marker 批量更新
+- 新增 `display.npz`
+- 完整轨迹改为分块写出
+
+### P2：下一轮治理
+
+- 拆分 `compute.py`
+- 全局变量状态化
+- 统一 `input.txt` 结构
+- 建立 pytest 测试
+- 清理废弃脚本与旧文档
+
+---
+
+## 八、最终评价
+
+三款工具的输出质量都不错，但擅长点不同。
+
+- **Claude** 最适合当前这个项目的“第一步”，因为它最能发现会直接影响运行正确性的缺陷。
+- **Codex** 最适合当前这个项目的“第二步”，因为它对性能瓶颈的层次判断最清楚，最容易产出高收益优化。
+- **Workbuddy** 最适合当前这个项目的“第三步”，因为它能把一套已经能跑的系统整理成一套更稳定、更易维护的工程。
+
+因此，综合最优解不是选其中一个，而是：
+
+**用 Claude 排雷，用 Codex 提速，用 Workbuddy 收尾治理。**

optimization/workbuddy.md

@@ -0,0 +1,303 @@
+# Dynamics 项目优化建议
+
+> 分析日期：2026-06-12
+> 版本：main (43 commits)
+
+---
+
+## 1. 架构优化
+
+### 1.1 compute.py 拆分（高优先级）
+
+**现状**：`compute.py` 1618 行，混合了以下职责：
+- 物理引擎（数值积分方法 ~10 个函数）
+- 文件 I/O（JSON/text 读写）
+- 外部引擎管理（subprocess 调用）
+- 参数加载与解析
+- 主入口 main()
+
+**建议**：按职责拆分为独立模块
+
+```
+compute/
+├── __init__.py          # 导出公共 API
+├── core.py              # 物理引擎: 数值积分方法, 力的计算 (~200 行)
+├── io.py                # 文件 I/O: save/load display/trajectory (~300 行)
+├── params.py            # 参数加载: load_coord_file, load_bond 等 (~150 行)
+├── runner.py            # 运行管理: run_simulation, run_engine, run_from_config (~200 行)
+├── engine_helper.py     # 外部引擎: 校准, subprocess, 进度监控 (~150 行)
+└── main.py              # 主入口: main(), if __name__ == "__main__" (~30 行)
+```
+
+**收益**：可单独测试物理算法、文件读写、引擎调用
+
+### 1.2 外部引擎代码重构（中优先级）
+
+**现状**：C/C++/Fortran 引擎各自实现了一套：
+- JSON 读取/写入（3 套独立实现）
+- 物理算法（复制粘贴）
+- display.txt / trajectory.txt 输出
+
+**建议**：
+- C 和 C++ 引擎共用头文件/函数库
+- 物理算法提取为公共库（`engines/common/`）
+- Fortran 更新支持 `save_trajectory` 兼容
+
+---
+
+## 2. 性能优化
+
+### 2.1 外部引擎校准耗时（高）
+
+**现状**：`run_engine()` 每次执行前先跑 1000 步校准测速（`_calib_nt = min(1000, max(100, total_steps // 10))`），对短时运行（~1000 步）校准占比高达 50%。
+
+**建议**：
+- 缓存校准结果：按 `(engine, NT, DT, NSTEP, n_atoms)` 哈希缓存
+- 或对短运行 (<5000 步) 跳过校准
+- 或只在校准目录为空时执行
+
+```python
+# 建议改动
+_calib_cache = os.path.join(script_dir, "engines", engine, "_calib_cache.json")
+if os.path.exists(_calib_cache):
+    with open(_calib_cache) as f:
+        _calib_data = json.load(f)
+    if _calib_data.get("key") == _calib_key:
+        _step_time = _calib_data["step_time"]
+        _overhead = _calib_data["overhead"]
+        skip_calibration()
+```
+
+### 2.2 display.txt 读取加速（已完成）
+
+**现状**：已从逐行 `split()+float()` 改为 `np.genfromtxt()` 批量解析，200帧×120原子 加载仅 0.087s（测试已通过）。
+
+### 2.3 Python 引擎性能（低）
+
+**现状**：Python 引擎比 C 引擎慢约 6-8 倍（43s vs 7s for 10 万步）。对 50 万步以上的长时间模拟不适用。
+
+**建议**：
+- 使用 Numba JIT 加速力的计算
+- 或使用 `numpy.vectorize` 替代 Python 循环（力的计算中已部分使用 numpy 数组操作）
+- 文档中注明 Python 引擎仅适用于小规模测试
+
+### 2.4 校准数据的进度条（已修复）
+
+校准后进度条使用实时帧号更新，不再使用时间估算，避免了 `0% → 100%` 跳变。但注意校准本身不显示进度。
+
+---
+
+## 3. 代码质量
+
+### 3.1 消除全局变量（高优先级）
+
+**现状**：`compute.py` 定义了 40+ 个模块级全局变量：
+
+```python
+box_a = None
+alpha = None
+ATOM_IDS = None
+...
+```
+
+所有算法函数直接依赖这些全局变量，导致：
+- 无法并发运行两个模拟
+- 无法单元测试单个函数
+- 函数签名无法自文档化
+
+**建议**：将全局状态封装为 `SimulationState` 类：
+
+```python
+@dataclass
+class SimulationState:
+    box_a: float = 10.0
+    atom_ids: np.ndarray = None
+    atom_masses: np.ndarray = None
+    ...
+
+def leapfrog(state: SimulationState, ...):
+    ...
+```
+
+### 3.2 draw.py 全局变量（中优先级）
+
+**现状**：类似的问题，`draw.py` 用了 ~30 个模块级变量：
+
+```python
+DISPLAY_X = ...
+DISP_ALL_X = ...
+N_FRAMES = ...
+frame_idx = 0
+```
+
+**建议**：封装为 `AnimationData` 和 `CameraState` 类
+
+### 3.3 废弃文件清理（低优先级）
+
+以下文件可能已废弃或重复：
+
+| 文件 | 行数 | 判断 |
+|------|------|------|
+| `sample.py` | 184 | 已不再被调用（抽帧已移至引擎内部） |
+| `build_release_zip.py` | - | 打包脚本，非核心功能 |
+| `export_web_data.py` | - | 导出脚本，非核心功能 |
+| `migrate_npz_outputs.py` | - | 迁移脚本，一次性使用 |
+| `tools/NumericalMethods.py` | - | 可能已过时 |
+| `build/` (CMake 相关) | - | 与独立引擎编译无关 |
+
+---
+
+## 4. 配置管理
+
+### 4.1 input.txt 格式统一（高优先级）
+
+**现状**：6 个案例的 `input.txt` 格式不统一：
+
+| 字段 | case01 | case06 |
+|------|--------|--------|
+| `save_trajectory` | ❌ 缺失 | ✅ 有 |
+| `camera_distance/elevation/azimuth` | ❌ 缺失 | ✅ 有 |
+| `move_camera` | ❌ 缺失 | ✅ 有 |
+| `step_sample` | 仍然写 1 | 已改为 0 |
+| `step_plot`配置 | 旧注释格式 | 新注释格式 |
+
+**建议**：
+
+```yaml
+# 统一模板
+flow:
+  step_simulate: 1
+  step_sample: 0           # 引擎已内置抽帧
+  step_animation: 1
+  save_trajectory: 0       # 默认不保留
+
+physics:
+  box_a: 10.0
+  G: [0, 0, -9.8]
+  method: leapfrog
+  T_total: 10.0
+  NSTEP: 100
+  DT: 0.001
+
+render:
+  use_marker: 0
+  alpha: [0.0]*6
+  ball_color: [0.9, 0.2, 0.2]
+  box_color: [0.8, 0.8, 0.85]
+  camera:
+    distance: 40.0
+    elevation: 0
+    azimuth: 0
+    move: false
+```
+
+### 4.2 YAML 结构优化（中优先级）
+
+**现状**：`ball_color_r/g/b` 三个独立的键，YAML 列表更自然：
+
+```yaml
+# 当前
+ball_color_r: 0.90
+ball_color_g: 0.20
+ball_color_b: 0.20
+
+# 建议
+ball_color: [0.90, 0.20, 0.20]
+```
+
+---
+
+## 5. 测试与质量
+
+### 5.1 缺少测试（高优先级）
+
+**现状**：无单元测试、无集成测试。每次修改只能手动跑案例。
+
+**建议**：
+- 为 `compute.py` 的物理算法（Leapfrog、Euler 等）添加 pytest 测试
+- 为文件 I/O（`save_display_txt` / `load_display_txt` 读写一致性）添加测试
+- 添加一个 `examples/case00` 最小验证案例（2 个原子，10 步）
+- CI：`pytest` + 简单的集成测试脚本
+
+### 5.2 缺少类型标注（中优先级）
+
+```python
+# 现状
+def leapfrog(x, y, z, vx, vy, vz, dt, m, g, b):
+
+# 建议
+def leapfrog(
+    x: np.ndarray, y: np.ndarray, z: np.ndarray,
+    vx: np.ndarray, vy: np.ndarray, vz: np.ndarray,
+    dt: float, m: np.ndarray, g: np.ndarray, b: np.ndarray
+) -> tuple[np.ndarray, np.ndarray, np.ndarray,
+           np.ndarray, np.ndarray, np.ndarray]:
+```
+
+---
+
+## 6. 引擎一致性
+
+### 6.1 Fortran 引擎更新（高）
+
+**现状**：Fortran 引擎**不支持** `save_trajectory=0` 模式，总是写 JSON 格式的 `trajectory.txt`。需要：
+
+1. 读取 `save_trajectory` 参数
+2. 计算时按 NSTEP 采样
+3. 写 display.txt 新格式
+4. 按条件跳过 trajectory.txt 输出
+
+参考 C 引擎的 `write_display_txt` 实现（约 40 行代码）。
+
+### 6.2 参数默认值一致性（中）
+
+各引擎的 `SimParams` 默认值不一致：
+
+| 参数 | C | C++ | Fortran |
+|------|---|---|---------|
+| `box_a` | 无默认 | 10.0 | 10.0 (代码中) |
+| `method` | leapfrog | leapfrog | leapfrog |
+| `save_trajectory` | 0 | 1 | 不支持 |
+
+---
+
+## 7. 用户体验
+
+### 7.1 进度条改进（低）
+
+- 外部引擎校准阶段不显示进度（建议加 "校准中…" 文本）
+- Fortran 引擎支持时，增加进度消息
+
+### 7.2 缓存一致性（已修复）
+
+已修复：`display.txt` 在 `save_trajectory=0` 时不再保留 `trajectory.txt`，外部引擎路径总是从 `trajectory.txt` 重新抽帧。
+
+### 7.3 move_camera.txt 实时生效（已修复）
+
+`draw.py` 直接读取 `move_camera.txt`，修改后重启动画即可生效，无需重新跑模拟。
+
+---
+
+## 8. 实施优先级
+
+| 优先级 | 内容 | 预估工作量 |
+|--------|------|-----------|
+| P0 | Fortran 引擎支持 save_trajectory=0 | 2-3h |
+| P0 | 全局变量封装为 SimulationState | 3-4h |
+| P1 | compute.py 模块拆分 | 4-6h |
+| P1 | 添加单元测试 (pytest) | 3-5h |
+| P1 | input.txt 格式统一 (6个案例) | 1-2h |
+| P2 | 校准缓存加速 | 1-2h |
+| P2 | 废弃文件清理 | 0.5h |
+| P3 | draw.py 封装 AnimationData | 2-3h |
+| P3 | 类型标注 | 2-3h |
+
+---
+
+## 9. 快速 wins（半小时内可完成）
+
+1. ✅ 更新 Fortran 引擎（参考 C 的 `write_display_txt` 迁移）
+2. ✅ 删除 `sample.py`（已废弃）
+3. ✅ 统一 6 个案例的 `input.txt` 格式
+4. ✅ 添加 `save_trajectory` 到 case01-05 的 default 配置
+5. ✅ 统一 `ball_color_r/g/b` 为 `ball_color` 列表

optimization/workbuddy_v1.md

@@ -0,0 +1,216 @@
+# Dynamics 项目优化方案 — 三方工具综合评价与整合
+
+> **分析日期**：2026-06-12  
+> **参与工具**：WorkBuddy（本助手）、Claude（Sonnet 4.6）、Codex（GPT-4o 代码分析）  
+> **分析范围**：`D:\Share\Data\aliyun-gitea\dynamics` 完整代码库  
+> **参考文档**：`optimization/workbuddy.md`、`optimization/claude.md`、`optimization/codex.md`
+
+---
+
+## 一、三方工具评价
+
+### 1.1 WorkBuddy（AI Agent，Senior Developer 角色）
+
+| 维度 | 评价 |
+|------|------|
+| **优势** | ✅ 实际运行和修改过整个代码库，已验证 display.txt 格式/渲染参数/运动相机等功能的正确性 — 不是静态分析，是实机验证 |
+| | ✅ 覆盖最广：架构、性能、代码质量、配置、测试、引擎一致性、UX 等 9 个维度 |
+| | ✅ 每个建议标注优先级(P0-P3)和预估工作量(小时级) |
+| | ✅ 发现了一些真正的坑：C++ 引擎 `save_trajectory` 默认值 1（导致行为不一致）、`alpha` 字段未传递到 display.txt 等 |
+| | ✅ 提供 input.txt 格式统一模板和 YAML 结构优化建议 |
+| **劣势** | ❌ 建议偏"宏观架构"（模块拆分、全局变量封装），缺乏具体的向量化代码 |
+| | ❌ 没有发现 B1（config 不存在作用域）、B2（绘图代码完全失效）等具体 Bug |
+| | ❌ 对 Python 引擎性能优化仅提到 Numba，没有给出弹簧力向量化等具体代码 |
+
+### 1.2 Claude（Sonnet 4.6）
+
+| 维度 | 评价 |
+|------|------|
+| **优势** | ✅ **Bug 挖掘能力极强**：发现了 `run_simulation` 内 `config` 未定义、`dynamics.py` 绘图块完全失效、`plot_wave.py` 格式不匹配等 5 个确切 Bug |
+| | ✅ **Python 引擎向量化代码极其具体**：弹簧力、引力 O(N²)、固定约束、驱动力的优化都给出了可直接替换的代码（含 `np.add.at`） |
+| | ✅ 引擎一致性分析全面：C/C++/Fortran 默认值差异、势能归属约定差异 |
+| | ✅ 标注了 `load_parameters` 已废弃、相机解析函数重复等代码质量问题 |
+| | ✅ 按投入产出比排序最高回报的 3 件事 |
+| **劣势** | ❌ 纯静态分析，**部分结论可能是误报**：B1 中 `config` 在 `run_simulation` 中不可用 → 但 `run_simulation` 是从 `run_from_config` 内部调用的，实际上 Python 引擎路径已验证可正常工作（后续对话中用户成功运行） |
+| | ❌ 没有分析 Fortran 引擎（可能只读了 C 引擎代码） |
+| | ❌ 没有发现 display.txt 格式最初加载过慢的问题（已修复） |
+| | ❌ 没有发现 `use_marker` 字段丢失导致 VisPy 卡顿的问题（已修复） |
+| | ❌ 没有分析 input.txt 格式统一性问题 |
+
+### 1.3 Codex（GPT-4o）
+
+| 维度 | 评价 |
+|------|------|
+| **优势** | ✅ **战略思维最好**：建议先加 profiling 再做优化（"没度量就没有优化"），反对盲目优化 |
+| | ✅ 提出 `display.npz` 二进制格式替代方案 — 所有工具中唯一想到这个的 |
+| | ✅ "分两层优化"理念清晰：计算性能(CPU) vs I/O 性能(磁盘)分开处理 |
+| | ✅ 强调 "最小一致性测试集" — 多引擎正确性验证的实用方案 |
+| | ✅ 推荐实施顺序最合理：1)计时 2)向量化 3)Marker 切片 4)二进制输出 5)修复绘图 6)架构重构 |
+| **劣势** | ❌ 代码细节最少，没有给出具体的向量化实现 |
+| | ❌ 没有发现任何确切的 Bug——描述的都是"可能的风险"而非"可复现的崩溃" |
+| | ❌ 架构建议偏泛（"减少全局变量""收敛到状态对象"），没有 WorkBuddy 的 SimulationState 代码示例 |
+| | ❌ 没有分析相机运动、渲染参数传递等最近新增功能 |
+
+---
+
+## 二、建议对比汇总
+
+### 2.1 Bug 发现对比
+
+| Bug 编号 | 描述 | WorkBuddy | Claude | Codex | 验证状态 |
+|----------|------|-----------|--------|-------|---------|
+| B1 | `run_simulation` 内 `config` 变量作用域 | ❌ 未发现 | ✅ 发现 | ❌ 未发现 | ⚠️ 疑似误报（Python 引擎实际可用） |
+| B2 | `dynamics.py` 绘图块死代码/变量未定义 | ❌ 未发现 | ✅ 发现 | ✅ 提及风险 | ✅ 确认（step_plot=1 时崩溃） |
+| B3 | `plot_wave.py` 旧格式加载新 display.txt | ❌ 未发现 | ✅ 发现 | ❌ 未发现 | ✅ 确认（step_plot_wave=1 时崩溃） |
+| B4 | case06 描述写 case01 | ❌ 未发现 | ✅ 发现 | ❌ 未发现 | ✅ 确认 |
+| B5 | `draw.py` 裸 `except` | ❌ 未发现 | ✅ 发现 | ❌ 未发现 | ✅ 确认 |
+| B6 | display.txt 旧格式加载过慢（已修复） | ✅ 发现修复 | ❌ 未发现 | ❌ 未发现 | ✅ 已修复 |
+| B7 | `use_marker` 丢失 → VisPy 卡顿（已修复） | ✅ 发现修复 | ❌ 未发现 | ❌ 未发现 | ✅ 已修复 |
+| B8 | `alpha` 未写入 display.txt header（已修复） | ✅ 发现修复 | ❌ 未发现 | ❌ 未发现 | ✅ 已修复 |
+| B9 | 运动相机数据缓存不刷新（已修复） | ✅ 发现修复 | ❌ 未发现 | ❌ 未发现 | ✅ 已修复 |
+
+### 2.2 优化建议对比
+
+| 优化项 | WorkBuddy | Claude | Codex | 综合评价 |
+|--------|-----------|--------|-------|---------|
+| 弹簧力向量化 | 仅提到 Numba | ✅ 完整代码含 `np.add.at` | ✅ 有思路但无代码 | **Claude 最佳** |
+| 固定约束优化 | 提到但无代码 | ✅ 有代码 | ✅ 有思路 | **Claude 最佳** |
+| display.npz 二进制格式 | ❌ 未想到 | ❌ 未想到 | ✅ 唯一想到 | **Codex 最佳** |
+| 全局变量封装 SimulationState | ✅ 有代码示例 | ✅ 提及 | ✅ 提及 | **WorkBuddy 最佳** |
+| compute.py 模块拆分 | ✅ 完整目录结构 | ❌ 未提及 | ❌ 未提及 | **WorkBuddy 最佳** |
+| input.txt 格式统一 | ✅ 完整模板 | ❌ 未提及 | ❌ 未提及 | **WorkBuddy 最佳** |
+| 先加 profiling 再优化 | ❌ 未提及 | ❌ 未提及 | ✅ 强烈建议 | **Codex 最佳** |
+| Fortran 引擎更新 | ✅ 提及 P0 | ✅ 提及 | ❌ 未提及 | **WorkBuddy+Claude** |
+| 外部引擎校准缓存 | ✅ 有代码示例 | ✅ 有方案 | ❌ 未提及 | **WorkBuddy+Claude** |
+| 多引擎一致性测试集 | ❌ 未提及 | ❌ 未提及 | ✅ 建议 | **Codex 最佳** |
+| Marker 更新切片化 | ❌ 未提及 | ❌ 未提及 | ✅ 建议 | **Codex 最佳** |
+| 驱动去 t_vec 分配 | ❌ 未提及 | ✅ 有代码 | ✅ 提及 | **Claude 最佳** |
+
+### 2.3 风格对比
+
+| 维度 | WorkBuddy | Claude | Codex |
+|------|-----------|--------|-------|
+| 粒度 | 宏观+微观 | 微观为主 | 宏观为主 |
+| 代码示例 | 中等（配置模板/架构） | 丰富（向量化/性能优化） | 最少 |
+| Bug 发现 | 运行时验证的 Bug | 静态分析发现的 Bug | 潜在风险 |
+| 战略层 | 中度 | 低 | 高 |
+| 实施顺序 | P0-P3 优先级 | 投入产出比排序 | 推荐实施顺序 |
+| 可信度 | 高（实际运行过） | 中（静态分析） | 中（静态分析） |
+
+---
+
+## 三、最终综合方案
+
+综合三方分析，推荐按以下 **6 个阶段** 实施，每阶段都有明确的可验证交付物。
+
+### 阶段一：立即修复已确认的 Bug（1 天）
+
+| 编号 | 内容 | 参考工具 | 工作量 |
+|------|------|---------|--------|
+| F1 | `plot_wave.py` 适配 `load_display_txt` 新格式 (B3) | Claude | 1-2h |
+| F2 | `dynamics.py` 绘图块添加保护，避免 `step_plot=1` 崩溃 (B2) | Claude | 30min |
+| F3 | 修复 `draw.py` 裸 `except` (B5) | Claude | 5min |
+| F4 | 修复 case06 描述文字 (B4) | Claude | 1min |
+| F5 | C++ 引擎 `save_trajectory` 默认值改为 0 (E2) | Claude | 5min |
+| F6 | **Fortran 引擎支持 `save_trajectory=0`**（参照 C 引擎实现） | WorkBuddy/Claude | 2-3h |
+| F7 | 相机解析函数去重：`dynamics.py` 导入 `compute._load_camera_motion` (Q1) | Claude | 15min |
+| F8 | 废弃 `load_parameters` + `main()` 删除或移入 tools/ (Q2) | Claude | 15min |
+
+**验证**：Python/C/C++/Fortran 全部 4 种引擎跑 case06，`step_plot: 1` + `step_plot_wave: 1` 不崩溃
+
+### 阶段二：Python 引擎性能优化（2 天）
+
+| 编号 | 内容 | 参考工具 | 工作量 |
+|------|------|---------|--------|
+| P1 | **弹簧力向量化**：`compute_force()` 中键循环改为 `np.add.at` 批量计算 | Claude（完整代码） | 1h |
+| P2 | 固定约束原地掩码写回，消除 `column_stack` | Claude | 30min |
+| P3 | `frame_indices` 列表改计数器 | Claude | 5min |
+| P4 | 驱动力去掉 `t_vec` 临时数组 | Claude/Codex | 5min |
+| P5 | `GRAVITY_INTERACTION` O(N²) 双重循环向量化（可选） | Claude | 1h |
+
+**验证**：`engine: python` 跑 case06 耗时缩短 5-15 倍（目标：从 43s → <5s）
+
+### 阶段三：I/O 与可视化优化（2 天）
+
+| 编号 | 内容 | 参考工具 | 工作量 |
+|------|------|---------|--------|
+| I1 | **新增 `display.npz` 二进制格式**：`save_display_npz` / `load_display_npz` | Codex（方案） | 1-2h |
+| I2 | `draw.py` 优先读 `.npz`，不存在时回退 `display.txt` | Codex | 1h |
+| I3 | `draw.py` Marker 更新改为切片赋值（消除 `for i in range(N_ATOMS)`） | Codex | 15min |
+| I4 | `save_trajectory=1` 时轨迹数据改为 `memmap` 或分块写入 | Codex | 2h |
+
+**验证**：读取 200 帧×120 原子数据 < 0.01s（现 0.087s），动画帧率 60fps
+
+### 阶段四：架构重构（3 天）
+
+| 编号 | 内容 | 参考工具 | 工作量 |
+|------|------|---------|--------|
+| A1 | **`compute.py` 模块拆分**：`core/io/params/runner/main` | WorkBuddy（目录结构） | 4-6h |
+| A2 | **全局变量封装为 `SimulationState` 类**，从 `compute_force()` 开始 | WorkBuddy | 3-4h |
+| A3 | `draw.py` 全局变量封装为 `AnimationData` + `CameraState` | WorkBuddy | 2-3h |
+| A4 | C 和 C++ 引擎共用公共头文件（`engines/common/`） | Workbuddy | 3h |
+
+**验证**：所有 6 个 case + 4 种引擎，结果与重构前一致
+
+### 阶段五：配置与文档统一（1 天）
+
+| 编号 | 内容 | 参考工具 | 工作量 |
+|------|------|---------|--------|
+| C1 | **统一 6 个案例的 input.txt 格式**（增加 save_trajectory/camera 等缺失字段） | WorkBuddy（模板） | 1-2h |
+| C2 | `ball_color_r/g/b` 改为 YAML 列表 `ball_color: [R,G,B]` | WorkBuddy | 30min |
+| C3 | 更新 README.md 匹配当前架构 | Codex | 1h |
+| C4 | 注释清理：边界条件分工、弹性势能归属约定 | Claude | 20min |
+
+**验证**：6 个 `run_dynamics.py` 全部可运行，README 描述与代码一致
+
+### 阶段六：测试与CI（2 天）
+
+| 编号 | 内容 | 参考工具 | 工作量 |
+|------|------|---------|--------|
+| T1 | **添加 pytest 单元测试**：物理算法（Leapfrog/Euler）、文件 I/O 读写一致性 | WorkBuddy | 3-5h |
+| T2 | **多引擎一致性测试**：4 种引擎跑 case07（2 原子 10 步最小案例），输出容差内一致 | Codex | 2h |
+| T3 | 添加基础 profiling 计时（总时间、力计算、I/O 分段） | Codex | 1h |
+| T4 | 添加外部引擎校准缓存 | WorkBuddy/Claude | 1-2h |
+
+**验证**：`pytest` 绿色通过，`python run_dynamics.py --engine python` 与 `--engine c` 结果一致
+
+---
+
+## 四、总体工作量估算
+
+| 阶段 | 内容 | 预估人天 | 依赖 |
+|------|------|---------|------|
+| 一 | Bug 修复 | 1 天 | — |
+| 二 | Python 引擎性能优化 | 2 天 | 阶段一 |
+| 三 | I/O 与可视化优化 | 2 天 | 阶段一 |
+| 四 | 架构重构 | 3 天 | 阶段一、二 |
+| 五 | 配置与文档统一 | 1 天 | 阶段四 |
+| 六 | 测试与 CI | 2 天 | 阶段二、四 |
+| **合计** | | **~11 人天** | |
+
+---
+
+## 五、如果只做 3 件事
+
+基于三方分析共识 + 投入产出比，最值得做的三件事：
+
+1. **弹簧力向量化**（Claude 给出完整代码，Python 引擎性能提升 5-15 倍）
+2. **Fortran 引擎支持 `save_trajectory=0`**（WorkBuddy + Claude 共识 P0，所有引擎行为一致的必要条件）
+3. **新增 `display.npz` 二进制格式**（Codex 独有见解，统一切换 I/O 性能瓶颈）
+
+---
+
+## 六、工具选择建议
+
+| 场景 | 推荐工具 | 原因 |
+|------|---------|------|
+| 找 Bug | **Claude** | Bug 发现能力最强（5 个确切 Bug） |
+| 性能优化（向量化） | **Claude** | 给出可直接替换的 Numpy 向量化代码 |
+| 架构重构 | **WorkBuddy** | 熟悉整体代码结构，有模块拆分/状态封装的具体方案 |
+| I/O 策略设计 | **Codex** | 战略思维好，提出二进制格式等创新方案 |
+| 实施顺序 | **Codex + WorkBuddy** | Codex 的 profiling-first 理念 + WorkBuddy 的 P0-P3 优先级 |
+| 确认 Bug 是否真实 | **实际运行测试** | 三方都是静态分析，最终需要运行确认 |
+
+---
+
+*本文档综合了 WorkBuddy（实机验证）、Claude Sonnet 4.6（静态 Bug 挖掘）、Codex/GPT-4o（战略思维）三方的分析。建议在实施每个阶段前，先用相关案例运行确认 Bug 现象和优化收益。*

optimization/workbuddy_v2.md

@@ -0,0 +1,241 @@
+# Dynamics 项目优化方案 — 三工具六版本综合分析 v2
+
+> **分析日期**：2026-06-12  
+> **参与工具**：WorkBuddy（Senior Developer 角色）、Claude（Sonnet 4.6）、Codex（GPT-4o）  
+> **分析范围**：`D:\Share\Data\aliyun-gitea\dynamics` 完整代码库  
+> **参考文档**：
+>   - 原始分析：`workbuddy.md`、`claude.md`、`codex.md`
+>   - 综合版本：`claude_v1.md`、`codex_v1.md`、`workbuddy_v1.md`
+
+---
+
+## 一、三工具角色定位
+
+三份原始分析报告和它们的 v1 版本揭示了三个工具最适合的角色，可以用一句话概括：
+
+> **Claude 是 QA / 审计员，Codex 是性能工程师，WorkBuddy 是架构负责人**
+
+这是一个高度互补的关系，而非互相竞争。
+
+### 1.1 角色矩阵
+
+| 维度 | Claude | Codex | WorkBuddy |
+|------|--------|-------|-----------|
+| **角色类比** | 代码审计员 | 性能工程师 | 架构负责人 |
+| **最擅长** | 找会崩溃的 Bug & 具体修复代码 | 定位性能瓶颈 & 战略优先级 | 模块拆分 & 工程治理 |
+| **最弱项** | 架构视野窄 | 代码细节少 | 微观 Bug 不敏感 |
+| **输出风格** | 逐条 Bug + 修复代码 | 分层分析 + 实施路线 | 全景目录 + 优先级表格 |
+| **适合阶段** | **第一步：排雷** | **第二步：提速** | **第三步：治理** |
+| **可执行性** | ⭐⭐⭐⭐⭐ 直接可改 | ⭐⭐⭐ 方案级 | ⭐⭐⭐ 框架级 |
+
+### 1.2 自我认知修正（v1 版本揭示的关键洞察）
+
+claude_v1.md 和 codex_v1.md 在分析 WorkBuddy 的原始报告时，指出了几项重要的自我修正：
+
+| 原始 WorkBuddy 结论 | v1 版本揭示的问题 | 修正 |
+|---------------------|-------------------|------|
+| 全局变量封装列为 P0（与 Fortran 同级） | Fortran 崩溃是功能阻断，重构是长期工作 | P0 → P2 ⬇️ |
+| 建议在 B1 修复中继续加全局变量 | 这与"封装全局变量"方向矛盾 | 短期补丁 + 长期封装需注明衔接 |
+| `compute.py` 拆模块建议偏理想化 | 未说明拆分顺序和风险控制 | 补充 `io.py` 优先，`runner.py` 最后 |
+
+---
+
+## 二、各工具 v1 版本评价
+
+### 2.1 claude_v1.md — Claude 的第二轮分析
+
+**性质**：Claude 对 workbuddy.md 和 claude.md 的综合
+
+**新增价值**：
+- 对 WorkBuddy 的优势/劣势做了 5 条以上具体分析（非泛泛批评）
+- 将 B1 修复方案做了更完善的设计（全局变量 + `global` 声明 + 作用域验证）
+- 给出了 B3 的完整分步修复代码（含 header 补充 `atom_masses` 字段）
+- **引擎评价表格**：对 C/C++/Fortran 三个引擎逐一评估功能/性能/一致性/代码质量
+- **"今天可以完成"清单**：5 个 1-5 分钟的小修复，非常实用
+
+**不足之处**：
+- 仍然基于静态分析，未发现 B1 实际上可能不是 Bug（已在实机验证中确认 Python 引擎可用）
+- 补丁式修复（继续加全局变量）与长期封装方向存在内在矛盾，虽已指出但未给出桥接方案
+
+### 2.2 codex_v1.md — Codex 的第二轮分析
+
+**性质**：Codex 对 workbuddy.md + claude.md + codex.md 三份的综合
+
+**新增价值**：
+- **"结论先行"式结构**：开篇即给出 "Claude=QA, Codex=Perf, WorkBuddy=Architect" 框架
+- **角色定位最清晰**：三段式建议（先 Claude 修正确性 → 再 Codex 做性能 → 最后 WorkBuddy 做治理）被 claude_v1.md 采纳，形成共识
+- 指出即使切到 C 引擎，文本 I/O 仍然是共享瓶颈 — 这是比仅关注 Python 引擎更深层的洞察
+- **性能分层最系统**：CPU 瓶颈（弹簧力）vs I/O 瓶颈（文本格式）vs 可视化瓶颈（逐原子更新）
+
+**不足之处**：
+- 对架构建议仍然偏泛，没有 WorkBuddy 或 claude_v1 的详细代码
+- 没有引擎一致性分析
+
+### 2.3 workbuddy_v1.md — 本工具的第一次综合
+
+**性质**：对三份原始报告的综合对比
+
+**自我评价**：
+| 维度 | 表现 | 与 v2 对比 |
+|------|------|-----------|
+| Bug 对比表格 | ✅ 全面（B1-B9，含验证状态） | 保持 |
+| 优化建议对比 | ✅ 12 项工具间对比 | 补充 v1 版本对比 |
+| 阶段划分 | 🟡 6 阶段偏细，v2 合并为 3 阶段 | 精简为 3 阶段 |
+| 工具选择建议 | ✅ 场景化推荐 | 保持并丰富 |
+| v1 版本覆盖 | ❌ 未覆盖 claude_v1/codex_v1 | 新增 ✅ |
+| 自我修正 | ❌ 未对自己做评价 | 新增 §1.2 ✅ |
+
+---
+
+## 三、Bug 核实与状态（B1-B9）
+
+### 3.1 已验证的真 Bug（可直接复现）
+
+| 编号 | 描述 | 触发条件 | 发现者 | 验证方式 | 优先级 |
+|------|------|---------|--------|---------|--------|
+| **B2** | `dynamics.py` 绘图块变量未定义 → NameError | `step_plot: 1` | Claude | 代码分析，当前案例设为 0 才未爆发 | 🔴 P0 |
+| **B3** | `plot_wave.py` 使用 `load_text_data` 读新格式 display.txt | `step_plot_wave: 1` | Claude | 代码分析，函数调用链不兼容 | 🔴 P0 |
+| **B4** | case06 描述文字写 case01 | 运行显示 | Claude | 肉眼可见 | 🟡 P1 |
+| **B5** | `draw.py` 裸 `except` 吞异常 | 任何运行 | Claude | 代码分析 | 🟡 P1 |
+| **E3** | Fortran 引擎不支持 `display.txt` | `engine: fortran` | 双方共识 | 代码分析，引擎无此功能 | 🔴 P0 |
+
+### 3.2 已修复的 Bug（当前版本已解决）
+
+| 编号 | 描述 | 发现者 | 修复阶段 |
+|------|------|--------|---------|
+| B6 | display.txt 加载过慢（逐行解析 60s+） | WorkBuddy | 已修复：改用 `np.genfromtxt` → 0.087s |
+| B7 | `use_marker` 未传递 → VisPy 卡顿 | WorkBuddy | 已修复：加入 header 传递 |
+| B8 | `alpha` 未写入 display.txt header | WorkBuddy | 已修复：C/C++/Python 引擎均已补全 |
+| B9 | 运动相机缓存不刷新 | WorkBuddy | 已修复：draw.py 直读 move_camera.txt |
+
+### 3.3 疑似误报的 Bug
+
+| 编号 | 描述 | 发现者 | 分析 |
+|------|------|--------|------|
+| **B1** | `run_simulation` 内 `config` 未定义 | Claude | ⚠️ Python 引擎已验证可用。`run_simulation` 虽不接收 `config` 参数，但在 `run_from_config` 内部调用时可能通过闭包/全局变量获得 `config`。虽建议修复（50% 概率是 Bug，50% 是侥幸可用），但优先级应降低 |
+
+---
+
+## 四、最终综合方案（3 阶段，约 8-10 人天）
+
+综合三工具六版本的全部建议，采用 **codex_v1 的三阶段框架** + **claude_v1 的详细修复步骤** + **workbuddy_v1 的对比表格**，形成以下方案。
+
+### 阶段 A：先修正确性（~2 天）— 以 Claude 为主
+
+> 目标：恢复所有功能到可用状态，禁止"带病优化"
+
+| 任务 | 参考 | 工作量 | 产出 |
+|------|------|--------|------|
+| A1 | `plot_wave.py` 适配 `load_display_txt` + header 补充 `atom_masses` | claude_v1 §1 | 1.5h | `step_plot_wave: 1` 可运行 |
+| A2 | `dynamics.py` 绘图块添加保护（短期） / 从 display.txt 重建（长期） | claude_v1 §3 | 30min | `step_plot: 1` 不崩溃 |
+| A3 | **Fortran 引擎写 display.txt**（参照 C 引擎 `write_display_txt`） | claude_v1 §4 | 2-3h | `engine: fortran` 可用 |
+| A4 | 小修复：case06 文字 / draw.py except / C++ 默认值 / 相机函数去重 / 废弃函数清理 | claude_v1 §5 | 30min | 代码整洁 |
+| A5 | B1 修复：`camera_distance/elevation/azimuth` 改为全局变量（短期补丁） | claude_v1 §2 | 15min | Python 引擎无潜在风险 |
+
+**验证标准**：4 种引擎跑 case06 均通过，`step_plot: 1` + `step_plot_wave: 1` 不崩溃
+
+### 阶段 B：再做性能优化（~3 天）— 以 Codex 框架 + Claude 代码为主
+
+> 目标：Python 引擎加速 5-15x，新增二进制格式消除 I/O 瓶颈
+
+| 任务 | 参考 | 工作量 | 预期收益 |
+|------|------|--------|---------|
+| B1 | **弹簧力向量化**：`np.add.at` 批量计算（链状体系可进一步优化为直接索引） | claude.md §P1 + claude_v1 §6 | 1h | Python 引擎 5-15x |
+| B2 | 固定约束原地掩码写回，消除 `column_stack` | claude.md §P2 | 30min | 每步减少 2 次临时数组分配 |
+| B3 | Marker 更新改为整列切片赋值 | codex_v1 阶段 B | 15min | 大原子数动画帧率提升 |
+| B4 | 驱动力去掉 `t_vec` / `frame_indices` 改计数器 | claude.md §P3-P4 | 10min | 微优化 |
+| B5 | **新增 `display.npz` 二进制格式**，draw.py 优先读二进制 | codex.md + codex_v1 | 2h | I/O 从文本级加速到二进制级 |
+| B6 | 外部引擎校准缓存 | claude_v1 §8 | 1h | 二次运行节省 10s+ |
+| B7 | `GRAVITY_INTERACTION` O(N²) 向量化（可选） | claude.md §P5 | 1h | 仅影响引力场景 |
+
+**验证标准**：`engine: python` 跑 case06 < 5s（当前 43s），`load_display_npz` < 0.01s
+
+### 阶段 C：最后做工程治理（~3-4 天）— 以 WorkBuddy 为主
+
+> 目标：从"能用"变成"好维护"
+
+| 任务 | 参考 | 工作量 | 说明 |
+|------|------|--------|------|
+| C1 | `compute.py` 模块拆分（`io.py` → `params.py` → `core.py` → `runner.py`） | workbuddy.md §1.1 | 4-6h | io.py 优先（与物理解耦），runner.py 最后 |
+| C2 | 全局变量 → `SimulationState` 类（是 A5 补丁的长期替代） | workbuddy.md §3.1 + claude_v1 §10 | 3-4h | 从 `compute_force` 开始，逐步替换 |
+| C3 | 统一 6 案例 input.txt 格式 + `ball_color` 列表化 | workbuddy.md §4.1-4.2 | 1-2h | 所有案例含 `save_trajectory`、`camera_*` 字段 |
+| C4 | pytest 单元测试（物理算法 + I/O 往返 + 参数边界） | workbuddy.md §5.1 | 3-5h | 简谐振子解析解验证 leapfrog |
+| C5 | 多引擎一致性测试 + profiling 计时 | codex_v1 + claude_v1 | 2h | 4 种引擎输出容差内一致 |
+| C6 | C/C++ 引擎共用公共头文件 | workbuddy.md §1.2 | 2-3h | `engines/common/` |
+| C7 | README 更新 / 废弃脚本清理 / 注释完善 | 综合 | 1h | 匹配当前架构 |
+
+**验证标准**：全部案例 + 全部引擎通过 pytest，README 与代码一致
+
+---
+
+## 五、三阶段总工作量与依赖
+
+```
+阶段 A（修复）   阶段 B（提速）   阶段 C（治理）
+   A1 1.5h ───── B1 1h ─────  C1 4-6h
+   A2 30min ──── B2 30min ─── C2 3-4h
+   A3 2-3h ───── B3 15min      C3 1-2h
+   A4 30min ──── B4 10min ─── C4 3-5h
+   A5 15min ──── B5 2h ────── C5 2h
+                  B6 1h         C6 2-3h
+                  B7 1h(可选)   C7 1h
+ ───────────── ───────────── ─────────────
+   ~5h           ~6h           ~15-20h
+   ~1天          ~1天          ~3-4天
+                               
+                总人天：约 8-10 天
+```
+
+**依赖规则**：
+- A 阶段无外部依赖，可直接开始
+- B 阶段依赖 A 阶段完成（否则优化可能建立在错误的代码上）
+- C 阶段依赖 B 阶段完成（架构重构应在性能热点明确后再做）
+
+---
+
+## 六、如果只做 5 件事
+
+基于三工具六版本全文分析共识，最高回报的 5 件事：
+
+| 排名 | 任务 | 工作量 | 来源 | 投入产出比 |
+|------|------|--------|------|-----------|
+| ⭐1 | **弹簧力向量化**（B1） | 1h | Claude 代码 | Python 引擎 5-15x，所有案例通用 |
+| ⭐2 | **Fortran 引擎写 display.txt**（A3） | 2-3h | 双方共识 P0 | 恢复第 4 种引擎，消除功能盲区 |
+| ⭐3 | **plot_wave.py 格式适配**（A1） | 1.5h | Claude B3 | 恢复波形动画功能（case05/06） |
+| ⭐4 | **display.npz 二进制格式**（B5） | 2h | Codex 独有 | 消除 I/O 瓶颈，所有引擎受益 |
+| ⭐5 | **校准缓存**（B6） | 1h | claude_v1 §8 | 二次运行节省 10s+，小投入大回报 |
+
+---
+
+## 七、工具选择速查表
+
+| 你需要 | 最佳工具 | 次选 | 为什么 |
+|--------|---------|------|--------|
+| 找会崩溃的 Bug | **Claude** | — | B1-B5 五个确认 Bug 全来自 Claude |
+| Python 引擎向量化 | **Claude** | Codex | 给出 `np.add.at` 可运行代码 |
+| 二进制格式策略 | **Codex** | — | 唯一想到 `display.npz` 的工具 |
+| 模块拆分方案 | **WorkBuddy** | — | 唯一给出完整目录结构的 |
+| 实施顺序规划 | **Codex** | WorkBuddy | "先修正确性→再提速→最后治理"逻辑最清晰 |
+| 多引擎一致性 | **Claude** | WorkBuddy | C/C++/Fortran 逐一评估 |
+| 配置统一 & 模板 | **WorkBuddy** | — | 提供了 6 案例统一 YAML 模板 |
+| 快速执行清单 | **codex_v1** | claude_v1 | 按"今天/本周/下阶段"列优先级 |
+| 对比多种方案 | **workbuddy_v1/v2** | claude_v1 | Bug 发现/优化项/风格三维对比表 |
+
+---
+
+## 八、附录：从 v1 到 v2 的改进
+
+| 维度 | workbuddy_v1 | workbuddy_v2 (当前) |
+|------|-------------|-------------------|
+| 参考文档数 | 3 份 | 6 份（含 3 份 v1 版本） |
+| 自我评价 | ❌ 无 | ✅ 新增 §1.2 自我修正 |
+| 阶段数 | 6 阶段 | 3 阶段（更精简） |
+| Bug 验证 | 🟡 部分标注 | ✅ 分类"真 Bug / 已修复 / 疑似误报" |
+| 工具角色定位 | 散落在各处 | ✅ 统一框架（QA/Perf/Architect） |
+| 快速执行 | 3 件事 | 5 件事（更全面） |
+| 工具选择速查 | 6 场景 | 10 场景（覆盖更全） |
+| 篇幅 | 216 行 | 约 250+ 行 |
+
+---
+
+*本文档综合了 WorkBuddy（实机验证 + 架构视野）、Claude Sonnet 4.6（Bug 挖掘 + 向量化代码）、Codex/GPT-4o（战略思维 + 二进制格式）三方共六份分析报告。建议按 A → B → C 三阶段顺序实施。*

plot_wave.py

@@ -14,17 +14,101 @@ import matplotlib.pyplot as plt
 from matplotlib.animation import FuncAnimation
 import os
 import sys
+import json
 
 sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
 import compute
 
 
 def load_disp_data(output_dir):
-    """加载 display.txt"""
-    disp_path = os.path.join(output_dir, "display.txt")
-    if not os.path.exists(disp_path):
-        raise FileNotFoundError(f"找不到 {disp_path}")
-    return compute.load_text_data(disp_path)
+    """加载 display.npz（优先）或 display.txt。"""
+    npz_path = os.path.join(output_dir, "display.npz")
+    txt_path = os.path.join(output_dir, "display.txt")
+    if os.path.exists(npz_path):
+        return compute.load_display_npz(npz_path)
+    if os.path.exists(txt_path):
+        return compute.load_display_txt(txt_path)
+    raise FileNotFoundError(f"找不到 display.npz 或 display.txt in {output_dir}")
+
+
+def _header_json(header_fields, key, default):
+    raw = header_fields.get(key, "")
+    if not raw:
+        return default
+    try:
+        return json.loads(raw)
+    except json.JSONDecodeError:
+        return default
+
+
+def _load_wave_dataset(output_dir):
+    """Load wave/energy plotting data from display metadata or sibling input files."""
+    disp_data = load_disp_data(output_dir)
+    header = disp_data["header_fields"]
+
+    x = disp_data["frames_x"]
+    y = disp_data["frames_y"]
+    z = disp_data["frames_z"]
+    vx = disp_data["frames_vx"]
+    vy = disp_data["frames_vy"]
+    vz = disp_data["frames_vz"]
+    atom_ids = np.array(disp_data["atom_ids"], dtype=np.int64)
+    n_frames = x.shape[0]
+    dt = float(header.get("DT", 0.001))
+    nstep = int(header.get("NSTEP", 1))
+    t = np.arange(n_frames, dtype=np.float64) * dt * nstep
+
+    masses = np.array(_header_json(header, "atom_masses", []), dtype=np.float64)
+    pos_0 = np.array(_header_json(header, "atom_positions", []), dtype=np.float64)
+    bond_pairs = np.array(_header_json(header, "bond_pairs", []), dtype=np.int64)
+    bond_stiffness = np.array(_header_json(header, "bond_stiffness", []), dtype=np.float64)
+    bond_rest_lengths = np.array(_header_json(header, "bond_rest_lengths", []), dtype=np.float64)
+    gravity_vec = _header_json(header, "G", [0.0, 0.0, 0.0])
+
+    # Backward-compatible fallback for older display.txt outputs.
+    if masses.size == 0 or pos_0.size == 0:
+        input_dir = os.path.join(os.path.dirname(output_dir), "input")
+        coord_path = os.path.join(input_dir, "coord.txt")
+        if os.path.exists(coord_path):
+            (_, masses_fb, _, positions_fb, _, _) = compute.load_coord_file(coord_path)
+            masses = np.array(masses_fb, dtype=np.float64)
+            pos_0 = np.array(positions_fb, dtype=np.float64)
+        else:
+            raise ValueError("display.txt 缺少 atom_masses/atom_positions 元数据，且未找到 input/coord.txt")
+
+    if bond_pairs.size == 0:
+        input_dir = os.path.join(os.path.dirname(output_dir), "input")
+        connection_path = os.path.join(input_dir, "connection.txt")
+        bond_path = os.path.join(input_dir, "bond.txt")
+        if os.path.exists(connection_path) and os.path.exists(bond_path):
+            bond_map = compute.load_bond_parameters(bond_path)
+            pairs_fb, _, stiffness_fb, rest_lengths_fb = compute.load_bond_connections(
+                connection_path, atom_ids, pos_0, bond_map)
+            bond_pairs = np.array(pairs_fb, dtype=np.int64)
+            bond_stiffness = np.array(stiffness_fb, dtype=np.float64)
+            bond_rest_lengths = np.array(rest_lengths_fb, dtype=np.float64)
+
+    return {
+        "n_frames": n_frames,
+        "t": t,
+        "x": x,
+        "y": y,
+        "z": z,
+        "vx": vx,
+        "vy": vy,
+        "vz": vz,
+        "pos_0": pos_0,
+        "masses": masses,
+        "atom_ids": atom_ids,
+        "bond_pairs": bond_pairs,
+        "bond_stiffness": bond_stiffness,
+        "bond_rest_lengths": bond_rest_lengths,
+        "gravity_field": int(header.get("gravity_field", 0)),
+        "gravity_interaction": int(header.get("gravity_interaction", 0)),
+        "gravity_strength": float(header.get("gravity_strength", 1.0)),
+        "G": gravity_vec,
+        "driving_force": int(header.get("driving_force", 0)),
+    }
 
 
 def compute_energy(x, y, z, vx, vy, vz, masses, mass_arr,
@@ -91,36 +175,36 @@ def plot_wave(output_dir, save_gif=False, save_mp4=False):
         save_gif:   是否保存 GIF
         save_mp4:   是否保存 MP4
     """
-    data = load_disp_data(output_dir)
+    data = _load_wave_dataset(output_dir)
 
     n_frames = int(data["n_frames"])
-    t = np.array(data["disp_t"])
+    t = np.array(data["t"])
 
     # 位置 / 速度
-    x  = np.array(data["disp_all_x"])
-    y  = np.array(data["disp_all_y"])
-    z  = np.array(data["disp_all_z"])
-    vx = np.array(data["disp_all_vx"])
-    vy = np.array(data["disp_all_vy"])
-    vz = np.array(data["disp_all_vz"])
+    x = np.array(data["x"])
+    y = np.array(data["y"])
+    z = np.array(data["z"])
+    vx = np.array(data["vx"])
+    vy = np.array(data["vy"])
+    vz = np.array(data["vz"])
 
     # 原子信息
-    pos_0    = np.array(data["atom_positions"])  # (n_atoms, 3)
-    masses   = np.array(data["atom_masses"])
+    pos_0 = np.array(data["pos_0"])
+    masses = np.array(data["masses"])
     atom_ids = np.array(data["atom_ids"])
-    n_atoms  = len(atom_ids)
+    n_atoms = len(atom_ids)
 
     # 成键
-    bond_pairs        = data.get("bond_pairs", [])
-    bond_stiffness    = np.array(data.get("bond_stiffness", []))
-    bond_rest_lengths = np.array(data.get("bond_rest_lengths", []))
+    bond_pairs = np.array(data.get("bond_pairs", []), dtype=np.int64)
+    bond_stiffness = np.array(data.get("bond_stiffness", []), dtype=np.float64)
+    bond_rest_lengths = np.array(data.get("bond_rest_lengths", []), dtype=np.float64)
 
     # 物理开关
-    gravity_field       = int(data.get("gravity_field", 0))
+    gravity_field = int(data.get("gravity_field", 0))
     gravity_interaction = int(data.get("gravity_interaction", 0))
-    G                   = data.get("G", [0, 0, 0])
-    gravity_strength    = float(data.get("gravity_strength", 1.0))
-    driving_force       = int(data.get("driving_force", 0))
+    G = data.get("G", [0, 0, 0])
+    gravity_strength = float(data.get("gravity_strength", 1.0))
+    driving_force = int(data.get("driving_force", 0))
 
     # ── 位移（偏离初始平衡位形）──
     dx = x - pos_0[np.newaxis, :, 0]   # 纵波（沿链方向 x）