feat: 新增波形能量动画系统 plot_wave.py

- 创建 plot_wave.py: 从 display.txt 读取原子位移数据
  绘制纵波(x) + 横波(y) + 横波(z) 波形随时间的动画
  同时绘制系统动能/弹性势能/总能量/输入功率(dE/dt)时变曲线
  输出 wave_animation.gif
- 所有 input.txt 新增 step_plot_wave: 0 开关
- case05 开启 step_plot_wave: 1
- dynamics.py disp_data 新增 bond_stiffness/bond_rest_lengths
- 更新案例文档

compute.py

@@ -697,6 +697,7 @@ def run_engine(engine, input_dir, output_dir, config):
         "elastic_force": int(config.get("elastic_force", 1)),
         "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)),
     }
     param_path = os.path.join(script_dir, "engines", engine, "param.json")
     os.makedirs(os.path.dirname(param_path), exist_ok=True)

draw.py

@@ -317,10 +317,11 @@ def update_camera_info(event=None):
     c = view.camera
     camera_info.text = (
         "Camera\n"
-        f"center = ({c.center[0]:.2f}, {c.center[1]:.2f}, {c.center[2]:.2f})\n"
+        f"center = ({c.center[0]:.1f}, {c.center[1]:.1f}, {c.center[2]:.1f})\n"
-        f"distance = {c.distance:.2f}\n"
+        f"distance = {c.distance:.1f}\n"
-        f"elevation = {c.elevation:.2f}\n"
+        f"elevation = {c.elevation:.1f}\n"
-        f"azimuth = {c.azimuth:.2f}"
+        f"azimuth = {c.azimuth:.1f}\n"
+        f"step = {PAN_SPEED:.1f}  |  {'透视' if _PERSPECTIVE else '正交'}"
     )
 
 
@@ -333,6 +334,7 @@ def update_ball_info(frame_idx, x, y, z, vx, vy, vz):
         f"t = {t:.2f} s  |  dt = {DT:.3f} s  |  nstep = {NSTEP}\n"
         f"Position: ({x:.2f}, {y:.2f}, {z:.2f})\n"
         f"Velocity: ({vx:.2f}, {vy:.2f}, {vz:.2f})\n"
+        f"W/S 沿Z轴  |  A/D 左右  |  Q/E 升降  |  C/X 步长  |  V 透视/正交"
     )
 
 
@@ -358,28 +360,99 @@ def reposition_camera_info(event=None):
     update_camera_info()
 
 
+# ── 平移速度 & 投影模式（全局变量）──
+PAN_SPEED = 1.0
+_PERSPECTIVE = True          # True=透视, False=正交
+
+
 def handle_view_interaction(event):
     update_camera_info()
 
 
-def rotate_about_screen_normal(angle):
-    if hasattr(view.camera, "roll"):
-        view.camera.roll = (view.camera.roll + angle) % 360
-    else:
-        view.camera.azimuth = (view.camera.azimuth + angle) % 360
-    update_camera_info()
-
-
 def handle_key_press(event):
+    global PAN_SPEED, _PERSPECTIVE
     key_name = ""
     if getattr(event, "text", None):
         key_name = event.text.lower()
     elif getattr(event, "key", None) is not None:
         key_name = str(event.key).lower()
-    if key_name == "q":
+
-        rotate_about_screen_normal(-90)
+    c = view.camera
-    elif key_name == "e":
+
-        rotate_about_screen_normal(90)
+    # ── 投影切换 ──
+    if key_name == "v":
+        _PERSPECTIVE = not _PERSPECTIVE
+        try:
+            if _PERSPECTIVE:
+                c.fov = 60.0
+            else:
+                c.fov = 0.0
+            print(f"[draw] 投影模式: {'透视' if _PERSPECTIVE else '正交'}")
+        except Exception:
+            print("[draw] 当前相机不支持 fov 切换")
+        update_camera_info()
+        return
+
+    # ── 步长控制 ──
+    if key_name == "c":
+        PAN_SPEED = min(PAN_SPEED * 1.5, 50.0)
+        print(f"[draw] 步长: {PAN_SPEED:.1f}")
+        update_camera_info()
+        return
+    elif key_name == "x":
+        PAN_SPEED = max(PAN_SPEED / 1.5, 0.05)
+        print(f"[draw] 步长: {PAN_SPEED:.1f}")
+        update_camera_info()
+        return
+
+    # ── 计算相机方向向量 ──
+    import math as _math
+    azim_rad = _math.radians(c.azimuth)
+    elev_rad = _math.radians(c.elevation)
+
+    # 视线方向（从 center 指向相机）
+    vd = np.array([
+        _math.cos(elev_rad) * _math.sin(azim_rad),
+        _math.sin(elev_rad),
+        _math.cos(elev_rad) * _math.cos(azim_rad),
+    ])
+    vd /= np.linalg.norm(vd)
+
+    # 屏幕右方向
+    world_up = np.array([0.0, 1.0, 0.0])
+    right = np.cross(vd, world_up)
+    rn = np.linalg.norm(right)
+    if rn > 1e-10:
+        right /= rn
+    else:
+        right = np.array([1.0, 0.0, 0.0])
+
+    # 屏幕上方向
+    up = np.cross(right, vd)
+    un = np.linalg.norm(up)
+    if un > 1e-10:
+        up /= un
+    else:
+        up = np.array([0.0, 1.0, 0.0])
+
+    pan = PAN_SPEED * 0.3
+
+    if key_name == "a":           # 右移（屏幕右方向）
+        c.center = tuple(np.array(c.center) + right * pan)
+    elif key_name == "d":         # 左移（屏幕右方向负向）
+        c.center = tuple(np.array(c.center) - right * pan)
+    elif key_name == "e":         # 上升（屏幕上方向，原 W 的功能）
+        c.center = tuple(np.array(c.center) + up * pan)
+    elif key_name == "q":         # 下降（屏幕上方向负向，原 S 的功能）
+        c.center = tuple(np.array(c.center) - up * pan)
+    elif key_name == "w":         # 相机沿 Z 轴上移（靠近场景）
+        c.center = (c.center[0], c.center[1], c.center[2] + pan)
+    elif key_name == "s":         # 相机沿 Z 轴下移（远离场景）
+        c.center = (c.center[0], c.center[1], c.center[2] - pan)
+    else:
+        return
+
+    update_camera_info()
 
 
 def reset_camera_view():

dynamics.py

@@ -125,7 +125,7 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
             print(f"[run] 同时指定了 T_total 和 NT，使用 NT={config['NT']}")
 
         # 显示步骤控制信息
-        steps_info = {k: config.get(k, 1) for k in ["step_simulate", "step_sample", "step_plot", "step_animation"]}
+        steps_info = {k: config.get(k, 1) for k in ["step_simulate", "step_sample", "step_plot", "step_plot_wave", "step_animation"]}
         step_flags = ", ".join(f"{k}={v}" for k, v in steps_info.items())
         print(f"[run] 步骤开关: {step_flags}")
 
@@ -295,6 +295,8 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
                 "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,
@@ -500,6 +502,15 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
         else:
             print("[run] 运行 python draw.py 查看动画。")
 
+        # 6. 波形能量动画（可选）
+        if config.get("step_plot_wave", 0):
+            try:
+                import plot_wave as pw
+                print("[run] 正在绘制波形与能量图…")
+                pw.plot_wave(str(output_dir_abs))
+            except Exception as e:
+                print(f"[run] 绘制波形图失败: {e}")
+
 
 def main():
     parser = argparse.ArgumentParser(description="物理模拟统一入口")

engines/c/main.c

@@ -37,6 +37,7 @@ typedef struct {
     int     elastic_force;       /* 弹簧键力开关 */
     int     damping_force;       /* 阻尼开关 */
     double  gravity_strength;    /* 万有引力强度 */
+    int     driving_force;       /* 驱动力开关 */
 } SimParams;
 
 /* ========================================================================
@@ -62,6 +63,139 @@ typedef struct {
     double *rest_lengths;
 } BondData;
 
+/* 前向声明 */
+static void *xmalloc(size_t sz);
+
+/* ========================================================================
+ * 驱动力数据
+ * ======================================================================== */
+typedef struct {
+    int    n_drivers;
+    int    *atom_idx;
+    double *amp_x, *amp_y, *amp_z;
+    double *freq_x, *freq_y, *freq_z;
+    double *phi_x, *phi_y, *phi_z;   /* radians */
+    int    *has_period;    /* 0=all, 1=limited cycles */
+    double *period_cycles; /* number of cycles */
+    double *freeze_x, *freeze_y, *freeze_z;
+} DriverData;
+
+/* 读取 driver.txt */
+static DriverData read_driver(const char *input_dir, const AtomData *atoms) {
+    DriverData d;
+    memset(&d, 0, sizeof(d));
+
+    char path[512];
+    snprintf(path, sizeof(path), "%s/driver.txt", input_dir);
+    FILE *f = fopen(path, "r");
+    if (!f) return d;
+
+    char line[1024];
+    if (!fgets(line, sizeof(line), f)) { fclose(f); return d; }
+
+    /* 第一遍：统计行数 */
+    int n_lines = 0;
+    while (fgets(line, sizeof(line), f)) {
+        char trimmed[1024];
+        int j = 0;
+        for (int i = 0; line[i]; i++) {
+            if (line[i] != ' ' && line[i] != '\t' && line[i] != '\n' && line[i] != '\r')
+                trimmed[j++] = line[i];
+        }
+        trimmed[j] = '\0';
+        if (strlen(trimmed) > 0 && trimmed[0] != '#') n_lines++;
+    }
+
+    if (n_lines == 0) { fclose(f); return d; }
+
+    /* 分配内存 */
+    d.n_drivers = n_lines;
+    d.atom_idx     = (int*)xmalloc(n_lines * sizeof(int));
+    d.amp_x = (double*)xmalloc(n_lines * sizeof(double));
+    d.amp_y = (double*)xmalloc(n_lines * sizeof(double));
+    d.amp_z = (double*)xmalloc(n_lines * sizeof(double));
+    d.freq_x = (double*)xmalloc(n_lines * sizeof(double));
+    d.freq_y = (double*)xmalloc(n_lines * sizeof(double));
+    d.freq_z = (double*)xmalloc(n_lines * sizeof(double));
+    d.phi_x = (double*)xmalloc(n_lines * sizeof(double));
+    d.phi_y = (double*)xmalloc(n_lines * sizeof(double));
+    d.phi_z = (double*)xmalloc(n_lines * sizeof(double));
+    d.has_period = (int*)xmalloc(n_lines * sizeof(int));
+    d.period_cycles = (double*)xmalloc(n_lines * sizeof(double));
+    d.freeze_x = (double*)xmalloc(n_lines * sizeof(double));
+    d.freeze_y = (double*)xmalloc(n_lines * sizeof(double));
+    d.freeze_z = (double*)xmalloc(n_lines * sizeof(double));
+
+    /* 初始化 freeze 数组 */
+    for (int i = 0; i < n_lines; i++) {
+        d.freeze_x[i] = d.freeze_y[i] = d.freeze_z[i] = 0.0;
+    }
+
+    /* 第二遍：解析 */
+    rewind(f);
+    fgets(line, sizeof(line), f); /* 跳过表头 */
+
+    int idx = 0;
+    while (idx < n_lines && fgets(line, sizeof(line), f)) {
+        char trimmed[1024];
+        int j = 0;
+        for (int i = 0; line[i]; i++) {
+            if (line[i] != ' ' && line[i] != '\t' && line[i] != '\n' && line[i] != '\r')
+                trimmed[j++] = line[i];
+        }
+        trimmed[j] = '\0';
+        if (strlen(trimmed) == 0 || trimmed[0] == '#') continue;
+
+        int atom_id;
+        double amp_x, amp_y, amp_z;
+        double freq_x, freq_y, freq_z;
+        double phi_x, phi_y, phi_z;
+        char period_str[256] = {0};
+
+        int n_parsed = sscanf(line,
+            "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %255s",
+            &atom_id,
+            &amp_x, &amp_y, &amp_z,
+            &freq_x, &freq_y, &freq_z,
+            &phi_x, &phi_y, &phi_z,
+            period_str);
+
+        if (n_parsed < 11) continue;
+
+        /* 通过原子 ID 匹配内部索引（线性搜索）*/
+        int ii = -1;
+        for (int k = 0; k < atoms->n_atoms; k++) {
+            if (atoms->atom_ids[k] == atom_id) { ii = k; break; }
+        }
+        if (ii < 0) continue;
+
+        d.atom_idx[idx] = ii;
+        d.amp_x[idx] = amp_x;
+        d.amp_y[idx] = amp_y;
+        d.amp_z[idx] = amp_z;
+        d.freq_x[idx] = freq_x;
+        d.freq_y[idx] = freq_y;
+        d.freq_z[idx] = freq_z;
+        /* 角度 → 弧度 */
+        d.phi_x[idx] = phi_x * M_PI / 180.0;
+        d.phi_y[idx] = phi_y * M_PI / 180.0;
+        d.phi_z[idx] = phi_z * M_PI / 180.0;
+
+        if (strcmp(period_str, "all") == 0 || strcmp(period_str, "-1") == 0) {
+            d.has_period[idx] = 0;
+            d.period_cycles[idx] = -1.0;
+        } else {
+            d.has_period[idx] = 1;
+            d.period_cycles[idx] = strtod(period_str, NULL);
+        }
+        idx++;
+    }
+    d.n_drivers = idx;
+
+    fclose(f);
+    return d;
+}
+
 /* ========================================================================
  * 轨迹缓冲区
  * ======================================================================== */
@@ -169,6 +303,7 @@ static SimParams read_params(const char *path) {
     p.elastic_force       = json_read_int(buf, "elastic_force");
     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");
     g_gravity_field       = p.gravity_field;
     g_gravity_interaction = p.gravity_interaction;
     g_elastic_force       = p.elastic_force;
@@ -529,6 +664,63 @@ static void leapfrog_step(
     }
 }
 
+/* ── 驱动力（与 Python apply_driving_force 一致）──────────────── */
+static void apply_driving_force(
+    int n, double *x, double *y, double *z,
+    double *vx, double *vy, double *vz,
+    double t, int step, double dt,
+    const DriverData *drivers)
+{
+    if (!drivers || drivers->n_drivers == 0) return;
+    for (int d = 0; d < drivers->n_drivers; d++) {
+        int idx = drivers->atom_idx[d];
+        /* 检查周期限制 */
+        if (drivers->has_period[d]) {
+            double max_freq = fmax(fabs(drivers->freq_x[d]),
+                                   fmax(fabs(drivers->freq_y[d]), fabs(drivers->freq_z[d])));
+            int period_steps = 0;
+            if (max_freq > 1e-12) {
+                period_steps = (int)(drivers->period_cycles[d] / max_freq / dt);
+            }
+            if (step > period_steps) {
+                /* 冻结 */
+                if (drivers->freeze_x) {
+                    x[idx] = drivers->freeze_x[d];
+                    y[idx] = drivers->freeze_y[d];
+                    z[idx] = drivers->freeze_z[d];
+                }
+                vx[idx] = vy[idx] = vz[idx] = 0.0;
+                continue;
+            }
+        }
+
+        double px = drivers->amp_x[d] * cos(2*M_PI*drivers->freq_x[d]*t + drivers->phi_x[d]);
+        double py = drivers->amp_y[d] * cos(2*M_PI*drivers->freq_y[d]*t + drivers->phi_y[d]);
+        double pz = drivers->amp_z[d] * cos(2*M_PI*drivers->freq_z[d]*t + drivers->phi_z[d]);
+        double vpx = -drivers->amp_x[d]*2*M_PI*drivers->freq_x[d]*sin(2*M_PI*drivers->freq_x[d]*t + drivers->phi_x[d]);
+        double vpy = -drivers->amp_y[d]*2*M_PI*drivers->freq_y[d]*sin(2*M_PI*drivers->freq_y[d]*t + drivers->phi_y[d]);
+        double vpz = -drivers->amp_z[d]*2*M_PI*drivers->freq_z[d]*sin(2*M_PI*drivers->freq_z[d]*t + drivers->phi_z[d]);
+
+        x[idx] = px; y[idx] = py; z[idx] = pz;
+        vx[idx] = vpx; vy[idx] = vpy; vz[idx] = vpz;
+
+        /* 记录冻结位置（周期结束时） */
+        if (drivers->has_period[d]) {
+            double max_freq = fmax(fabs(drivers->freq_x[d]),
+                                   fmax(fabs(drivers->freq_y[d]), fabs(drivers->freq_z[d])));
+            int period_steps = 0;
+            if (max_freq > 1e-12) {
+                period_steps = (int)(drivers->period_cycles[d] / max_freq / dt);
+            }
+            if (step == period_steps) {
+                drivers->freeze_x[d] = px;
+                drivers->freeze_y[d] = py;
+                drivers->freeze_z[d] = pz;
+            }
+        }
+    }
+}
+
 /* ── 分发器：调用对应积分方法 + 边界条件 + 自由度约束（与 Python 一致）── */
 static void apply_step(
     const char *method,
@@ -643,7 +835,8 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
         if (b > 0) fputc(',', f);
         fprintf(f, "%.15g", bonds->rest_lengths[b]);
     }
-    fprintf(f, "]\n");
+    fprintf(f, "],\n");
+    fprintf(f, "  \"driving_force\": %d\n", params->driving_force);
 
     fprintf(f, "}\n");
     fclose(f);
@@ -668,8 +861,14 @@ int main(int argc, char **argv) {
     AtomData  atoms  = read_coord(input_dir);
     BondData  bonds  = read_bonds(input_dir, &atoms);
 
-    printf("[C-engine] 原子数=%d, 键数=%d, NT=%d, DT=%.6g, method=%s\n",
+    DriverData drivers;
-           atoms.n_atoms, bonds.n_bonds, params.NT, params.DT, params.method);
+    drivers.n_drivers = 0;
+    if (params.driving_force) {
+        drivers = read_driver(input_dir, &atoms);
+    }
+
+    printf("[C-engine] 原子数=%d, 键数=%d, 驱动=%d, NT=%d, DT=%.6g, method=%s\n",
+           atoms.n_atoms, bonds.n_bonds, drivers.n_drivers, params.NT, params.DT, params.method);
 
     int n = atoms.n_atoms;
     double *x  = (double*)xmalloc(n * sizeof(double));
@@ -701,6 +900,8 @@ int main(int argc, char **argv) {
 
     /* 预热 */
     for (int s = 0; s < params.warmup_steps; s++) {
+        double tw = (s + 1) * params.DT;
+        if (params.driving_force) apply_driving_force(n, x, y, z, vx, vy, vz, tw, s, params.DT, &drivers);
         apply_step(params.method, n, x, y, z, vx, vy, vz,
                    atoms.masses, params.G, params.B, &bonds, atoms.fixed,
                    atoms.pos_0,
@@ -709,6 +910,8 @@ int main(int argc, char **argv) {
 
     /* 记录 */
     for (int s = 0; s < record_steps; s++) {
+        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];

engines/cpp/main.cpp

@@ -42,6 +42,7 @@ struct SimParams {
     int elastic_force = 1;
     int damping_force = 0;
     double gravity_strength = 1.0;
+    int driving_force = 0;
 };
 
 // ========================================================================
@@ -65,6 +66,20 @@ struct BondData {
     std::vector<double> rest_lengths;
 };
 
+// ========================================================================
+// 驱动力数据
+// ========================================================================
+struct DriverData {
+    int n_drivers = 0;
+    std::vector<int> atom_idx;       // internal atom indices
+    std::vector<double> amp_x, amp_y, amp_z;
+    std::vector<double> freq_x, freq_y, freq_z;
+    std::vector<double> phi_x, phi_y, phi_z;  // radians
+    std::vector<int> has_period;     // 0=all, 1=limited
+    std::vector<double> period_cycles;
+    std::vector<double> freeze_x, freeze_y, freeze_z;
+};
+
 // ========================================================================
 // 辅助函数
 // ========================================================================
@@ -149,6 +164,7 @@ static SimParams read_params(const std::string &path) {
     p.elastic_force       = json_read_int(buf, "elastic_force");
     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");
     return p;
 }
 
@@ -233,6 +249,57 @@ static BondData read_bonds(const std::string &input_dir) {
     return b;
 }
 
+/* 读取 driver.txt */
+static DriverData read_driver(const std::string &input_dir, const AtomData &atoms) {
+    DriverData d;
+    std::string path = input_dir + "/driver.txt";
+    std::ifstream f(path);
+    if (!f) { std::cerr << "[C++-engine] 警告: 无法打开 " << path << std::endl; return d; }
+
+    std::string header;
+    std::getline(f, header);  // skip header
+
+    int n;
+    double ax, ay, az, fx, fy, fz, px, py, pz;
+    std::string period_str;
+
+    while (f >> n >> ax >> ay >> az >> fx >> fy >> fz >> px >> py >> pz >> period_str) {
+        // Find atom index by id
+        int idx = -1;
+        for (size_t i = 0; i < atoms.ids.size(); i++) {
+            if (atoms.ids[i] == n) { idx = i; break; }
+        }
+        if (idx < 0) {
+            std::cerr << "[C++-engine] 警告: driver.txt 原子 " << n << " 不在 coord.txt 中" << std::endl;
+            continue;
+        }
+        d.atom_idx.push_back(idx);
+        d.amp_x.push_back(ax); d.amp_y.push_back(ay); d.amp_z.push_back(az);
+        d.freq_x.push_back(fx); d.freq_y.push_back(fy); d.freq_z.push_back(fz);
+        // Convert degrees to radians
+        const double DEG2RAD = M_PI / 180.0;
+        d.phi_x.push_back(px * DEG2RAD);
+        d.phi_y.push_back(py * DEG2RAD);
+        d.phi_z.push_back(pz * DEG2RAD);
+
+        if (period_str == "all") {
+            d.has_period.push_back(0);
+            d.period_cycles.push_back(-1.0);
+        } else {
+            d.has_period.push_back(1);
+            d.period_cycles.push_back(std::stod(period_str));
+        }
+        d.freeze_x.push_back(0.0);
+        d.freeze_y.push_back(0.0);
+        d.freeze_z.push_back(0.0);
+        d.n_drivers++;
+    }
+
+    if (d.n_drivers > 0)
+        std::cout << "[C++-engine] 已加载驱动力: " << d.n_drivers << " 条定义" << std::endl;
+    return d;
+}
+
 // ========================================================================
 // 物理核心
 // ========================================================================
@@ -539,6 +606,68 @@ static void apply_step(
     }
 }
 
+// ========================================================================
+// 驱动力应用
+// ========================================================================
+
+static void apply_driving_force(
+    int n, double *x, double *y, double *z,
+    double *vx, double *vy, double *vz,
+    double t, int step, double dt,
+    DriverData &drivers)
+{
+    if (drivers.n_drivers == 0) return;
+    for (int d = 0; d < drivers.n_drivers; d++) {
+        int idx = drivers.atom_idx[d];
+
+        // Check period limits
+        if (drivers.has_period[d]) {
+            double max_freq = std::max({std::fabs(drivers.freq_x[d]),
+                                        std::fabs(drivers.freq_y[d]),
+                                        std::fabs(drivers.freq_z[d])});
+            int period_steps = 0;
+            if (max_freq > 1e-12) {
+                period_steps = static_cast<int>(drivers.period_cycles[d] / max_freq / dt);
+            }
+            if (step > period_steps) {
+                // Frozen: keep last position, zero velocity
+                x[idx] = drivers.freeze_x[d];
+                y[idx] = drivers.freeze_y[d];
+                z[idx] = drivers.freeze_z[d];
+                vx[idx] = vy[idx] = vz[idx] = 0.0;
+                continue;
+            }
+        }
+
+        const double TWO_PI = 2.0 * M_PI;
+        double px = drivers.amp_x[d] * std::cos(TWO_PI * drivers.freq_x[d] * t + drivers.phi_x[d]);
+        double py = drivers.amp_y[d] * std::cos(TWO_PI * drivers.freq_y[d] * t + drivers.phi_y[d]);
+        double pz = drivers.amp_z[d] * std::cos(TWO_PI * drivers.freq_z[d] * t + drivers.phi_z[d]);
+        double vpx = -drivers.amp_x[d] * TWO_PI * drivers.freq_x[d] * std::sin(TWO_PI * drivers.freq_x[d] * t + drivers.phi_x[d]);
+        double vpy = -drivers.amp_y[d] * TWO_PI * drivers.freq_y[d] * std::sin(TWO_PI * drivers.freq_y[d] * t + drivers.phi_y[d]);
+        double vpz = -drivers.amp_z[d] * TWO_PI * drivers.freq_z[d] * std::sin(TWO_PI * drivers.freq_z[d] * t + drivers.phi_z[d]);
+
+        x[idx] = px; y[idx] = py; z[idx] = pz;
+        vx[idx] = vpx; vy[idx] = vpy; vz[idx] = vpz;
+
+        // Record freeze position at the last driving step
+        if (drivers.has_period[d]) {
+            double max_freq = std::max({std::fabs(drivers.freq_x[d]),
+                                        std::fabs(drivers.freq_y[d]),
+                                        std::fabs(drivers.freq_z[d])});
+            int period_steps = 0;
+            if (max_freq > 1e-12) {
+                period_steps = static_cast<int>(drivers.period_cycles[d] / max_freq / dt);
+            }
+            if (step == period_steps) {
+                drivers.freeze_x[d] = px;
+                drivers.freeze_y[d] = py;
+                drivers.freeze_z[d] = pz;
+            }
+        }
+    }
+}
+
 // ========================================================================
 // JSON 输出
 // ========================================================================
@@ -621,7 +750,9 @@ static void write_trajectory_json(
         if (b > 0) f << ',';
         f << bonds.rest_lengths[b];
     }
-    f << "]\n";
+    f << "],\n";
+
+    f << "  \"driving_force\": " << params.driving_force << "\n";
 
     f << "}\n";
 }
@@ -647,6 +778,11 @@ int main(int argc, char **argv) {
     AtomData  atoms  = read_coord(input_dir);
     BondData  bonds  = read_bonds(input_dir);
 
+    DriverData drivers;
+    if (params.driving_force) {
+        drivers = read_driver(input_dir, atoms);
+    }
+
     std::cout << "[C++-engine] 原子数=" << atoms.ids.size()
               << ", 键数=" << bonds.stiffness.size()
               << ", NT=" << params.NT << ", DT=" << params.DT
@@ -676,6 +812,9 @@ int main(int argc, char **argv) {
 
     // 预热
     for (int s = 0; s < params.warmup_steps; s++) {
+        double tw = (s + 1) * params.DT;
+        if (params.driving_force)
+            apply_driving_force(n, x.data(), y.data(), z.data(), vx.data(), vy.data(), vz.data(), tw, s, params.DT, drivers);
         apply_step(params.method, n, x.data(), y.data(), z.data(),
                    vx.data(), vy.data(), vz.data(),
                    atoms.masses.data(), params.G, params.B,
@@ -688,6 +827,9 @@ int main(int argc, char **argv) {
 
     // 记录
     for (int s = 0; s < record_steps; s++) {
+        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];

engines/fortran/main.f90

@@ -18,9 +18,10 @@ program dynamics_f90
     double precision :: DT, box_a
     double precision :: G(3), B(3)
     integer :: gravity_field, gravity_interaction, elastic_force, damping_force
+    integer :: driving_force
     double precision :: gravity_strength
     character(len=32) :: method
-    double precision :: t0, t1, elapsed
+    double precision :: t0, t1, elapsed, tw
 
     ! 原子数据
     integer, allocatable :: atom_ids(:)
@@ -32,6 +33,16 @@ program dynamics_f90
     integer, allocatable :: bond_pairs(:, :)
     double precision, allocatable :: bond_stiffness(:), bond_rest_lengths(:)
 
+    ! 驱动力数据
+    integer :: n_drivers
+    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(:)
+    double precision, allocatable :: drv_phi_x(:), drv_phi_y(:), drv_phi_z(:)
+    integer, allocatable :: drv_has_period(:)
+    double precision, allocatable :: drv_period_cycles(:)
+    double precision, allocatable :: drv_freeze_x(:), drv_freeze_y(:), drv_freeze_z(:)
+
     ! 运行时位置/速度
     double precision, allocatable :: x(:), y(:), z(:)
     double precision, allocatable :: vx(:), vy(:), vz(:)
@@ -59,7 +70,7 @@ program dynamics_f90
     ! 读取 param.json
     call read_params(param_path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
                      gravity_field, gravity_interaction, &
-                     elastic_force, damping_force, gravity_strength)
+                     elastic_force, damping_force, gravity_strength, driving_force)
 
     ! 读取 coord.txt
     call read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed)
@@ -68,6 +79,17 @@ program dynamics_f90
     call read_bonds(input_dir, n_atoms, atom_ids, pos_0, n_bonds, &
                     bond_pairs, bond_stiffness, bond_rest_lengths)
 
+    ! 读取驱动力
+    n_drivers = 0
+    if (driving_force /= 0) then
+        call read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
+                         drv_atom_idx, drv_amp_x, drv_amp_y, drv_amp_z, &
+                         drv_freq_x, drv_freq_y, drv_freq_z, &
+                         drv_phi_x, drv_phi_y, drv_phi_z, &
+                         drv_has_period, drv_period_cycles, &
+                         drv_freeze_x, drv_freeze_y, drv_freeze_z)
+    end if
+
     write(*, '("[Fortran-engine] 原子数=", i0, ", 键数=", i0, &
            &", NT=", i0, ", DT=", f0.6, ", method=", a)') &
         n_atoms, n_bonds, NT, DT, trim(method)
@@ -87,6 +109,16 @@ program dynamics_f90
 
     ! 预热
     do s = 1, warmup_steps
+        if (driving_force /= 0 .and. n_drivers > 0) then
+            tw = (s * 1.0d0) * DT
+            call apply_driving(n, x, y, z, vx, vy, vz, tw, s-1, DT, &
+                               n_drivers, drv_atom_idx, &
+                               drv_amp_x, drv_amp_y, drv_amp_z, &
+                               drv_freq_x, drv_freq_y, drv_freq_z, &
+                               drv_phi_x, drv_phi_y, drv_phi_z, &
+                               drv_has_period, drv_period_cycles, &
+                               drv_freeze_x, drv_freeze_y, drv_freeze_z)
+        end if
         call apply_step(method, n, x, y, z, vx, vy, vz, masses, G, B, &
                         n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths, &
                         fixed, box_a, DT, &
@@ -97,6 +129,16 @@ program dynamics_f90
 
     ! 记录
     do s = 1, record_steps
+        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, &
+                               n_drivers, drv_atom_idx, &
+                               drv_amp_x, drv_amp_y, drv_amp_z, &
+                               drv_freq_x, drv_freq_y, drv_freq_z, &
+                               drv_phi_x, drv_phi_y, drv_phi_z, &
+                               drv_has_period, drv_period_cycles, &
+                               drv_freeze_x, drv_freeze_y, drv_freeze_z)
+        end if
         traj_x(s, :)  = x;  traj_y(s, :)  = y;  traj_z(s, :)  = z
         traj_vx(s, :) = vx; traj_vy(s, :) = vy; traj_vz(s, :) = vz
         call apply_step(method, n, x, y, z, vx, vy, vz, masses, G, B, &
@@ -111,7 +153,8 @@ program dynamics_f90
     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, &
-                    n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths)
+                    n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths, &
+                    driving_force)
 
     call cpu_time(t1)
     elapsed = t1 - t0
@@ -121,6 +164,14 @@ program dynamics_f90
     deallocate(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz)
     deallocate(atom_ids, masses, radii, pos_0, vel_0, fixed)
     if (allocated(bond_pairs)) deallocate(bond_pairs, bond_stiffness, bond_rest_lengths)
+    if (allocated(drv_atom_idx)) then
+        deallocate(drv_atom_idx)
+        deallocate(drv_amp_x, drv_amp_y, drv_amp_z)
+        deallocate(drv_freq_x, drv_freq_y, drv_freq_z)
+        deallocate(drv_phi_x, drv_phi_y, drv_phi_z)
+        deallocate(drv_has_period, drv_period_cycles)
+        deallocate(drv_freeze_x, drv_freeze_y, drv_freeze_z)
+    end if
 
 contains
 
@@ -129,11 +180,11 @@ contains
 ! ========================================================================
 subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
                        gravity_field, gravity_interaction, &
-                       elastic_force, damping_force, gravity_strength)
+                       elastic_force, damping_force, gravity_strength, driving_force)
     character(len=*), intent(in) :: path
     double precision, intent(out) :: box_a, DT, G(3), B(3), gravity_strength
     integer, intent(out) :: NT, NSTEP, warmup_steps
-    integer, intent(out) :: gravity_field, gravity_interaction, elastic_force, damping_force
+    integer, intent(out) :: gravity_field, gravity_interaction, elastic_force, damping_force, driving_force
     character(len=*), intent(out) :: method
     character(len=8096) :: buf
     character(len=256) :: line
@@ -144,6 +195,7 @@ subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
     G = (/ 0.0d0, 0.0d0, -9.8d0 /); B = 0.0d0
     gravity_field = 1; gravity_interaction = 0
     elastic_force = 1; damping_force = 0; gravity_strength = 1.0d0
+    driving_force = 0
 
     open(newunit=u, file=trim(path), status='old', action='read', iostat=ios)
     if (ios /= 0) then
@@ -173,6 +225,7 @@ subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
     elastic_force       = json_get_int(buf, 'elastic_force', 1)
     damping_force       = json_get_int(buf, 'damping_force', 0)
     gravity_strength    = json_get_double(buf, 'gravity_strength', 1.0d0)
+    driving_force       = json_get_int(buf, 'driving_force', 0)
 end subroutine read_params
 
 ! ========================================================================
@@ -685,6 +738,172 @@ subroutine apply_step(method, n, x, y, z, vx, vy, vz, m, G, B, &
     end do
 end subroutine apply_step
 
+! ========================================================================
+! 读取驱动力 driver.txt
+! ========================================================================
+subroutine read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
+                       drv_atom_idx, drv_amp_x, drv_amp_y, drv_amp_z, &
+                       drv_freq_x, drv_freq_y, drv_freq_z, &
+                       drv_phi_x, drv_phi_y, drv_phi_z, &
+                       drv_has_period, drv_period_cycles, &
+                       drv_freeze_x, drv_freeze_y, drv_freeze_z)
+    character(len=*), intent(in) :: input_dir
+    integer, intent(in) :: n_atoms, atom_ids(n_atoms)
+    integer, intent(out) :: n_drivers
+    integer, allocatable, intent(out) :: drv_atom_idx(:)
+    double precision, allocatable, intent(out) :: drv_amp_x(:), drv_amp_y(:), drv_amp_z(:)
+    double precision, allocatable, intent(out) :: drv_freq_x(:), drv_freq_y(:), drv_freq_z(:)
+    double precision, allocatable, intent(out) :: drv_phi_x(:), drv_phi_y(:), drv_phi_z(:)
+    integer, allocatable, intent(out) :: drv_has_period(:)
+    double precision, allocatable, intent(out) :: drv_period_cycles(:)
+    double precision, allocatable, intent(out) :: drv_freeze_x(:), drv_freeze_y(:), drv_freeze_z(:)
+
+    character(len=512) :: path, line, period_str
+    integer :: u, ios, i, n, idx, n_cap
+    double precision :: ax, ay, az, fx, fy, fz, px, py, pz
+    integer, parameter :: MX = 256
+    integer :: idx_tmp(MX), per_tmp(MX)
+    double precision :: ax_tmp(MX), ay_tmp(MX), az_tmp(MX)
+    double precision :: fxx_tmp(MX), fyy_tmp(MX), fzz_tmp(MX)
+    double precision :: pxx_tmp(MX), pyy_tmp(MX), pzz_tmp(MX)
+    double precision :: pc_tmp(MX)
+    double precision :: fzx_tmp(MX), fzy_tmp(MX), fzz_tmp2(MX)
+
+    n_drivers = 0
+    path = trim(input_dir) // '/driver.txt'
+    open(newunit=u, file=trim(path), status='old', action='read', iostat=ios)
+    if (ios /= 0) then
+        write(*, '("[Fortran-engine] 警告: 无法打开 ", a)') trim(path)
+        return
+    end if
+
+    read(u, '(a)', iostat=ios)  ! skip header
+
+    do
+        read(u, *, iostat=ios) n, ax, ay, az, fx, fy, fz, px, py, pz, period_str
+        if (ios /= 0) exit
+        n_drivers = n_drivers + 1
+        if (n_drivers > MX) exit
+
+        ! Find atom index by id
+        idx = -1
+        do i = 1, n_atoms
+            if (atom_ids(i) == n) then
+                idx = i
+                exit
+            end if
+        end do
+        if (idx < 0) cycle
+
+        idx_tmp(n_drivers) = idx - 1  ! 0-based index
+        ax_tmp(n_drivers) = ax; ay_tmp(n_drivers) = ay; az_tmp(n_drivers) = az
+        fxx_tmp(n_drivers) = fx; fyy_tmp(n_drivers) = fy; fzz_tmp(n_drivers) = fz
+        ! degrees to radians
+        pxx_tmp(n_drivers) = px * 3.141592653589793d0 / 180.0d0
+        pyy_tmp(n_drivers) = py * 3.141592653589793d0 / 180.0d0
+        pzz_tmp(n_drivers) = pz * 3.141592653589793d0 / 180.0d0
+
+        if (trim(period_str) == 'all') then
+            per_tmp(n_drivers) = 0
+            pc_tmp(n_drivers) = -1.0d0
+        else
+            per_tmp(n_drivers) = 1
+            read(period_str, *) pc_tmp(n_drivers)
+        end if
+        fzx_tmp(n_drivers) = 0.0d0
+        fzy_tmp(n_drivers) = 0.0d0
+        fzz_tmp2(n_drivers) = 0.0d0
+    end do
+    close(u)
+
+    if (n_drivers <= 0) return
+
+    allocate(drv_atom_idx(n_drivers))
+    allocate(drv_amp_x(n_drivers), drv_amp_y(n_drivers), drv_amp_z(n_drivers))
+    allocate(drv_freq_x(n_drivers), drv_freq_y(n_drivers), drv_freq_z(n_drivers))
+    allocate(drv_phi_x(n_drivers), drv_phi_y(n_drivers), drv_phi_z(n_drivers))
+    allocate(drv_has_period(n_drivers), drv_period_cycles(n_drivers))
+    allocate(drv_freeze_x(n_drivers), drv_freeze_y(n_drivers), drv_freeze_z(n_drivers))
+
+    drv_atom_idx = idx_tmp(1:n_drivers)
+    drv_amp_x = ax_tmp(1:n_drivers); drv_amp_y = ay_tmp(1:n_drivers); drv_amp_z = az_tmp(1:n_drivers)
+    drv_freq_x = fxx_tmp(1:n_drivers); drv_freq_y = fyy_tmp(1:n_drivers); drv_freq_z = fzz_tmp(1:n_drivers)
+    drv_phi_x = pxx_tmp(1:n_drivers); drv_phi_y = pyy_tmp(1:n_drivers); drv_phi_z = pzz_tmp(1:n_drivers)
+    drv_has_period = per_tmp(1:n_drivers)
+    drv_period_cycles = pc_tmp(1:n_drivers)
+    drv_freeze_x = fzx_tmp(1:n_drivers); drv_freeze_y = fzy_tmp(1:n_drivers); drv_freeze_z = fzz_tmp2(1:n_drivers)
+
+    write(*, '("[Fortran-engine] 已加载驱动力: ", i0, " 条定义")') n_drivers
+end subroutine read_driver
+
+! 施加驱动力
+subroutine apply_driving(n, x, y, z, vx, vy, vz, t, step, dt, &
+                         n_drivers, drv_atom_idx, &
+                         drv_amp_x, drv_amp_y, drv_amp_z, &
+                         drv_freq_x, drv_freq_y, drv_freq_z, &
+                         drv_phi_x, drv_phi_y, drv_phi_z, &
+                         drv_has_period, drv_period_cycles, &
+                         drv_freeze_x, drv_freeze_y, drv_freeze_z)
+    integer, intent(in) :: n, step, n_drivers
+    integer, intent(in) :: drv_atom_idx(n_drivers), drv_has_period(n_drivers)
+    double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
+    double precision, intent(in) :: t, dt
+    double precision, intent(in) :: drv_amp_x(n_drivers), drv_amp_y(n_drivers), drv_amp_z(n_drivers)
+    double precision, intent(in) :: drv_freq_x(n_drivers), drv_freq_y(n_drivers), drv_freq_z(n_drivers)
+    double precision, intent(in) :: drv_phi_x(n_drivers), drv_phi_y(n_drivers), drv_phi_z(n_drivers)
+    double precision, intent(in) :: drv_period_cycles(n_drivers)
+    double precision, intent(inout) :: drv_freeze_x(n_drivers), drv_freeze_y(n_drivers), drv_freeze_z(n_drivers)
+
+    integer :: d, idx, period_steps
+    double precision :: max_freq, px, py, pz, vpx, vpy, vpz, TWO_PI
+
+    TWO_PI = 2.0d0 * 3.141592653589793d0
+
+    do d = 1, n_drivers
+        idx = drv_atom_idx(d) + 1  ! Fortran 1-based indexing
+
+        ! Check period
+        if (drv_has_period(d) /= 0) then
+            max_freq = max(abs(drv_freq_x(d)), abs(drv_freq_y(d)), abs(drv_freq_z(d)))
+            period_steps = 0
+            if (max_freq > 1.0d-12) then
+                period_steps = int(drv_period_cycles(d) / max_freq / dt)
+            end if
+            if (step > period_steps) then
+                x(idx) = drv_freeze_x(d)
+                y(idx) = drv_freeze_y(d)
+                z(idx) = drv_freeze_z(d)
+                vx(idx) = 0.0d0; vy(idx) = 0.0d0; vz(idx) = 0.0d0
+                cycle
+            end if
+        end if
+
+        px = drv_amp_x(d) * cos(TWO_PI * drv_freq_x(d) * t + drv_phi_x(d))
+        py = drv_amp_y(d) * cos(TWO_PI * drv_freq_y(d) * t + drv_phi_y(d))
+        pz = drv_amp_z(d) * cos(TWO_PI * drv_freq_z(d) * t + drv_phi_z(d))
+        vpx = -drv_amp_x(d) * TWO_PI * drv_freq_x(d) * sin(TWO_PI * drv_freq_x(d) * t + drv_phi_x(d))
+        vpy = -drv_amp_y(d) * TWO_PI * drv_freq_y(d) * sin(TWO_PI * drv_freq_y(d) * t + drv_phi_y(d))
+        vpz = -drv_amp_z(d) * TWO_PI * drv_freq_z(d) * sin(TWO_PI * drv_freq_z(d) * t + drv_phi_z(d))
+
+        x(idx) = px; y(idx) = py; z(idx) = pz
+        vx(idx) = vpx; vy(idx) = vpy; vz(idx) = vpz
+
+        ! Record freeze position
+        if (drv_has_period(d) /= 0) then
+            max_freq = max(abs(drv_freq_x(d)), abs(drv_freq_y(d)), abs(drv_freq_z(d)))
+            period_steps = 0
+            if (max_freq > 1.0d-12) then
+                period_steps = int(drv_period_cycles(d) / max_freq / dt)
+            end if
+            if (step == period_steps) then
+                drv_freeze_x(d) = px
+                drv_freeze_y(d) = py
+                drv_freeze_z(d) = pz
+            end if
+        end if
+    end do
+end subroutine apply_driving
+
 ! ========================================================================
 ! JSON 输出
 ! ========================================================================
@@ -692,9 +911,9 @@ end subroutine apply_step
 subroutine write_json(outdir, tx, ty, tz, tvx, tvy, tvz, &
                       nsteps, nat, aid, amass, &
                       NT, DT, NSTEP, warmup, method, G, B, &
-                      nb, bp, bk, br)
+                      nb, bp, bk, br, driving_force)
     character(len=*), intent(in) :: outdir, method
-    integer, intent(in) :: nsteps, nat, NT, NSTEP, warmup, nb, bp(nb, 2), aid(nat)
+    integer, intent(in) :: nsteps, nat, NT, NSTEP, warmup, nb, bp(nb, 2), aid(nat), driving_force
     double precision, intent(in) :: tx(nsteps, nat), ty(nsteps, nat), tz(nsteps, nat)
     double precision, intent(in) :: tvx(nsteps, nat), tvy(nsteps, nat), tvz(nsteps, nat)
     double precision, intent(in) :: DT, G(3), B(3), bk(nb), br(nb), amass(nat)
@@ -798,6 +1017,9 @@ subroutine write_json(outdir, tx, ty, tz, tvx, tvy, tvz, &
         write(u, '(a)') '  "bond_rest_lengths": []'
     end if
 
+    write(buf, '(a, i0)') '  "driving_force": ', driving_force
+    write(u, '(a)') trim(buf)
+
     write(u, '(a)') '}'
     close(u)
 end subroutine write_json

examples/case01/input/input.txt

@@ -8,6 +8,7 @@
 step_simulate:  1   # 运行物理模拟 → output/trajectory.txt
 step_sample:    1   # 抽帧 → output/display.txt
 step_plot:      1   # 绘制轨迹/能量图 → output/trajectory_plots.png
+step_plot_wave: 0   # 绘制波形能量动画 → output/wave_animation.gif
 step_animation: 1   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
 force_calc: 0       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
 

examples/case02/input/input.txt

@@ -8,6 +8,7 @@
 step_simulate:  1   # 运行物理模拟 → output/trajectory.txt
 step_sample:    1   # 抽帧 → output/display.txt
 step_plot:      1   # 绘制轨迹/能量图 → output/trajectory_plots.png
+step_plot_wave: 0   # 绘制波形能量动画 → output/wave_animation.gif
 step_animation: 1   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
 force_calc: 0       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
 

examples/case03/input/input.txt

@@ -8,6 +8,7 @@
 step_simulate:  1   # 运行物理模拟 → output/trajectory.txt
 step_sample:    1   # 抽帧 → output/display.txt
 step_plot:      0   # 绘制轨迹/能量图 → output/trajectory_plots.png
+step_plot_wave: 0   # 绘制波形能量动画 → output/wave_animation.gif
 step_animation: 1   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
 force_calc: 0       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
 

examples/case04/input/input.txt

@@ -8,6 +8,7 @@
 step_simulate:  1   # 运行物理模拟 → output/trajectory.txt
 step_sample:    1   # 抽帧 → output/display.txt
 step_plot:      0   # 绘制轨迹/能量图 → output/trajectory_plots.png
+step_plot_wave: 0   # 绘制波形能量动画 → output/wave_animation.gif
 step_animation: 1   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
 force_calc: 0       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
 

examples/case05/doc/index.html

@@ -356,10 +356,14 @@ step_animation: 1   # 播放动画</pre>
       <tr><th>操作</th><th>效果</th></tr>
       <tr><td>鼠标拖动</td><td>旋转视角</td></tr>
       <tr><td>滚轮</td><td>缩放</td></tr>
+      <tr><td>W / S 键</td><td>相机沿 Z 轴向前 / 向后移动（靠近/远离场景）</td></tr>
+      <tr><td>A / D 键</td><td>视角向右 / 向左平移</td></tr>
+      <tr><td>E / Q 键</td><td>视角上升 / 下降（屏幕方向）</td></tr>
+      <tr><td>C / X 键</td><td>增大 / 减小步长</td></tr>
+      <tr><td>V 键</td><td>切换透视 / 正交投影</td></tr>
       <tr><td>左上角 <strong>reset</strong> 按钮</td><td>复位视角到初始位置</td></tr>
       <tr><td>左上角 <strong>info</strong> 按钮</td><td>切换信息面板显示/隐藏</td></tr>
       <tr><td>左上角 <strong>axes</strong> 按钮</td><td>切换坐标轴显示/隐藏</td></tr>
-      <tr><td>Q / E 键</td><td>画面绕视向旋转 -90° / +90°</td></tr>
     </table>
   </div>
 </section>
@@ -395,6 +399,7 @@ step_animation: 1   # 播放动画</pre>
       <tr><td>step_simulate</td><td>跳过模拟（加载已有轨迹）</td><td>运行物理模拟</td></tr>
       <tr><td>step_sample</td><td>跳过抽帧</td><td>从轨迹抽取显示帧</td></tr>
       <tr><td>step_plot</td><td>不生成图表</td><td>生成轨迹/能量图</td></tr>
+      <tr><td><strong>step_plot_wave</strong></td><td>不生成波形图</td><td>生成波形能量动画 GIF</td></tr>
       <tr><td>step_animation</td><td>不启动动画</td><td>自动打开 VisPy 3D 窗口</td></tr>
       <tr><td>force_calc</td><td>自动检测缓存</td><td>强制重新计算</td></tr>
     </table>
@@ -418,7 +423,8 @@ step_animation: 1   # 播放动画</pre>
 │   ├── trajectory.txt     # 全量轨迹数据（10000 步 × 60 原子）
 │   ├── display.txt        # 抽帧后的动画数据（200 帧 × 60 原子）
 │   ├── dynamics.log       # 计算日志
-│   └── animation.log      # 动画启动日志（闪退时排查用）
+│   ├── animation.log      # 动画启动日志（闪退时排查用）
+│   └── wave_animation.gif # 波形能量动画（step_plot_wave=1 时生成）
 ├── doc/
 │   └── index.html         # <span class="cm">本文档</span>
 ├── Readme.md              # 案例简介

examples/case05/input/bond.txt

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

examples/case05/input/coord.txt

@@ -1,61 +1,61 @@
 n mass radius x y z vx vy vz fix_x fix_y fix_z
-1 1 0.1 0 0 1 0 0 0 1 1 0
+1 1 0.1 0 0 1 0 0 0   0 1 1
-2 1 0.1 1 0 0 0 0 0 1 1 0
+2 1 0.1 1 0 0 0 0 0   0 1 1
-3 1 0.1 2 0 0 0 0 0 1 1 0
+3 1 0.1 2 0 0 0 0 0   0 1 1
-4 1 0.1 3 0 0 0 0 0 1 1 0
+4 1 0.1 3 0 0 0 0 0   0 1 1
-5 1 0.1 4 0 0 0 0 0 1 1 0
+5 1 0.1 4 0 0 0 0 0   0 1 1
-6 1 0.1 5 0 0 0 0 0 1 1 0
+6 1 0.1 5 0 0 0 0 0   0 1 1
-7 1 0.1 6 0 0 0 0 0 1 1 0
+7 1 0.1 6 0 0 0 0 0   0 1 1
-8 1 0.1 7 0 0 0 0 0 1 1 0
+8 1 0.1 7 0 0 0 0 0   0 1 1
-9 1 0.1 8 0 0 0 0 0 1 1 0
+9 1 0.1 8 0 0 0 0 0   0 1 1
-10 1 0.1 9 0 0 0 0 0 1 1 0
+10 1 0.1 9 0 0 0 0 0   0 1 1
-11 1 0.1 10 0 0 0 0 0 1 1 0
+11 1 0.1 10 0 0 0 0 0   0 1 1
-12 1 0.1 11 0 0 0 0 0 1 1 0
+12 1 0.1 11 0 0 0 0 0   0 1 1
-13 1 0.1 12 0 0 0 0 0 1 1 0
+13 1 0.1 12 0 0 0 0 0   0 1 1
-14 1 0.1 13 0 0 0 0 0 1 1 0
+14 1 0.1 13 0 0 0 0 0   0 1 1
-15 1 0.1 14 0 0 0 0 0 1 1 0
+15 1 0.1 14 0 0 0 0 0   0 1 1
-16 1 0.1 15 0 0 0 0 0 1 1 0
+16 1 0.1 15 0 0 0 0 0   0 1 1
-17 1 0.1 16 0 0 0 0 0 1 1 0
+17 1 0.1 16 0 0 0 0 0   0 1 1
-18 1 0.1 17 0 0 0 0 0 1 1 0
+18 1 0.1 17 0 0 0 0 0   0 1 1
-19 1 0.1 18 0 0 0 0 0 1 1 0
+19 1 0.1 18 0 0 0 0 0   0 1 1
-20 1 0.1 19 0 0 0 0 0 1 1 0
+20 1 0.1 19 0 0 0 0 0   0 1 1
-21 1 0.1 20 0 0 0 0 0 1 1 0
+21 1 0.1 20 0 0 0 0 0   0 1 1
-22 1 0.1 21 0 0 0 0 0 1 1 0
+22 1 0.1 21 0 0 0 0 0   0 1 1
-23 1 0.1 22 0 0 0 0 0 1 1 0
+23 1 0.1 22 0 0 0 0 0   0 1 1
-24 1 0.1 23 0 0 0 0 0 1 1 0
+24 1 0.1 23 0 0 0 0 0   0 1 1
-25 1 0.1 24 0 0 0 0 0 1 1 0
+25 1 0.1 24 0 0 0 0 0   0 1 1
-26 1 0.1 25 0 0 0 0 0 1 1 0
+26 1 0.1 25 0 0 0 0 0   0 1 1
-27 1 0.1 26 0 0 0 0 0 1 1 0
+27 1 0.1 26 0 0 0 0 0   0 1 1
-28 1 0.1 27 0 0 0 0 0 1 1 0
+28 1 0.1 27 0 0 0 0 0   0 1 1
-29 1 0.1 28 0 0 0 0 0 1 1 0
+29 1 0.1 28 0 0 0 0 0   0 1 1
-30 1 0.1 29 0 0 0 0 0 1 1 0
+30 1 0.1 29 0 0 0 0 0   0 1 1
-31 1 0.1 30 0 0 0 0 0 1 1 0
+31 1 0.1 30 0 0 0 0 0   0 1 1
-32 1 0.1 31 0 0 0 0 0 1 1 0
+32 1 0.1 31 0 0 0 0 0   0 1 1
-33 1 0.1 32 0 0 0 0 0 1 1 0
+33 1 0.1 32 0 0 0 0 0   0 1 1
-34 1 0.1 33 0 0 0 0 0 1 1 0
+34 1 0.1 33 0 0 0 0 0   0 1 1
-35 1 0.1 34 0 0 0 0 0 1 1 0
+35 1 0.1 34 0 0 0 0 0   0 1 1
-36 1 0.1 35 0 0 0 0 0 1 1 0
+36 1 0.1 35 0 0 0 0 0   0 1 1
-37 1 0.1 36 0 0 0 0 0 1 1 0
+37 1 0.1 36 0 0 0 0 0   0 1 1
-38 1 0.1 37 0 0 0 0 0 1 1 0
+38 1 0.1 37 0 0 0 0 0   0 1 1
-39 1 0.1 38 0 0 0 0 0 1 1 0
+39 1 0.1 38 0 0 0 0 0   0 1 1
-40 1 0.1 39 0 0 0 0 0 1 1 0
+40 1 0.1 39 0 0 0 0 0   0 1 1
-41 1 0.1 40 0 0 0 0 0 1 1 0
+41 1 0.1 40 0 0 0 0 0   0 1 1
-42 1 0.1 41 0 0 0 0 0 1 1 0
+42 1 0.1 41 0 0 0 0 0   0 1 1
-43 1 0.1 42 0 0 0 0 0 1 1 0
+43 1 0.1 42 0 0 0 0 0   0 1 1
-44 1 0.1 43 0 0 0 0 0 1 1 0
+44 1 0.1 43 0 0 0 0 0   0 1 1
-45 1 0.1 44 0 0 0 0 0 1 1 0
+45 1 0.1 44 0 0 0 0 0   0 1 1
-46 1 0.1 45 0 0 0 0 0 1 1 0
+46 1 0.1 45 0 0 0 0 0   0 1 1
-47 1 0.1 46 0 0 0 0 0 1 1 0
+47 1 0.1 46 0 0 0 0 0   0 1 1
-48 1 0.1 47 0 0 0 0 0 1 1 0
+48 1 0.1 47 0 0 0 0 0   0 1 1
-49 1 0.1 48 0 0 0 0 0 1 1 0
+49 1 0.1 48 0 0 0 0 0   0 1 1
-50 1 0.1 49 0 0 0 0 0 1 1 0
+50 1 0.1 49 0 0 0 0 0   0 1 1
-51 1 0.1 50 0 0 0 0 0 1 1 0
+51 1 0.1 50 0 0 0 0 0   0 1 1
-52 1 0.1 51 0 0 0 0 0 1 1 0
+52 1 0.1 51 0 0 0 0 0   0 1 1
-53 1 0.1 52 0 0 0 0 0 1 1 0
+53 1 0.1 52 0 0 0 0 0   0 1 1
-54 1 0.1 53 0 0 0 0 0 1 1 0
+54 1 0.1 53 0 0 0 0 0   0 1 1
-55 1 0.1 54 0 0 0 0 0 1 1 0
+55 1 0.1 54 0 0 0 0 0   0 1 1
-56 1 0.1 55 0 0 0 0 0 1 1 0
+56 1 0.1 55 0 0 0 0 0   0 1 1
-57 1 0.1 56 0 0 0 0 0 1 1 0
+57 1 0.1 56 0 0 0 0 0   0 1 1
-58 1 0.1 57 0 0 0 0 0 1 1 0
+58 1 0.1 57 0 0 0 0 0   0 1 1
-59 1 0.1 58 0 0 0 0 0 1 1 0
+59 1 0.1 58 0 0 0 0 0   0 1 1
-60 1 0.1 59 0 0 0 0 0 1 1 0
+60 1 0.1 59 0 0 0 0 0   1 1 1

examples/case05/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     1      0      0    0.1     0     0    90    all
+1   0.5     0     0      1      0      0     0     0     0    all

examples/case05/input/input.txt

@@ -5,15 +5,16 @@
 # ── 流程控制 ──────────────────────────────────
 # 每步用 0/1 单独开关，1=执行，0=跳过
 # 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧
-step_simulate:  1   # 运行物理模拟 → output/trajectory.txt
+step_simulate:  0   # 运行物理模拟 → output/trajectory.txt
-step_sample:    1   # 抽帧 → output/display.txt
+step_sample:    0   # 抽帧 → output/display.txt
 step_plot:      0   # 绘制轨迹/能量图 → output/trajectory_plots.png
-step_animation: 1   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
+step_animation: 0   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
+step_plot_wave: 1   # 绘制波形图
 force_calc: 0       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
 
 # ── 计算引擎 ──────────────────────────────────
 # 可选: python, c, cpp, fortran, java
-engine: python        # 默认使用 Python 引擎
+engine: python      # 默认使用 python 引擎
 
 # ── 盒子 ──────────────────────────────────────
 box_a: 80.0         # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
@@ -60,7 +61,7 @@ warmup_steps: 0           # 默认 0（立即开始记录）
 
 # 总模拟时间（秒），程序自动计算 NT = T_total / DT
 # 如果同时指定了 NT，以 NT 为准
-T_total: 100.0
+T_total: 10.0
 
 # 抽帧间隔（每 NSTEP 步取一帧用于动画）
 NSTEP: 50

examples/case06/Readme.md

@@ -0,0 +1,42 @@
+# case05: 一维原子链驱动力学模拟
+
+60 个原子沿 x 轴排列，相邻原子用弹簧连接。原子 1 受驱动力作用。
+
+## 物理设定
+
+| 参数 | 值 |
+|---|---|
+| 原子数 | 60 |
+| 排列 | 沿 x 轴等间距排列，间距为 1 |
+| 约束 | 原子**只沿 z 方向**振动（fix_x=1, fix_y=1, fix_z=0） |
+| 弹簧 | 劲度系数 k=1.0，原长 L₀=1.0 |
+| 重力 | 无 |
+| 万有引力 | 无 |
+| 阻尼 | 无 |
+| 驱动力 | 原子 1（详见 driver.txt） |
+| 算法 | leapfrog（蛙跳法，能量守恒） |
+
+## 驱动力
+
+原子 1 的位置由 `input/driver.txt` 中的驱动力公式决定：
+
+```math
+z(t) = A_z \cdot \cos(2\pi f_z t + \phi_z)
+```
+
+当前参数：A_z = 5.0, f_z = 1.0 Hz, φ_z = 90°（全程驱动）。
+
+受驱原子完全忽略 coord.txt 中的初始坐标和 fix 约束，位置/速度由驱动力解析确定。
+
+## 动力学行为
+
+原子 1 的受迫振动通过弹簧逐次传递给相邻原子，形成沿链传播的横波。由于横向振动的几何非线性（弹簧张力主要在 x 方向），z 方向有效刚度低，波速较慢，呈现 FPU 型非线性动力学特征。
+
+## 使用方法
+
+```bash
+cd examples/case05
+python run_dynamics.py
+```
+
+配置参数详见 `input/input.txt`，驱动力定义见 `input/driver.txt`，完整文档见 `doc/index.html`。

examples/case06/doc/index.html

@@ -0,0 +1,473 @@
+<!DOCTYPE html>
+<html lang="zh-CN">
+<head>
+<meta charset="UTF-8">
+<meta name="viewport" content="width=device-width, initial-scale=1.0">
+<title>case05 — 一维原子链驱动力学模拟 | 物理原理 &amp; 使用文档</title>
+<style>
+  :root {
+    --bg: #f8f9fa;
+    --card: #fff;
+    --text: #1a1a2e;
+    --accent: #2563eb;
+    --accent-light: #dbeafe;
+    --code-bg: #1e293b;
+    --code-text: #e2e8f0;
+    --border: #e2e8f0;
+    --muted: #64748b;
+  }
+  * { margin: 0; padding: 0; box-sizing: border-box; }
+  body {
+    font-family: -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Noto Sans SC", sans-serif;
+    background: var(--bg);
+    color: var(--text);
+    line-height: 1.7;
+  }
+
+  /* ── Header ── */
+  .hero {
+    background: linear-gradient(135deg, #1e293b 0%, #334155 100%);
+    color: #fff;
+    padding: 56px 24px 48px;
+    text-align: center;
+  }
+  .hero h1 { font-size: 2rem; font-weight: 700; letter-spacing: -0.02em; }
+  .hero .subtitle {
+    margin-top: 10px;
+    font-size: 1.05rem;
+    opacity: 0.8;
+  }
+  .hero .badge {
+    display: inline-block;
+    margin-top: 14px;
+    padding: 4px 14px;
+    border-radius: 999px;
+    background: rgba(255,255,255,0.12);
+    font-size: 0.82rem;
+  }
+
+  /* ── Layout ── */
+  .container { max-width: 820px; margin: 0 auto; padding: 32px 20px; }
+
+  section { margin-bottom: 44px; }
+  h2 {
+    font-size: 1.35rem;
+    font-weight: 600;
+    margin-bottom: 16px;
+    padding-bottom: 8px;
+    border-bottom: 2px solid var(--accent);
+    display: inline-block;
+  }
+  h3 {
+    font-size: 1.05rem;
+    font-weight: 600;
+    margin: 20px 0 10px;
+  }
+
+  p, li { margin-bottom: 10px; }
+  ul, ol { padding-left: 22px; }
+  strong { color: var(--accent); }
+
+  /* ── Cards ── */
+  .card {
+    background: var(--card);
+    border-radius: 12px;
+    padding: 20px 24px;
+    margin-bottom: 16px;
+    border: 1px solid var(--border);
+    box-shadow: 0 1px 3px rgba(0,0,0,0.04);
+  }
+
+  /* ── Formula / Code blocks ── */
+  .formula {
+    background: var(--card);
+    border-left: 4px solid var(--accent);
+    padding: 14px 20px;
+    margin: 14px 0;
+    font-family: "Times New Roman", "STIX", serif;
+    font-size: 1.05rem;
+    overflow-x: auto;
+    border-radius: 0 8px 8px 0;
+  }
+  code {
+    background: var(--accent-light);
+    padding: 2px 7px;
+    border-radius: 4px;
+    font-family: "JetBrains Mono", "Fira Code", monospace;
+    font-size: 0.88em;
+  }
+  pre {
+    background: var(--code-bg);
+    color: var(--code-text);
+    padding: 16px 20px;
+    border-radius: 10px;
+    overflow-x: auto;
+    font-size: 0.85rem;
+    line-height: 1.5;
+    margin: 14px 0;
+  }
+  pre .cm { color: #94a3b8; font-style: italic; }   /* comment */
+
+  /* ── Table ── */
+  table {
+    width: 100%;
+    border-collapse: collapse;
+    margin: 14px 0;
+    font-size: 0.92rem;
+  }
+  th, td {
+    padding: 8px 12px;
+    text-align: left;
+    border-bottom: 1px solid var(--border);
+  }
+  th { background: var(--accent-light); font-weight: 600; }
+
+  /* ── TOC ── */
+  .toc { counter-reset: toc; }
+  .toc li { counter-increment: toc; list-style: none; margin-bottom: 6px; }
+  .toc li::before { content: counter(toc) ". "; font-weight: 600; color: var(--accent); }
+  .toc a { color: var(--accent); text-decoration: none; }
+  .toc a:hover { text-decoration: underline; }
+
+  /* ── Flow diagram ── */
+  .flow { display: flex; flex-wrap: wrap; gap: 8px; align-items: center; justify-content: center; margin: 16px 0; }
+  .flow-step {
+    background: var(--accent-light);
+    border: 1px solid var(--accent);
+    border-radius: 8px;
+    padding: 8px 16px;
+    font-size: 0.88rem;
+    font-weight: 500;
+  }
+  .flow-arrow { color: var(--muted); font-size: 1.2rem; }
+
+  @media (max-width: 600px) {
+    .hero h1 { font-size: 1.5rem; }
+    .flow { flex-direction: column; }
+    .flow-arrow { transform: rotate(90deg); }
+  }
+</style>
+</head>
+<body>
+
+<!-- ============================================================ -->
+<!--  Header                                                       -->
+<!-- ============================================================ -->
+<header class="hero">
+  <h1>一维原子链驱动力学模拟</h1>
+  <p class="subtitle">60 个原子沿 x 轴排列 · 弹簧连接 · z 方向受迫振动</p>
+  <span class="badge">case05  ·  examples/case05</span>
+</header>
+
+<div class="container">
+
+<!-- ============================================================ -->
+<!--  TOC                                                          -->
+<!-- ============================================================ -->
+<section>
+  <h2>目录</h2>
+  <ol class="toc">
+    <li><a href="#physics">物理原理</a></li>
+    <li><a href="#algorithm">数值算法</a></li>
+    <li><a href="#driver">驱动力模型</a></li>
+    <li><a href="#usage">使用方法</a></li>
+    <li><a href="#params">参数参考</a></li>
+    <li><a href="#files">文件结构</a></li>
+    <li><a href="#troubleshoot">常见问题</a></li>
+  </ol>
+</section>
+
+<!-- ============================================================ -->
+<!--  1. Physics                                                   -->
+<!-- ============================================================ -->
+<section id="physics">
+  <h2>一、物理原理</h2>
+
+  <div class="card">
+    <h3>1.1 一维原子链</h3>
+    <p>60 个原子沿 <strong>x 轴</strong> 等间距排列，原子间距为 1。相邻原子之间用 <strong>理想弹簧</strong> 连接，弹簧的劲度系数 <em>k</em> = 1.0，原长 <em>L</em>₀ = 1.0（与原子间距一致，初始状态弹簧无拉伸）。</p>
+    <p>每个原子被限制在 <strong>z 方向</strong> 自由振动，x 和 y 方向锁定（<code>fix_x=1, fix_y=1, fix_z=0</code>）。</p>
+  </div>
+
+  <div class="card">
+    <h3>1.2 弹簧力（胡克定律）</h3>
+    <p>当原子 <em>i</em> 和 <em>j</em> 之间有弹簧连接时，原子 <em>i</em> 受到的弹簧力为：</p>
+    <div class="formula">
+      <strong>F</strong> = −<em>k</em> · (<em>d</em> − <em>L</em>₀) · <strong>u</strong><sub><em>ij</em></sub>
+    </div>
+    <p>其中 <em>d</em> = |<strong>r</strong><sub><em>j</em></sub> − <strong>r</strong><sub><em>i</em></sub>| 为两原子间距离，<strong>u</strong><sub><em>ij</em></sub> 为从 <em>i</em> 指向 <em>j</em> 的单位向量。由于原子只在 z 方向振动，弹簧在 z 方向的分量是 <strong>几何非线性</strong> 的——对于小振幅近似，z 方向等效于一个三次方恢复力（FPU 型非线性）。</p>
+  </div>
+
+  <div class="card">
+    <h3>1.3 运动方程</h3>
+    <p>对于第 <em>i</em> 个自由原子（非受驱），牛顿第二定律给出：</p>
+    <div class="formula">
+      <em>m</em> · <strong>a</strong><sub><em>i</em></sub> = <strong>F</strong><sub><em>i</em></sub><sup>spring</sup> + <strong>F</strong><sub><em>i</em></sub><sup>driving</sup>
+    </div>
+    <p>本案例中 <strong>唯一的外力</strong> 来自驱动力（仅施加于原子 1）。无重力、无万有引力、无阻尼，系统总能量守恒。</p>
+  </div>
+
+  <div class="card">
+    <h3>1.4 波传播</h3>
+    <p>原子 1 的受迫振动通过弹簧逐次传递给相邻原子，形成沿链传播的 <strong>横波</strong>。由于横向振动的几何非线性（弹簧大部分张力在 x 方向，z 方向的有效刚度远小于 1），波的传播速度较慢，且高阶频率成分会在链中产生复杂的非线性动力学行为（类似 FPU 回波现象）。</p>
+  </div>
+</section>
+
+<!-- ============================================================ -->
+<!--  2. Algorithm                                                 -->
+<!-- ============================================================ -->
+<section id="algorithm">
+  <h2>二、数值算法</h2>
+
+  <div class="card">
+    <h3>2.1 蛙跳法（Leapfrog / Velocity-Verlet）</h3>
+    <p>采用能量守恒特性优异的 <strong>蛙跳法</strong>（二阶辛积分器），更新公式为：</p>
+    <div class="formula">
+      <strong>v</strong>(<em>t</em> + ½Δ<em>t</em>) = <strong>v</strong>(<em>t</em>) + ½ <strong>a</strong>(<em>t</em>) · Δ<em>t</em><br>
+      <strong>r</strong>(<em>t</em> + Δ<em>t</em>) = <strong>r</strong>(<em>t</em>) + <strong>v</strong>(<em>t</em> + ½Δ<em>t</em>) · Δ<em>t</em><br>
+      <strong>a</strong>(<em>t</em> + Δ<em>t</em>) = <strong>F</strong>(<strong>r</strong>(<em>t</em> + Δ<em>t</em>), <strong>v</strong>(<em>t</em> + ½Δ<em>t</em>)) / <em>m</em><br>
+      <strong>v</strong>(<em>t</em> + Δ<em>t</em>) = <strong>v</strong>(<em>t</em> + ½Δ<em>t</em>) + ½ <strong>a</strong>(<em>t</em> + Δ<em>t</em>) · Δ<em>t</em>
+    </div>
+    <p>蛙跳法在长时间模拟中能量漂移极小（本案例验证 <strong>&lt; 0.004%</strong>），适合无阻尼的保守系统。</p>
+  </div>
+
+  <div class="card">
+    <h3>2.2 时间步长与采样</h3>
+    <table>
+      <tr><th>参数</th><th>值</th><th>说明</th></tr>
+      <tr><td>DT</td><td>0.01 s</td><td>积分步长（远小于 1/ω ≈ 0.16 s，满足稳定性条件）</td></tr>
+      <tr><td>T_total</td><td>100 s</td><td>总模拟时间 → NT = 10000 步</td></tr>
+      <tr><td>NSTEP</td><td>50</td><td>每 NSTEP 步取一帧用于动画 → 200 帧</td></tr>
+      <tr><td>method</td><td>leapfrog</td><td>蛙跳法（Velocity-Verlet）</td></tr>
+    </table>
+  </div>
+
+  <div class="card">
+    <h3>2.3 计算流程</h3>
+    <div class="flow">
+      <span class="flow-step">读入 coord.txt<br>connection.txt<br>bond.txt</span>
+      <span class="flow-arrow">→</span>
+      <span class="flow-step">施加驱动力<br>（驱动原子 1）</span>
+      <span class="flow-arrow">→</span>
+      <span class="flow-step">记录轨迹</span>
+      <span class="flow-arrow">→</span>
+      <span class="flow-step">蛙跳法<br>更新位置/速度</span>
+      <span class="flow-arrow">→</span>
+      <span class="flow-step">固定约束<br>（x, y 锁定）</span>
+      <span class="flow-arrow">→</span>
+      <span class="flow-step" style="background:#fef3c7;border-color:#f59e0b;">循环<br>NT 次</span>
+    </div>
+    <p style="margin-top:12px;">注意：驱动力在 <strong>每次积分前</strong> 施加，确保受驱原子的位置正确传递给弹簧力计算。</p>
+  </div>
+</section>
+
+<!-- ============================================================ -->
+<!--  3. Driving Force                                             -->
+<!-- ============================================================ -->
+<section id="driver">
+  <h2>三、驱动力模型</h2>
+
+  <div class="card">
+    <h3>3.1 定义文件</h3>
+    <p>驱动力由 <code>input/driver.txt</code> 定义，格式如下：</p>
+    <pre>n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period
+1     0     0     5      0      0      1     0     0    90    all</pre>
+  </div>
+
+  <div class="card">
+    <h3>3.2 数学公式</h3>
+    <p>受驱原子的位置由下式决定（<strong>完全替换</strong> coord.txt 中的初始坐标和固定约束）：</p>
+    <div class="formula">
+      <strong>r</strong>(<em>t</em>) = <strong>A</strong> · cos(2π<em>f</em> · <em>t</em> + <strong>φ</strong>)
+    </div>
+    <p>速度由解析导数给出：</p>
+    <div class="formula">
+      <strong>v</strong>(<em>t</em>) = −<strong>A</strong> · 2π<em>f</em> · sin(2π<em>f</em> · <em>t</em> + <strong>φ</strong>)
+    </div>
+    <p>其中 <strong>A</strong> = (amp_x, amp_y, amp_z)，<strong>f</strong> = (freq_x, freq_y, freq_z) 为不同方向的驱动频率，<strong>φ</strong> = (phi_x, phi_y, phi_z) 为相位（<strong>角度制</strong>，代码自动转换为弧度）。</p>
+  </div>
+
+  <div class="card">
+    <h3>3.3 本案例驱动参数</h3>
+    <table>
+      <tr><th>参数</th><th>值</th><th>含义</th></tr>
+      <tr><td>amp_z</td><td>5.0</td><td>z 方向驱动振幅</td></tr>
+      <tr><td>freq_z</td><td>1.0 Hz</td><td>驱动频率（周期 1 s）</td></tr>
+      <tr><td>phi_z</td><td>90°</td><td>驱动相位 → z(0) = 5·cos(90°) = 0</td></tr>
+      <tr><td>period</td><td>all</td><td>全程驱动，永不停止</td></tr>
+    </table>
+    <div class="formula">
+      <em>z</em>(<em>t</em>) = 5.0 · cos(2π · 1.0 · <em>t</em> + 90°)
+    </div>
+  </div>
+
+  <div class="card">
+    <h3>3.4 有限周期驱动</h3>
+    <p><code>period</code> 参数支持三种模式：</p>
+    <ul>
+      <li><strong>all</strong> — 全程驱动</li>
+      <li><strong>数值</strong> — 驱动指定周期数后 <strong>静止</strong>（冻结在最终位置，速度归零）。例如 <code>period: 1</code> 表示驱动 1 个完整周期后停止。</li>
+    </ul>
+  </div>
+
+  <div class="card">
+    <h3>3.5 驱动与固定约束的关系</h3>
+    <p>对于受驱原子（<code>driver.txt</code> 中 <code>n</code> 指定的原子），其在 <code>coord.txt</code> 中的初始坐标和 <code>fix_x/fix_y/fix_z</code> 约束被 <strong>完全忽略</strong>。原子的位置和速度完全由驱动力公式决定。</p>
+  </div>
+</section>
+
+<!-- ============================================================ -->
+<!--  4. Usage                                                     -->
+<!-- ============================================================ -->
+<section id="usage">
+  <h2>四、使用方法</h2>
+
+  <div class="card">
+    <h3>4.1 完整运行（模拟 + 动画）</h3>
+    <pre>cd examples/case05
+python run_dynamics.py</pre>
+    <p>这步会依次执行：物理模拟 → 抽帧 → 打开 VisPy 3D 动画窗口。</p>
+  </div>
+
+  <div class="card">
+    <h3>4.2 仅查看已有结果</h3>
+    <p>如果已经跑完模拟且生成了 <code>output/display.txt</code>，可以通过修改 <code>input.txt</code> 跳过计算，只开动画：</p>
+    <pre>step_simulate:  0   # 跳过模拟
+step_sample:    0   # 跳过抽帧
+step_animation: 1   # 播放动画</pre>
+    <p>然后运行：<code>python run_dynamics.py</code></p>
+  </div>
+
+  <div class="card">
+    <h3>4.3 手动 3D 动画</h3>
+    <p>也可以单独启动 VisPy 窗口：</p>
+    <pre>python ../../draw.py output/</pre>
+  </div>
+
+  <div class="card">
+    <h3>4.4 强制重新计算</h3>
+    <p>修改参数后需要重新运行模拟时，设置：</p>
+    <pre>force_calc: 1        # 忽略缓存，强制重新计算</pre>
+  </div>
+
+  <div class="card">
+    <h3>4.5 动画交互</h3>
+    <table>
+      <tr><th>操作</th><th>效果</th></tr>
+      <tr><td>鼠标拖动</td><td>旋转视角</td></tr>
+      <tr><td>滚轮</td><td>缩放</td></tr>
+      <tr><td>A / D 键</td><td>视角向左 / 向右平移</td></tr>
+      <tr><td>W / S 键</td><td>视角向上 / 向下平移</td></tr>
+      <tr><td>Q / E 键</td><td>降低 / 提高平移速度</td></tr>
+      <tr><td>左上角 <strong>reset</strong> 按钮</td><td>复位视角到初始位置</td></tr>
+      <tr><td>左上角 <strong>info</strong> 按钮</td><td>切换信息面板显示/隐藏</td></tr>
+      <tr><td>左上角 <strong>axes</strong> 按钮</td><td>切换坐标轴显示/隐藏</td></tr>
+    </table>
+  </div>
+</section>
+
+<!-- ============================================================ -->
+<!--  5. Parameters                                                -->
+<!-- ============================================================ -->
+<section id="params">
+  <h2>五、参数参考</h2>
+
+  <div class="card">
+    <h3>5.1 input.txt 关键参数</h3>
+    <table>
+      <tr><th>参数</th><th>默认值</th><th>说明</th></tr>
+      <tr><td>gravity_field</td><td>0</td><td>均匀重力场（已关闭）</td></tr>
+      <tr><td>gravity_interaction</td><td>0</td><td>原子间万有引力（已关闭）</td></tr>
+      <tr><td>elastic_force</td><td>1</td><td>弹簧键力（已开启）</td></tr>
+      <tr><td>damping_force</td><td>0</td><td>阻尼（已关闭）</td></tr>
+      <tr><td><strong>driving_force</strong></td><td><strong>1</strong></td><td>驱动力开关（1=开启，需 driver.txt）</td></tr>
+      <tr><td>method</td><td>leapfrog</td><td>数值积分方法</td></tr>
+      <tr><td>DT</td><td>0.01</td><td>积分步长 (s)</td></tr>
+      <tr><td>T_total</td><td>100.0</td><td>总模拟时间 (s)</td></tr>
+      <tr><td>NSTEP</td><td>50</td><td>抽帧步数间隔</td></tr>
+      <tr><td>engine</td><td>python</td><td>计算引擎（python / c / cpp / fortran）</td></tr>
+      <tr><td>use_marker</td><td>1</td><td>渲染模式（0=Sphere 网格, 1=Marker GPU 实例化）</td></tr>
+    </table>
+  </div>
+
+  <div class="card">
+    <h3>5.2 流程控制参数</h3>
+    <table>
+      <tr><th>参数</th><th>0</th><th>1</th></tr>
+      <tr><td>step_simulate</td><td>跳过模拟（加载已有轨迹）</td><td>运行物理模拟</td></tr>
+      <tr><td>step_sample</td><td>跳过抽帧</td><td>从轨迹抽取显示帧</td></tr>
+      <tr><td>step_plot</td><td>不生成图表</td><td>生成轨迹/能量图</td></tr>
+      <tr><td>step_animation</td><td>不启动动画</td><td>自动打开 VisPy 3D 窗口</td></tr>
+      <tr><td>force_calc</td><td>自动检测缓存</td><td>强制重新计算</td></tr>
+    </table>
+  </div>
+</section>
+
+<!-- ============================================================ -->
+<!--  6. File Structure                                            -->
+<!-- ============================================================ -->
+<section id="files">
+  <h2>六、文件结构</h2>
+
+  <pre>case05/
+├── input/
+│   ├── input.txt          # 主配置文件（YAML 格式）
+│   ├── coord.txt          # 原子坐标（60 个原子）
+│   ├── connection.txt     # 弹簧连接关系（59 条键）
+│   ├── bond.txt           # 弹簧参数（k=1.0, L₀=1.0）
+│   └── <strong>driver.txt</strong>       # <span class="cm">驱动力定义（本案例新增）</span>
+├── output/
+│   ├── trajectory.txt     # 全量轨迹数据（10000 步 × 60 原子）
+│   ├── display.txt        # 抽帧后的动画数据（200 帧 × 60 原子）
+│   ├── dynamics.log       # 计算日志
+│   └── animation.log      # 动画启动日志（闪退时排查用）
+├── doc/
+│   └── index.html         # <span class="cm">本文档</span>
+├── Readme.md              # 案例简介
+└── run_dynamics.py        # 案例运行入口</pre>
+</section>
+
+<!-- ============================================================ -->
+<!--  7. Troubleshooting                                           -->
+<!-- ============================================================ -->
+<section id="troubleshoot">
+  <h2>七、常见问题</h2>
+
+  <div class="card">
+    <h3>7.1 动画窗口闪退</h3>
+    <p>如果 VisPy 窗口一闪就消失，请检查：</p>
+    <ul>
+      <li><code>output/animation.log</code> 中是否有错误信息</li>
+      <li><code>output/display.txt</code> 是否存在（需先跑 <code>step_sample: 1</code>）</li>
+    </ul>
+  </div>
+
+  <div class="card">
+    <h3>7.2 原子不振动</h3>
+    <p>可能原因：</p>
+    <ul>
+      <li><strong>NSTEP 过大</strong>：抽帧间隔大于驱动周期的一半时，动画会丢失振动细节。建议 NSTEP ≤ 1/(freq · DT · 10)</li>
+      <li><strong>相位 φ 使采样点落在零值</strong>：试试 <code>phi_z: 0</code> 让原子在 t=0 处于振幅峰值</li>
+      <li>确认 <code>driving_force: 1</code> 且 <code>driver.txt</code> 中 amp_z 不为 0</li>
+    </ul>
+  </div>
+
+  <div class="card">
+    <h3>7.3 渲染性能慢</h3>
+    <p>原子数多时动画卡顿：</p>
+    <ul>
+      <li>设置 <code>use_marker: 1</code>（使用 GPU 实例化渲染替代独立网格球体）</li>
+      <li>增大 <code>NSTEP</code> 减少动画帧数</li>
+    </ul>
+  </div>
+</section>
+
+<hr style="border:none;border-top:1px solid var(--border);margin:40px 0;">
+
+<footer style="text-align:center;color:var(--muted);font-size:0.85rem;margin-bottom:40px;">
+  Dynamics Simulation Framework &nbsp;·&nbsp; 生成于 2026-06-10
+</footer>
+
+</div>
+</body>
+</html>

examples/case06/input/bond.txt

@@ -0,0 +1,2 @@
+bond_name k      rest_length
+k1        100.0  1.0

examples/case06/input/connection.txt

@@ -0,0 +1,60 @@
+n1 n2 bond_name
+1 2 k1
+2 3 k1
+3 4 k1
+4 5 k1
+5 6 k1
+6 7 k1
+7 8 k1
+8 9 k1
+9 10 k1
+10 11 k1
+11 12 k1
+12 13 k1
+13 14 k1
+14 15 k1
+15 16 k1
+16 17 k1
+17 18 k1
+18 19 k1
+19 20 k1
+20 21 k1
+21 22 k1
+22 23 k1
+23 24 k1
+24 25 k1
+25 26 k1
+26 27 k1
+27 28 k1
+28 29 k1
+29 30 k1
+30 31 k1
+31 32 k1
+32 33 k1
+33 34 k1
+34 35 k1
+35 36 k1
+36 37 k1
+37 38 k1
+38 39 k1
+39 40 k1
+40 41 k1
+41 42 k1
+42 43 k1
+43 44 k1
+44 45 k1
+45 46 k1
+46 47 k1
+47 48 k1
+48 49 k1
+49 50 k1
+50 51 k1
+51 52 k1
+52 53 k1
+53 54 k1
+54 55 k1
+55 56 k1
+56 57 k1
+57 58 k1
+58 59 k1
+59 60 k1

examples/case06/input/coord.txt

@@ -0,0 +1,61 @@
+n mass radius x y z vx vy vz fix_x fix_y fix_z
+1 1 0.1 0 0 1 0 0 0   0 1 1
+2 1 0.1 1 0 0 0 0 0   0 1 1
+3 1 0.1 2 0 0 0 0 0   0 1 1
+4 1 0.1 3 0 0 0 0 0   0 1 1
+5 1 0.1 4 0 0 0 0 0   0 1 1
+6 1 0.1 5 0 0 0 0 0   0 1 1
+7 1 0.1 6 0 0 0 0 0   0 1 1
+8 1 0.1 7 0 0 0 0 0   0 1 1
+9 1 0.1 8 0 0 0 0 0   0 1 1
+10 1 0.1 9 0 0 0 0 0   0 1 1
+11 1 0.1 10 0 0 0 0 0   0 1 1
+12 1 0.1 11 0 0 0 0 0   0 1 1
+13 1 0.1 12 0 0 0 0 0   0 1 1
+14 1 0.1 13 0 0 0 0 0   0 1 1
+15 1 0.1 14 0 0 0 0 0   0 1 1
+16 1 0.1 15 0 0 0 0 0   0 1 1
+17 1 0.1 16 0 0 0 0 0   0 1 1
+18 1 0.1 17 0 0 0 0 0   0 1 1
+19 1 0.1 18 0 0 0 0 0   0 1 1
+20 1 0.1 19 0 0 0 0 0   0 1 1
+21 1 0.1 20 0 0 0 0 0   0 1 1
+22 1 0.1 21 0 0 0 0 0   0 1 1
+23 1 0.1 22 0 0 0 0 0   0 1 1
+24 1 0.1 23 0 0 0 0 0   0 1 1
+25 1 0.1 24 0 0 0 0 0   0 1 1
+26 1 0.1 25 0 0 0 0 0   0 1 1
+27 1 0.1 26 0 0 0 0 0   0 1 1
+28 1 0.1 27 0 0 0 0 0   0 1 1
+29 1 0.1 28 0 0 0 0 0   0 1 1
+30 1 0.1 29 0 0 0 0 0   0 1 1
+31 1 0.1 30 0 0 0 0 0   0 1 1
+32 1 0.1 31 0 0 0 0 0   0 1 1
+33 1 0.1 32 0 0 0 0 0   0 1 1
+34 1 0.1 33 0 0 0 0 0   0 1 1
+35 1 0.1 34 0 0 0 0 0   0 1 1
+36 1 0.1 35 0 0 0 0 0   0 1 1
+37 1 0.1 36 0 0 0 0 0   0 1 1
+38 1 0.1 37 0 0 0 0 0   0 1 1
+39 1 0.1 38 0 0 0 0 0   0 1 1
+40 1 0.1 39 0 0 0 0 0   0 1 1
+41 1 0.1 40 0 0 0 0 0   0 1 1
+42 1 0.1 41 0 0 0 0 0   0 1 1
+43 1 0.1 42 0 0 0 0 0   0 1 1
+44 1 0.1 43 0 0 0 0 0   0 1 1
+45 1 0.1 44 0 0 0 0 0   0 1 1
+46 1 0.1 45 0 0 0 0 0   0 1 1
+47 1 0.1 46 0 0 0 0 0   0 1 1
+48 1 0.1 47 0 0 0 0 0   0 1 1
+49 1 0.1 48 0 0 0 0 0   0 1 1
+50 1 0.1 49 0 0 0 0 0   0 1 1
+51 1 0.1 50 0 0 0 0 0   0 1 1
+52 1 0.1 51 0 0 0 0 0   0 1 1
+53 1 0.1 52 0 0 0 0 0   0 1 1
+54 1 0.1 53 0 0 0 0 0   0 1 1
+55 1 0.1 54 0 0 0 0 0   0 1 1
+56 1 0.1 55 0 0 0 0 0   0 1 1
+57 1 0.1 56 0 0 0 0 0   0 1 1
+58 1 0.1 57 0 0 0 0 0   0 1 1
+59 1 0.1 58 0 0 0 0 0   0 1 1
+60 1 0.1 59 0 0 0 0 0   1 1 1

examples/case06/input/driver.txt

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

examples/case06/input/input.txt

@@ -0,0 +1,96 @@
+# 物理模拟参数配置
+# 格式：YAML
+# 用法：python run_dynamics.py
+
+# ── 流程控制 ──────────────────────────────────
+# 每步用 0/1 单独开关，1=执行，0=跳过
+# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧
+step_simulate:  1   # 运行物理模拟 → output/trajectory.txt
+step_sample:    1   # 抽帧 → output/display.txt
+step_plot:      0   # 绘制轨迹/能量图 → output/trajectory_plots.png
+step_animation: 1   # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
+force_calc: 0       # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
+
+# ── 计算引擎 ──────────────────────────────────
+# 可选: python, c, cpp, fortran, java
+engine: python      # 默认使用 python 引擎
+
+# ── 盒子 ──────────────────────────────────────
+box_a: 80.0         # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
+
+# ── 初始构型 ──────────────────────────────────
+# 坐标文件格式：
+#   第一行：n mass radius x y z vx vy vz fix_x fix_y fix_z
+#   后续行：原子序号 质量 半径 x y z vx vy vz fix_x fix_y fix_z
+coord_file: input/coord.txt
+connection_file: input/connection.txt
+bond_file: input/bond.txt
+driver_file: input/driver.txt     # 驱动力定义文件（driving_force=1 时生效）
+
+# 绘图/动画展示的原子序号（对应 coord_file 第一列 n）
+plot_atom: 1
+
+# ── 物理参数 ──────────────────────────────────
+# 三个方向分量分别对应 x, y, z
+G: [0.0, 0.0, 0.0]   # 重力场分量 (m/s²)
+B: [0.0, 0.0, 0.0]   # 阻尼分量
+
+# ── 力开关（0=关闭, 1=开启）──────────────────
+gravity_field: 0         # 均匀重力场 (G)
+gravity_interaction: 0   # 原子间万有引力
+elastic_force: 1         # 弹簧键力
+damping_force: 0         # 阻尼 (B)
+driving_force: 1         # 驱动力（需 driver_file 定义）
+#
+gravity_strength: 1.0    # 万有引力强度（仅 gravity_interaction=1 时有效）
+
+# ── 数值算法 ──────────────────────────────────
+# 可选：
+#   explicit_euler  显式欧拉法
+#   implicit_euler  隐式欧拉法
+#   midpoint        中点法
+#   leapfrog        蛙跳法
+method: leapfrog
+
+# ── 步骤控制 ──────────────────────────────────
+# 以下参数控制哪些步骤被执行和保存
+
+# 预热步数：模拟开始时跳过不保存的步数（用于稳定初始状态）
+warmup_steps: 0           # 默认 0（立即开始记录）
+
+# 总模拟时间（秒），程序自动计算 NT = T_total / DT
+# 如果同时指定了 NT，以 NT 为准
+T_total: 50.0
+
+# 抽帧间隔（每 NSTEP 步取一帧用于动画）
+NSTEP: 50
+
+# ── 时间步长 ──────────────────────────────────
+DT:    0.001         # 时间步长 (s)
+
+# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧
+sample_start: null        # null 表示从头开始（帧索引从 0 起）
+sample_end: null          # null 表示到末尾
+
+
+
+# ── 渲染方式 ──────────────────────────────────
+# 3D 动画中原子渲染方式：
+#   0 = Sphere （网格球体，效果精细，原子数少时推荐）
+#   1 = Marker （GPU 实例化点，原子数多时性能更佳）
+use_marker: 1
+
+# ── 显示参数 ──────────────────────────────────
+# 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序
+alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5]
+
+# 小球颜色
+# 小球半径从 coord_file 的 radius 列读取
+ball_color_r: 0.20  # R 分量 (0~1)
+ball_color_g: 0.60  # G 分量
+ball_color_b: 0.90  # B 分量
+
+# 盒子面颜色
+box_color_r:  0.80
+box_color_g:  0.80
+box_color_b:  0.85

examples/case06/run_dynamics.py

@@ -0,0 +1,54 @@
+"""
+Case runner for Dynamics case05 — 1D atomic chain.
+
+This script keeps program and data separated:
+  - program: ../../dynamics.py
+  - input:   ./input
+  - output:  ./output
+"""
+
+from __future__ import annotations
+
+import argparse
+import importlib.util
+from pathlib import Path
+
+
+CASE_DIR = Path(__file__).resolve().parent
+DYNAMICS_PATH = Path("..") / ".." / "dynamics.py"
+INPUT_DIR = Path("input")
+OUTPUT_DIR = Path("output")
+CONFIG_FILE = INPUT_DIR / "input.txt"
+
+
+def load_dynamics_module(module_path: Path):
+    spec = importlib.util.spec_from_file_location("dynamics_module", module_path)
+    if spec is None or spec.loader is None:
+        raise ImportError(f"无法加载 dynamics.py: {module_path}")
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+def main():
+    parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01")
+    parser.add_argument("--no-plot", action="store_true", help="跳过 matplotlib 绘图")
+    args = parser.parse_args()
+
+    dynamics_path = (CASE_DIR / DYNAMICS_PATH).resolve()
+    input_dir = (CASE_DIR / INPUT_DIR).resolve()
+    output_dir = (CASE_DIR / OUTPUT_DIR).resolve()
+    config_path = (CASE_DIR / CONFIG_FILE).resolve()
+
+    module = load_dynamics_module(dynamics_path)
+    module.run_case(
+        config_path=config_path,
+        runtime_base=CASE_DIR,
+        input_dir=input_dir,
+        output_dir=output_dir,
+        no_plot=args.no_plot,
+    )
+
+
+if __name__ == "__main__":
+    main()

plot_wave.py

@@ -0,0 +1,260 @@
+"""
+plot_wave.py
+============
+波形与能量动态图：读取 display.txt，绘制原子位移波形
+（纵波 + 2 个横波）和系统能量/输入功率随时间变化的二维动画。
+
+用法：
+    python plot_wave.py                          # 使用 dynamics 根目录下 output/
+    python plot_wave.py examples/case05/output    # 指定案例输出目录
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+import os
+import sys
+
+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)
+
+
+def compute_energy(x, y, z, vx, vy, vz, masses, mass_arr,
+                   bond_pairs, bond_stiffness, bond_rest_lengths,
+                   gravity_field, G, gravity_interaction, gravity_strength):
+    """计算系统各能量分量。
+
+    Returns:
+        ek_sys:  系统动能 (n_frames,)
+        us_sys:  系统弹性势能 (n_frames,)
+        ug_sys:  系统重力势能 (n_frames,)
+        ugr_sys: 系统万有引力势能 (n_frames,)
+    """
+    n_frames = x.shape[0]
+    masses_2d = masses[np.newaxis, :]  # (1, n_atoms)
+
+    # 动能 Ek = ½ m v²
+    ek = 0.5 * masses_2d * (vx**2 + vy**2 + vz**2)
+    ek_sys = np.sum(ek, axis=1)
+
+    # 弹性势能 Us = ½ k (d - d₀)²
+    us_sys = np.zeros(n_frames)
+    if 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 = x[:, j] - x[:, i]
+            dy = y[:, j] - y[:, i]
+            dz = z[:, j] - z[:, i]
+            dist = np.sqrt(dx**2 + dy**2 + dz**2)
+            stretch = dist - bond_rest_lengths[b_idx]
+            us_sys += 0.5 * bond_stiffness[b_idx] * stretch**2
+
+    # 均匀重力场势能 Ug = -m G·r
+    ug_sys = np.zeros(n_frames)
+    if gravity_field:
+        G_vec = np.array(G)
+        ug_sys = -masses_2d * (G_vec[0] * x + G_vec[1] * y + G_vec[2] * z)
+        ug_sys = np.sum(ug_sys, axis=1)
+
+    # 万有引力势能 Ug_grav = -G_grav Σ m_i m_j / r
+    ugr_sys = np.zeros(n_frames)
+    if gravity_interaction:
+        n_atoms = len(masses)
+        # 为避免巨大计算量，仅当原子数较少时计算
+        if n_atoms <= 200:
+            for i in range(n_atoms):
+                for j in range(i + 1, n_atoms):
+                    dx = x[:, j] - x[:, i]
+                    dy = y[:, j] - y[:, i]
+                    dz = z[:, j] - 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
+                    ugr_sys += pair_pe
+
+    return ek_sys, us_sys, ug_sys, ugr_sys
+
+
+def plot_wave(output_dir):
+    """主绘图函数：读取 display.txt 并生成波形+能量动画。"""
+    data = load_disp_data(output_dir)
+
+    n_frames = int(data["n_frames"])
+    t = np.array(data["disp_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"])
+
+    # 原子信息
+    pos_0    = np.array(data["atom_positions"])  # (n_atoms, 3)
+    masses   = np.array(data["atom_masses"])
+    atom_ids = np.array(data["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", []))
+
+    # 物理开关
+    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))
+
+    # ── 位移（偏离初始平衡位形）──
+    dx = x - pos_0[np.newaxis, :, 0]   # 纵波（沿链方向 x）
+    dy = y - pos_0[np.newaxis, :, 1]   # 横波 1（y 方向）
+    dz = z - pos_0[np.newaxis, :, 2]   # 横波 2（z 方向）
+
+    # ── 能量 ──
+    ek, us, ug, ugr = compute_energy(
+        x, y, z, vx, vy, vz, masses, masses,
+        bond_pairs, bond_stiffness, bond_rest_lengths,
+        gravity_field, G, gravity_interaction, gravity_strength)
+    e_total = ek + us + ug + ugr
+    power = np.gradient(e_total, t)  # 输入功率 = dE/dt
+
+    # ── 波形纵轴范围（全局统一，避免抖动）──
+    def get_ylim(disp_data):
+        vmax = np.max(np.abs(disp_data))
+        if vmax < 1e-10:
+            return -1.0, 1.0
+        margin = vmax * 0.2
+        return -vmax - margin, vmax + margin
+
+    ylims = [get_ylim(dx), get_ylim(dy), get_ylim(dz)]
+    e_max = max(np.max(e_total), 0.01) * 1.3
+    p_abs = np.max(np.abs(power))
+    p_max = max(p_abs * 1.3, 0.01)
+
+    atom_idx = np.arange(n_atoms)
+
+    # ── 图形设置 ──
+    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei', 'DejaVu Sans']
+    plt.rcParams['axes.unicode_minus'] = False
+
+    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
+    fig.suptitle("波形与能量分析", fontsize=16)
+
+    # ── 3 个波形子图 ──
+    titles = [
+        "纵波 (x 方向位移)",
+        "横波 (y 方向位移)",
+        "横波 (z 方向位移)",
+    ]
+    colors = ["#2563eb", "#ea580c", "#16a34a"]
+    wave_configs = list(zip(
+        [axes[0, 0], axes[0, 1], axes[1, 0]],
+        [dx, dy, dz], titles, colors, ylims
+    ))
+
+    wave_lines = []
+    time_texts = []
+    for ax, disp, title, color, yl in wave_configs:
+        ax.set_xlim(0, n_atoms - 1)
+        ax.set_ylim(yl)
+        ax.set_xlabel("原子序号")
+        ax.set_ylabel("位移")
+        ax.set_title(title)
+        ax.grid(True, alpha=0.3)
+        line, = ax.plot([], [], color=color, linewidth=1.5)
+        wave_lines.append(line)
+        tt = ax.text(0.02, 0.95, "", transform=ax.transAxes,
+                     fontsize=11, verticalalignment="top")
+        time_texts.append(tt)
+
+    # ── 能量+功率子图 ──
+    ax_ep = axes[1, 1]
+    ax_ep.set_xlim(t[0], t[-1])
+    ep_ylow = min(-0.1 * max(e_max, p_max), -p_max * 0.1)
+    ep_yhigh = max(e_max, p_max)
+    ax_ep.set_ylim(ep_ylow, ep_yhigh)
+    ax_ep.set_xlabel("时间 (s)")
+    ax_ep.set_ylabel("能量 / 功率")
+    ax_ep.set_title("系统能量与输入功率")
+    ax_ep.grid(True, alpha=0.3)
+
+    line_ek, = ax_ep.plot([], [], "b-",  lw=1.5, label="动能")
+    line_us, = ax_ep.plot([], [], "orange", lw=1.5, label="弹性势能")
+    line_et, = ax_ep.plot([], [], "r--", lw=1.5, label="总能量")
+    line_pw, = ax_ep.plot([], [], "g-",  lw=1.5, alpha=0.7, label="输入功率 (dE/dt)")
+    if gravity_field:
+        line_ug, = ax_ep.plot([], [], "purple", lw=1.0, alpha=0.5, label="重力势能")
+    else:
+        line_ug = None
+    if gravity_interaction and n_atoms <= 200:
+        line_ugr, = ax_ep.plot([], [], "brown", lw=1.0, alpha=0.5, label="万有引力势能")
+    else:
+        line_ugr = None
+
+    ax_ep.legend(loc="upper right", fontsize=9)
+
+    plt.tight_layout(rect=[0, 0, 1, 0.95])
+
+    # ── 动画更新 ──
+    def update(frame):
+        # 波形
+        for i in range(3):
+            wave_lines[i].set_data(atom_idx, [dx, dy, dz][i][frame])
+            time_texts[i].set_text(f"t = {t[frame]:.2f} s  |  帧 {frame+1}/{n_frames}")
+
+        # 能量（累计到当前帧）
+        cur_t = t[:frame + 1]
+        line_ek.set_data(cur_t, ek[:frame + 1])
+        line_us.set_data(cur_t, us[:frame + 1])
+        line_et.set_data(cur_t, e_total[:frame + 1])
+        line_pw.set_data(cur_t, power[:frame + 1])
+        if line_ug:
+            line_ug.set_data(cur_t, ug[:frame + 1])
+        if line_ugr:
+            line_ugr.set_data(cur_t, ugr[:frame + 1])
+
+        # 能量图 x 轴动态扩展
+        ax_ep.set_xlim(t[0], max(t[frame] + max(t[-1] * 0.05, 1), t[-1]))
+
+        artists = wave_lines + time_texts + [line_ek, line_us, line_et, line_pw]
+        if line_ug:  artists.append(line_ug)
+        if line_ugr: artists.append(line_ugr)
+        return artists
+
+    ani = FuncAnimation(fig, update, frames=n_frames, interval=50, blit=True)
+
+    # 保存 GIF
+    gif_path = os.path.join(output_dir, "wave_animation.gif")
+    ani.save(gif_path, writer="pillow", fps=min(20, max(1, n_frames // 5)))
+    print(f"[plot_wave] GIF 已保存: {gif_path}")
+
+    # 尝试保存 MP4
+    try:
+        mp4_path = os.path.join(output_dir, "wave_animation.mp4")
+        ani.save(mp4_path, writer="ffmpeg", fps=min(20, max(1, n_frames // 5)))
+        print(f"[plot_wave] MP4 已保存: {mp4_path}")
+    except Exception:
+        pass
+
+    plt.show()
+    return gif_path
+
+
+if __name__ == "__main__":
+    script_dir = os.path.dirname(os.path.abspath(__file__))
+    if len(sys.argv) > 1:
+        output_dir = os.path.abspath(sys.argv[1])
+    else:
+        output_dir = compute.get_output_dir(script_dir)
+    plot_wave(output_dir)