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
2026-06-12 18:36:37 +08:00
parent e1ade53fff
commit e62e536cee
12 changed files with 627 additions and 355 deletions
@@ -533,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; }
@@ -650,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,
@@ -659,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 一致）──────────────── */
@@ -895,12 +901,13 @@ static void write_display_txt(const char *path, const Trajectory *traj,
    FILE *f = fopen(path, "w");
    if (!f) die("无法写入 display.txt");

-    int n_frames    = traj->n_steps;
+    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;

-    fprintf(f, "number of frames: %d\n", n_frames);
+    /* 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);
@@ -1015,6 +1022,20 @@ int main(int argc, char **argv) {
        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 一致）*/
    if (params.driving_force) apply_driving_force(n, x, y, z, vx, vy, vz, 0.0, 0, params.DT, &drivers);
@@ -420,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; }
@@ -543,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,
@@ -552,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 一致）── */
@@ -737,15 +744,17 @@ static void write_display_txt(
    if (!f) die("无法写入 " + path);
    std::cout << "[Cpp-engine] 正在写入显示数据…" << std::endl;

-    double T_total = params.NT * params.DT;
+    int dynamic_steps = params.NT - params.warmup_steps;
+    double T_total = dynamic_steps * params.DT;

-    f << "number of frames: " << n_steps << "\n";
+    /* 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: " << (params.NT - 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";
@@ -925,6 +934,21 @@ int main(int argc, char **argv) {
    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 一致）
    if (params.driving_force)
@@ -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
@@ -524,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
@@ -652,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, &
@@ -662,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

 ! ── 分发器 + 边界条件 + 自由度约束 ──