From a16d2239a1050e48cd2ae20c36d97a655f5abdaa Mon Sep 17 00:00:00 2001 From: Ying-Li Niu <64801511@qq.com> Date: Wed, 17 Jun 2026 15:48:30 +0800 Subject: [PATCH] =?UTF-8?q?refactor:=20=E5=BC=95=E6=93=8E=E6=BA=90?= =?UTF-8?q?=E7=A0=81=E7=A7=BB=E8=87=B3=20src/=EF=BC=8C=E7=BC=96=E8=AF=91?= =?UTF-8?q?=E4=BA=A7=E7=89=A9=E7=BB=9F=E4=B8=80=E8=87=B3=20release/?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 源码: engines/c/ → engines/src/c/ - 源码: engines/cpp/ → engines/src/cpp/ - 源码: engines/fortran/ → engines/src/fortran/ - 编译产物: engines/release/ (静态编译) - dynamics_c.exe (516KB, C 引擎) - dynamics_cpp.exe (3.3MB, C++ 引擎) - dynamics_f90.exe (988KB, Fortran 引擎) - compute.py: engine_map 指向 release/ - compute.py: param.json 和校准缓存移至 release/ - engine_dll.py: DLL 搜索路径指向 release/ - release/ 目录纳入 git 版本管理，用户克隆后直接可用 --- compute.py | 18 +- engines/engine_dll.py | 2 +- engines/release/_calib_c.json | 1 + engines/release/c.json | 51 ++ engines/src/c/Makefile | 82 ++ engines/src/c/dynamics_lib.c | 554 +++++++++++++ engines/src/c/main.c | 1114 ++++++++++++++++++++++++++ engines/src/cpp/Makefile | 49 ++ engines/src/cpp/dynamics_lib.cpp | 450 +++++++++++ engines/src/cpp/main.cpp | 1025 ++++++++++++++++++++++++ engines/src/fortran/Makefile | 47 ++ engines/src/fortran/dynamics_lib.f90 | 483 +++++++++++ engines/src/fortran/main.f90 | 1071 +++++++++++++++++++++++++ 13 files changed, 4937 insertions(+), 10 deletions(-) create mode 100644 engines/release/_calib_c.json create mode 100644 engines/release/c.json create mode 100644 engines/src/c/Makefile create mode 100644 engines/src/c/dynamics_lib.c create mode 100644 engines/src/c/main.c create mode 100644 engines/src/cpp/Makefile create mode 100644 engines/src/cpp/dynamics_lib.cpp create mode 100644 engines/src/cpp/main.cpp create mode 100644 engines/src/fortran/Makefile create mode 100644 engines/src/fortran/dynamics_lib.f90 create mode 100644 engines/src/fortran/main.f90 diff --git a/compute.py b/compute.py index 1aae868..ceabc8c 100644 --- a/compute.py +++ b/compute.py @@ -1134,11 +1134,11 @@ def run_engine(engine, input_dir, output_dir, config): script_dir = os.path.dirname(os.path.abspath(__file__)) system = platform.system().lower() engine_map = { - "c": "engines/c/build/dynamics_c", - "cpp": "engines/cpp/build/dynamics_cpp", - "c++": "engines/cpp/build/dynamics_cpp", - "fortran": "engines/fortran/build/dynamics_f90", - "f90": "engines/fortran/build/dynamics_f90", + "c": "engines/release/dynamics_c", + "cpp": "engines/release/dynamics_cpp", + "c++": "engines/release/dynamics_cpp", + "fortran": "engines/release/dynamics_f90", + "f90": "engines/release/dynamics_f90", "python": None, # 特殊处理：用 sys.executable 调用 main.py } if engine not in engine_map: @@ -1208,7 +1208,7 @@ def run_engine(engine, input_dir, output_dir, config): "camera_center_y": float(config.get("camera_center_y", 0.0)), "camera_center_z": float(config.get("camera_center_z", 0.0)), } - param_path = os.path.join(script_dir, "engines", engine, "param.json") + param_path = os.path.join(script_dir, "engines", "release", f"{engine}.json") os.makedirs(os.path.dirname(param_path), exist_ok=True) with open(param_path, "w", encoding="utf-8") as f: json.dump(param_json, f, indent=2) @@ -1226,7 +1226,7 @@ def run_engine(engine, input_dir, output_dir, config): n_atoms_calib = len(ATOM_IDS) if ATOM_IDS is not None else 0 # 尝试读取缓存；当 n_atoms 相同且 NT 在 50% 范围内时视为有效 - _cache_path = os.path.join(script_dir, "engines", engine, "_calib_cache.json") + _cache_path = os.path.join(script_dir, "engines", "release", f"_calib_{engine}.json") _step_time = None try: with open(_cache_path, encoding="utf-8") as _cf: @@ -1241,10 +1241,10 @@ def run_engine(engine, input_dir, output_dir, config): if _step_time is None: _calib_param = dict(param_json) _calib_param["NT"] = _calib_nt - _calib_path = os.path.join(script_dir, "engines", engine, "_calib.json") + _calib_path = os.path.join(script_dir, "engines", "release", f"_calib_{engine}.json") with open(_calib_path, "w", encoding="utf-8") as _cf: json.dump(_calib_param, _cf, indent=2) - _calib_outdir = os.path.join(script_dir, "engines", engine, "_calib_out") + _calib_outdir = os.path.join(script_dir, "engines", "release", f"_calib_{engine}_out") os.makedirs(_calib_outdir, exist_ok=True) _ct0 = time.time() subprocess.run( diff --git a/engines/engine_dll.py b/engines/engine_dll.py index cb32d34..3e3249c 100644 --- a/engines/engine_dll.py +++ b/engines/engine_dll.py @@ -67,7 +67,7 @@ def _dll_candidates(engine: str) -> list[str]: ext = _SUFFIX.get(sys, ".so") eng_dir = _ENGINE_DIR.get(engine, engine) name = _DLL_NAME.get(engine, f"dynamics_{engine}") - base = os.path.join(_HERE, eng_dir, "build", name) + base = os.path.join(_HERE, "release", name) return [ base + ext, base + ".dll", diff --git a/engines/release/_calib_c.json b/engines/release/_calib_c.json new file mode 100644 index 0000000..480f0d2 --- /dev/null +++ b/engines/release/_calib_c.json @@ -0,0 +1 @@ +{"n_atoms": 120, "nt": 1000000, "step_time": 0.0009682028412818909} \ No newline at end of file diff --git a/engines/release/c.json b/engines/release/c.json new file mode 100644 index 0000000..8e2bdd8 --- /dev/null +++ b/engines/release/c.json @@ -0,0 +1,51 @@ +{ + "box_a": 300.0, + "NT": 1000000, + "DT": 0.001, + "NSTEP": 500, + "warmup_steps": 0, + "method": "leapfrog", + "G": [ + 0.0, + 0.0, + 0.0 + ], + "B": [ + 0.01, + 0.0, + 0.01 + ], + "gravity_field": 0, + "gravity_interaction": 0, + "elastic_force": 1, + "damping_force": 1, + "gravity_strength": 1.0, + "driving_force": 1, + "save_trajectory": 0, + "alpha": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "ball_radius": 0.5, + "ball_color": [ + 0.2, + 0.6, + 0.9 + ], + "box_color": [ + 0.8, + 0.8, + 0.85 + ], + "use_marker": 1, + "camera_distance": 120.0, + "camera_elevation": 0.0, + "camera_azimuth": 0.0, + "camera_center_x": 60.0, + "camera_center_y": 0.0, + "camera_center_z": 0.0 +} \ No newline at end of file diff --git a/engines/src/c/Makefile b/engines/src/c/Makefile new file mode 100644 index 0000000..03e99ce --- /dev/null +++ b/engines/src/c/Makefile @@ -0,0 +1,82 @@ +# engines/c/Makefile +# 跨平台编译：make → 本地系统编译 +# make linux → Linux 交叉编译（需 x86_64-linux-gnu-gcc） +# make windows → Windows 交叉编译（需 x86_64-w64-mingw32-gcc） +# make macos → macOS 交叉编译（需 osxcross 工具链） + +CC = gcc +CFLAGS = -O3 -march=native -Wall -Wextra +LDFLAGS = -lm +SRCS = main.c +LIB_SRC = dynamics_lib.c + +# 自动检测系统 +UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) + +# 目标文件名：统一使用 .exe 后缀（方便 Python 跨平台调用） +TARGET = build/dynamics_c.exe + +# DLL 目标（平台自动选择后缀） +ifeq ($(UNAME_S),Linux) + DLL_TARGET = build/dynamics_c.so + 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 linux windows macos + +all: $(TARGET) + +dll: $(DLL_TARGET) + +$(TARGET): $(SRCS) | build + $(CC) $(CFLAGS) -o $@ $(SRCS) $(LDFLAGS) + @echo " === C engine built: $@ ===" + +$(DLL_TARGET): $(LIB_SRC) | build + $(CC) $(CFLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) $(LDFLAGS) + @echo " === C DLL built: $@ ===" + +build: + mkdir -p build + +# ── 交叉编译 ───────────────────────────────── +# Linux → Linux (x86_64) +linux: CROSS_PREFIX = x86_64-linux-gnu- +linux: CC = $(CROSS_PREFIX)gcc +linux: CFLAGS = -O3 -march=x86-64 -Wall -Wextra +linux: $(SRCS) | build + $(CC) $(CFLAGS) -o build/dynamics_c_linux.exe $(SRCS) $(LDFLAGS) + @echo " === Linux binary: build/dynamics_c_linux.exe ===" + +# 任意平台 → Windows (x86_64) +# 需要安装 MinGW 交叉编译器： +# apt install mingw-w64 (Debian/Ubuntu) +# brew install mingw-w64 (macOS) +windows: CROSS_PREFIX = x86_64-w64-mingw32- +windows: CC = $(CROSS_PREFIX)gcc +windows: CFLAGS = -O3 -march=x86-64 -Wall -Wextra +windows: $(SRCS) | build + $(CC) $(CFLAGS) -o build/dynamics_c_win.exe $(SRCS) $(LDFLAGS) + @echo " === Windows binary: build/dynamics_c_win.exe ===" + +# 任意平台 → macOS (x86_64) +# 需要安装 osxcross 工具链 +macos: CROSS_PREFIX = x86_64-apple-darwin- +macos: CC = $(CROSS_PREFIX)gcc +macos: CFLAGS = -O3 -march=x86-64 -Wall -Wextra +macos: $(SRCS) | build + $(CC) $(CFLAGS) -o build/dynamics_c_mac.exe $(SRCS) $(LDFLAGS) + @echo " === macOS binary: build/dynamics_c_mac.exe ===" + +# ── 编译所有平台 ────────────────────────────── +all-platforms: linux windows macos + +clean: + rm -rf build *.o diff --git a/engines/src/c/dynamics_lib.c b/engines/src/c/dynamics_lib.c new file mode 100644 index 0000000..c83d984 --- /dev/null +++ b/engines/src/c/dynamics_lib.c @@ -0,0 +1,554 @@ +/** + * 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/src/c/main.c b/engines/src/c/main.c new file mode 100644 index 0000000..661c915 --- /dev/null +++ b/engines/src/c/main.c @@ -0,0 +1,1114 @@ +/** + * 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 new file mode 100644 index 0000000..93cefc7 --- /dev/null +++ b/engines/src/cpp/Makefile @@ -0,0 +1,49 @@ +# engines/cpp/Makefile + +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 版本冲突 +ifeq ($(UNAME_S),Windows) + STATIC_FLAGS = -static-libgcc -static-libstdc++ +else + STATIC_FLAGS = +endif + +TARGET = build/dynamics_cpp.exe + +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: $(TARGET) + +dll: $(DLL_TARGET) + +$(TARGET): $(SRCS) | build + $(CXX) $(CXXFLAGS) $(STATIC_FLAGS) -o $@ $(SRCS) + @echo " === C++ engine built: $@ ===" + +$(DLL_TARGET): $(LIB_SRC) | build + $(CXX) $(CXXFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) + @echo " === C++ DLL built: $@ ===" + +build: + mkdir -p build + +clean: + rm -rf build *.o diff --git a/engines/src/cpp/dynamics_lib.cpp b/engines/src/cpp/dynamics_lib.cpp new file mode 100644 index 0000000..bf5bcdc --- /dev/null +++ b/engines/src/cpp/dynamics_lib.cpp @@ -0,0 +1,450 @@ +/** + * 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/src/cpp/main.cpp b/engines/src/cpp/main.cpp new file mode 100644 index 0000000..338677f --- /dev/null +++ b/engines/src/cpp/main.cpp @@ -0,0 +1,1025 @@ +/** + * 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 new file mode 100644 index 0000000..5c2b6a3 --- /dev/null +++ b/engines/src/fortran/Makefile @@ -0,0 +1,47 @@ +# engines/fortran/Makefile + +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 +else + STATIC_FLAGS = +endif + +TARGET = build/dynamics_f90.exe + +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: $(TARGET) + +dll: $(DLL_TARGET) + +$(TARGET): $(SRCS) | build + $(FC) $(FFLAGS) $(STATIC_FLAGS) -o $@ $(SRCS) + @echo " === Fortran engine built: $@ ===" + +$(DLL_TARGET): $(LIB_SRC) | build + $(FC) $(FFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) + @echo " === Fortran DLL built: $@ ===" + +build: + mkdir -p build + +clean: + rm -rf build *.o *.mod diff --git a/engines/src/fortran/dynamics_lib.f90 b/engines/src/fortran/dynamics_lib.f90 new file mode 100644 index 0000000..565f0c4 --- /dev/null +++ b/engines/src/fortran/dynamics_lib.f90 @@ -0,0 +1,483 @@ +! 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/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