diff --git a/engines/c/main.c b/engines/c/main.c
index 1f6d21f..441fffa 100644
--- a/engines/c/main.c
+++ b/engines/c/main.c
@@ -507,6 +507,12 @@ static void limit_in_box(double *pos, double *vel, double lo, double hi) {
     if (*pos < lo)  { *pos = lo;  *vel = -*vel; }
 }
 
+/* 周期边界回绕：超出边界的位置绕到对侧（与 Python wrap_position 一致）*/
+static void wrap_position(double *pos, double lo, double hi) {
+    if (*pos > hi)  *pos = lo;
+    if (*pos < lo)  *pos = hi;
+}
+
 /* ── 显式欧拉法 ──────────── */
 static void explicit_euler_step(
     int n, double *x, double *y, double *z,
@@ -752,6 +758,13 @@ static void apply_step(
         limit_in_box(&z[i],  &vz[i], -box_a, box_a);
     }
 
+    /* 周期边界回绕（与 Python wrap_position 一致）*/
+    for (int i = 0; i < n; i++) {
+        wrap_position(&x[i], -box_a, box_a);
+        wrap_position(&y[i], -box_a, box_a);
+        wrap_position(&z[i], -box_a, box_a);
+    }
+
     /* 逐自由度固定约束（与 Python apply_fixed_constraints 一致） */
     for (int i = 0; i < n; i++) {
         if (fixed[i*3+0]) { x[i] = pos_0[i*3+0]; vx[i] = 0.0; }
@@ -899,6 +912,9 @@ int main(int argc, char **argv) {
     traj.vz = traj.vy + record_steps * n;
 
     /* 预热 */
+    /* 初始时刻 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);
+
     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);
diff --git a/engines/cpp/main.cpp b/engines/cpp/main.cpp
index 53b8b97..e00e4f6 100644
--- a/engines/cpp/main.cpp
+++ b/engines/cpp/main.cpp
@@ -386,6 +386,12 @@ static void limit_in_box(double &pos, double &vel, double lo, double hi) {
     if (pos < lo)  { pos = lo;  vel = -vel; }
 }
 
+/* 周期边界回绕（与 Python wrap_position 一致）*/
+static void wrap_position(double &pos, double lo, double hi) {
+    if (pos > hi) pos = lo;
+    if (pos < lo) pos = hi;
+}
+
 // ========================================================================
 // 四种积分方法（只做位置/速度更新，不含边界条件）
 // 与 Python: Explicit_Euler_Method / Implicit_Euler_Method /
@@ -598,6 +604,13 @@ static void apply_step(
         limit_in_box(z[i],  vz[i], -box_a, box_a);
     }
 
+    // 周期边界回绕（与 Python wrap_position 一致）
+    for (int i = 0; i < n; i++) {
+        wrap_position(x[i], -box_a, box_a);
+        wrap_position(y[i], -box_a, box_a);
+        wrap_position(z[i], -box_a, box_a);
+    }
+
     // 逐自由度固定约束（与 Python apply_fixed_constraints 一致）
     for (int i = 0; i < n; i++) {
         if (fixed[i*3+0]) { x[i] = pos_0[i*3];   vx[i] = 0.0; }
@@ -811,6 +824,10 @@ int main(int argc, char **argv) {
     std::vector<double> traj_vz(record_steps * n);
 
     // 预热
+    // 初始时刻 t=0 驱动力（与 Python run_simulation 一致）
+    if (params.driving_force)
+        apply_driving_force(n, x.data(), y.data(), z.data(), vx.data(), vy.data(), vz.data(), 0.0, 0, params.DT, drivers);
+
     for (int s = 0; s < params.warmup_steps; s++) {
         double tw = (s + 1) * params.DT;
         if (params.driving_force)
diff --git a/engines/fortran/main.f90 b/engines/fortran/main.f90
index 60362e5..0dcaef4 100644
--- a/engines/fortran/main.f90
+++ b/engines/fortran/main.f90
@@ -108,6 +108,17 @@ program dynamics_f90
     allocate(traj_vx(record_steps, n), traj_vy(record_steps, n), traj_vz(record_steps, n))
 
     ! 预热
+    ! 初始时刻 t=0 驱动力（与 Python run_simulation 一致）
+    if (driving_force /= 0 .and. n_drivers > 0) then
+        call apply_driving(n, x, y, z, vx, vy, vz, 0.0d0, 0, 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
+
     do s = 1, warmup_steps
         if (driving_force /= 0 .and. n_drivers > 0) then
             tw = (s * 1.0d0) * DT
@@ -512,6 +523,14 @@ subroutine limit_in_box(pos, vel, lo, hi)
     end if
 end subroutine limit_in_box
 
+! 周期边界回绕（与 Python wrap_position 一致）
+subroutine wrap_position(pos, lo, hi)
+    double precision, intent(inout) :: pos
+    double precision, intent(in) :: lo, hi
+    if (pos > hi) pos = lo
+    if (pos < lo) pos = hi
+end subroutine wrap_position
+
 ! ── 显式欧拉法 ────────────
 subroutine explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, &
                                nb, bp, bk, br, fixed, dt, &
@@ -730,6 +749,13 @@ subroutine apply_step(method, n, x, y, z, vx, vy, vz, m, G, B, &
         call limit_in_box(z(i),  vz(i), -box_a, box_a)
     end do
 
+    ! 周期边界回绕（与 Python wrap_position 一致）
+    do i = 1, n
+        call wrap_position(x(i), -box_a, box_a)
+        call wrap_position(y(i), -box_a, box_a)
+        call wrap_position(z(i), -box_a, box_a)
+    end do
+
     ! 逐自由度固定约束（与 Python apply_fixed_constraints 一致）
     do i = 1, n
         if (fixed(i,1) /= 0) then; x(i) = pos_0(i, 1); vx(i) = 0.0d0; end if