diff --git a/.gitignore b/.gitignore index 972facf..952390b 100644 --- a/.gitignore +++ b/.gitignore @@ -18,10 +18,9 @@ pip-wheel-metadata/ venv/ ENV/ -# ── C / C++ 编译产物 ───────────────────────────────────────── -# Makefile 构建输出（engines/c/build/） -engines/c/build/ -engines/cpp/build/ +# ── C / C++ / Fortran 编译产物 ──────────────────────── +# 源码目录 engines/src/*/ 中可能产生的构建输出 +engines/src/*/build/ # CMake 构建目录（根目录或自定义 build 目录） CMakeCache.txt @@ -46,7 +45,7 @@ build_*/ # 可执行文件（保留源码，排除编译出的二进制） # 注意：Windows 下 .exe 后缀的可执行文件 *.exe -# 但 engines/c/Makefile 里指定了 build/ 目录，已由上面覆盖 +# 可在 engines/src/*/ 中用 make dll 编译引擎 DLL # 运行时生成的引擎参数文件（每次运行都会覆盖） engines/*/param.json diff --git a/engines/c/Makefile b/engines/c/Makefile deleted file mode 100644 index 495f8d9..0000000 --- a/engines/c/Makefile +++ /dev/null @@ -1,41 +0,0 @@ -# engines/c/Makefile -# 编译 DLL（主程序通过 ctypes 直接调用） -# make dll → 本地系统编译 -# make linux → Linux 交叉编译（需 x86_64-linux-gnu-gcc） -# make windows → Windows 交叉编译（需 x86_64-w64-mingw32-gcc） - -CC = gcc -CFLAGS = -O3 -march=native -Wall -Wextra -LDFLAGS = -lm -LIB_SRC = dynamics_lib.c - -# 自动检测系统 -UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) - -# DLL 目标（平台自动选择后缀） -ifeq ($(UNAME_S),Linux) - DLL_TARGET = build/dynamics_c.so - DLL_FLAGS = -shared -fPIC -else ifeq ($(UNAME_S),Darwin) - DLL_TARGET = build/dynamics_c.dylib - DLL_FLAGS = -dynamiclib -else - DLL_TARGET = build/dynamics_c.dll - DLL_FLAGS = -shared -endif - -.PHONY: all dll clean - -all: dll - -dll: $(DLL_TARGET) - -$(DLL_TARGET): $(LIB_SRC) | build - $(CC) $(CFLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) $(LDFLAGS) - @echo " === C DLL built: $@ ===" - -build: - mkdir -p build - -clean: - rm -rf build *.o diff --git a/engines/c/dynamics_lib.c b/engines/c/dynamics_lib.c deleted file mode 100644 index c83d984..0000000 --- a/engines/c/dynamics_lib.c +++ /dev/null @@ -1,554 +0,0 @@ -/** - * engines/c/dynamics_lib.c - * ------------------------- - * 纯计算 DLL：无文件 I/O，所有数据由 Python 以 NumPy 数组传入， - * 结果直接写入 Python 预分配的输出数组。 - * 算法与 main.c 和 compute.py 保持完全一致。 - * - * 编译（Windows DLL）: - * gcc -O3 -march=native -shared -o build/dynamics_c.dll dynamics_lib.c -lm - * 编译（Linux .so）: - * gcc -O3 -march=native -shared -fPIC -o build/dynamics_c.so dynamics_lib.c -lm - * 编译（macOS .dylib）: - * gcc -O3 -march=native -dynamiclib -o build/dynamics_c.dylib dynamics_lib.c -lm - */ - -#ifdef _WIN32 -# define EXPORT __declspec(dllexport) -#else -# define EXPORT __attribute__((visibility("default"))) -#endif - -#include -#include -#include -#include - -/* ── 驱动力结构体 ─────────────────────────────────────────── */ -typedef struct { - int n_drivers; - const int *idx; /* [n_drivers] 0-based local atom index */ - const double *amp; /* [n_drivers*3] (ax,ay,az) interleaved */ - const double *freq; /* [n_drivers*3] */ - const double *phi; /* [n_drivers*3] radians */ - const double *eq; /* [n_drivers*3] equilibrium positions */ - const double *ncycles; /* [n_drivers] 0=unlimited */ - const int *has_period; /* [n_drivers] */ - /* mutable freeze positions (allocated internally) */ - double *freeze; /* [n_drivers*3] */ -} Drivers; - -/* ── 加速度：保守力（弹簧键 + 均匀重力场）────────────────── */ -static void accel_conservative( - int n, const double *x, const double *y, const double *z, - const double *m, - double Gx, double Gy, double Gz, - int gravity_field, int elastic_force, - int n_bonds, const int *bond_pairs, - const double *bond_k, const double *bond_r0, - double *ax, double *ay, double *az) -{ - for (int i = 0; i < n; i++) { - ax[i] = gravity_field ? Gx : 0.0; - ay[i] = gravity_field ? Gy : 0.0; - az[i] = gravity_field ? Gz : 0.0; - } - - if (!elastic_force || n_bonds == 0) return; - - for (int b = 0; b < n_bonds; b++) { - int ii = bond_pairs[b*2]; - int jj = bond_pairs[b*2+1]; - double dx = x[jj] - x[ii]; - double dy = y[jj] - y[ii]; - double dz = z[jj] - z[ii]; - double dist = sqrt(dx*dx + dy*dy + dz*dz); - if (dist < 1e-12) continue; - double k = bond_k[b]; - double r0 = bond_r0[b]; - double fac = k * (dist - r0) / dist; - double fx = fac * dx, fy = fac * dy, fz_b = fac * dz; - ax[ii] += fx / m[ii]; ay[ii] += fy / m[ii]; az[ii] += fz_b / m[ii]; - ax[jj] -= fx / m[jj]; ay[jj] -= fy / m[jj]; az[jj] -= fz_b / m[jj]; - } -} - -/* ── 完整加速度（含阻尼）────────────────────────────────── */ -static void accel_full( - int n, const double *x, const double *y, const double *z, - const double *vx, const double *vy, const double *vz, - const double *m, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bond_pairs, - const double *bond_k, const double *bond_r0, - double *ax, double *ay, double *az) -{ - accel_conservative(n, x, y, z, m, Gx, Gy, Gz, - gravity_field, elastic_force, - n_bonds, bond_pairs, bond_k, bond_r0, - ax, ay, az); - if (damping_force) { - for (int i = 0; i < n; i++) { - ax[i] -= Bx * vx[i] / m[i]; - ay[i] -= By * vy[i] / m[i]; - az[i] -= Bz * vz[i] / m[i]; - } - } -} - -/* ── 边界：反弹（与 main.c limit_in_box 一致）────────────── */ -static inline void _limit1(double *p, double *v, double lo, double hi) { - if (*p > hi) { *p = hi; *v = -fabs(*v); } - if (*p < lo) { *p = lo; *v = fabs(*v); } -} - -/* ── 边界：回绕（与 main.c wrap_position 一致）──────────── */ -static inline void _wrap1(double *p, double lo, double hi) { - if (*p > hi) *p = lo; - if (*p < lo) *p = hi; -} - -/* ── 边界 + 固定约束（与 main.c apply_step 末尾一致）──────── */ -static void apply_boundary_and_constraints( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const int *fixed, const double *pos_init, - double box_a) -{ - double lo = -box_a, hi = box_a; - - /* 反弹 */ - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - _limit1(&x[i], &vx[i], lo, hi); - _limit1(&y[i], &vy[i], lo, hi); - _limit1(&z[i], &vz[i], lo, hi); - } - - /* 回绕 */ - for (int i = 0; i < n; i++) { - _wrap1(&x[i], lo, hi); - _wrap1(&y[i], lo, hi); - _wrap1(&z[i], lo, hi); - } - - /* 逐自由度固定约束：与 main.c 和 Python apply_fixed_constraints 一致 */ - for (int i = 0; i < n; i++) { - if (fixed[i*3+0]) { x[i] = pos_init[i*3+0]; vx[i] = 0.0; } - if (fixed[i*3+1]) { y[i] = pos_init[i*3+1]; vy[i] = 0.0; } - if (fixed[i*3+2]) { z[i] = pos_init[i*3+2]; vz[i] = 0.0; } - } -} - -/* ══════════════════════════════════════════════════════════ - * 蛙跳法（与 main.c leapfrog_step 完全一致） - * x(t), v(t-dt/2) → x(t+dt), v(t+dt/2) - * 无阻尼：纯辛积分。有阻尼：半隐式处理 α = 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 int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, - double dt) -{ - double *ax = (double*)alloca(n*sizeof(double)*3); - double *ay = ax+n; double *az = ay+n; - - accel_conservative(n, x, y, z, m, Gx, Gy, Gz, - gravity_field, elastic_force, - n_bonds, bp, bk, br0, ax, ay, az); - - int has_damp = damping_force && (Bx != 0.0 || By != 0.0 || Bz != 0.0); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - if (has_damp) { - double ax_ = Bx*dt/(2.0*m[i]); - double ay_ = By*dt/(2.0*m[i]); - double az_ = Bz*dt/(2.0*m[i]); - vx[i] = (vx[i]*(1.0-ax_) + ax[i]*dt) / (1.0+ax_); - vy[i] = (vy[i]*(1.0-ay_) + ay[i]*dt) / (1.0+ay_); - vz[i] = (vz[i]*(1.0-az_) + az[i]*dt) / (1.0+az_); - } 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; - } -} - -/* ══════════════════════════════════════════════════════════ - * 显式欧拉法（与 main.c explicit_euler_step 一致） - * ══════════════════════════════════════════════════════════ */ -static void euler_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, - double dt) -{ - double *ax = (double*)alloca(n*sizeof(double)*3); - double *ay = ax+n; double *az = ay+n; - accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, ax, ay, az); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && 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; - } -} - -/* ══════════════════════════════════════════════════════════ - * 隐式欧拉法（与 main.c implicit_euler_step 完全一致） - * - * main.c 逻辑： - * 1. 用 v_next ≈ (v + G·dt)/(1 + γ·dt) 预测（只含重力+阻尼，不含弹簧） - * 2. 用 (x, v_next) 计算完整加速度 a_next - * 3. v += a_next·dt; x += v·dt - * ══════════════════════════════════════════════════════════ */ -static void implicit_euler_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, - double dt) -{ - double *vxn = (double*)alloca(n*sizeof(double)*3); - double *vyn = vxn+n; double *vzn = vyn+n; - - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) { - vxn[i] = vyn[i] = vzn[i] = 0.0; continue; - } - double gx = Bx / m[i], gy = By / m[i], gz = Bz / m[i]; - vxn[i] = (vx[i] + Gx*dt) / (1.0 + gx*dt); - vyn[i] = (vy[i] + Gy*dt) / (1.0 + gy*dt); - vzn[i] = (vz[i] + Gz*dt) / (1.0 + gz*dt); - } - - double *ax = (double*)alloca(n*sizeof(double)*3); - double *ay = ax+n; double *az = ay+n; - accel_full(n, x, y, z, vxn, vyn, vzn, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, ax, ay, az); - - for (int i = 0; i < n; i++) { - if (fixed[i*3] && 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; - } -} - -/* ══════════════════════════════════════════════════════════ - * 中点法（与 main.c midpoint_step 完全一致） - * - * main.c 逻辑： - * 1. a = accel(x, v) - * 2. xm = x + 0.5·v·dt; vm = v + 0.5·a·dt - * 3. x = x + vm·dt (位置更新用 vm，即中点速度) - * 4. am = accel(xm, vm) - * 5. v = v + am·dt - * ══════════════════════════════════════════════════════════ */ -static void midpoint_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, - double dt) -{ - /* Allocate in one block for cache locality */ - double *buf = (double*)alloca(n*sizeof(double)*9); - double *ax = buf; - double *ay = ax+n; double *az = ay+n; - double *xm = az+n; double *ym = xm+n; double *zm = ym+n; - double *vxm = zm+n; double *vym = vxm+n; double *vzm = vym+n; - - accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, ax, ay, az); - - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) { - xm[i]=x[i]; ym[i]=y[i]; zm[i]=z[i]; - vxm[i]=vym[i]=vzm[i]=0.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; - /* position updated with midpoint velocity (same as main.c) */ - x[i] = x[i] + vxm[i]*dt; - y[i] = y[i] + vym[i]*dt; - z[i] = z[i] + vzm[i]*dt; - } - - double *axm = (double*)alloca(n*sizeof(double)*3); - double *aym = axm+n; double *azm = aym+n; - accel_full(n, xm, ym, zm, vxm, vym, vzm, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, axm, aym, azm); - - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - vx[i] += axm[i]*dt; - vy[i] += aym[i]*dt; - vz[i] += azm[i]*dt; - } -} - -/* ── 驱动力（与 main.c apply_driving_force 一致）────────── */ -static void apply_driving( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - double t, int step, double dt, Drivers *drv) -{ - (void)n; - if (!drv || drv->n_drivers == 0) return; - const double TWO_PI = 2.0 * 3.14159265358979323846; - - for (int d = 0; d < drv->n_drivers; d++) { - int idx = drv->idx[d]; - double fx = drv->freq[d*3+0]; - double fy = drv->freq[d*3+1]; - double fz = drv->freq[d*3+2]; - - if (drv->has_period[d]) { - double mf = fabs(fx) > fabs(fy) ? fabs(fx) : fabs(fy); - if (fabs(fz) > mf) mf = fabs(fz); - int period_steps = 0; - if (mf > 1e-12) - period_steps = (int)(drv->ncycles[d] / mf / dt); - if (step > period_steps) { - x[idx] = drv->freeze[d*3+0]; - y[idx] = drv->freeze[d*3+1]; - z[idx] = drv->freeze[d*3+2]; - vx[idx] = vy[idx] = vz[idx] = 0.0; - continue; - } - - double px = drv->eq[d*3+0] + drv->amp[d*3+0]*cos(TWO_PI*fx*t + drv->phi[d*3+0]); - double py = drv->eq[d*3+1] + drv->amp[d*3+1]*cos(TWO_PI*fy*t + drv->phi[d*3+1]); - double pz = drv->eq[d*3+2] + drv->amp[d*3+2]*cos(TWO_PI*fz*t + drv->phi[d*3+2]); - if (step == period_steps) { - drv->freeze[d*3+0] = px; - drv->freeze[d*3+1] = py; - drv->freeze[d*3+2] = pz; - } - } - - x[idx] = drv->eq[d*3+0] + drv->amp[d*3+0]*cos(TWO_PI*fx*t + drv->phi[d*3+0]); - y[idx] = drv->eq[d*3+1] + drv->amp[d*3+1]*cos(TWO_PI*fy*t + drv->phi[d*3+1]); - z[idx] = drv->eq[d*3+2] + drv->amp[d*3+2]*cos(TWO_PI*fz*t + drv->phi[d*3+2]); - vx[idx] = -drv->amp[d*3+0]*TWO_PI*fx*sin(TWO_PI*fx*t + drv->phi[d*3+0]); - vy[idx] = -drv->amp[d*3+1]*TWO_PI*fy*sin(TWO_PI*fy*t + drv->phi[d*3+1]); - vz[idx] = -drv->amp[d*3+2]*TWO_PI*fz*sin(TWO_PI*fz*t + drv->phi[d*3+2]); - } -} - -/* ══════════════════════════════════════════════════════════ - * 导出函数：run_dynamics - * - * 与 main.c 的计算顺序完全一致： - * 1. leapfrog 初始化 v(-dt/2) - * 2. 初始驱动 t=0 - * 3. 预热循环（不记录） - * 4. 记录循环：drive → record → step → boundary → constraints - * - * 参数说明（所有数组均为 C-contiguous 行优先 float64/int32）： - * n_atoms 原子数 - * pos_init 初始位置 [n_atoms*3] x0,y0,z0, x1,y1,z1, ... - * vel_init 初始速度 [n_atoms*3] - * masses 质量 [n_atoms] - * fixed 自由度约束 [n_atoms*3] int32, 1=固定 - * n_bonds 键数 - * bond_pairs 键对 [n_bonds*2] int32, 0-based local index - * bond_k 刚度 [n_bonds] - * bond_r0 平衡键长 [n_bonds] - * box_a 盒子半边长 - * dt 时间步长 - * NT 总步数（含预热） - * NSTEP 抽帧间隔 - * warmup_steps 预热步数 - * method_id 0=euler 1=implicit 2=midpoint 3=leapfrog - * Gx/Gy/Gz 均匀重力场加速度分量 - * Bx/By/Bz 阻尼系数分量 - * gravity_field / elastic_force / damping_force 力开关 - * gravity_strength 原子间引力强度（暂未实现，留接口） - * n_drivers 驱动原子数 - * drv_idx 驱动原子局部索引 [n_drivers] int32 - * drv_amp 振幅 [n_drivers*3] - * drv_freq 频率 [n_drivers*3] - * drv_phi 初相（弧度）[n_drivers*3] - * drv_eq 平衡位置 [n_drivers*3] - * drv_ncycles 周期数 [n_drivers] 0=不限 - * drv_has_period [n_drivers] int32 - * n_frames 输出帧数（Python 预计算：(NT-warmup)/NSTEP 向上取整） - * out_x/y/z/vx/vy/vz 输出数组 [n_frames*n_atoms] 由 Python 预分配 - * progress_cb 进度回调（可为 NULL） - * - * 返回：0=成功，负数=错误 - * ══════════════════════════════════════════════════════════ */ -EXPORT int run_dynamics( - int n_atoms, - const double *pos_init, - const double *vel_init, - const double *masses, - const int *fixed, - int n_bonds, - const int *bond_pairs, - const double *bond_k, - const double *bond_r0, - double box_a, double dt, - int NT, int NSTEP, int warmup_steps, int method_id, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - double gravity_strength, - int n_drivers, - const int *drv_idx, - const double *drv_amp, - const double *drv_freq, - const double *drv_phi, - const double *drv_eq, - const double *drv_ncycles, - const int *drv_has_period, - int n_frames, - double *out_x, double *out_y, double *out_z, - double *out_vx, double *out_vy, double *out_vz, - void (*progress_cb)(int step, int total)) -{ - (void)gravity_strength; /* 原子间引力暂未实现 */ - - int n = n_atoms; - - /* ── 工作数组 ── */ - double *x = (double*)malloc(n*sizeof(double)); - double *y = (double*)malloc(n*sizeof(double)); - double *z = (double*)malloc(n*sizeof(double)); - double *vx = (double*)malloc(n*sizeof(double)); - double *vy = (double*)malloc(n*sizeof(double)); - double *vz = (double*)malloc(n*sizeof(double)); - if (!x||!y||!z||!vx||!vy||!vz) return -1; - - for (int i = 0; i < n; i++) { - x[i]=pos_init[i*3+0]; y[i]=pos_init[i*3+1]; z[i]=pos_init[i*3+2]; - vx[i]=vel_init[i*3+0]; vy[i]=vel_init[i*3+1]; vz[i]=vel_init[i*3+2]; - } - - /* ── 驱动结构 ── */ - Drivers drv; - drv.n_drivers = n_drivers; - drv.idx = drv_idx; - drv.amp = drv_amp; - drv.freq = drv_freq; - drv.phi = drv_phi; - drv.eq = drv_eq; - drv.ncycles = drv_ncycles; - drv.has_period = drv_has_period; - drv.freeze = NULL; - if (n_drivers > 0) { - drv.freeze = (double*)calloc(n_drivers*3, sizeof(double)); - if (!drv.freeze) { free(x);free(y);free(z);free(vx);free(vy);free(vz); return -2; } - } - - /* ── 内联步进宏 ── */ -#define DO_STEP() do { \ - switch (method_id) { \ - case 0: euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - case 1: implicit_euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - case 2: midpoint_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - default: leapfrog_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - } \ - apply_boundary_and_constraints(n,x,y,z,vx,vy,vz,fixed,pos_init,box_a); \ -} while(0) - - /* ── 蛙跳法：初始化 v(-dt/2) = v(0) - 0.5·a_c(0)·dt ── */ - if (method_id == 3) { - double *ax0 = (double*)alloca(n*sizeof(double)*3); - double *ay0 = ax0+n; double *az0 = ay0+n; - accel_conservative(n, x, y, z, masses, Gx, Gy, Gz, - gravity_field, elastic_force, - n_bonds, bond_pairs, bond_k, bond_r0, - ax0, ay0, az0); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - vx[i] -= 0.5*ax0[i]*dt; - vy[i] -= 0.5*ay0[i]*dt; - vz[i] -= 0.5*az0[i]*dt; - } - } - - /* ── 初始驱动 t=0（与 main.c 一致：leapfrog init 之后施加）── */ - if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, 0.0, 0, dt, &drv); - - /* ── 预热（不记录）── */ - for (int s = 0; s < warmup_steps; s++) { - double tw = (s + 1) * dt; - if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, tw, s, dt, &drv); - DO_STEP(); - } - - /* ── 记录循环 ── */ - int record_steps = NT - warmup_steps; - int prog_interval = record_steps / 100; - if (prog_interval < 1) prog_interval = 1; - int frame_idx = 0; - - for (int s = 0; s < record_steps; s++) { - if (progress_cb && s % prog_interval == 0 && s > 0) - progress_cb(s, record_steps); - - double t = (s + warmup_steps) * dt; - if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, t, s, dt, &drv); - - /* 抽帧记录（drive 之后，step 之前，与 main.c 一致）*/ - if (s % NSTEP == 0 && frame_idx < n_frames) { - int base = frame_idx * n; - for (int i = 0; i < n; i++) { - out_x [base+i] = x[i]; out_y [base+i] = y[i]; out_z [base+i] = z[i]; - out_vx[base+i] = vx[i]; out_vy[base+i] = vy[i]; out_vz[base+i] = vz[i]; - } - frame_idx++; - } - DO_STEP(); - } - -#undef DO_STEP - - free(x); free(y); free(z); - free(vx); free(vy); free(vz); - if (drv.freeze) free(drv.freeze); - return 0; -} diff --git a/engines/cpp/Makefile b/engines/cpp/Makefile deleted file mode 100644 index fa307d8..0000000 --- a/engines/cpp/Makefile +++ /dev/null @@ -1,43 +0,0 @@ -# engines/cpp/Makefile -# 编译 DLL（主程序通过 ctypes 直接调用） - -CXX = g++ -LIB_SRC = dynamics_lib.cpp - -UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) - -CXXFLAGS = -O3 -march=native -std=c++17 -Wall -Wextra -D_USE_MATH_DEFINES - -# Windows 下静态链接运行时，避免 libstdc++-6.dll / libgcc_s_seh-1.dll 版本冲突 -ifeq ($(UNAME_S),Windows) - STATIC_FLAGS = -static-libgcc -static-libstdc++ -else - STATIC_FLAGS = -endif - -ifeq ($(UNAME_S),Linux) - DLL_TARGET = build/dynamics_cpp.so - DLL_FLAGS = -shared -fPIC -else ifeq ($(UNAME_S),Darwin) - DLL_TARGET = build/dynamics_cpp.dylib - DLL_FLAGS = -dynamiclib -else - DLL_TARGET = build/dynamics_cpp.dll - DLL_FLAGS = -shared -endif - -.PHONY: all dll clean - -all: dll - -dll: $(DLL_TARGET) - -$(DLL_TARGET): $(LIB_SRC) | build - $(CXX) $(CXXFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) - @echo " === C++ DLL built: $@ ===" - -build: - mkdir -p build - -clean: - rm -rf build *.o diff --git a/engines/cpp/dynamics_lib.cpp b/engines/cpp/dynamics_lib.cpp deleted file mode 100644 index bf5bcdc..0000000 --- a/engines/cpp/dynamics_lib.cpp +++ /dev/null @@ -1,450 +0,0 @@ -/** - * engines/cpp/dynamics_lib.cpp - * ----------------------------- - * 纯计算 DLL（C++ 版）：无文件 I/O，所有数据由 Python 以 NumPy 数组传入。 - * 算法与 main.cpp / compute.py 保持完全一致。 - * - * 编译（Windows）: - * g++ -O3 -march=native -std=c++17 -shared -o build/dynamics_cpp.dll dynamics_lib.cpp - * 编译（Linux）: - * g++ -O3 -march=native -std=c++17 -shared -fPIC -o build/dynamics_cpp.so dynamics_lib.cpp - * 编译（macOS）: - * g++ -O3 -march=native -std=c++17 -dynamiclib -o build/dynamics_cpp.dylib dynamics_lib.cpp - */ - -#ifdef _WIN32 -# define EXPORT extern "C" __declspec(dllexport) -#else -# define EXPORT extern "C" __attribute__((visibility("default"))) -#endif - -#include -#include -#include -#include - -/* ── 驱动力结构体 ─────────────────────────────────────────── */ -struct Drivers { - int n_drivers = 0; - const int *idx = nullptr; - const double *amp = nullptr; - const double *freq = nullptr; - const double *phi = nullptr; - const double *eq = nullptr; - const double *ncycles = nullptr; - const int *has_period = nullptr; - std::vector freeze; /* [n_drivers*3] 冻结位置（period 结束时锁定）*/ -}; - -/* ── 加速度：保守力（弹簧键 + 均匀重力场）────────────────── */ -static void accel_conservative( - int n, const double *x, const double *y, const double *z, - const double *m, - double Gx, double Gy, double Gz, - int gravity_field, int elastic_force, - int n_bonds, const int *bond_pairs, - const double *bond_k, const double *bond_r0, - double *ax, double *ay, double *az) -{ - for (int i = 0; i < n; i++) { - ax[i] = gravity_field ? Gx : 0.0; - ay[i] = gravity_field ? Gy : 0.0; - az[i] = gravity_field ? Gz : 0.0; - } - if (!elastic_force || n_bonds == 0) return; - for (int b = 0; b < n_bonds; b++) { - int ii = bond_pairs[b*2]; - int jj = bond_pairs[b*2+1]; - double dx = x[jj]-x[ii], dy = y[jj]-y[ii], dz = z[jj]-z[ii]; - double dist = std::sqrt(dx*dx + dy*dy + dz*dz); - if (dist < 1e-12) continue; - double fac = bond_k[b] * (dist - bond_r0[b]) / dist; - double fx = fac*dx, fy = fac*dy, fz_b = fac*dz; - ax[ii] += fx/m[ii]; ay[ii] += fy/m[ii]; az[ii] += fz_b/m[ii]; - ax[jj] -= fx/m[jj]; ay[jj] -= fy/m[jj]; az[jj] -= fz_b/m[jj]; - } -} - -/* ── 完整加速度（含阻尼）────────────────────────────────── */ -static void accel_full( - int n, const double *x, const double *y, const double *z, - const double *vx, const double *vy, const double *vz, - const double *m, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bond_pairs, - const double *bond_k, const double *bond_r0, - double *ax, double *ay, double *az) -{ - accel_conservative(n, x, y, z, m, Gx, Gy, Gz, - gravity_field, elastic_force, - n_bonds, bond_pairs, bond_k, bond_r0, - ax, ay, az); - if (damping_force) { - for (int i = 0; i < n; i++) { - ax[i] -= Bx * vx[i] / m[i]; - ay[i] -= By * vy[i] / m[i]; - az[i] -= Bz * vz[i] / m[i]; - } - } -} - -/* ── 边界：反弹 ──────────────────────────────────────────── */ -static inline void _limit1(double &p, double &v, double lo, double hi) { - if (p > hi) { p = hi; v = -std::fabs(v); } - if (p < lo) { p = lo; v = std::fabs(v); } -} - -/* ── 边界：回绕 ──────────────────────────────────────────── */ -static inline void _wrap1(double &p, double lo, double hi) { - if (p > hi) p = lo; - if (p < lo) p = hi; -} - -/* ── 边界 + 固定约束 ────────────────────────────────────── */ -static void apply_boundary_and_constraints( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const int *fixed, const double *pos_init, double box_a) -{ - double lo = -box_a, hi = box_a; - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - _limit1(x[i], vx[i], lo, hi); - _limit1(y[i], vy[i], lo, hi); - _limit1(z[i], vz[i], lo, hi); - } - for (int i = 0; i < n; i++) { - _wrap1(x[i], lo, hi); - _wrap1(y[i], lo, hi); - _wrap1(z[i], lo, hi); - } - for (int i = 0; i < n; i++) { - if (fixed[i*3+0]) { x[i] = pos_init[i*3+0]; vx[i] = 0.0; } - if (fixed[i*3+1]) { y[i] = pos_init[i*3+1]; vy[i] = 0.0; } - if (fixed[i*3+2]) { z[i] = pos_init[i*3+2]; vz[i] = 0.0; } - } -} - -/* ══════════════════════════════════════════════════════════ - * 蛙跳法（与 main.cpp leapfrog_step 完全一致） - * ══════════════════════════════════════════════════════════ */ -static void leapfrog_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, double dt) -{ - std::vector buf(n * 3); - double *ax = buf.data(), *ay = ax+n, *az = ay+n; - accel_conservative(n, x, y, z, m, Gx, Gy, Gz, - gravity_field, elastic_force, - n_bonds, bp, bk, br0, ax, ay, az); - bool has_damp = damping_force && (Bx != 0.0 || By != 0.0 || Bz != 0.0); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - if (has_damp) { - double ax_ = Bx*dt/(2.0*m[i]); - double ay_ = By*dt/(2.0*m[i]); - double az_ = Bz*dt/(2.0*m[i]); - vx[i] = (vx[i]*(1.0-ax_) + ax[i]*dt) / (1.0+ax_); - vy[i] = (vy[i]*(1.0-ay_) + ay[i]*dt) / (1.0+ay_); - vz[i] = (vz[i]*(1.0-az_) + az[i]*dt) / (1.0+az_); - } 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; - } -} - -/* ══════════════════════════════════════════════════════════ - * 显式欧拉法 - * ══════════════════════════════════════════════════════════ */ -static void euler_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, double dt) -{ - std::vector buf(n * 3); - double *ax = buf.data(), *ay = ax+n, *az = ay+n; - accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, ax, ay, az); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && 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; - } -} - -/* ══════════════════════════════════════════════════════════ - * 隐式欧拉法（与 main.cpp implicit_euler_step 完全一致） - * ══════════════════════════════════════════════════════════ */ -static void implicit_euler_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, double dt) -{ - std::vector vbuf(n * 3), abuf(n * 3); - double *vxn = vbuf.data(), *vyn = vxn+n, *vzn = vyn+n; - double *ax = abuf.data(), *ay = ax+n, *az = ay+n; - - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) { - vxn[i] = vyn[i] = vzn[i] = 0.0; continue; - } - double gx = Bx/m[i], gy = By/m[i], gz = Bz/m[i]; - vxn[i] = (vx[i] + Gx*dt) / (1.0 + gx*dt); - vyn[i] = (vy[i] + Gy*dt) / (1.0 + gy*dt); - vzn[i] = (vz[i] + Gz*dt) / (1.0 + gz*dt); - } - accel_full(n, x, y, z, vxn, vyn, vzn, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, ax, ay, az); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && 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; - } -} - -/* ══════════════════════════════════════════════════════════ - * 中点法（与 main.cpp midpoint_step 完全一致） - * ══════════════════════════════════════════════════════════ */ -static void midpoint_step( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - const double *m, const int *fixed, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - int n_bonds, const int *bp, const double *bk, const double *br0, double dt) -{ - std::vector buf(n * 9); - double *ax = buf.data(); - double *ay = ax+n; double *az = ay+n; - double *xm = az+n; double *ym = xm+n; double *zm = ym+n; - double *vxm = zm+n; double *vym = vxm+n; double *vzm = vym+n; - - accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, ax, ay, az); - - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) { - xm[i]=x[i]; ym[i]=y[i]; zm[i]=z[i]; - vxm[i]=vym[i]=vzm[i]=0.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; - } - - std::vector abuf(n * 3); - double *axm = abuf.data(), *aym = axm+n, *azm = aym+n; - accel_full(n, xm, ym, zm, vxm, vym, vzm, m, Gx, Gy, Gz, Bx, By, Bz, - gravity_field, elastic_force, damping_force, - n_bonds, bp, bk, br0, axm, aym, azm); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - vx[i] += axm[i]*dt; - vy[i] += aym[i]*dt; - vz[i] += azm[i]*dt; - } -} - -/* ── 驱动力 ─────────────────────────────────────────────── */ -static void apply_driving( - int n, double *x, double *y, double *z, - double *vx, double *vy, double *vz, - double t, int step, double dt, Drivers &drv) -{ - (void)n; - if (drv.n_drivers == 0) return; - constexpr double TWO_PI = 2.0 * 3.14159265358979323846; - - for (int d = 0; d < drv.n_drivers; d++) { - int idx = drv.idx[d]; - double fx = drv.freq[d*3+0]; - double fy = drv.freq[d*3+1]; - double fz = drv.freq[d*3+2]; - - if (drv.has_period[d]) { - double mf = std::fabs(fx) > std::fabs(fy) ? std::fabs(fx) : std::fabs(fy); - if (std::fabs(fz) > mf) mf = std::fabs(fz); - int period_steps = 0; - if (mf > 1e-12) - period_steps = (int)(drv.ncycles[d] / mf / dt); - if (step > period_steps) { - x[idx] = drv.freeze[d*3+0]; - y[idx] = drv.freeze[d*3+1]; - z[idx] = drv.freeze[d*3+2]; - vx[idx] = vy[idx] = vz[idx] = 0.0; - continue; - } - double px = drv.eq[d*3+0] + drv.amp[d*3+0]*std::cos(TWO_PI*fx*t + drv.phi[d*3+0]); - double py = drv.eq[d*3+1] + drv.amp[d*3+1]*std::cos(TWO_PI*fy*t + drv.phi[d*3+1]); - double pz = drv.eq[d*3+2] + drv.amp[d*3+2]*std::cos(TWO_PI*fz*t + drv.phi[d*3+2]); - if (step == period_steps) { - drv.freeze[d*3+0] = px; - drv.freeze[d*3+1] = py; - drv.freeze[d*3+2] = pz; - } - } - x[idx] = drv.eq[d*3+0] + drv.amp[d*3+0]*std::cos(TWO_PI*fx*t + drv.phi[d*3+0]); - y[idx] = drv.eq[d*3+1] + drv.amp[d*3+1]*std::cos(TWO_PI*fy*t + drv.phi[d*3+1]); - z[idx] = drv.eq[d*3+2] + drv.amp[d*3+2]*std::cos(TWO_PI*fz*t + drv.phi[d*3+2]); - vx[idx] = -drv.amp[d*3+0]*TWO_PI*fx*std::sin(TWO_PI*fx*t + drv.phi[d*3+0]); - vy[idx] = -drv.amp[d*3+1]*TWO_PI*fy*std::sin(TWO_PI*fy*t + drv.phi[d*3+1]); - vz[idx] = -drv.amp[d*3+2]*TWO_PI*fz*std::sin(TWO_PI*fz*t + drv.phi[d*3+2]); - } -} - -/* ══════════════════════════════════════════════════════════ - * 导出函数：run_dynamics（接口与 C 版完全相同） - * ══════════════════════════════════════════════════════════ */ -EXPORT int run_dynamics( - int n_atoms, - const double *pos_init, - const double *vel_init, - const double *masses, - const int *fixed, - int n_bonds, - const int *bond_pairs, - const double *bond_k, - const double *bond_r0, - double box_a, double dt, - int NT, int NSTEP, int warmup_steps, int method_id, - double Gx, double Gy, double Gz, - double Bx, double By, double Bz, - int gravity_field, int elastic_force, int damping_force, - double gravity_strength, - int n_drivers, - const int *drv_idx, - const double *drv_amp, - const double *drv_freq, - const double *drv_phi, - const double *drv_eq, - const double *drv_ncycles, - const int *drv_has_period, - int n_frames, - double *out_x, double *out_y, double *out_z, - double *out_vx, double *out_vy, double *out_vz, - void (*progress_cb)(int step, int total)) -{ - (void)gravity_strength; - int n = n_atoms; - - std::vector xv(n), yv(n), zv(n); - std::vector vxv(n), vyv(n), vzv(n); - for (int i = 0; i < n; i++) { - xv[i]=pos_init[i*3+0]; yv[i]=pos_init[i*3+1]; zv[i]=pos_init[i*3+2]; - vxv[i]=vel_init[i*3+0]; vyv[i]=vel_init[i*3+1]; vzv[i]=vel_init[i*3+2]; - } - double *x=xv.data(), *y=yv.data(), *z=zv.data(); - double *vx=vxv.data(), *vy=vyv.data(), *vz=vzv.data(); - - Drivers drv; - drv.n_drivers = n_drivers; - drv.idx = drv_idx; - drv.amp = drv_amp; - drv.freq = drv_freq; - drv.phi = drv_phi; - drv.eq = drv_eq; - drv.ncycles = drv_ncycles; - drv.has_period = drv_has_period; - if (n_drivers > 0) - drv.freeze.assign(n_drivers * 3, 0.0); - -#define DO_STEP() do { \ - switch (method_id) { \ - case 0: euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - case 1: implicit_euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - case 2: midpoint_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - default: leapfrog_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \ - gravity_field,elastic_force,damping_force, \ - n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \ - } \ - apply_boundary_and_constraints(n,x,y,z,vx,vy,vz,fixed,pos_init,box_a); \ -} while(0) - - /* 蛙跳法：初始化 v(-dt/2) */ - if (method_id == 3) { - std::vector ibuf(n * 3); - double *ax0=ibuf.data(), *ay0=ax0+n, *az0=ay0+n; - accel_conservative(n, x, y, z, masses, Gx, Gy, Gz, - gravity_field, elastic_force, - n_bonds, bond_pairs, bond_k, bond_r0, - ax0, ay0, az0); - for (int i = 0; i < n; i++) { - if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue; - vx[i] -= 0.5*ax0[i]*dt; - vy[i] -= 0.5*ay0[i]*dt; - vz[i] -= 0.5*az0[i]*dt; - } - } - - /* 初始驱动 t=0 */ - if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, 0.0, 0, dt, drv); - - /* 预热 */ - for (int s = 0; s < warmup_steps; s++) { - double tw = (s + 1) * dt; - if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, tw, s, dt, drv); - DO_STEP(); - } - - /* 记录循环 */ - int record_steps = NT - warmup_steps; - int prog_interval = std::max(1, record_steps / 100); - int frame_idx = 0; - - for (int s = 0; s < record_steps; s++) { - if (progress_cb && s % prog_interval == 0 && s > 0) - progress_cb(s, record_steps); - - double t = (s + warmup_steps) * dt; - if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, t, s, dt, drv); - - if (s % NSTEP == 0 && frame_idx < n_frames) { - int base = frame_idx * n; - for (int i = 0; i < n; i++) { - out_x [base+i] = x[i]; out_y [base+i] = y[i]; out_z [base+i] = z[i]; - out_vx[base+i] = vx[i]; out_vy[base+i] = vy[i]; out_vz[base+i] = vz[i]; - } - frame_idx++; - } - DO_STEP(); - } - -#undef DO_STEP - return 0; -} diff --git a/engines/engine_dll.py b/engines/engine_dll.py index 3e3249c..3fdae70 100644 --- a/engines/engine_dll.py +++ b/engines/engine_dll.py @@ -12,9 +12,11 @@ Python ctypes 包装器：加载 C/C++/Fortran 动态链接库并调用 run_dyna # arrays: dict with keys x, y, z, vx, vy, vz shape=(n_frames, n_atoms) DLL 编译（C 版本）： - Windows: gcc -O3 -shared -o engines/c/build/dynamics_c.dll engines/c/dynamics_lib.c -lm - Linux: gcc -O3 -shared -fPIC -o engines/c/build/dynamics_c.so engines/c/dynamics_lib.c -lm - macOS: gcc -O3 -dynamiclib -o engines/c/build/dynamics_c.dylib engines/c/dynamics_lib.c -lm + Windows: gcc -O3 -shared -o engines/release/dynamics_c.dll engines/src/c/dynamics_lib.c -lm + Linux: gcc -O3 -shared -fPIC -o engines/release/dynamics_c.so engines/src/c/dynamics_lib.c -lm + macOS: gcc -O3 -dynamiclib -o engines/release/dynamics_c.dylib engines/src/c/dynamics_lib.c -lm +或用 make dll 一键编译： + cd engines/src/c && make dll """ import ctypes diff --git a/engines/fortran/Makefile b/engines/fortran/Makefile deleted file mode 100644 index 723c91b..0000000 --- a/engines/fortran/Makefile +++ /dev/null @@ -1,41 +0,0 @@ -# engines/fortran/Makefile -# 编译 DLL（主程序通过 ctypes 直接调用） - -FC = gfortran -FFLAGS = -O3 -march=native -Wall -Wextra -LIB_SRC = dynamics_lib.f90 - -UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) - -ifeq ($(UNAME_S),Windows) - STATIC_FLAGS = -static -else - STATIC_FLAGS = -endif - -ifeq ($(UNAME_S),Linux) - DLL_TARGET = build/dynamics_f90.so - DLL_FLAGS = -shared -fPIC -else ifeq ($(UNAME_S),Darwin) - DLL_TARGET = build/dynamics_f90.dylib - DLL_FLAGS = -dynamiclib -else - DLL_TARGET = build/dynamics_f90.dll - DLL_FLAGS = -shared -fPIC -endif - -.PHONY: all dll clean - -all: dll - -dll: $(DLL_TARGET) - -$(DLL_TARGET): $(LIB_SRC) | build - $(FC) $(FFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) - @echo " === Fortran DLL built: $@ ===" - -build: - mkdir -p build - -clean: - rm -rf build *.o *.mod diff --git a/engines/fortran/dynamics_lib.f90 b/engines/fortran/dynamics_lib.f90 deleted file mode 100644 index 565f0c4..0000000 --- a/engines/fortran/dynamics_lib.f90 +++ /dev/null @@ -1,483 +0,0 @@ -! engines/fortran/dynamics_lib.f90 -! --------------------------------- -! 纯计算 DLL（Fortran 版）：无文件 I/O，由 Python ctypes 调用。 -! 算法与 main.f90 / compute.py 完全一致。 -! 使用 iso_c_binding 导出 C 兼容接口。 -! -! 编译（Windows）: -! gfortran -O3 -march=native -shared -fPIC -o build/dynamics_f90.dll dynamics_lib.f90 -! 编译（Linux）: -! gfortran -O3 -march=native -shared -fPIC -o build/dynamics_f90.so dynamics_lib.f90 -! 编译（macOS）: -! gfortran -O3 -march=native -dynamiclib -o build/dynamics_f90.dylib dynamics_lib.f90 - -module dynamics_dll - use iso_c_binding, only: c_int, c_double, c_funptr, c_f_procpointer, c_associated - implicit none - private - - real(c_double), parameter :: TWO_PI = 2.0d0 * 3.14159265358979323846d0 - - public :: run_dynamics - -contains - -! ── 保守加速度 ─────────────────────────────────────────────── -subroutine accel_conservative(n, x, y, z, m, Gx, Gy, Gz, & - gravity_field, elastic_force, & - n_bonds, bond_pairs, bond_k, bond_r0, & - ax, ay, az) - integer, intent(in) :: n, gravity_field, elastic_force, n_bonds - real(c_double), intent(in) :: x(n), y(n), z(n), m(n) - real(c_double), intent(in) :: Gx, Gy, Gz - integer, intent(in) :: bond_pairs(2, n_bonds) - real(c_double), intent(in) :: bond_k(n_bonds), bond_r0(n_bonds) - real(c_double), intent(out) :: ax(n), ay(n), az(n) - - integer :: b, ii, jj - real(c_double) :: dx, dy, dz, dist, fac, fx, fy, fz_b - - if (gravity_field /= 0) then - ax = Gx; ay = Gy; az = Gz - else - ax = 0.0d0; ay = 0.0d0; az = 0.0d0 - end if - - if (elastic_force == 0 .or. n_bonds == 0) return - - do b = 1, n_bonds - ii = bond_pairs(1, b) + 1 ! 0-based → 1-based - jj = bond_pairs(2, b) + 1 - dx = x(jj)-x(ii); dy = y(jj)-y(ii); dz = z(jj)-z(ii) - dist = sqrt(dx*dx + dy*dy + dz*dz) - if (dist < 1.0d-12) cycle - fac = bond_k(b) * (dist - bond_r0(b)) / dist - fx = fac*dx; fy = fac*dy; fz_b = fac*dz - ax(ii) = ax(ii) + fx/m(ii); ay(ii) = ay(ii) + fy/m(ii); az(ii) = az(ii) + fz_b/m(ii) - ax(jj) = ax(jj) - fx/m(jj); ay(jj) = ay(jj) - fy/m(jj); az(jj) = az(jj) - fz_b/m(jj) - end do -end subroutine - -! ── 完整加速度（含阻尼）────────────────────────────────────── -subroutine accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az) - integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds - real(c_double), intent(in) :: x(n), y(n), z(n), vx(n), vy(n), vz(n), m(n) - real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz - integer, intent(in) :: bond_pairs(2, n_bonds) - real(c_double), intent(in) :: bond_k(n_bonds), bond_r0(n_bonds) - real(c_double), intent(out) :: ax(n), ay(n), az(n) - - integer :: i - - call accel_conservative(n, x, y, z, m, Gx, Gy, Gz, & - gravity_field, elastic_force, & - n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az) - if (damping_force /= 0) then - do i = 1, n - ax(i) = ax(i) - Bx*vx(i)/m(i) - ay(i) = ay(i) - By*vy(i)/m(i) - az(i) = az(i) - Bz*vz(i)/m(i) - end do - end if -end subroutine - -! ── 边界 + 固定约束 ────────────────────────────────────────── -subroutine apply_bc(n, x, y, z, vx, vy, vz, fixed, pos_init, box_a) - integer, intent(in) :: n - real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n) - integer, intent(in) :: fixed(3, n) - real(c_double), intent(in) :: pos_init(3, n), box_a - - integer :: i - real(c_double) :: lo, hi - - lo = -box_a; hi = box_a - - ! 反弹 - do i = 1, n - if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle - if (x(i)>hi) then; x(i)=hi; vx(i)=-abs(vx(i)); end if - if (x(i)hi) then; y(i)=hi; vy(i)=-abs(vy(i)); end if - if (y(i)hi) then; z(i)=hi; vz(i)=-abs(vz(i)); end if - if (z(i)hi) x(i)=lo; if (x(i)hi) y(i)=lo; if (y(i)hi) z(i)=lo; if (z(i) 0.0d0) - do i = 1, n - if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle - if (has_damp) then - ax_ = Bx*dt/(2.0d0*m(i)); ay_ = By*dt/(2.0d0*m(i)); az_ = Bz*dt/(2.0d0*m(i)) - vx(i) = (vx(i)*(1.0d0-ax_) + ax(i)*dt)/(1.0d0+ax_) - vy(i) = (vy(i)*(1.0d0-ay_) + ay(i)*dt)/(1.0d0+ay_) - vz(i) = (vz(i)*(1.0d0-az_) + az(i)*dt)/(1.0d0+az_) - 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 - -! ── 显式欧拉法 ─────────────────────────────────────────────── -subroutine euler_step(n, x, y, z, vx, vy, vz, m, fixed, & - Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, dt) - integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds - real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n) - real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds) - integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds) - real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt - - real(c_double) :: ax(n), ay(n), az(n) - integer :: i - - call accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az) - do i = 1, n - if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=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 - -! ── 隐式欧拉法 ─────────────────────────────────────────────── -subroutine implicit_euler_step(n, x, y, z, vx, vy, vz, m, fixed, & - Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, dt) - integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds - real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n) - real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds) - integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds) - real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt - - real(c_double) :: vxn(n), vyn(n), vzn(n), ax(n), ay(n), az(n) - real(c_double) :: gamma_x, gamma_y, gamma_z - integer :: i - - do i = 1, n - if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) then - vxn(i)=0.0d0; vyn(i)=0.0d0; vzn(i)=0.0d0; cycle - end if - gamma_x = Bx/m(i); gamma_y = By/m(i); gamma_z = Bz/m(i) - vxn(i) = (vx(i)+Gx*dt)/(1.0d0+gamma_x*dt) - vyn(i) = (vy(i)+Gy*dt)/(1.0d0+gamma_y*dt) - vzn(i) = (vz(i)+Gz*dt)/(1.0d0+gamma_z*dt) - end do - call accel_full(n, x, y, z, vxn, vyn, vzn, m, Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az) - do i = 1, n - if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=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 - -! ── 中点法 ─────────────────────────────────────────────────── -subroutine midpoint_step(n, x, y, z, vx, vy, vz, m, fixed, & - Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, dt) - integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds - real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n) - real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds) - integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds) - real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt - - real(c_double) :: ax(n), ay(n), az(n) - real(c_double) :: xm(n), ym(n), zm(n), vxm(n), vym(n), vzm(n) - real(c_double) :: axm(n), aym(n), azm(n) - integer :: i - - call accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az) - do i = 1, n - if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) then - xm(i)=x(i); ym(i)=y(i); zm(i)=z(i) - vxm(i)=0.0d0; vym(i)=0.0d0; vzm(i)=0.0d0; 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_full(n, xm, ym, zm, vxm, vym, vzm, m, Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, axm, aym, azm) - do i = 1, n - if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle - vx(i)=vx(i)+axm(i)*dt; vy(i)=vy(i)+aym(i)*dt; vz(i)=vz(i)+azm(i)*dt - end do -end subroutine - -! ── 驱动力施加 ─────────────────────────────────────────────── -subroutine apply_drive(n, x, y, z, vx, vy, vz, t, step, dt, & - nd, drv_idx, drv_amp, drv_freq, drv_phi, & - drv_eq, drv_ncycles, drv_has_period, freeze) - integer, intent(in) :: n, nd, step - real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n) - real(c_double), intent(in) :: t, dt - integer, intent(in) :: drv_idx(nd), drv_has_period(nd) - real(c_double), intent(in) :: drv_amp(3,nd), drv_freq(3,nd) - real(c_double), intent(in) :: drv_phi(3,nd), drv_eq(3,nd) - real(c_double), intent(in) :: drv_ncycles(nd) - real(c_double), intent(inout) :: freeze(3,nd) - - integer :: d, idx, ps - real(c_double) :: fx, fy, fz, mf, px, py, pz - - do d = 1, nd - idx = drv_idx(d) + 1 ! 0-based → 1-based - fx = drv_freq(1,d); fy = drv_freq(2,d); fz = drv_freq(3,d) - - if (drv_has_period(d) /= 0) then - mf = max(abs(fx), max(abs(fy), abs(fz))) - ps = 0 - if (mf > 1.0d-12) ps = int(drv_ncycles(d)/mf/dt) - if (step > ps) then - x(idx)=freeze(1,d); y(idx)=freeze(2,d); z(idx)=freeze(3,d) - vx(idx)=0.0d0; vy(idx)=0.0d0; vz(idx)=0.0d0 - cycle - end if - px = drv_eq(1,d)+drv_amp(1,d)*cos(TWO_PI*fx*t+drv_phi(1,d)) - py = drv_eq(2,d)+drv_amp(2,d)*cos(TWO_PI*fy*t+drv_phi(2,d)) - pz = drv_eq(3,d)+drv_amp(3,d)*cos(TWO_PI*fz*t+drv_phi(3,d)) - if (step == ps) then - freeze(1,d)=px; freeze(2,d)=py; freeze(3,d)=pz - end if - end if - x(idx) = drv_eq(1,d)+drv_amp(1,d)*cos(TWO_PI*fx*t+drv_phi(1,d)) - y(idx) = drv_eq(2,d)+drv_amp(2,d)*cos(TWO_PI*fy*t+drv_phi(2,d)) - z(idx) = drv_eq(3,d)+drv_amp(3,d)*cos(TWO_PI*fz*t+drv_phi(3,d)) - vx(idx) = -drv_amp(1,d)*TWO_PI*fx*sin(TWO_PI*fx*t+drv_phi(1,d)) - vy(idx) = -drv_amp(2,d)*TWO_PI*fy*sin(TWO_PI*fy*t+drv_phi(2,d)) - vz(idx) = -drv_amp(3,d)*TWO_PI*fz*sin(TWO_PI*fz*t+drv_phi(3,d)) - end do -end subroutine - -! ══════════════════════════════════════════════════════════════ -! 导出函数：run_dynamics（C 兼容接口，bind(C)） -! 接口与 C/C++ DLL 完全相同（扁平 C-contiguous 数组）。 -! ══════════════════════════════════════════════════════════════ -integer(c_int) function run_dynamics( & - n_atoms, pos_init, vel_init, masses, fixed, & - n_bonds, bond_pairs, bond_k, bond_r0, & - box_a, dt, NT, NSTEP, warmup_steps, method_id, & - Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, gravity_strength, & - n_drivers, drv_idx, drv_amp, drv_freq, drv_phi, drv_eq, & - drv_ncycles, drv_has_period, & - n_frames, out_x, out_y, out_z, out_vx, out_vy, out_vz, & - progress_cb) & - bind(C, name="run_dynamics") - - integer(c_int), value, intent(in) :: n_atoms, n_bonds, NT, NSTEP - integer(c_int), value, intent(in) :: warmup_steps, method_id - integer(c_int), value, intent(in) :: gravity_field, elastic_force, damping_force - integer(c_int), value, intent(in) :: n_drivers, n_frames - real(c_double), value, intent(in) :: box_a, dt - real(c_double), value, intent(in) :: Gx, Gy, Gz, Bx, By, Bz - real(c_double), value, intent(in) :: gravity_strength - - ! 扁平数组：Python 传入 C-contiguous int32/float64 - ! Fortran 以列优先解释，维度反转：(3,n) 对应 C 的 n×3 - real(c_double), intent(in) :: pos_init(3, n_atoms) - real(c_double), intent(in) :: vel_init(3, n_atoms) - real(c_double), intent(in) :: masses(n_atoms) - integer(c_int), intent(in) :: fixed(3, n_atoms) - integer(c_int), intent(in) :: bond_pairs(2, n_bonds) - real(c_double), intent(in) :: bond_k(n_bonds), bond_r0(n_bonds) - integer(c_int), intent(in) :: drv_idx(n_drivers) - real(c_double), intent(in) :: drv_amp(3, n_drivers) - real(c_double), intent(in) :: drv_freq(3, n_drivers) - real(c_double), intent(in) :: drv_phi(3, n_drivers) - real(c_double), intent(in) :: drv_eq(3, n_drivers) - real(c_double), intent(in) :: drv_ncycles(n_drivers) - integer(c_int), intent(in) :: drv_has_period(n_drivers) - - real(c_double), intent(out) :: out_x(n_atoms, n_frames) - real(c_double), intent(out) :: out_y(n_atoms, n_frames) - real(c_double), intent(out) :: out_z(n_atoms, n_frames) - real(c_double), intent(out) :: out_vx(n_atoms, n_frames) - real(c_double), intent(out) :: out_vy(n_atoms, n_frames) - real(c_double), intent(out) :: out_vz(n_atoms, n_frames) - - type(c_funptr), value, intent(in) :: progress_cb - - ! 进度回调接口 - abstract interface - subroutine cb_iface(step, total) bind(C) - use iso_c_binding - integer(c_int), value :: step, total - end subroutine - end interface - procedure(cb_iface), pointer :: cb_ptr - - integer :: n, s, frame_idx, record_steps, prog_interval, nd - real(c_double) :: t, tw - real(c_double), allocatable :: x(:), y(:), z(:), vx(:), vy(:), vz(:) - real(c_double), allocatable :: ax0(:), ay0(:), az0(:) - real(c_double), allocatable :: freeze(:,:) - logical :: has_cb - - n = n_atoms - nd = n_drivers - - allocate(x(n), y(n), z(n), vx(n), vy(n), vz(n)) - do s = 1, n - x(s) = pos_init(1,s); y(s) = pos_init(2,s); z(s) = pos_init(3,s) - vx(s) = vel_init(1,s); vy(s) = vel_init(2,s); vz(s) = vel_init(3,s) - end do - - allocate(freeze(3, max(nd,1))) - freeze = 0.0d0 - - has_cb = c_associated(progress_cb) - if (has_cb) call c_f_procpointer(progress_cb, cb_ptr) - - ! ── 蛙跳法：初始化 v(-dt/2) ───────────────────────────── - if (method_id == 3) then - allocate(ax0(n), ay0(n), az0(n)) - call accel_conservative(n, x, y, z, masses, Gx, Gy, Gz, & - gravity_field, elastic_force, & - n_bonds, bond_pairs, bond_k, bond_r0, ax0, ay0, az0) - do s = 1, n - if (fixed(1,s)/=0 .and. fixed(2,s)/=0 .and. fixed(3,s)/=0) cycle - vx(s)=vx(s)-0.5d0*ax0(s)*dt - vy(s)=vy(s)-0.5d0*ay0(s)*dt - vz(s)=vz(s)-0.5d0*az0(s)*dt - end do - deallocate(ax0, ay0, az0) - end if - - ! ── 初始驱动 t=0 ───────────────────────────────────────── - if (nd > 0) call apply_drive(n, x, y, z, vx, vy, vz, 0.0d0, 0, dt, & - nd, drv_idx, drv_amp, drv_freq, drv_phi, & - drv_eq, drv_ncycles, drv_has_period, freeze) - - ! ── 预热 ───────────────────────────────────────────────── - do s = 0, warmup_steps-1 - tw = (s+1)*dt - if (nd>0) call apply_drive(n, x, y, z, vx, vy, vz, tw, s, dt, & - nd, drv_idx, drv_amp, drv_freq, drv_phi, & - drv_eq, drv_ncycles, drv_has_period, freeze) - call do_step(n, x, y, z, vx, vy, vz, masses, fixed, & - Gx, Gy, Gz, Bx, By, Bz, & - gravity_field, elastic_force, damping_force, & - n_bonds, bond_pairs, bond_k, bond_r0, dt, method_id, & - pos_init, box_a) - end do - - ! ── 记录循环 ───────────────────────────────────────────── - record_steps = NT - warmup_steps - prog_interval = max(1, record_steps/100) - frame_idx = 0 - - do s = 0, record_steps-1 - if (has_cb .and. mod(s, prog_interval)==0 .and. s>0) call cb_ptr(s, record_steps) - - t = (s+warmup_steps)*dt - if (nd>0) call apply_drive(n, x, y, z, vx, vy, vz, t, s, dt, & - nd, drv_idx, drv_amp, drv_freq, drv_phi, & - drv_eq, drv_ncycles, drv_has_period, freeze) - - if (mod(s, NSTEP)==0 .and. frame_idx/dev/null || echo Windows) -# 目标文件名：统一使用 .exe 后缀（方便 Python 跨平台调用） -TARGET = build/dynamics_c.exe +# DLL 输出到 engines/release/ +DLL_DIR = ../../release -# DLL 目标（平台自动选择后缀） ifeq ($(UNAME_S),Linux) - DLL_TARGET = build/dynamics_c.so + DLL_TARGET = $(DLL_DIR)/dynamics_c.so DLL_FLAGS = -shared -fPIC else ifeq ($(UNAME_S),Darwin) - DLL_TARGET = build/dynamics_c.dylib + DLL_TARGET = $(DLL_DIR)/dynamics_c.dylib DLL_FLAGS = -dynamiclib else - DLL_TARGET = build/dynamics_c.dll + DLL_TARGET = $(DLL_DIR)/dynamics_c.dll 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 +$(DLL_TARGET): $(LIB_SRC) | $(DLL_DIR) $(CC) $(CFLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) $(LDFLAGS) @echo " === C DLL built: $@ ===" -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 +$(DLL_DIR): + mkdir -p $(DLL_DIR) clean: - rm -rf build *.o + rm -f $(DLL_TARGET) diff --git a/engines/src/c/main.c b/engines/src/c/main.c deleted file mode 100644 index 661c915..0000000 --- a/engines/src/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/src/cpp/Makefile b/engines/src/cpp/Makefile index 93cefc7..08e071f 100644 --- a/engines/src/cpp/Makefile +++ b/engines/src/cpp/Makefile @@ -1,49 +1,46 @@ -# engines/cpp/Makefile +# engines/src/cpp/Makefile +# 编译 DLL 到 engines/release/（主程序通过 ctypes 直接调用） CXX = g++ -SRCS = main.cpp LIB_SRC = dynamics_lib.cpp UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) CXXFLAGS = -O3 -march=native -std=c++17 -Wall -Wextra -D_USE_MATH_DEFINES -# Windows 下静态链接运行时，避免 libstdc++-6.dll / libgcc_s_seh-1.dll 版本冲突 +# Windows 下静态链接运行时 ifeq ($(UNAME_S),Windows) STATIC_FLAGS = -static-libgcc -static-libstdc++ else STATIC_FLAGS = endif -TARGET = build/dynamics_cpp.exe +# DLL 输出到 engines/release/ +DLL_DIR = ../../release ifeq ($(UNAME_S),Linux) - DLL_TARGET = build/dynamics_cpp.so + DLL_TARGET = $(DLL_DIR)/dynamics_cpp.so DLL_FLAGS = -shared -fPIC else ifeq ($(UNAME_S),Darwin) - DLL_TARGET = build/dynamics_cpp.dylib + DLL_TARGET = $(DLL_DIR)/dynamics_cpp.dylib DLL_FLAGS = -dynamiclib else - DLL_TARGET = build/dynamics_cpp.dll + DLL_TARGET = $(DLL_DIR)/dynamics_cpp.dll DLL_FLAGS = -shared 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 +$(DLL_TARGET): $(LIB_SRC) | $(DLL_DIR) $(CXX) $(CXXFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) @echo " === C++ DLL built: $@ ===" -build: - mkdir -p build +$(DLL_DIR): + mkdir -p $(DLL_DIR) clean: - rm -rf build *.o + rm -f $(DLL_TARGET) diff --git a/engines/src/cpp/main.cpp b/engines/src/cpp/main.cpp deleted file mode 100644 index 338677f..0000000 --- a/engines/src/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/src/fortran/Makefile b/engines/src/fortran/Makefile index 5c2b6a3..a505886 100644 --- a/engines/src/fortran/Makefile +++ b/engines/src/fortran/Makefile @@ -1,47 +1,44 @@ -# engines/fortran/Makefile +# engines/src/fortran/Makefile +# 编译 DLL 到 engines/release/（主程序通过 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 +# DLL 输出到 engines/release/ +DLL_DIR = ../../release ifeq ($(UNAME_S),Linux) - DLL_TARGET = build/dynamics_f90.so + DLL_TARGET = $(DLL_DIR)/dynamics_f90.so DLL_FLAGS = -shared -fPIC else ifeq ($(UNAME_S),Darwin) - DLL_TARGET = build/dynamics_f90.dylib + DLL_TARGET = $(DLL_DIR)/dynamics_f90.dylib DLL_FLAGS = -dynamiclib else - DLL_TARGET = build/dynamics_f90.dll + DLL_TARGET = $(DLL_DIR)/dynamics_f90.dll DLL_FLAGS = -shared -fPIC 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 +$(DLL_TARGET): $(LIB_SRC) | $(DLL_DIR) $(FC) $(FFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) @echo " === Fortran DLL built: $@ ===" -build: - mkdir -p build +$(DLL_DIR): + mkdir -p $(DLL_DIR) clean: - rm -rf build *.o *.mod + rm -f $(DLL_TARGET) diff --git a/engines/src/fortran/main.f90 b/engines/src/fortran/main.f90 deleted file mode 100644 index 2e90c05..0000000 --- a/engines/src/fortran/main.f90 +++ /dev/null @@ -1,1071 +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 - 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) - - ! 读取 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) - - 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) - 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 - character(len=*), intent(out) :: method - character(len=8096) :: buf - character(len=256) :: line - integer :: u, ios - - 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 - - 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) -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 - -! 从 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) - character(len=*), intent(in) :: outdir, method - integer, intent(in) :: n_frames, nat, NT, NSTEP, warmup, nb - integer, intent(in) :: gravity_field, elastic_force, damping_force, driving_force - 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) - - ! ── 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