feat: 为 C/C++/Fortran 引擎增加驱动力(driving_force)支持

- param.json 新增 driving_force 开关
- C 引擎: 新增 DriverData 结构体、read_driver()、apply_driving_force()
- C++ 引擎: 同上（C++ 风格实现）
- Fortran 引擎: 同上（Fortran 90 风格实现）
- 修复 JSON 输出末尾逗号导致加载失败的问题
- 编译通过并验证 C 引擎运行正常（100000步/6.6s）

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)

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