From 685234c84f4bdfc91ffdf3ebbc2bbfdb1144c28f Mon Sep 17 00:00:00 2001 From: Ying-Li Niu <64801511@qq.com> Date: Thu, 11 Jun 2026 09:19:34 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E4=B8=BA=20C/C++/Fortran=20=E5=BC=95?= =?UTF-8?q?=E6=93=8E=E5=A2=9E=E5=8A=A0=E9=A9=B1=E5=8A=A8=E5=8A=9B(driving?= =?UTF-8?q?=5Fforce)=E6=94=AF=E6=8C=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - param.json 新增 driving_force 开关 - C 引擎: 新增 DriverData 结构体、read_driver()、apply_driving_force() - C++ 引擎: 同上（C++ 风格实现） - Fortran 引擎: 同上（Fortran 90 风格实现） - 修复 JSON 输出末尾逗号导致加载失败的问题 - 编译通过并验证 C 引擎运行正常（100000步/6.6s） --- compute.py | 1 + engines/c/main.c | 209 +++++++++++++++++++++++++++++++++- engines/cpp/main.cpp | 144 +++++++++++++++++++++++- engines/fortran/main.f90 | 236 +++++++++++++++++++++++++++++++++++++-- 4 files changed, 579 insertions(+), 11 deletions(-) diff --git a/compute.py b/compute.py index 81f9cc1..2ba9e23 100644 --- a/compute.py +++ b/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) diff --git a/engines/c/main.c b/engines/c/main.c index 7c2bc19..1f6d21f 100644 --- a/engines/c/main.c +++ b/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, + &_x, &_y, &_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", - atoms.n_atoms, bonds.n_bonds, params.NT, params.DT, params.method); + DriverData drivers; + 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]; diff --git a/engines/cpp/main.cpp b/engines/cpp/main.cpp index 65eeaf3..53b8b97 100644 --- a/engines/cpp/main.cpp +++ b/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 rest_lengths; }; +// ======================================================================== +// 驱动力数据 +// ======================================================================== +struct DriverData { + int n_drivers = 0; + std::vector atom_idx; // internal atom indices + std::vector amp_x, amp_y, amp_z; + std::vector freq_x, freq_y, freq_z; + std::vector phi_x, phi_y, phi_z; // radians + std::vector has_period; // 0=all, 1=limited + std::vector period_cycles; + std::vector 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(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(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]; diff --git a/engines/fortran/main.f90 b/engines/fortran/main.f90 index 6edfdb9..60362e5 100644 --- a/engines/fortran/main.f90 +++ b/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