From e371fa8db162432f51f7c9e9bd9c9f4ff72f0079 Mon Sep 17 00:00:00 2001 From: Ying-Li Niu <64801511@qq.com> Date: Sun, 21 Jun 2026 05:33:52 +0800 Subject: [PATCH] =?UTF-8?q?refactor:=20=E7=A7=BB=E9=99=A4=E6=97=A7?= =?UTF-8?q?=E5=AD=90=E8=BF=9B=E7=A8=8B=E6=A8=A1=E5=BC=8F=E7=9A=84=20main.*?= =?UTF-8?q?=20=E6=BA=90=E6=96=87=E4=BB=B6=EF=BC=8CMakefile=20=E5=8F=AA?= =?UTF-8?q?=E7=BC=96=E8=AF=91=20DLL?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 删除 engines/c/main.c, engines/cpp/main.cpp, engines/fortran/main.f90 - Makefile 移除 all/exe 构建目标，默认只编译 DLL --- engines/c/Makefile | 51 +- engines/c/main.c | 1114 ------------------------------------- engines/cpp/Makefile | 10 +- engines/cpp/main.cpp | 1025 ---------------------------------- engines/fortran/Makefile | 12 +- engines/fortran/main.f90 | 1134 -------------------------------------- 6 files changed, 10 insertions(+), 3336 deletions(-) delete mode 100644 engines/c/main.c delete mode 100644 engines/cpp/main.cpp delete mode 100644 engines/fortran/main.f90 diff --git a/engines/c/Makefile b/engines/c/Makefile index 03e99ce..495f8d9 100644 --- a/engines/c/Makefile +++ b/engines/c/Makefile @@ -1,21 +1,17 @@ # engines/c/Makefile -# 跨平台编译：make → 本地系统编译 -# make linux → Linux 交叉编译（需 x86_64-linux-gnu-gcc） +# 编译 DLL（主程序通过 ctypes 直接调用） +# make dll → 本地系统编译 +# make linux → Linux 交叉编译（需 x86_64-linux-gnu-gcc） # make windows → Windows 交叉编译（需 x86_64-w64-mingw32-gcc） -# make macos → macOS 交叉编译（需 osxcross 工具链） CC = gcc CFLAGS = -O3 -march=native -Wall -Wextra LDFLAGS = -lm -SRCS = main.c LIB_SRC = dynamics_lib.c # 自动检测系统 UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) -# 目标文件名：统一使用 .exe 后缀（方便 Python 跨平台调用） -TARGET = build/dynamics_c.exe - # DLL 目标（平台自动选择后缀） ifeq ($(UNAME_S),Linux) DLL_TARGET = build/dynamics_c.so @@ -28,17 +24,12 @@ else DLL_FLAGS = -shared endif -# ── 本地编译 ───────────────────────────────── -.PHONY: all dll clean linux windows macos +.PHONY: all dll clean -all: $(TARGET) +all: dll dll: $(DLL_TARGET) -$(TARGET): $(SRCS) | build - $(CC) $(CFLAGS) -o $@ $(SRCS) $(LDFLAGS) - @echo " === C engine built: $@ ===" - $(DLL_TARGET): $(LIB_SRC) | build $(CC) $(CFLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) $(LDFLAGS) @echo " === C DLL built: $@ ===" @@ -46,37 +37,5 @@ $(DLL_TARGET): $(LIB_SRC) | build build: mkdir -p build -# ── 交叉编译 ───────────────────────────────── -# Linux → Linux (x86_64) -linux: CROSS_PREFIX = x86_64-linux-gnu- -linux: CC = $(CROSS_PREFIX)gcc -linux: CFLAGS = -O3 -march=x86-64 -Wall -Wextra -linux: $(SRCS) | build - $(CC) $(CFLAGS) -o build/dynamics_c_linux.exe $(SRCS) $(LDFLAGS) - @echo " === Linux binary: build/dynamics_c_linux.exe ===" - -# 任意平台 → Windows (x86_64) -# 需要安装 MinGW 交叉编译器： -# apt install mingw-w64 (Debian/Ubuntu) -# brew install mingw-w64 (macOS) -windows: CROSS_PREFIX = x86_64-w64-mingw32- -windows: CC = $(CROSS_PREFIX)gcc -windows: CFLAGS = -O3 -march=x86-64 -Wall -Wextra -windows: $(SRCS) | build - $(CC) $(CFLAGS) -o build/dynamics_c_win.exe $(SRCS) $(LDFLAGS) - @echo " === Windows binary: build/dynamics_c_win.exe ===" - -# 任意平台 → macOS (x86_64) -# 需要安装 osxcross 工具链 -macos: CROSS_PREFIX = x86_64-apple-darwin- -macos: CC = $(CROSS_PREFIX)gcc -macos: CFLAGS = -O3 -march=x86-64 -Wall -Wextra -macos: $(SRCS) | build - $(CC) $(CFLAGS) -o build/dynamics_c_mac.exe $(SRCS) $(LDFLAGS) - @echo " === macOS binary: build/dynamics_c_mac.exe ===" - -# ── 编译所有平台 ────────────────────────────── -all-platforms: linux windows macos - clean: rm -rf build *.o diff --git a/engines/c/main.c b/engines/c/main.c deleted file mode 100644 index 661c915..0000000 --- a/engines/c/main.c +++ /dev/null @@ -1,1114 +0,0 @@ -/** - * engines/c/main.c - * ----------------- - * C 语言动力学模拟引擎。 - * 与 Python 版 (compute.py) 算法保持一致。 - * - * 输入: param.json 数值参数（Python 从 YAML 转换得来） - * /coord.txt - * /connection.txt - * /bond.txt - * 输出: /trajectory.txt (JSON 格式，与 Python 版兼容) - * - * 编译: cmake --build build --target dynamics_c - * 用法: ./build/dynamics_c - */ - -#include -#include -#include -#include -#include - -/* ======================================================================== - * 配置参数（从 param.json 读取） - * ======================================================================== */ -typedef struct { - double box_a; /* 盒子半边长 */ - int NT; /* 总步数 */ - double DT; /* 时间步长 */ - int NSTEP; /* 抽帧间隔 */ - int warmup_steps; /* 预热步数 */ - char method[32]; /* 算法名称 */ - double G[3]; /* 重力分量 */ - double B[3]; /* 阻尼分量 */ - int gravity_field; /* 均匀重力场开关 */ - int gravity_interaction; /* 原子间万有引力开关 */ - int elastic_force; /* 弹簧键力开关 */ - int damping_force; /* 阻尼开关 */ - double gravity_strength; /* 万有引力强度 */ - int driving_force; /* 驱动力开关 */ - int save_trajectory; /* 是否保存完整轨迹文件 */ - double alpha[6]; /* 盒子透明度 */ - double ball_radius; - double ball_color[3]; - double box_color[3]; - int use_marker; - double camera_distance, camera_elevation, camera_azimuth; -} SimParams; - -/* ======================================================================== - * 原子数据 - * ======================================================================== */ -typedef struct { - int n_atoms; - int *atom_ids; - double *masses; - double *radii; - double *pos_0; /* 初始位置 (n_atoms*3) */ - double *vel_0; /* 初始速度 (n_atoms*3) */ - int *fixed; /* 固定约束 (n_atoms*3) */ -} AtomData; - -/* ======================================================================== - * 成键数据 - * ======================================================================== */ -typedef struct { - int n_bonds; - int *pairs; /* (n_bonds*2) */ - double *stiffness; - 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 *eq_x, *eq_y, *eq_z; /* 平衡位置（初始坐标） */ - 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.eq_x = (double*)xmalloc(n_lines * sizeof(double)); - d.eq_y = (double*)xmalloc(n_lines * sizeof(double)); - d.eq_z = (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/eq 数组 */ - for (int i = 0; i < n_lines; i++) { - d.eq_x[i] = d.eq_y[i] = d.eq_z[i] = 0.0; - 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.eq_x[idx] = atoms->pos_0[ii*3+0]; - d.eq_y[idx] = atoms->pos_0[ii*3+1]; - d.eq_z[idx] = atoms->pos_0[ii*3+2]; - 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; -} - -/* ======================================================================== - * 轨迹缓冲区 - * ======================================================================== */ -typedef struct { - int n_steps; - int n_atoms; - double *x, *y, *z; - double *vx, *vy, *vz; -} Trajectory; - -/* ======================================================================== - * 辅助函数 - * ======================================================================== */ - -static void die(const char *msg) { - fprintf(stderr, "[C-engine] 错误: %s\n", msg); - exit(1); -} - -static void *xmalloc(size_t sz) { - void *p = malloc(sz); - if (!p) die("内存分配失败"); - return p; -} - -/* 从 JSON 中读取 double 值 */ -static double json_read_double(const char *json, const char *key) { - char search[256]; - snprintf(search, sizeof(search), "\"%s\"", key); - const char *p = strstr(json, search); - if (!p) return 0.0; - p = strchr(p, ':'); - if (!p) return 0.0; - p++; - while (*p == ' ' || *p == '\t' || *p == '\n') p++; - return strtod(p, NULL); -} - -static int json_read_int(const char *json, const char *key) { - return (int)json_read_double(json, key); -} - -/* 从 JSON 中读取字符串值（写入 dst，最多 dst_sz 字节） */ -static void json_read_string(const char *json, const char *key, char *dst, int dst_sz) { - char search[256]; - snprintf(search, sizeof(search), "\"%s\"", key); - const char *p = strstr(json, search); - if (!p) { dst[0] = '\0'; return; } - p = strchr(p, ':'); - if (!p) { dst[0] = '\0'; return; } - p++; - while (*p == ' ' || *p == '\t' || *p == '\n') p++; - if (*p != '"') { dst[0] = '\0'; return; } - p++; - int i = 0; - while (*p && *p != '"' && i < dst_sz - 1) { dst[i++] = *p++; } - dst[i] = '\0'; -} - -/* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */ -static void json_read_double3(const char *json, const char *key, double out[3]) { - char search[256]; - snprintf(search, sizeof(search), "\"%s\"", key); - const char *p = strstr(json, search); - if (!p) { out[0]=out[1]=out[2]=0; return; } - p = strchr(p, '['); - if (!p) { out[0]=out[1]=out[2]=0; return; } - p++; - for (int i = 0; i < 3; i++) { - while (*p == ' ' || *p == '\t' || *p == '\n' || *p == ',') p++; - out[i] = strtod(p, (char**)&p); - } -} - -static void json_read_double6(const char *json, const char *key, double out[6]) { - char search[256]; - snprintf(search, sizeof(search), "\"%s\"", key); - const char *p = strstr(json, search); - if (!p) { for (int i=0;i<6;i++) out[i]=0; return; } - p = strchr(p, '['); - if (!p) { for (int i=0;i<6;i++) out[i]=0; return; } - p++; - for (int i = 0; i < 6; i++) { - while (*p == ' ' || *p == '\t' || *p == '\n' || *p == ',' || *p == ']') p++; - out[i] = strtod(p, (char**)&p); - } -} - -/* 读取 param.json */ -static int g_gravity_field = 1; -static int g_gravity_interaction = 0; -static int g_elastic_force = 1; -static int g_damping_force = 0; -static double g_gravity_strength = 1.0; - -static SimParams read_params(const char *path) { - FILE *f = fopen(path, "rb"); - if (!f) die("无法打开 param.json"); - fseek(f, 0, SEEK_END); - long sz = ftell(f); - fseek(f, 0, SEEK_SET); - char *buf = (char*)xmalloc(sz + 1); - fread(buf, 1, sz, f); - buf[sz] = '\0'; - fclose(f); - - SimParams p; - p.box_a = json_read_double(buf, "box_a"); - p.NT = json_read_int(buf, "NT"); - p.DT = json_read_double(buf, "DT"); - p.NSTEP = json_read_int(buf, "NSTEP"); - p.warmup_steps = json_read_int(buf, "warmup_steps"); - strcpy(p.method, "leapfrog"); /* 默认 */ - json_read_string(buf, "method", p.method, sizeof(p.method)); - json_read_double3(buf, "G", p.G); - json_read_double3(buf, "B", p.B); - p.gravity_field = json_read_int(buf, "gravity_field"); - p.gravity_interaction = json_read_int(buf, "gravity_interaction"); - 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"); - p.save_trajectory = json_read_int(buf, "save_trajectory"); - /* 渲染参数 */ - json_read_double6(buf, "alpha", p.alpha); - p.ball_radius = json_read_double(buf, "ball_radius"); - json_read_double3(buf, "ball_color", p.ball_color); - json_read_double3(buf, "box_color", p.box_color); - p.use_marker = json_read_int(buf, "use_marker"); - p.camera_distance = json_read_double(buf, "camera_distance"); - p.camera_elevation = json_read_double(buf, "camera_elevation"); - p.camera_azimuth = json_read_double(buf, "camera_azimuth"); - g_gravity_field = p.gravity_field; - g_gravity_interaction = p.gravity_interaction; - g_elastic_force = p.elastic_force; - g_damping_force = p.damping_force; - g_gravity_strength = p.gravity_strength; - - free(buf); - return p; -} - -/* 读取 coord.txt */ -static AtomData read_coord(const char *input_dir) { - char path[512]; - snprintf(path, sizeof(path), "%s/coord.txt", input_dir); - FILE *f = fopen(path, "r"); - if (!f) die("无法打开 coord.txt"); - - /* 跳过第一行表头 */ - char line[1024]; - if (!fgets(line, sizeof(line), f)) die("coord.txt 为空"); - - int capacity = 16; - AtomData a; - a.n_atoms = 0; - a.atom_ids = (int*)xmalloc(capacity * sizeof(int)); - a.masses = (double*)xmalloc(capacity * sizeof(double)); - a.radii = (double*)xmalloc(capacity * sizeof(double)); - a.pos_0 = (double*)xmalloc(capacity * 3 * sizeof(double)); - a.vel_0 = (double*)xmalloc(capacity * 3 * sizeof(double)); - a.fixed = (int*)xmalloc(capacity * 3 * sizeof(int)); - - while (fgets(line, sizeof(line), f)) { - if (a.n_atoms >= capacity) { - capacity *= 2; - a.atom_ids = realloc(a.atom_ids, capacity * sizeof(int)); - a.masses = realloc(a.masses, capacity * sizeof(double)); - a.radii = realloc(a.radii, capacity * sizeof(double)); - a.pos_0 = realloc(a.pos_0, capacity * 3 * sizeof(double)); - a.vel_0 = realloc(a.vel_0, capacity * 3 * sizeof(double)); - a.fixed = realloc(a.fixed, capacity * 3 * sizeof(int)); - } - int id, fx, fy, fz; - double mass, rad, px, py, pz, vx, vy, vz; - int n_parsed = sscanf(line, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %d", - &id, &mass, &rad, &px, &py, &pz, &vx, &vy, &vz, &fx, &fy, &fz); - if (n_parsed == 9) { - fx = fy = fz = 0; - } else if (n_parsed != 12) { - continue; - } - int i = a.n_atoms; - a.atom_ids[i] = id; - a.masses[i] = mass; - a.radii[i] = rad; - a.pos_0[i*3+0] = px; a.pos_0[i*3+1] = py; a.pos_0[i*3+2] = pz; - a.vel_0[i*3+0] = vx; a.vel_0[i*3+1] = vy; a.vel_0[i*3+2] = vz; - a.fixed[i*3+0] = fx; a.fixed[i*3+1] = fy; a.fixed[i*3+2] = fz; - a.n_atoms++; - } - fclose(f); - - if (a.n_atoms <= 0) die("coord.txt 原子数无效"); - return a; -} - -/* 读取 connection.txt */ -static BondData read_bonds(const char *input_dir, const AtomData *atoms) { - char path[512]; - BondData b; - b.n_bonds = 0; - b.pairs = NULL; - b.stiffness = NULL; - b.rest_lengths = NULL; - - snprintf(path, sizeof(path), "%s/connection.txt", input_dir); - FILE *f = fopen(path, "r"); - if (!f) return b; - - char line[256]; - if (!fgets(line, sizeof(line), f)) { fclose(f); return b; } - - int n_lines = 0, tmp_a, tmp_b; - char bond_name[256]; - while (fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name) == 3) n_lines++; - rewind(f); - fgets(line, sizeof(line), f); // 再次跳过表头 - - if (n_lines == 0) { fclose(f); return b; } - - b.n_bonds = n_lines; - b.pairs = (int*)xmalloc(n_lines * 2 * sizeof(int)); - b.stiffness = (double*)xmalloc(n_lines * sizeof(double)); - b.rest_lengths = (double*)xmalloc(n_lines * sizeof(double)); - - char bond_path[512]; - snprintf(bond_path, sizeof(bond_path), "%s/bond.txt", input_dir); - FILE *fb = fopen(bond_path, "r"); - - for (int i = 0; i < n_lines; i++) { - fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name); - b.pairs[i*2+0] = tmp_a - 1; - b.pairs[i*2+1] = tmp_b - 1; - b.stiffness[i] = 1.0; - b.rest_lengths[i] = 2.0; - if (fb) { - char name[256], header[256]; - double k, r0; - rewind(fb); - fgets(header, sizeof(header), fb); // 跳过表头行 - while (fscanf(fb, "%s %lf %lf", name, &k, &r0) == 3) { - if (strcmp(name, bond_name) == 0) { - b.stiffness[i] = k; - b.rest_lengths[i] = r0; - break; - } - } - } - } - fclose(f); - if (fb) fclose(fb); - return b; -} - -/* ======================================================================== - * 物理核心（与 Python compute.py 对应） - * ======================================================================== */ - -/* 加速度计算（各力独立开关控制） */ -static void compute_acceleration( - int n, const double *x, const double *y, const double *z, - const double *vx, const double *vy, const double *vz, - const double *m, const double G[3], const double B[3], - const BondData *bonds, - double *ax, double *ay, double *az) -{ - /* 先清零 */ - for (int i = 0; i < n; i++) { - ax[i] = 0.0; ay[i] = 0.0; az[i] = 0.0; - } - - /* 均匀重力场 */ - if (g_gravity_field) { - for (int i = 0; i < n; i++) { - ax[i] += G[0]; - ay[i] += G[1]; - az[i] += G[2]; - } - } - - /* 阻尼 */ - if (g_damping_force) { - for (int i = 0; i < n; i++) { - ax[i] -= B[0] * vx[i] / m[i]; - ay[i] -= B[1] * vy[i] / m[i]; - az[i] -= B[2] * vz[i] / m[i]; - } - } - - /* 弹簧键力 */ - if (g_elastic_force) { - for (int b = 0; b < bonds->n_bonds; b++) { - int i = bonds->pairs[b*2+0]; - int j = bonds->pairs[b*2+1]; - double dx = x[j] - x[i]; - double dy = y[j] - y[i]; - double dz = z[j] - z[i]; - double dist = sqrt(dx*dx + dy*dy + dz*dz); - if (dist < 1e-12) continue; - double stretch = dist - bonds->rest_lengths[b]; - double fmag = bonds->stiffness[b] * stretch; - double ux = dx / dist, uy = dy / dist, uz = dz / dist; - double fx = fmag * ux, fy = fmag * uy, fz = fmag * uz; - ax[i] += fx / m[i]; ay[i] += fy / m[i]; az[i] += fz / m[i]; - ax[j] -= fx / m[j]; ay[j] -= fy / m[j]; az[j] -= fz / m[j]; - } - } - - /* 万有引力（所有原子对之间） */ - if (g_gravity_interaction) { - for (int i = 0; i < n; i++) { - for (int j = i + 1; j < n; j++) { - double dx = x[j] - x[i]; - double dy = y[j] - y[i]; - double dz = z[j] - z[i]; - double r2 = dx*dx + dy*dy + dz*dz; - if (r2 <= 1e-12) continue; - double r = sqrt(r2); - double f_mag = g_gravity_strength * m[i] * m[j] / r2; - double fx_g = f_mag * dx / r; - double fy_g = f_mag * dy / r; - double fz_g = f_mag * dz / r; - ax[i] += fx_g / m[i]; ay[i] += fy_g / m[i]; az[i] += fz_g / m[i]; - ax[j] -= fx_g / m[j]; ay[j] -= fy_g / m[j]; az[j] -= fz_g / m[j]; - } - } - } -} - -/* 保守力加速度（不含阻尼），供真蛙跳法专用。 - 通过传入零速度调用 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; } - 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, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData *bonds, const int *fixed, double dt) -{ - 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); - - 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; - y[i] += vy[i] * dt; - z[i] += vz[i] * dt; - vx[i] += ax[i] * dt; - vy[i] += ay[i] * dt; - vz[i] += az[i] * dt; - } -} - -/* ── 隐式欧拉法 ──────────── */ -static void implicit_euler_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData *bonds, const int *fixed, double dt) -{ - double *vxn = (double*)alloca(n * sizeof(double)); - double *vyn = (double*)alloca(n * sizeof(double)); - double *vzn = (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]) { - vxn[i] = 0; vyn[i] = 0; vzn[i] = 0; continue; - } - double gx = B[0] / m[i], gy = B[1] / m[i], gz = B[2] / m[i]; - vxn[i] = (vx[i] + G[0] * dt) / (1.0 + gx * dt); - vyn[i] = (vy[i] + G[1] * dt) / (1.0 + gy * dt); - vzn[i] = (vz[i] + G[2] * dt) / (1.0 + gz * dt); - } - - 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, vxn, vyn, vzn, m, G, B, bonds, ax, ay, az); - - 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; - 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; - } -} - -/* ── 中点法 ──────────── */ -static void midpoint_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData *bonds, const int *fixed, double dt) -{ - 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); - - double *xm = (double*)alloca(n * sizeof(double)); - double *ym = (double*)alloca(n * sizeof(double)); - double *zm = (double*)alloca(n * sizeof(double)); - double *vxm = (double*)alloca(n * sizeof(double)); - double *vym = (double*)alloca(n * sizeof(double)); - double *vzm = (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]) { - xm[i]=ym[i]=zm[i]=vxm[i]=vym[i]=vzm[i]=0; continue; - } - xm[i] = x[i] + 0.5 * vx[i] * dt; - ym[i] = y[i] + 0.5 * vy[i] * dt; - zm[i] = z[i] + 0.5 * vz[i] * dt; - vxm[i] = vx[i] + 0.5 * ax[i] * dt; - vym[i] = vy[i] + 0.5 * ay[i] * dt; - vzm[i] = vz[i] + 0.5 * az[i] * dt; - x[i] = x[i] + vxm[i] * dt; - y[i] = y[i] + vym[i] * dt; - z[i] = z[i] + vzm[i] * dt; - } - - compute_acceleration(n, xm, ym, zm, vxm, vym, vzm, m, G, B, bonds, ax, ay, az); - - 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; - vy[i] += ay[i] * dt; - vz[i] += az[i] * dt; - } -} - -/* ── 蛙跳法（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, - const double *m, const double G[3], const double B[3], - const BondData *bonds, const int *fixed, double dt) -{ - double *ax = (double*)alloca(n * sizeof(double)); - double *ay = (double*)alloca(n * sizeof(double)); - double *az = (double*)alloca(n * sizeof(double)); - - /* 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; - 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; - } -} - -/* ── 驱动力（与 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->eq_x[d] + drivers->amp_x[d] * cos(2*M_PI*drivers->freq_x[d]*t + drivers->phi_x[d]); - double py = drivers->eq_y[d] + drivers->amp_y[d] * cos(2*M_PI*drivers->freq_y[d]*t + drivers->phi_y[d]); - double pz = drivers->eq_z[d] + 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, - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData *bonds, const int *fixed, - const double *pos_0, - double box_a, double dt) -{ - if (strcmp(method, "explicit_euler") == 0) { - explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); - } else if (strcmp(method, "implicit_euler") == 0) { - implicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); - } else if (strcmp(method, "midpoint") == 0) { - midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); - } else if (strcmp(method, "leapfrog") == 0) { - leapfrog_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); - } else { - fprintf(stderr, "[C-engine] 未知算法: %s\n", method); - exit(1); - } - - /* 边界条件（与 Python Limit_in_box 一致） */ - for (int i = 0; i < n; i++) { - if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; - limit_in_box(&x[i], &vx[i], -box_a, box_a); - limit_in_box(&y[i], &vy[i], -box_a, box_a); - 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; } - if (fixed[i*3+1]) { y[i] = pos_0[i*3+1]; vy[i] = 0.0; } - if (fixed[i*3+2]) { z[i] = pos_0[i*3+2]; vz[i] = 0.0; } - } -} - -// ======================================================================== -// JSON 输出 -// ======================================================================== - -static void write_trajectory_json(const char *path, const Trajectory *traj, - const SimParams *params, const AtomData *atoms, - const BondData *bonds) -{ - FILE *f = fopen(path, "w"); - if (!f) die("无法写入 trajectory.txt"); - - fprintf(f, "{\n"); - - const char *names[] = {"traj_x","traj_y","traj_z","traj_vx","traj_vy","traj_vz"}; - double *arrs[] = {traj->x, traj->y, traj->z, traj->vx, traj->vy, traj->vz}; - - printf("[C-engine] 正在写入轨迹数据…\n"); - fflush(stdout); - for (int a = 0; a < 6; a++) { - fprintf(f, " \"%s\": [\n", names[a]); - for (int t = 0; t < traj->n_steps; t++) { - fprintf(f, " ["); - for (int i = 0; i < traj->n_atoms; i++) { - fprintf(f, "%.8g", arrs[a][t * traj->n_atoms + i]); - if (i < traj->n_atoms - 1) fputc(',', f); - } - fprintf(f, "]"); - if (t < traj->n_steps - 1) fputc(',', f); - fputc('\n', f); - } - fprintf(f, " ]"); - fputc(',', f); - fputc('\n', f); - } - - /* 标量参数 */ - fprintf(f, " \"NT\": %d,\n", params->NT); - fprintf(f, " \"DT\": %.8g,\n", params->DT); - fprintf(f, " \"NSTEP\": %d,\n", params->NSTEP); - fprintf(f, " \"method\": \"%s\",\n", params->method); - fprintf(f, " \"warmup_steps\": %d,\n", params->warmup_steps); - fprintf(f, " \"G\": [%.8g, %.8g, %.8g],\n", params->G[0], params->G[1], params->G[2]); - fprintf(f, " \"B\": [%.8g, %.8g, %.8g],\n", params->B[0], params->B[1], params->B[2]); - - fprintf(f, " \"atom_ids\": ["); - for (int i = 0; i < atoms->n_atoms; i++) { - if (i > 0) fputc(',', f); - fprintf(f, "%d", atoms->atom_ids[i]); - } - fprintf(f, "],\n"); - - fprintf(f, " \"atom_masses\": ["); - for (int i = 0; i < atoms->n_atoms; i++) { - if (i > 0) fputc(',', f); - fprintf(f, "%.8g", atoms->masses[i]); - } - fprintf(f, "],\n"); - - fprintf(f, " \"bond_pairs\": ["); - for (int b = 0; b < bonds->n_bonds; b++) { - if (b > 0) fputc(',', f); - fprintf(f, "[%d, %d]", bonds->pairs[b*2], bonds->pairs[b*2+1]); - } - fprintf(f, "],\n"); - - fprintf(f, " \"bond_stiffness\": ["); - for (int b = 0; b < bonds->n_bonds; b++) { - if (b > 0) fputc(',', f); - fprintf(f, "%.8g", bonds->stiffness[b]); - } - fprintf(f, "],\n"); - - fprintf(f, " \"bond_rest_lengths\": ["); - for (int b = 0; b < bonds->n_bonds; b++) { - if (b > 0) fputc(',', f); - fprintf(f, "%.8g", bonds->rest_lengths[b]); - } - fprintf(f, "],\n"); - fprintf(f, " \"driving_force\": %d\n", params->driving_force); - - fprintf(f, "}\n"); - fclose(f); -} - -static void write_display_txt(const char *path, const Trajectory *traj, - const SimParams *params, const AtomData *atoms) -{ - FILE *f = fopen(path, "w"); - if (!f) die("无法写入 display.txt"); - - 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; - - /* 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); - fprintf(f, "method: %s\n", params->method); - fprintf(f, "warmup_steps: %d\n", params->warmup_steps); - fprintf(f, "dynamic_steps: %d\n", dynamic_steps); - fprintf(f, "T_total: %.16g\n", T_total); - fprintf(f, "box_a: %.16g\n", params->box_a); - fprintf(f, "alpha: %.16g,%.16g,%.16g,%.16g,%.16g,%.16g\n", - params->alpha[0], params->alpha[1], params->alpha[2], - params->alpha[3], params->alpha[4], params->alpha[5]); - fprintf(f, "ball_radius: %.16g\n", params->ball_radius); - fprintf(f, "ball_color_r: %.16g\n", params->ball_color[0]); - fprintf(f, "ball_color_g: %.16g\n", params->ball_color[1]); - fprintf(f, "ball_color_b: %.16g\n", params->ball_color[2]); - fprintf(f, "box_color_r: %.16g\n", params->box_color[0]); - fprintf(f, "box_color_g: %.16g\n", params->box_color[1]); - fprintf(f, "box_color_b: %.16g\n", params->box_color[2]); - fprintf(f, "use_marker: %d\n", params->use_marker); - fprintf(f, "camera_distance: %.16g\n", params->camera_distance); - fprintf(f, "camera_elevation: %.16g\n", params->camera_elevation); - fprintf(f, "camera_azimuth: %.16g\n", params->camera_azimuth); - fprintf(f, "\n"); - - if (params->driving_force) { - fprintf(f, "driving_force: 1\n"); - } else { - fprintf(f, "driving_force: 0\n"); - } - fprintf(f, "\n"); - - for (int t = 0; t < n_frames; t++) { - fprintf(f, "frame: %3d\n", t + 1); - fprintf(f, "n x y z vx vy vz\n"); - for (int i = 0; i < n_particles; i++) { - int idx = t * n_particles + i; - fprintf(f, "%4d %12.6f %12.6f %12.6f %10.6f %10.6f %10.6f\n", - atoms->atom_ids[i], - traj->x[idx], traj->y[idx], traj->z[idx], - traj->vx[idx], traj->vy[idx], traj->vz[idx]); - } - } - - fclose(f); -} - -// ======================================================================== -// 主函数 -// ======================================================================== - -int main(int argc, char **argv) { - if (argc < 4) { - fprintf(stderr, "用法: %s \n", argv[0]); - return 1; - } - const char *input_dir = argv[1]; - const char *output_dir = argv[2]; - const char *param_path = argv[3]; - - clock_t t0 = clock(); - - SimParams params = read_params(param_path); - AtomData atoms = read_coord(input_dir); - BondData bonds = read_bonds(input_dir, &atoms); - - 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)); - double *y = (double*)xmalloc(n * sizeof(double)); - double *z = (double*)xmalloc(n * sizeof(double)); - double *vx = (double*)xmalloc(n * sizeof(double)); - double *vy = (double*)xmalloc(n * sizeof(double)); - double *vz = (double*)xmalloc(n * sizeof(double)); - for (int i = 0; i < n; i++) { - x[i] = atoms.pos_0[i*3+0]; - y[i] = atoms.pos_0[i*3+1]; - z[i] = atoms.pos_0[i*3+2]; - vx[i] = atoms.vel_0[i*3+0]; - vy[i] = atoms.vel_0[i*3+1]; - vz[i] = atoms.vel_0[i*3+2]; - } - - /* 分配轨迹缓冲区 */ - int record_steps = params.NT - params.warmup_steps; - Trajectory traj; - traj.n_atoms = n; - if (params.save_trajectory) { - traj.n_steps = record_steps; - traj.x = (double*)xmalloc(record_steps * n * sizeof(double) * 6); - traj.y = traj.x + record_steps * n; - traj.z = traj.y + record_steps * n; - traj.vx = traj.z + record_steps * n; - traj.vy = traj.vx + record_steps * n; - traj.vz = traj.vy + record_steps * n; - } else { - int sampled_steps = (record_steps + params.NSTEP - 1) / params.NSTEP; - if (sampled_steps < 1) sampled_steps = 1; - traj.n_steps = sampled_steps; - traj.x = (double*)xmalloc(sampled_steps * n * sizeof(double) * 6); - traj.y = traj.x + sampled_steps * n; - traj.z = traj.y + sampled_steps * n; - traj.vx = traj.z + sampled_steps * n; - traj.vy = traj.vx + sampled_steps * n; - 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); - - 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, - params.box_a, params.DT); - } - - /* 记录 */ - int _prog_interval = record_steps / 100; - if (_prog_interval < 1) _prog_interval = 1; - int sample_idx = 0; - for (int s = 0; s < record_steps; s++) { - if (s % _prog_interval == 0 && s > 0) { - printf("[C-engine] progress: %d/%d\n", s, record_steps); - fflush(stdout); - } - 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); - int do_record = params.save_trajectory || (s % params.NSTEP == 0); - if (do_record) { - int idx = params.save_trajectory ? s : sample_idx; - for (int i = 0; i < n; i++) { - traj.x[ idx * n + i] = x[i]; - traj.y[ idx * n + i] = y[i]; - traj.z[ idx * n + i] = z[i]; - traj.vx[idx * n + i] = vx[i]; - traj.vy[idx * n + i] = vy[i]; - traj.vz[idx * n + i] = vz[i]; - } - if (!params.save_trajectory) sample_idx++; - } - apply_step(params.method, n, x, y, z, vx, vy, vz, - atoms.masses, params.G, params.B, &bonds, atoms.fixed, - atoms.pos_0, - params.box_a, params.DT); - } - - char out_path[512]; - if (params.save_trajectory) { - snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir); - write_trajectory_json(out_path, &traj, ¶ms, &atoms, &bonds); - } else { - snprintf(out_path, sizeof(out_path), "%s/display.txt", output_dir); - write_display_txt(out_path, &traj, ¶ms, &atoms); - } - - clock_t t1 = clock(); - double elapsed = (double)(t1 - t0) / CLOCKS_PER_SEC; - printf("[C-engine] 计算完成: %d 步, %.3f s\n", record_steps, elapsed); - - /* 清理 */ - free(traj.x); - free(x); free(y); free(z); - free(vx); free(vy); free(vz); - free(atoms.atom_ids); - free(atoms.masses); free(atoms.radii); - free(atoms.pos_0); free(atoms.vel_0); - free(atoms.fixed); - if (bonds.pairs) { free(bonds.pairs); free(bonds.stiffness); free(bonds.rest_lengths); } - - return 0; -} diff --git a/engines/cpp/Makefile b/engines/cpp/Makefile index 93cefc7..fa307d8 100644 --- a/engines/cpp/Makefile +++ b/engines/cpp/Makefile @@ -1,7 +1,7 @@ # engines/cpp/Makefile +# 编译 DLL（主程序通过 ctypes 直接调用） CXX = g++ -SRCS = main.cpp LIB_SRC = dynamics_lib.cpp UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) @@ -15,8 +15,6 @@ else STATIC_FLAGS = endif -TARGET = build/dynamics_cpp.exe - ifeq ($(UNAME_S),Linux) DLL_TARGET = build/dynamics_cpp.so DLL_FLAGS = -shared -fPIC @@ -30,14 +28,10 @@ endif .PHONY: all dll clean -all: $(TARGET) +all: dll dll: $(DLL_TARGET) -$(TARGET): $(SRCS) | build - $(CXX) $(CXXFLAGS) $(STATIC_FLAGS) -o $@ $(SRCS) - @echo " === C++ engine built: $@ ===" - $(DLL_TARGET): $(LIB_SRC) | build $(CXX) $(CXXFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) @echo " === C++ DLL built: $@ ===" diff --git a/engines/cpp/main.cpp b/engines/cpp/main.cpp deleted file mode 100644 index 338677f..0000000 --- a/engines/cpp/main.cpp +++ /dev/null @@ -1,1025 +0,0 @@ -/** - * engines/cpp/main.cpp - * -------------------- - * C++ 动力学模拟引擎。 - * 与 Python 版 (compute.py) 算法保持一致。 - * - * 输入: param.json, /coord.txt, connection.txt, bond.txt - * 输出: /trajectory.txt (JSON, 与 Python 版兼容) - * - * 编译: - * g++ -O3 -march=native -std=c++17 -o build/dynamics_cpp main.cpp - * - * 用法: - * ./build/dynamics_cpp - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// ======================================================================== -// 配置参数（从 param.json 读取） -// ======================================================================== -struct SimParams { - double box_a = 10.0; - int NT = 10000; - double DT = 0.001; - int NSTEP = 100; - int warmup_steps = 0; - std::string method = "leapfrog"; - double G[3] = {0, 0, -9.8}; - double B[3] = {0, 0, 0}; - int gravity_field = 1; - int gravity_interaction = 0; - int elastic_force = 1; - int damping_force = 0; - double gravity_strength = 1.0; - int driving_force = 0; - int save_trajectory = 1; - double alpha[6] = {0,0,0,0,0,0}; - double ball_radius = 0.5; - double ball_color[3] = {0.9, 0.2, 0.2}; - double box_color[3] = {0.8, 0.8, 0.85}; - int use_marker = 0; - double camera_distance = 40.0; - double camera_elevation = 0; - double camera_azimuth = 0; -}; - -// ======================================================================== -// 原子数据 -// ======================================================================== -struct AtomData { - std::vector ids; - std::vector masses; - std::vector radii; - std::vector pos_0; // (n_atoms * 3) - std::vector vel_0; // (n_atoms * 3) - std::vector fixed; // (n_atoms * 3), 0/1 flags -}; - -// ======================================================================== -// 成键数据 -// ======================================================================== -struct BondData { - std::vector pairs; // (n_bonds * 2) - std::vector stiffness; - 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 eq_x, eq_y, eq_z; // 平衡位置（初始坐标） - std::vector freeze_x, freeze_y, freeze_z; -}; - -// ======================================================================== -// 辅助函数 -// ======================================================================== - -static void die(const std::string &msg) { - std::cerr << "[C++-engine] 错误: " << msg << std::endl; - exit(1); -} - -/* 读整个文件为字符串 */ -static std::string read_file(const std::string &path) { - std::ifstream f(path, std::ios::binary); - if (!f) die("无法打开 " + path); - std::ostringstream ss; - ss << f.rdbuf(); - return ss.str(); -} - -/* 从 JSON 中查找 key，返回冒号后的数值 */ -static double json_read_double(const std::string &json, const std::string &key) { - auto pos = json.find("\"" + key + "\""); - if (pos == std::string::npos) return 0.0; - pos = json.find(':', pos); - if (pos == std::string::npos) return 0.0; - while (pos < json.size() && (json[pos] == ':' || json[pos] == ' ' || json[pos] == '\t' || json[pos] == '\n')) pos++; - return std::stod(json.substr(pos)); -} - -static int json_read_int(const std::string &json, const std::string &key) { - return static_cast(json_read_double(json, key)); -} - -/* 从 JSON 中读取字符串值 */ -static std::string json_read_string(const std::string &json, const std::string &key) { - auto pos = json.find("\"" + key + "\""); - if (pos == std::string::npos) return ""; - pos = json.find(':', pos); - if (pos == std::string::npos) return ""; - while (pos < json.size() && (json[pos] == ':' || json[pos] == ' ' || json[pos] == '\t' || json[pos] == '\n')) pos++; - if (pos >= json.size()) return ""; - // 找到引号 - if (json[pos] != '"') return ""; - pos++; - std::string result; - while (pos < json.size() && json[pos] != '"') { - result += json[pos]; - pos++; - } - return result; -} - -/* 检查 JSON 中是否存在某个 key */ -static bool json_has_key(const std::string &json, const std::string &key) { - return json.find("\"" + key + "\"") != std::string::npos; -} - -/* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */ -static void json_read_double3(const std::string &json, const std::string &key, double out[3]) { - auto pos = json.find("\"" + key + "\""); - if (pos == std::string::npos) { out[0] = out[1] = out[2] = 0; return; } - pos = json.find('[', pos); - if (pos == std::string::npos) { out[0] = out[1] = out[2] = 0; return; } - pos++; - for (int i = 0; i < 3; i++) { - while (pos < json.size() && (json[pos] == ' ' || json[pos] == '\t' || json[pos] == '\n' || json[pos] == ',')) pos++; - char *end; - out[i] = std::strtod(json.c_str() + pos, &end); - pos = end - json.c_str(); - } -} - -/* 读取 JSON 数组到 double[6] */ -static void json_read_double6(const std::string &json, const std::string &key, double out[6]) { - auto pos = json.find("\"" + key + "\""); - if (pos == std::string::npos) { for (int i=0;i<6;i++) out[i]=0; return; } - pos = json.find('[', pos); - if (pos == std::string::npos) { for (int i=0;i<6;i++) out[i]=0; return; } - pos++; - for (int i = 0; i < 6; i++) { - while (pos < json.size() && (json[pos]==' '||json[pos]=='\t'||json[pos]=='\n'||json[pos]==','||json[pos]==']')) pos++; - char *end; - out[i] = std::strtod(json.c_str() + pos, &end); - pos = end - json.c_str(); - } -} - -/* 解析 param.json */ -static SimParams read_params(const std::string &path) { - std::string buf = read_file(path); - SimParams p; - p.box_a = json_read_double(buf, "box_a"); - p.NT = json_read_int(buf, "NT"); - p.DT = json_read_double(buf, "DT"); - p.NSTEP = json_read_int(buf, "NSTEP"); - p.warmup_steps = json_read_int(buf, "warmup_steps"); - std::string m = json_read_string(buf, "method"); - if (!m.empty()) p.method = m; - json_read_double3(buf, "G", p.G); - json_read_double3(buf, "B", p.B); - p.gravity_field = json_read_int(buf, "gravity_field"); - p.gravity_interaction = json_read_int(buf, "gravity_interaction"); - 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"); - // save_trajectory 默认 1（全量），仅在 JSON 中存在该 key 时覆盖 - if (json_has_key(buf, "save_trajectory")) - p.save_trajectory = json_read_int(buf, "save_trajectory"); - json_read_double6(buf, "alpha", p.alpha); - p.ball_radius = json_read_double(buf, "ball_radius"); - json_read_double3(buf, "ball_color", p.ball_color); - json_read_double3(buf, "box_color", p.box_color); - p.use_marker = json_read_int(buf, "use_marker"); - p.camera_distance = json_read_double(buf, "camera_distance"); - p.camera_elevation = json_read_double(buf, "camera_elevation"); - p.camera_azimuth = json_read_double(buf, "camera_azimuth"); - return p; -} - -/* 读取 coord.txt */ -static AtomData read_coord(const std::string &input_dir) { - std::string path = input_dir + "/coord.txt"; - std::ifstream f(path); - if (!f) die("无法打开 " + path); - - std::string header; - std::getline(f, header); // 跳过表头 - - AtomData a; - int id, fx, fy, fz; - double mass, rad, px, py, pz, vx, vy, vz; - std::string line; - - while (std::getline(f, line)) { - if (line.empty() || line[0] == '#') continue; - int n_parsed = std::sscanf(line.c_str(), "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %d", - &id, &mass, &rad, &px, &py, &pz, &vx, &vy, &vz, &fx, &fy, &fz); - if (n_parsed == 9) { - fx = fy = fz = 0; - } else if (n_parsed != 12) { - continue; - } - a.ids.push_back(id); - a.masses.push_back(mass); - a.radii.push_back(rad); - a.pos_0.push_back(px); a.pos_0.push_back(py); a.pos_0.push_back(pz); - a.vel_0.push_back(vx); a.vel_0.push_back(vy); a.vel_0.push_back(vz); - a.fixed.push_back(fx); a.fixed.push_back(fy); a.fixed.push_back(fz); - } - - if (a.ids.empty()) die("coord.txt 中没有原子数据"); - return a; -} - -/* 读取 connection.txt 和 bond.txt */ -static BondData read_bonds(const std::string &input_dir) { - BondData b; - std::string conn_path = input_dir + "/connection.txt"; - std::ifstream f(conn_path); - if (!f) return b; // 无成键 - - std::string header; - std::getline(f, header); // 跳过表头 - - // 先读取 bond.txt 获得键参数映射 - std::string bond_path = input_dir + "/bond.txt"; - std::ifstream fb(bond_path); - - int a1, a2; - std::string bond_name; - std::vector> conn_lines; - while (f >> a1 >> a2 >> bond_name) { - conn_lines.emplace_back(a1 - 1, a2 - 1, bond_name); - } - - for (auto &[i, j, name] : conn_lines) { - double k = 1.0, r0 = 2.0; - if (fb) { - fb.clear(); - fb.seekg(0); - std::string bn, header; - double bk, br; - std::getline(fb, header); // 跳过表头行 - while (fb >> bn >> bk >> br) { - if (bn == name) { - k = bk; - r0 = br; - break; - } - } - } - b.pairs.push_back(i); - b.pairs.push_back(j); - b.stiffness.push_back(k); - b.rest_lengths.push_back(r0); - } - - 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.eq_x.push_back(atoms.pos_0[idx*3+0]); - d.eq_y.push_back(atoms.pos_0[idx*3+1]); - d.eq_z.push_back(atoms.pos_0[idx*3+2]); - 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; -} - -// ======================================================================== -// 物理核心 -// ======================================================================== - -/* 加速度计算（各力独立开关控制）——与 Python compute_acceleration 一致 */ -static void compute_acceleration( - int n, - const double *x, const double *y, const double *z, - const double *vx, const double *vy, const double *vz, - const double *m, const double G[3], const double B[3], - const BondData &bonds, - int gravity_field, int gravity_interaction, - int elastic_force, int damping_force, - double gravity_strength, - double *ax, double *ay, double *az) -{ - // 清零 - std::fill(ax, ax + n, 0.0); - std::fill(ay, ay + n, 0.0); - std::fill(az, az + n, 0.0); - - // 均匀重力场 - if (gravity_field) { - for (int i = 0; i < n; i++) { - ax[i] += G[0]; - ay[i] += G[1]; - az[i] += G[2]; - } - } - - // 阻尼 - if (damping_force) { - for (int i = 0; i < n; i++) { - ax[i] -= B[0] * vx[i] / m[i]; - ay[i] -= B[1] * vy[i] / m[i]; - az[i] -= B[2] * vz[i] / m[i]; - } - } - - // 弹簧力 - if (elastic_force) { - int nb = static_cast(bonds.stiffness.size()); - for (int b = 0; b < nb; b++) { - int i = bonds.pairs[b * 2]; - int j = bonds.pairs[b * 2 + 1]; - double dx = x[j] - x[i]; - double dy = y[j] - y[i]; - double dz = z[j] - z[i]; - double dist = std::sqrt(dx * dx + dy * dy + dz * dz); - if (dist < 1e-12) continue; - double stretch = dist - bonds.rest_lengths[b]; - double fmag = bonds.stiffness[b] * stretch; - double ux = dx / dist, uy = dy / dist, uz = dz / dist; - double fx = fmag * ux, fy = fmag * uy, fz = fmag * uz; - ax[i] += fx / m[i]; ay[i] += fy / m[i]; az[i] += fz / m[i]; - ax[j] -= fx / m[j]; ay[j] -= fy / m[j]; az[j] -= fz / m[j]; - } - } - - // 万有引力（所有原子对之间） - if (gravity_interaction) { - for (int i = 0; i < n; i++) { - for (int j = i + 1; j < n; j++) { - double dx = x[j] - x[i]; - double dy = y[j] - y[i]; - double dz = z[j] - z[i]; - double r2 = dx * dx + dy * dy + dz * dz; - if (r2 <= 1e-12) continue; - double r = std::sqrt(r2); - double f_mag = gravity_strength * m[i] * m[j] / r2; - double fx_g = f_mag * dx / r; - double fy_g = f_mag * dy / r; - double fz_g = f_mag * dz / r; - ax[i] += fx_g / m[i]; ay[i] += fy_g / m[i]; az[i] += fz_g / m[i]; - ax[j] -= fx_g / m[j]; ay[j] -= fy_g / m[j]; az[j] -= fz_g / m[j]; - } - } - } -} - -/* 保守力加速度（不含阻尼），供真蛙跳法专用。 - 传入 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 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; } - 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 / -// Midpoint_Method / Leapfrog_Method 保持一致 -// ======================================================================== - -/* ── 显式欧拉法 ──────────── */ -static void explicit_euler_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData &bonds, const int *fixed, double dt, - int gravity_field, int gravity_interaction, - int elastic_force, int damping_force, - double gravity_strength) -{ - std::vector 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()); - 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; - y[i] += vy[i] * dt; - z[i] += vz[i] * dt; - vx[i] += ax[i] * dt; - vy[i] += ay[i] * dt; - vz[i] += az[i] * dt; - } -} - -/* ── 隐式欧拉法 ──────────── */ -static void implicit_euler_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData &bonds, const int *fixed, double dt, - int gravity_field, int gravity_interaction, - int elastic_force, int damping_force, - double gravity_strength) -{ - std::vector ax(n), ay(n), az(n); - for (int i = 0; i < n; i++) { - if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) { - ax[i] = ay[i] = az[i] = 0; - continue; - } - double gamma_x = B[0] / m[i]; - double gamma_y = B[1] / m[i]; - double gamma_z = B[2] / m[i]; - // 隐式更新速度（重力 + 阻尼） - double vxn = (vx[i] + G[0] * dt) / (1.0 + gamma_x * dt); - double vyn = (vy[i] + G[1] * dt) / (1.0 + gamma_y * dt); - double vzn = (vz[i] + G[2] * dt) / (1.0 + gamma_z * dt); - // 用隐式速度 + 当前位置算加速度（包含各力开关） - // 注意：Python 中 compute_acceleration(x, y, z, vx_next, ...) 用新速度+旧位置 - double tpx = x[i], tpy = y[i], tpz = z[i]; - compute_acceleration(1, &tpx, &tpy, &tpz, &vxn, &vyn, &vzn, - &m[i], G, B, bonds, - gravity_field, gravity_interaction, - elastic_force, damping_force, gravity_strength, - &ax[i], &ay[i], &az[i]); - 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; - } -} - -/* ── 中点法 ──────────── */ -static void midpoint_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData &bonds, const int *fixed, double dt, - int gravity_field, int gravity_interaction, - int elastic_force, int damping_force, - double gravity_strength) -{ - std::vector 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()); - for (int i = 0; i < n; i++) { - if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; - double x_mid = x[i] + 0.5 * vx[i] * dt; - double y_mid = y[i] + 0.5 * vy[i] * dt; - double z_mid = z[i] + 0.5 * vz[i] * dt; - double vx_mid = vx[i] + 0.5 * ax[i] * dt; - double vy_mid = vy[i] + 0.5 * ay[i] * dt; - double vz_mid = vz[i] + 0.5 * az[i] * dt; - x[i] += vx_mid * dt; - y[i] += vy_mid * dt; - z[i] += vz_mid * dt; - double ax_mid, ay_mid, az_mid; - compute_acceleration(1, &x_mid, &y_mid, &z_mid, &vx_mid, &vy_mid, &vz_mid, - &m[i], G, B, bonds, - gravity_field, gravity_interaction, - elastic_force, damping_force, gravity_strength, - &ax_mid, &ay_mid, &az_mid); - vx[i] += ax_mid * dt; - vy[i] += ay_mid * dt; - vz[i] += az_mid * dt; - } -} - -/* ── 蛙跳法（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, - const double *m, const double G[3], const double B[3], - const BondData &bonds, const int *fixed, double dt, - int gravity_field, int gravity_interaction, - int elastic_force, int damping_force, - double gravity_strength) -{ - std::vector ax(n), ay(n), az(n); - - // 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; - 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; - } -} - -/* ── 分发器：调用对应积分方法 + 边界条件（与 Python apply_motion_update 一致）── */ -static void apply_step( - const std::string &method, - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const double G[3], const double B[3], - const BondData &bonds, const int *fixed, - const double *pos_0, - double box_a, double dt, - int gravity_field, int gravity_interaction, - int elastic_force, int damping_force, - double gravity_strength) -{ - // 积分 - if (method == "explicit_euler") { - explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt, - gravity_field, gravity_interaction, - elastic_force, damping_force, gravity_strength); - } else if (method == "implicit_euler") { - implicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt, - gravity_field, gravity_interaction, - elastic_force, damping_force, gravity_strength); - } else if (method == "midpoint") { - midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt, - gravity_field, gravity_interaction, - elastic_force, damping_force, gravity_strength); - } else if (method == "leapfrog") { - leapfrog_full_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt, - gravity_field, gravity_interaction, - elastic_force, damping_force, gravity_strength); - } else { - die("未知算法: " + method); - } - - // 边界条件（与 Python Limit_in_box 一致） - for (int i = 0; i < n; i++) { - if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; - limit_in_box(x[i], vx[i], -box_a, box_a); - limit_in_box(y[i], vy[i], -box_a, box_a); - 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; } - if (fixed[i*3+1]) { y[i] = pos_0[i*3+1]; vy[i] = 0.0; } - if (fixed[i*3+2]) { z[i] = pos_0[i*3+2]; vz[i] = 0.0; } - } -} - -// ======================================================================== -// 驱动力应用 -// ======================================================================== - -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.eq_x[d] + drivers.amp_x[d] * std::cos(TWO_PI * drivers.freq_x[d] * t + drivers.phi_x[d]); - double py = drivers.eq_y[d] + drivers.amp_y[d] * std::cos(TWO_PI * drivers.freq_y[d] * t + drivers.phi_y[d]); - double pz = drivers.eq_z[d] + 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; - } - } - } -} - -// ======================================================================== -// display.txt 输出（save_trajectory=0 时使用） -// ======================================================================== - -static void write_display_txt( - const std::string &path, - const std::vector &x, const std::vector &y, - const std::vector &z, const std::vector &vx, - const std::vector &vy, const std::vector &vz, - int n_steps, int n_atoms, - const SimParams ¶ms, const AtomData &atoms) -{ - std::ofstream f(path); - if (!f) die("无法写入 " + path); - std::cout << "[Cpp-engine] 正在写入显示数据…" << std::endl; - - int dynamic_steps = params.NT - params.warmup_steps; - double T_total = dynamic_steps * params.DT; - - /* 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: " << dynamic_steps << "\n"; - f << "T_total: " << T_total << "\n"; - f << "box_a: " << params.box_a << "\n"; - f << "driving_force: " << params.driving_force << "\n"; - f << "alpha: " << params.alpha[0] << "," << params.alpha[1] << "," << params.alpha[2] << "," - << params.alpha[3] << "," << params.alpha[4] << "," << params.alpha[5] << "\n"; - f << "ball_radius: " << params.ball_radius << "\n"; - f << "ball_color_r: " << params.ball_color[0] << "\n"; - f << "ball_color_g: " << params.ball_color[1] << "\n"; - f << "ball_color_b: " << params.ball_color[2] << "\n"; - f << "box_color_r: " << params.box_color[0] << "\n"; - f << "box_color_g: " << params.box_color[1] << "\n"; - f << "box_color_b: " << params.box_color[2] << "\n"; - f << "use_marker: " << params.use_marker << "\n"; - f << "camera_distance: " << params.camera_distance << "\n"; - f << "camera_elevation: " << params.camera_elevation << "\n"; - f << "camera_azimuth: " << params.camera_azimuth << "\n\n"; - - f << std::fixed << std::setprecision(6); - for (int t = 0; t < n_steps; t++) { - f << "frame: " << (t + 1) << "\n"; - f << "n x y z vx vy vz\n"; - for (int i = 0; i < n_atoms; i++) { - int base = t * n_atoms + i; - f << std::setw(3) << atoms.ids[i] - << std::setw(12) << x[base] - << std::setw(12) << y[base] - << std::setw(12) << z[base] - << std::setw(12) << vx[base] - << std::setw(12) << vy[base] - << std::setw(12) << vz[base] << "\n"; - } - } -} - -// ======================================================================== -// JSON 输出 -// ======================================================================== - -static void write_trajectory_json( - const std::string &path, - const std::vector &x, const std::vector &y, - const std::vector &z, const std::vector &vx, - const std::vector &vy, const std::vector &vz, - int n_steps, int n_atoms, - const SimParams ¶ms, const AtomData &atoms, const BondData &bonds) -{ - std::ofstream f(path); - if (!f) die("无法写入 " + path); - std::cout << "[Cpp-engine] 正在写入轨迹数据…" << std::endl; - f << std::setprecision(8); - - f << "{\n"; - - // 轨迹数组 - const std::string names[] = {"traj_x","traj_y","traj_z","traj_vx","traj_vy","traj_vz"}; - const std::vector *arrs[] = {&x, &y, &z, &vx, &vy, &vz}; - - for (int a = 0; a < 6; a++) { - f << " \"" << names[a] << "\": [\n"; - const auto &data = *arrs[a]; - for (int t = 0; t < n_steps; t++) { - f << " ["; - for (int i = 0; i < n_atoms; i++) { - f << data[t * n_atoms + i]; - if (i < n_atoms - 1) f << ','; - } - f << ']'; - if (t < n_steps - 1) f << ','; - f << '\n'; - } - f << " ],\n"; - } - - // 标量参数 - f << " \"NT\": " << params.NT << ",\n"; - f << " \"DT\": " << params.DT << ",\n"; - f << " \"NSTEP\": " << params.NSTEP << ",\n"; - f << " \"method\": \"" << params.method << "\",\n"; - f << " \"warmup_steps\": " << params.warmup_steps << ",\n"; - f << " \"G\": [" << params.G[0] << ", " << params.G[1] << ", " << params.G[2] << "],\n"; - f << " \"B\": [" << params.B[0] << ", " << params.B[1] << ", " << params.B[2] << "],\n"; - - // 原子信息 - f << " \"atom_ids\": ["; - for (size_t i = 0; i < atoms.ids.size(); i++) { - if (i > 0) f << ','; - f << atoms.ids[i]; - } - f << "],\n"; - - f << " \"atom_masses\": ["; - for (size_t i = 0; i < atoms.masses.size(); i++) { - if (i > 0) f << ','; - f << atoms.masses[i]; - } - f << "],\n"; - - // 成键 - f << " \"bond_pairs\": ["; - for (size_t b = 0; b < bonds.stiffness.size(); b++) { - if (b > 0) f << ','; - f << "[" << bonds.pairs[b * 2] << ", " << bonds.pairs[b * 2 + 1] << "]"; - } - f << "],\n"; - - f << " \"bond_stiffness\": ["; - for (size_t b = 0; b < bonds.stiffness.size(); b++) { - if (b > 0) f << ','; - f << bonds.stiffness[b]; - } - f << "],\n"; - - f << " \"bond_rest_lengths\": ["; - for (size_t b = 0; b < bonds.rest_lengths.size(); b++) { - if (b > 0) f << ','; - f << bonds.rest_lengths[b]; - } - f << "],\n"; - - f << " \"driving_force\": " << params.driving_force << "\n"; - - f << "}\n"; -} - -// ======================================================================== -// 主函数 -// ======================================================================== - -int main(int argc, char **argv) { - if (argc < 4) { - std::cerr << "用法: " << argv[0] << " " << std::endl; - return 1; - } - - std::string input_dir = argv[1]; - std::string output_dir = argv[2]; - std::string param_path = argv[3]; - - auto t0 = std::chrono::high_resolution_clock::now(); - - // 读取参数和输入 - SimParams params = read_params(param_path); - 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 - << ", method=" << params.method << std::endl; - - int n = static_cast(atoms.ids.size()); - - // 初始化位置/速度 - std::vector x(n), y(n), z(n), vx(n), vy(n), vz(n); - for (int i = 0; i < n; i++) { - x[i] = atoms.pos_0[i * 3]; - y[i] = atoms.pos_0[i * 3 + 1]; - z[i] = atoms.pos_0[i * 3 + 2]; - vx[i] = atoms.vel_0[i * 3]; - vy[i] = atoms.vel_0[i * 3 + 1]; - vz[i] = atoms.vel_0[i * 3 + 2]; - } - - // 分配轨迹缓冲区 - int record_steps = params.NT - params.warmup_steps; - int nstep_sampling = (params.save_trajectory == 0) ? params.NSTEP : 1; - int buf_steps = (params.save_trajectory == 0) ? (record_steps / nstep_sampling) : record_steps; - std::vector traj_x(buf_steps * n); - std::vector traj_y(buf_steps * n); - std::vector traj_z(buf_steps * n); - std::vector traj_vx(buf_steps * n); - std::vector traj_vy(buf_steps * n); - std::vector traj_vz(buf_steps * n); - - // 真蛙跳初始化：v(0) 反推 v(-dt/2) = v(0) - 0.5*a_c(0)*dt - if (params.method == "leapfrog") { - std::vector 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) - 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) - 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, - bonds, atoms.fixed.data(), - atoms.pos_0.data(), - params.box_a, params.DT, - params.gravity_field, params.gravity_interaction, - params.elastic_force, params.damping_force, params.gravity_strength); - } - - // 记录 - int _prog_int = record_steps / 100; - if (_prog_int < 1) _prog_int = 1; - int si = 0; // 采样帧索引 - for (int s = 0; s < record_steps; s++) { - if (s % _prog_int == 0 && s > 0) { - std::cout << "[Cpp-engine] progress: " << s << "/" << record_steps << std::endl; - } - 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); - // 保存当前帧（采样模式仅每 NSTEP 步保存一次） - if (s % nstep_sampling == 0) { - for (int i = 0; i < n; i++) { - traj_x[si * n + i] = x[i]; - traj_y[si * n + i] = y[i]; - traj_z[si * n + i] = z[i]; - traj_vx[si * n + i] = vx[i]; - traj_vy[si * n + i] = vy[i]; - traj_vz[si * n + i] = vz[i]; - } - si++; - } - - apply_step(params.method, n, x.data(), y.data(), z.data(), - vx.data(), vy.data(), vz.data(), - atoms.masses.data(), params.G, params.B, - bonds, atoms.fixed.data(), - atoms.pos_0.data(), - params.box_a, params.DT, - params.gravity_field, params.gravity_interaction, - params.elastic_force, params.damping_force, params.gravity_strength); - } - - // 输出轨迹 - if (params.save_trajectory == 0) { - std::string out_path = output_dir + "/display.txt"; - write_display_txt(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, - si, n, params, atoms); - } else { - std::string out_path = output_dir + "/trajectory.txt"; - write_trajectory_json(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, - record_steps, n, params, atoms, bonds); - } - - auto t1 = std::chrono::high_resolution_clock::now(); - double elapsed = std::chrono::duration(t1 - t0).count(); - std::cout << "[C++-engine] 计算完成: " << record_steps << " 步, " << elapsed << " s" << std::endl; - - return 0; -} diff --git a/engines/fortran/Makefile b/engines/fortran/Makefile index 5c2b6a3..723c91b 100644 --- a/engines/fortran/Makefile +++ b/engines/fortran/Makefile @@ -1,20 +1,18 @@ # engines/fortran/Makefile +# 编译 DLL（主程序通过 ctypes 直接调用） FC = gfortran FFLAGS = -O3 -march=native -Wall -Wextra -SRCS = main.f90 LIB_SRC = dynamics_lib.f90 UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) ifeq ($(UNAME_S),Windows) - STATIC_FLAGS = -static-libgcc -static-libgfortran -static-libquadmath + STATIC_FLAGS = -static else STATIC_FLAGS = endif -TARGET = build/dynamics_f90.exe - ifeq ($(UNAME_S),Linux) DLL_TARGET = build/dynamics_f90.so DLL_FLAGS = -shared -fPIC @@ -28,14 +26,10 @@ endif .PHONY: all dll clean -all: $(TARGET) +all: dll dll: $(DLL_TARGET) -$(TARGET): $(SRCS) | build - $(FC) $(FFLAGS) $(STATIC_FLAGS) -o $@ $(SRCS) - @echo " === Fortran engine built: $@ ===" - $(DLL_TARGET): $(LIB_SRC) | build $(FC) $(FFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) @echo " === Fortran DLL built: $@ ===" diff --git a/engines/fortran/main.f90 b/engines/fortran/main.f90 deleted file mode 100644 index 06a8afc..0000000 --- a/engines/fortran/main.f90 +++ /dev/null @@ -1,1134 +0,0 @@ -! engines/fortran/main.f90 -! ------------------------ -! Fortran 动力学模拟引擎。 -! 与 Python 版 (compute.py) 算法保持一致。 -! -! 输入: param.json, /coord.txt, connection.txt, bond.txt -! 输出: /trajectory.txt (JSON) -! -! 编译: cmake --build build --target dynamics_f90 -! 用法: ./build/dynamics_f90 - -program dynamics_f90 - implicit none - - character(len=256) :: input_dir, output_dir, param_path - integer :: narg, i, n, s - integer :: NT, NSTEP, warmup_steps, n_atoms, n_bonds - 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 - integer :: use_marker - character(len=512) :: alpha_str - double precision :: t0, t1, elapsed, tw - - ! 原子数据 - integer, allocatable :: atom_ids(:) - double precision, allocatable :: masses(:), radii(:) - double precision, allocatable :: pos_0(:, :), vel_0(:, :) - integer, allocatable :: fixed(:, :) - - ! 成键数据 - integer, allocatable :: bond_pairs(:, :) - double precision, allocatable :: bond_stiffness(:), bond_rest_lengths(:) - - ! 驱动力数据 - integer :: n_drivers, prog_step - 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_eq_x(:), drv_eq_y(:), drv_eq_z(:) - double precision, allocatable :: drv_freeze_x(:), drv_freeze_y(:), drv_freeze_z(:) - - ! 运行时位置/速度 - double precision, allocatable :: x(:), y(:), z(:) - double precision, allocatable :: vx(:), vy(:), vz(:) - - ! 轨迹缓冲区 - integer :: record_steps, n_frames, frame_idx - double precision, allocatable :: traj_x(:, :), traj_y(:, :), traj_z(:, :) - double precision, allocatable :: traj_vx(:, :), traj_vy(:, :), traj_vz(:, :) - - ! ======================================================================== - ! 主流程 - ! ======================================================================== - narg = command_argument_count() - if (narg < 3) then - write(*,*) "用法: dynamics_f90 " - stop - end if - - call get_command_argument(1, input_dir) - call get_command_argument(2, output_dir) - call get_command_argument(3, param_path) - - call cpu_time(t0) - - ! 读取 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, driving_force, & - use_marker, alpha_str) - - ! 读取 coord.txt - call read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed) - - ! 读取成键信息 - 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, pos_0, 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_eq_x, drv_eq_y, drv_eq_z, & - 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) - - ! 初始化位置/速度 - n = n_atoms - allocate(x(n), y(n), z(n), vx(n), vy(n), vz(n)) - do i = 1, n - x(i) = pos_0(i, 1); y(i) = pos_0(i, 2); z(i) = pos_0(i, 3) - vx(i) = vel_0(i, 1); vy(i) = vel_0(i, 2); vz(i) = vel_0(i, 3) - end do - - ! 分配轨迹缓冲区（只保存采样帧，不保存每一步） - record_steps = NT - warmup_steps - n_frames = max(1, record_steps / max(1, NSTEP)) - allocate(traj_x(n_frames, n), traj_y(n_frames, n), traj_z(n_frames, n)) - allocate(traj_vx(n_frames, n), traj_vy(n_frames, n), traj_vz(n_frames, 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 - 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_eq_x, drv_eq_y, drv_eq_z, & - 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 - 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_eq_x, drv_eq_y, drv_eq_z, & - 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, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, & - pos_0) - end do - - ! 记录（每 NSTEP 步采一帧） - prog_step = max(1, record_steps / 100) - frame_idx = 0 - do s = 1, record_steps - if (mod(s, prog_step) == 0) then - write(*, '("[Fortran-engine] progress: ", i0, "/", i0)') s, record_steps - flush(6) - end if - ! 采帧：在每个 NSTEP 区间的起始时刻记录 - if (mod(s-1, max(1, NSTEP)) == 0 .and. frame_idx < n_frames) then - frame_idx = frame_idx + 1 - traj_x(frame_idx, :) = x; traj_y(frame_idx, :) = y; traj_z(frame_idx, :) = z - traj_vx(frame_idx, :) = vx; traj_vy(frame_idx, :) = vy; traj_vz(frame_idx, :) = vz - end if - 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_eq_x, drv_eq_y, drv_eq_z, & - 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, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, & - pos_0) - end do - - ! 输出 display.txt - write(*, '("[Fortran-engine] 正在写入 display.txt (", i0, " 帧)…")') n_frames - flush(6) - call write_display_txt(output_dir, n_frames, n_atoms, atom_ids, & - traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, & - NT, DT, NSTEP, warmup_steps, method, G, B, & - n_bonds, gravity_field, elastic_force, damping_force, & - driving_force, box_a, gravity_strength, & - use_marker, alpha_str) - - call cpu_time(t1) - elapsed = t1 - t0 - write(*, '("[Fortran-engine] 计算完成: ", i0, " 步, ", f0.6, " s")') record_steps, elapsed - - deallocate(x, y, z, vx, vy, vz) - 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_eq_x, drv_eq_y, drv_eq_z) - deallocate(drv_freeze_x, drv_freeze_y, drv_freeze_z) - end if - -contains - -! ======================================================================== -! 读取 param.json -! ======================================================================== -subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, driving_force, & - use_marker, alpha_str) - 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, driving_force - integer, intent(out) :: use_marker - character(len=*), intent(out) :: method, alpha_str - character(len=8096) :: buf - character(len=256) :: line - integer :: u, ios, i - - box_a = 10.0d0; NT = 10000; DT = 0.001d0; NSTEP = 100; warmup_steps = 0 - method = 'leapfrog' - 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 - use_marker = 0 - alpha_str = '0.2' - - 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 - buf = '' - do - read(u, '(a)', iostat=ios) line - if (ios /= 0) exit - buf = trim(buf) // trim(line) // char(10) - end do - close(u) - if (ios > 0) return - - box_a = json_get_double(buf, 'box_a', box_a) - NT = json_get_int(buf, 'NT', NT) - DT = json_get_double(buf, 'DT', DT) - NSTEP = json_get_int(buf, 'NSTEP', NSTEP) - warmup_steps = json_get_int(buf, 'warmup_steps', warmup_steps) - call json_get_string(buf, 'method', method) - if (len_trim(method) == 0) method = 'leapfrog' - call json_get_double3(buf, 'G', G) - call json_get_double3(buf, 'B', B) - gravity_field = json_get_int(buf, 'gravity_field', 1) - gravity_interaction = json_get_int(buf, 'gravity_interaction', 0) - 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) - use_marker = json_get_int(buf, 'use_marker', 0) - call json_get_alpha(buf, 'alpha', alpha_str) -end subroutine read_params - -! ======================================================================== -! JSON 辅助解析 -! ======================================================================== -function json_get_double(buf, key, default) result(val) - character(len=*), intent(in) :: buf, key - double precision, intent(in) :: default - double precision :: val - integer :: p, ios - val = default - p = index(buf, '"' // key // '"') - if (p == 0) return - p = index(buf(p:), ':') + p - do while (p <= len(buf) .and. (buf(p:p) == ':' .or. buf(p:p) == ' ' & - .or. buf(p:p) == char(9) .or. buf(p:p) == char(10))) - p = p + 1 - end do - read(buf(p:), *, iostat=ios) val -end function json_get_double - -function json_get_int(buf, key, default) result(val) - character(len=*), intent(in) :: buf, key - integer, intent(in) :: default - integer :: val - val = nint(json_get_double(buf, key, dble(default))) -end function json_get_int - -subroutine json_get_double3(buf, key, vals) - character(len=*), intent(in) :: buf, key - double precision, intent(out) :: vals(3) - integer :: p, i, ios - vals = 0.0d0 - p = index(buf, '"' // key // '"') - if (p == 0) return - p = index(buf(p:), '[') + p - do i = 1, 3 - do while (buf(p:p) == ' ' .or. buf(p:p) == ',' .or. & - buf(p:p) == char(9) .or. buf(p:p) == char(10)) - p = p + 1 - end do - read(buf(p:), *, iostat=ios) vals(i) - if (ios /= 0) exit - do while (p <= len(buf) .and. buf(p:p) /= ',' .and. buf(p:p) /= ']') - p = p + 1 - end do - end do -end subroutine json_get_double3 - -! 读取 alpha：兼容单值 (0.2) 或 6 元数组 [0.0,0.0,0.0,0.0,0.0,0.0] -! 输出为逗号分隔字符串，与 draw.py 的解析格式一致 -subroutine json_get_alpha(buf, key, val) - character(len=*), intent(in) :: buf, key - character(len=*), intent(out) :: val - double precision :: d6(6) - integer :: p, i, ios - character(len=32) :: tmp - val = '0.2' - p = index(buf, '"' // key // '"') - if (p == 0) return - p = index(buf(p:), ':') + p - ! skip whitespace - do while (p <= len(buf) .and. (buf(p:p) == ':' .or. buf(p:p) == ' ' & - .or. buf(p:p) == char(9) .or. buf(p:p) == char(10))) - p = p + 1 - end do - if (buf(p:p) == '[') then - ! 6 元数组 - d6 = 0.2d0 - p = p + 1 - do i = 1, 6 - do while (p <= len(buf) .and. (buf(p:p) == ' ' .or. buf(p:p) == ',' & - .or. buf(p:p) == char(9) .or. buf(p:p) == char(10))) - p = p + 1 - end do - if (p > len(buf) .or. buf(p:p) == ']') exit - read(buf(p:), *, iostat=ios) d6(i) - if (ios /= 0) exit - do while (p <= len(buf) .and. buf(p:p) /= ',' .and. buf(p:p) /= ']') - p = p + 1 - end do - end do - val = '' - do i = 1, 6 - write(tmp, '(f0.6)') d6(i) - if (i > 1) val = trim(val) // ',' - val = trim(val) // trim(tmp) - end do - else - ! 单值标量 - read(buf(p:), *, iostat=ios) d6(1) - if (ios == 0) then - write(tmp, '(f0.6)') d6(1) - val = trim(tmp) - end if - end if -end subroutine json_get_alpha - -! 从 JSON 中读取字符串值 -subroutine json_get_string(buf, key, val) - character(len=*), intent(in) :: buf, key - character(len=*), intent(out) :: val - integer :: p, i, j - val = '' - p = index(buf, '"' // key // '"') - if (p == 0) return - p = index(buf(p:), ':') + p - do while (p <= len(buf) .and. (buf(p:p) == ':' .or. buf(p:p) == ' ' & - .or. buf(p:p) == char(9) .or. buf(p:p) == char(10))) - p = p + 1 - end do - if (buf(p:p) /= '"') return - p = p + 1 - i = 1 - do while (p <= len(buf) .and. buf(p:p) /= '"' .and. i <= len(val)) - val(i:i) = buf(p:p) - i = i + 1; p = p + 1 - end do -end subroutine json_get_string - -! ======================================================================== -! 读取 coord.txt -! ======================================================================== -subroutine read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed) - character(len=*), intent(in) :: input_dir - integer, intent(out) :: n_atoms - integer, allocatable, intent(out) :: atom_ids(:) - double precision, allocatable, intent(out) :: masses(:), radii(:) - double precision, allocatable, intent(out) :: pos_0(:, :), vel_0(:, :) - integer, allocatable, intent(out) :: fixed(:, :) - - character(len=512) :: path, line - integer :: u, ios, ncols, n_cap, id, fx, fy, fz, i - double precision :: mass, rad, px, py, pz, vx, vy, vz - integer, parameter :: MX = 4096 - integer :: ids_tmp(MX), fix_tmp(MX, 3) - double precision :: m_tmp(MX), r_tmp(MX), p_tmp(MX, 3), v_tmp(MX, 3) - - n_atoms = 0 - path = trim(input_dir) // '/coord.txt' - open(newunit=u, file=trim(path), status='old', action='read', iostat=ios) - if (ios /= 0) then; write(*, '("[Fortran-engine] 错误: 无法打开 ", a)') trim(path); stop; end if - - read(u, '(a)', iostat=ios) line - ncols = 0 - do i = 1, len_trim(line) - if (line(i:i) /= ' ') then - if (i == 1) then - ncols = ncols + 1 - else if (line(i-1:i-1) == ' ') then - ncols = ncols + 1 - end if - end if - end do - - do - read(u, '(a)', iostat=ios) line - if (ios /= 0) exit - if (len_trim(line) == 0 .or. line(1:1) == '#') cycle - n_atoms = n_atoms + 1 - if (n_atoms > MX) exit - if (ncols >= 12) then - read(line, *) id, mass, rad, px, py, pz, vx, vy, vz, fx, fy, fz - else - read(line, *) id, mass, rad, px, py, pz, vx, vy, vz - fx = 0; fy = 0; fz = 0 - end if - ids_tmp(n_atoms) = id - m_tmp(n_atoms) = mass; r_tmp(n_atoms) = rad - p_tmp(n_atoms, :) = (/ px, py, pz /) - v_tmp(n_atoms, :) = (/ vx, vy, vz /) - fix_tmp(n_atoms, :) = (/ fx, fy, fz /) - end do - close(u) - if (n_atoms <= 0) then; write(*, '("[Fortran-engine] 错误: coord.txt 为空")') ; stop; end if - - allocate(atom_ids(n_atoms), masses(n_atoms), radii(n_atoms)) - allocate(pos_0(n_atoms, 3), vel_0(n_atoms, 3), fixed(n_atoms, 3)) - atom_ids = ids_tmp(1:n_atoms) - masses = m_tmp(1:n_atoms); radii = r_tmp(1:n_atoms) - pos_0 = p_tmp(1:n_atoms, :); vel_0 = v_tmp(1:n_atoms, :) - fixed = fix_tmp(1:n_atoms, :) -end subroutine read_coord - -! ======================================================================== -! 读取成键信息 (connection.txt + bond.txt) -! ======================================================================== -subroutine read_bonds(input_dir, n_atoms, atom_ids, pos_0, n_bonds, & - bond_pairs, bond_stiffness, bond_rest_lengths) - character(len=*), intent(in) :: input_dir - integer, intent(in) :: n_atoms, atom_ids(n_atoms) - double precision, intent(in) :: pos_0(n_atoms, 3) - integer, intent(out) :: n_bonds - integer, allocatable, intent(out) :: bond_pairs(:, :) - double precision, allocatable, intent(out) :: bond_stiffness(:), bond_rest_lengths(:) - - character(len=512) :: conn_path - integer :: u, ios, a1, a2, idx1, idx2, i - integer, parameter :: MX = 4096 - integer :: ptmp(MX, 2), idx_map(9999) - double precision :: ktmp(MX), r0tmp(MX), dx, dy, dz - character(len=64) :: name - - n_bonds = 0; idx_map = -1 - do i = 1, n_atoms - if (atom_ids(i) > 0 .and. atom_ids(i) <= 9998) idx_map(atom_ids(i)) = i - end do - - conn_path = trim(input_dir) // '/connection.txt' - open(newunit=u, file=trim(conn_path), status='old', action='read', iostat=ios) - if (ios /= 0) return - - read(u, *, iostat=ios) - do - read(u, *, iostat=ios) a1, a2, name - if (ios /= 0) exit - call find_bond(trim(input_dir), trim(name), ktmp(n_bonds+1), r0tmp(n_bonds+1)) - idx1 = idx_map(a1); idx2 = idx_map(a2) - if (idx1 < 1 .or. idx2 < 1) cycle - if (r0tmp(n_bonds+1) <= 0.0d0) then - dx = pos_0(idx2, 1) - pos_0(idx1, 1) - dy = pos_0(idx2, 2) - pos_0(idx1, 2) - dz = pos_0(idx2, 3) - pos_0(idx1, 3) - r0tmp(n_bonds+1) = sqrt(dx*dx + dy*dy + dz*dz) - end if - n_bonds = n_bonds + 1 - ptmp(n_bonds, :) = (/ idx1 - 1, idx2 - 1 /) - if (n_bonds >= MX) exit - end do - close(u) - - if (n_bonds <= 0) return - allocate(bond_pairs(n_bonds, 2), bond_stiffness(n_bonds), bond_rest_lengths(n_bonds)) - bond_pairs(:, :) = ptmp(1:n_bonds, :) - bond_stiffness = ktmp(1:n_bonds) - bond_rest_lengths = r0tmp(1:n_bonds) -end subroutine read_bonds - -subroutine find_bond(bond_dir, bond_name, k, r0) - character(len=*), intent(in) :: bond_dir, bond_name - double precision, intent(out) :: k, r0 - character(len=512) :: path, bn - integer :: u, ios - double precision :: bk, br - k = 1.0d0; r0 = -1.0d0 - path = trim(bond_dir) // '/bond.txt' - open(newunit=u, file=trim(path), status='old', action='read', iostat=ios) - if (ios /= 0) return - read(u, *, iostat=ios) - do - read(u, *, iostat=ios) bn, bk, br - if (ios /= 0) exit - if (trim(bn) == trim(bond_name)) then - k = bk - if (br > 0.0d0) r0 = br - exit - end if - end do - close(u) -end subroutine find_bond - -! ======================================================================== -! 物理核心 -! ======================================================================== - -! 加速度计算（各力独立开关控制） -pure subroutine 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) - integer, intent(in) :: n, nb, bp(nb, 2) - integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force - double precision, intent(in) :: 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), gravity_strength - double precision, intent(out) :: ax(n), ay(n), az(n) - integer :: i, j, ib - double precision :: dx, dy, dz, dist, s, fm, ux, uy, uz, r2, fg - - ! 清零 - ax = 0.0d0; ay = 0.0d0; az = 0.0d0 - - ! 均匀重力场 - if (gravity_field /= 0) then - do i = 1, n - ax(i) = G(1); ay(i) = G(2); az(i) = G(3) - end do - end if - - ! 阻尼 - if (damping_force /= 0) then - do i = 1, n - ax(i) = ax(i) - B(1) * vx(i) / m(i) - ay(i) = ay(i) - B(2) * vy(i) / m(i) - az(i) = az(i) - B(3) * vz(i) / m(i) - end do - end if - - ! 弹簧键力 - if (elastic_force /= 0) then - do ib = 1, nb - i = bp(ib, 1) + 1; j = bp(ib, 2) + 1 - dx = x(j) - x(i); dy = y(j) - y(i); dz = z(j) - z(i) - dist = sqrt(dx*dx + dy*dy + dz*dz) - if (dist < 1.0d-12) cycle - s = (dist - br(ib)) / dist - fm = bk(ib) * s - ux = fm * dx; uy = fm * dy; uz = fm * dz - ax(i) = ax(i) + ux / m(i); ay(i) = ay(i) + uy / m(i); az(i) = az(i) + uz / m(i) - ax(j) = ax(j) - ux / m(j); ay(j) = ay(j) - uy / m(j); az(j) = az(j) - uz / m(j) - end do - end if - - ! 万有引力（所有原子对之间） - if (gravity_interaction /= 0) then - do i = 1, n - 1 - do j = i + 1, n - dx = x(j) - x(i); dy = y(j) - y(i); dz = z(j) - z(i) - r2 = dx*dx + dy*dy + dz*dz - if (r2 <= 1.0d-12) cycle - dist = sqrt(r2) - fg = gravity_strength * m(i) * m(j) / r2 - ux = fg * dx / dist; uy = fg * dy / dist; uz = fg * dz / dist - ax(i) = ax(i) + ux / m(i); ay(i) = ay(i) + uy / m(i); az(i) = az(i) + uz / m(i) - ax(j) = ax(j) - ux / m(j); ay(j) = ay(j) - uy / m(j); az(j) = az(j) - uz / m(j) - end do - end do - 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 - double precision, intent(in) :: lo, hi - if (pos > hi) then - pos = hi; vel = -vel - else if (pos < lo) then - pos = lo; vel = -vel - 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, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength) - integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3) - integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force - 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) - 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) - do i = 1, n - if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle - x(i) = x(i) + vx(i) * dt - y(i) = y(i) + vy(i) * dt - z(i) = z(i) + vz(i) * dt - vx(i) = vx(i) + ax(i) * dt - vy(i) = vy(i) + ay(i) * dt - vz(i) = vz(i) + az(i) * dt - end do -end subroutine explicit_euler_step - -! ── 隐式欧拉法 ──────────── -subroutine implicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, & - nb, bp, bk, br, fixed, dt, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength) - integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3) - integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force - 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 :: vxn(n), vyn(n), vzn(n) - double precision :: gx, gy, gz - integer :: i - - ! 隐式速度（重力 + 阻尼） - do i = 1, n - if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then - vxn(i) = 0; vyn(i) = 0; vzn(i) = 0; cycle - end if - gx = B(1) / m(i); gy = B(2) / m(i); gz = B(3) / m(i) - vxn(i) = (vx(i) + G(1) * dt) / (1.0d0 + gx * dt) - vyn(i) = (vy(i) + G(2) * dt) / (1.0d0 + gy * dt) - vzn(i) = (vz(i) + G(3) * dt) / (1.0d0 + gz * dt) - end do - - ! 用当前位 + 隐式速度算完整加速度 - call accel(n, x, y, z, vxn, vyn, vzn, m, G, B, nb, bp, bk, br, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, ax, ay, az) - - 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 - vy(i) = vy(i) + ay(i) * dt - vz(i) = vz(i) + az(i) * dt - x(i) = x(i) + vx(i) * dt - y(i) = y(i) + vy(i) * dt - z(i) = z(i) + vz(i) * dt - end do -end subroutine implicit_euler_step - -! ── 中点法 ──────────── -subroutine midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, & - nb, bp, bk, br, fixed, dt, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength) - integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3) - integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force - 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 :: xm(n), ym(n), zm(n), vxm(n), vym(n), vzm(n) - 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) - do i = 1, n - if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then - xm(i)=0; ym(i)=0; zm(i)=0; vxm(i)=0; vym(i)=0; vzm(i)=0; cycle - end if - xm(i) = x(i) + 0.5d0 * vx(i) * dt - ym(i) = y(i) + 0.5d0 * vy(i) * dt - zm(i) = z(i) + 0.5d0 * vz(i) * dt - vxm(i) = vx(i) + 0.5d0 * ax(i) * dt - vym(i) = vy(i) + 0.5d0 * ay(i) * dt - vzm(i) = vz(i) + 0.5d0 * az(i) * dt - x(i) = x(i) + vxm(i) * dt - y(i) = y(i) + vym(i) * dt - z(i) = z(i) + vzm(i) * dt - end do - - call accel(n, xm, ym, zm, vxm, vym, vzm, m, G, B, nb, bp, bk, br, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, ax, ay, az) - 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 - vy(i) = vy(i) + ay(i) * dt - vz(i) = vz(i) + az(i) * dt - end do -end subroutine midpoint_step - -! 真蛙跳一步：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, & - elastic_force, damping_force, gravity_strength) - integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3) - integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force - 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 :: alphax, alphay, alphaz - logical :: has_damping - integer :: i - - ! 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 - 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 -end subroutine leapfrog_full - -! ── 分发器 + 边界条件 + 自由度约束 ── -subroutine apply_step(method, n, x, y, z, vx, vy, vz, m, G, B, & - nb, bp, bk, br, fixed, box_a, dt, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, & - pos_0) - character(len=*), intent(in) :: method - integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3) - integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force - 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), box_a, dt, gravity_strength - double precision, intent(in) :: pos_0(n, 3) - double precision :: mm(n) - integer :: i - - mm = m - - if (trim(method) == 'explicit_euler') then - call explicit_euler_step(n, x, y, z, vx, vy, vz, mm, G, B, & - nb, bp, bk, br, fixed, dt, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength) - else if (trim(method) == 'implicit_euler') then - call implicit_euler_step(n, x, y, z, vx, vy, vz, mm, G, B, & - nb, bp, bk, br, fixed, dt, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength) - else if (trim(method) == 'midpoint') then - call midpoint_step(n, x, y, z, vx, vy, vz, mm, G, B, & - nb, bp, bk, br, fixed, dt, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength) - else if (trim(method) == 'leapfrog') then - call leapfrog_full(n, x, y, z, vx, vy, vz, mm, G, B, & - nb, bp, bk, br, fixed, dt, & - gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength) - else - write(*, '("[Fortran-engine] 未知算法: ", a)') trim(method) - stop - end if - - ! 边界条件 - do i = 1, n - if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle - call limit_in_box(x(i), vx(i), -box_a, box_a) - call limit_in_box(y(i), vy(i), -box_a, box_a) - 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 - if (fixed(i,2) /= 0) then; y(i) = pos_0(i, 2); vy(i) = 0.0d0; end if - if (fixed(i,3) /= 0) then; z(i) = pos_0(i, 3); vz(i) = 0.0d0; end if - end do -end subroutine apply_step - -! ======================================================================== -! 读取驱动力 driver.txt -! ======================================================================== -subroutine read_driver(input_dir, n_atoms, atom_ids, pos_0, 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_eq_x, drv_eq_y, drv_eq_z, & - drv_freeze_x, drv_freeze_y, drv_freeze_z) - character(len=*), intent(in) :: input_dir - integer, intent(in) :: n_atoms, atom_ids(n_atoms) - double precision, intent(in) :: pos_0(n_atoms, 3) - 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_eq_x(:), drv_eq_y(:), drv_eq_z(:) - 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) - double precision :: eqx_tmp(MX), eqy_tmp(MX), eqz_tmp(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 - eqx_tmp(n_drivers) = pos_0(idx, 1) - eqy_tmp(n_drivers) = pos_0(idx, 2) - eqz_tmp(n_drivers) = pos_0(idx, 3) - 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_eq_x(n_drivers), drv_eq_y(n_drivers), drv_eq_z(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_eq_x = eqx_tmp(1:n_drivers); drv_eq_y = eqy_tmp(1:n_drivers); drv_eq_z = eqz_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_eq_x, drv_eq_y, drv_eq_z, & - 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(in) :: drv_eq_x(n_drivers), drv_eq_y(n_drivers), drv_eq_z(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_eq_x(d) + drv_amp_x(d) * cos(TWO_PI * drv_freq_x(d) * t + drv_phi_x(d)) - py = drv_eq_y(d) + drv_amp_y(d) * cos(TWO_PI * drv_freq_y(d) * t + drv_phi_y(d)) - pz = drv_eq_z(d) + 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 - -! ======================================================================== -! display.txt 输出（与 compute.py save_display_txt 格式一致） -! ======================================================================== - -subroutine write_display_txt(outdir, n_frames, nat, aid, & - tx, ty, tz, tvx, tvy, tvz, & - NT, DT, NSTEP, warmup, method, G, B, & - nb, gravity_field, elastic_force, damping_force, & - driving_force, box_a, gravity_strength, & - use_marker, alpha_str) - character(len=*), intent(in) :: outdir, method, alpha_str - integer, intent(in) :: n_frames, nat, NT, NSTEP, warmup, nb - integer, intent(in) :: gravity_field, elastic_force, damping_force, driving_force - integer, intent(in) :: use_marker - integer, intent(in) :: aid(nat) - double precision, intent(in) :: tx(n_frames, nat), ty(n_frames, nat), tz(n_frames, nat) - double precision, intent(in) :: tvx(n_frames, nat), tvy(n_frames, nat), tvz(n_frames, nat) - double precision, intent(in) :: DT, G(3), B(3), box_a, gravity_strength - - character(len=512) :: path, buf - integer :: u, f, a, ios - integer :: dynamic_steps - double precision :: T_total - - dynamic_steps = NT - warmup - T_total = NT * DT - - path = trim(outdir) // '/display.txt' - open(newunit=u, file=trim(path), status='replace', action='write', iostat=ios) - if (ios /= 0) then - write(*, '("[Fortran-engine] 错误: 无法写入 ", a)') trim(path) - return - end if - - ! ── header ──────────────────────────────────────────────────────────── - write(u, '("number of frames: ", i0)') n_frames - write(u, '("number of particles: ", i0)') nat - write(u, '("DT: ", g0)') DT - write(u, '("NSTEP: ", i0)') NSTEP - write(u, '("method: ", a)') trim(method) - write(u, '("NT: ", i0)') NT - write(u, '("warmup_steps: ", i0)') warmup - write(u, '("dynamic_steps: ", i0)') dynamic_steps - write(u, '("T_total: ", g0)') T_total - write(u, '("box_a: ", g0)') box_a - write(u, '("gravity_field: ", i0)') gravity_field - write(u, '("elastic_force: ", i0)') elastic_force - write(u, '("damping_force: ", i0)') damping_force - write(u, '("driving_force: ", i0)') driving_force - write(u, '("gravity_strength: ", g0)') gravity_strength - write(buf, '("G: [", g0, ", ", g0, ", ", g0, "]")') G(1), G(2), G(3) - write(u, '(a)') trim(buf) - write(buf, '("B: [", g0, ", ", g0, ", ", g0, "]")') B(1), B(2), B(3) - write(u, '(a)') trim(buf) - write(u, '("number_of_frames: ", i0)') n_frames - write(u, '("number_of_particles: ", i0)') nat - write(buf, '("X_MIN: ", g0)') -box_a; write(u, '(a)') trim(buf) - write(buf, '("X_MAX: ", g0)') box_a; write(u, '(a)') trim(buf) - write(buf, '("Y_MIN: ", g0)') -box_a; write(u, '(a)') trim(buf) - write(buf, '("Y_MAX: ", g0)') box_a; write(u, '(a)') trim(buf) - write(buf, '("Z_MIN: ", g0)') -box_a; write(u, '(a)') trim(buf) - write(buf, '("Z_MAX: ", g0)') box_a; write(u, '(a)') trim(buf) - write(u, '("use_marker: ", i0)') use_marker - write(u, '("alpha: ", a)') trim(alpha_str) - - ! ── frame data ──────────────────────────────────────────────────────── - do f = 1, n_frames - write(u, '()') ! 空行 - write(u, '("frame: ", i0)') f - write(u, '("n x y z vx vy vz")') - do a = 1, nat - write(u, '(i0, 6(f13.6))') aid(a), & - tx(f,a), ty(f,a), tz(f,a), tvx(f,a), tvy(f,a), tvz(f,a) - end do - end do - - close(u) - write(*, '("[Fortran-engine] display.txt 已保存: ", a)') trim(path) - flush(6) -end subroutine write_display_txt - -end program dynamics_f90