diff --git a/CMakeLists.txt b/CMakeLists.txt index a1354ad..bcf17e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,9 +53,11 @@ function(add_engine_target name lang src) add_executable(${name} ${src}) endif() - # 输出到 engines//build/（Python 调度器期望的位置） + # 从源码路径推导输出目录：engines/c/main.c → engines/c/build/ + # 这样 Python 的 run_engine() 能在固定路径找到可执行文件 + string(REGEX REPLACE "/[^/]+$" "/build" _out_subdir "${src}") set_target_properties(${name} PROPERTIES - RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/engines/${name}/build" + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/${_out_subdir}" ) # macOS 上禁止生成 .app 包（命令行工具不需要） set_target_properties(${name} PROPERTIES diff --git a/INSTALL.md b/INSTALL.md index 4d237f5..161a15e 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -461,7 +461,7 @@ py -3 examples/case01/run_dynamics.py ### C 引擎验证 ```bash -# 修改 parameters.yaml 中 engine: python → engine: c +# 修改 input.txt 中 engine: python → engine: c # 或用命令行指定 py -3 examples/case01/run_dynamics.py ``` @@ -543,7 +543,7 @@ pip install -i https://pypi.tuna.tsinghua.edu.cn/simple numpy pyyaml tqdm matplo ### Q: 如何使用 C 引擎？ ```yaml -# 1. 修改 examples/case01/input/parameters.yaml +# 1. 修改 examples/case01/input/input.txt engine: c # 2. 先编译 diff --git a/README.md b/README.md index db022f5..bfc87f9 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ 多语言动力学数值积分与轨迹可视化框架。 支持 Python、C、C++、Fortran 多种语言实现同一套物理引擎，共享同一套可视化管线 -（抽帧、绘图、3D 动画），直接在 `parameters.yaml` 中切换引擎即可对比性能。 +（抽帧、绘图、3D 动画），直接在 `input.txt` 中切换引擎即可对比性能。 ## 快速开始 @@ -46,7 +46,7 @@ Dynamics/ ├── examples/ │ └── case01/ │ ├── input/ -│ │ ├── parameters.yaml ← 唯一配置文件 +│ │ ├── input.txt ← 唯一配置文件（YAML 格式） │ │ ├── coord.txt │ │ ├── connection.txt │ │ └── bond.txt @@ -59,7 +59,7 @@ Dynamics/ └── .workbuddy/memory/ ← 工作记忆（AI 辅助开发用） ``` -## 配置文件 (`parameters.yaml`) +## 配置文件 (`input.txt` — YAML 格式) 所有控制集中在一个文件里： @@ -91,7 +91,7 @@ alpha: [0.0, 0.0, 0.2, 0.2, 0.0, 0.0] ### 一键切换 ```yaml -# parameters.yaml 中改一行 +# input.txt 中改一行 engine: c # 切换到 C 引擎 engine: python # 切回 Python ``` @@ -187,7 +187,7 @@ pip install numpy pyyaml tqdm matplotlib vispy PyQt5 ## 工作原理 ``` -parameters.yaml +input.txt │ ▼ dynamics.py ──→ run_from_config() 加载全局参数 diff --git a/build_release_zip.py b/build_release_zip.py index f7775b7..2bf396a 100644 --- a/build_release_zip.py +++ b/build_release_zip.py @@ -19,7 +19,7 @@ INCLUDE_FILES = [ "sample.py", "draw.py", "plot_trajectory.py", - "examples\\case01\\input\\parameters.yaml", + "examples\\case01\\input\\input.txt", "examples\\case01\\input\\coord.txt", ] diff --git a/compute.py b/compute.py index f297418..68d7fa8 100644 --- a/compute.py +++ b/compute.py @@ -64,6 +64,13 @@ box_color_r = None box_color_g = None box_color_b = None +# 力开关 +GRAVITY_FIELD = 1 # 均匀重力场 +GRAVITY_INTERACTION = 0 # 原子间万有引力 +ELASTIC_FORCE = 1 # 弹簧键力 +DAMPING_FORCE = 0 # 阻尼 +GRAVITY_STRENGTH = 1.0 + # 派生边界（根据 box_a 计算） X_MIN = None X_MAX = None @@ -156,9 +163,11 @@ def load_coord_file(coord_path): header = f.readline().strip().split() expected = ["n", "mass", "radius", "x", "y", "z", "vx", "vy", "vz"] legacy = expected + ["fixed_x", "fixed_y", "fixed_z"] - if header not in (expected, legacy): + legacy_new = expected + ["fix_x", "fix_y", "fix_z"] + if header not in (expected, legacy, legacy_new): raise ValueError( - f"坐标文件表头应为: {' '.join(expected)}，实际为: {' '.join(header)}") + f"坐标文件表头应为: {' '.join(expected)} 或加三列 fix_x fix_y fix_z，" + f"实际为: {' '.join(header)}") for line_no, line in enumerate(f, start=2): stripped = line.strip() @@ -167,7 +176,7 @@ def load_coord_file(coord_path): parts = stripped.split() if header == expected and len(parts) != 9: raise ValueError(f"{coord_path}:{line_no} 应有 9 列，实际为 {len(parts)} 列") - if header == legacy and len(parts) != 12: + if header in (legacy, legacy_new) and len(parts) != 12: raise ValueError(f"{coord_path}:{line_no} 应有 12 列，实际为 {len(parts)} 列") rows.append(parts) @@ -182,7 +191,7 @@ def load_coord_file(coord_path): velocities = np.array([[float(row[6]), float(row[7]), float(row[8])] for row in rows], dtype=np.float64) fixed = np.zeros((len(rows), 3), dtype=np.int64) - if header == legacy: + if header in (legacy, legacy_new): fixed = np.array([[int(row[9]), int(row[10]), int(row[11])] for row in rows], dtype=np.int64) @@ -369,6 +378,7 @@ def run_from_config(config, out_dir=None): global ball_radius, ball_color_r, ball_color_g, ball_color_b global box_color_r, box_color_g, box_color_b global warmup_steps, sample_start, sample_end + global GRAVITY_FIELD, GRAVITY_INTERACTION, ELASTIC_FORCE, DAMPING_FORCE, GRAVITY_STRENGTH box_a = float(config.get("box_a", 10.0)) raw_alpha = config.get("alpha", 0.2) @@ -427,6 +437,14 @@ def run_from_config(config, out_dir=None): box_color_g = float(config.get("box_color_g", 0.8)) box_color_b = float(config.get("box_color_b", 0.85)) + # 力开关 + global GRAVITY_FIELD, GRAVITY_INTERACTION, ELASTIC_FORCE, DAMPING_FORCE, GRAVITY_STRENGTH + GRAVITY_FIELD = int(config.get("gravity_field", 1)) + GRAVITY_INTERACTION = int(config.get("gravity_interaction", 0)) + ELASTIC_FORCE = int(config.get("elastic_force", 1)) + DAMPING_FORCE = int(config.get("damping_force", 0)) + GRAVITY_STRENGTH = float(config.get("gravity_strength", 1.0)) + print(f"[compute] 使用算法: {METHOD}") print(f"[compute] 已加载成键信息: {len(BOND_PAIRS)} 条键") if config.get("_skip_run", False): @@ -470,7 +488,8 @@ def run_from_config(config, out_dir=None): f.write(f" mass: {M}\n") f.write("=" * 50 + "\n") print(f"[compute] 日志已保存至: {log_path}") - print(f"[compute] 计算耗时: {elapsed:.3f} s") + record = len(traj_x) + print(f"[compute] 计算完成: {record} 步, {elapsed:.3f} s") return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz @@ -525,8 +544,14 @@ def run_engine(engine, input_dir, output_dir, config): "DT": float(config["DT"]), "NSTEP": int(config.get("NSTEP", 1)), "warmup_steps": int(config.get("warmup_steps", 0)), + "method": normalize_method_name(config.get("method", "leapfrog")), "G": [float(v) for v in G], "B": [float(v) for v in B], + "gravity_field": int(config.get("gravity_field", 1)), + "gravity_interaction": int(config.get("gravity_interaction", 0)), + "elastic_force": int(config.get("elastic_force", 1)), + "damping_force": int(config.get("damping_force", 0)), + "gravity_strength": float(config.get("gravity_strength", 1.0)), } param_path = os.path.join(script_dir, "engines", engine, "param.json") os.makedirs(os.path.dirname(param_path), exist_ok=True) @@ -540,24 +565,87 @@ def run_engine(engine, input_dir, output_dir, config): print(f"[compute] input: {input_dir}") print(f"[compute] output: {output_dir}") + # ── 预校准：跑 NT=1000 步测速 ──────────────────────────────── + total_steps = int(config["NT"]) - int(config.get("warmup_steps", 0)) + _calib_nt = min(1000, max(100, total_steps // 10)) + _calib_param = dict(param_json) + _calib_param["NT"] = _calib_nt + _calib_path = os.path.join(script_dir, "engines", engine, "_calib.json") + with open(_calib_path, "w", encoding="utf-8") as _cf: + json.dump(_calib_param, _cf, indent=2) + _ct0 = time.time() + subprocess.run( + [engine_path, os.path.abspath(input_dir), os.devnull, _calib_path], + capture_output=True, timeout=60) + os.remove(_calib_path) + _calib_elapsed = max(time.time() - _ct0, 0.001) + _overhead = _calib_elapsed * 0.15 + _step_time = max(_calib_elapsed - _overhead, 0.0001) / _calib_nt + _est_total = max(_calib_elapsed, _overhead + _step_time * total_steps) + t_start = time.time() t_start_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") - result = subprocess.run( - [engine_path, os.path.abspath(input_dir), os.path.abspath(output_dir), param_path], - capture_output=True, text=True, timeout=600) - t_end = time.time() - elapsed = t_end - t_start - t_end_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") - # 打印引擎输出 - if result.stdout: - for line in result.stdout.strip().split("\n"): - line = line.strip() - if line: - print(f" {line}") - if result.returncode != 0: - print(f"[compute] 引擎错误:\n{result.stderr}") - raise RuntimeError(f"引擎 {engine} 返回错误码 {result.returncode}") + _p = subprocess.Popen( + [engine_path, os.path.abspath(input_dir), os.path.abspath(output_dir), param_path], + stdout=subprocess.PIPE, stderr=subprocess.PIPE, + text=True, encoding='utf-8', errors='replace') + _engine_lines = [] + + try: + from tqdm import tqdm as _tqdm + _pbar = _tqdm(total=total_steps, desc=f"[compute] 引擎 {engine}", + unit="步", bar_format='{l_bar}{bar}| {n_fmt}/{total_fmt} [{elapsed}<{remaining}]') + except ImportError: + _pbar = None + + try: + while True: + _line = _p.stdout.readline() if _p.stdout else '' + if _line: + _line = _line.strip() + if _line: + _engine_lines.append(_line) + if _p.poll() is not None: + if _p.stdout: + for _r in _p.stdout: + _r = _r.strip() + if _r: + _engine_lines.append(_r) + break + if _pbar is not None and _est_total > 0: + _pbar.n = int(min((time.time() - t_start) / _est_total, 0.99) * total_steps) + _pbar.refresh() + time.sleep(0.2) + finally: + if _pbar is not None: + _pbar.n = total_steps + _pbar.refresh() + _pbar.close() + if _p.poll() is None: + _p.kill() + _p.wait(timeout=5) + + _rc = _p.returncode + t_end = time.time() + t_end_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + elapsed = t_end - t_start + + if _rc != 0: + _err = _p.stderr.read() if _p.stderr else '' + if _err: + print(f"[compute] 引擎错误:\n{_err}") + raise RuntimeError(f"引擎 {engine} 返回错误码 {_rc}") + + if _engine_lines: + _log_path = os.path.join(output_dir, "dynamics.log") + try: + with open(_log_path, "a", encoding="utf-8") as _lf: + _lf.write("\n--- 引擎输出 ---\n") + for _ln in _engine_lines: + _lf.write(_ln + "\n") + except OSError: + pass # 加载输出的 trajectory.txt traj_path = os.path.join(os.path.abspath(output_dir), "trajectory.txt") @@ -636,6 +724,11 @@ def save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, out_d "box_color_r": box_color_r, "box_color_g": box_color_g, "box_color_b": box_color_b, + "gravity_field": GRAVITY_FIELD, + "gravity_interaction": GRAVITY_INTERACTION, + "elastic_force": ELASTIC_FORCE, + "damping_force": DAMPING_FORCE, + "gravity_strength": GRAVITY_STRENGTH, } save_text_data(out_path, payload) print(f"[compute] 轨迹数据已保存至: {out_path}") @@ -752,16 +845,27 @@ def normalize_method_name(method): def compute_force(x, y, z, vx, vy, vz, m, g, b): """Compute total force from the current state. - Current model: - - gravity: F = m * g_vec - - linear drag: F_drag = -B_vec * v - - spring bonds: Hooke force based on bonded pair distance + Each force type is independently controlled by a global switch: + - GRAVITY_FIELD: uniform gravity F = m * g_vec + - DAMPING_FORCE: linear drag F_drag = -B_vec * v + - ELASTIC_FORCE: spring bonds based on bonded pair distance + - GRAVITY_INTERACTION: Newtonian gravity between atom pairs """ - fx = m * g[0] - b[0] * vx - fy = m * g[1] - b[1] * vy - fz = m * g[2] - b[2] * vz + fx = np.zeros_like(x) + fy = np.zeros_like(y) + fz = np.zeros_like(z) - if BOND_PAIRS is not None and len(BOND_PAIRS) > 0: + if GRAVITY_FIELD: + fx += m * g[0] + fy += m * g[1] + fz += m * g[2] + + if DAMPING_FORCE: + fx -= b[0] * vx + fy -= b[1] * vy + fz -= b[2] * vz + + if ELASTIC_FORCE and BOND_PAIRS is not None and len(BOND_PAIRS) > 0: for bond_idx, (idx_1, idx_2) in enumerate(BOND_PAIRS): dx = x[idx_2] - x[idx_1] dy = y[idx_2] - y[idx_1] @@ -787,6 +891,25 @@ def compute_force(x, y, z, vx, vy, vz, m, g, b): fy[idx_2] -= fy_bond fz[idx_2] -= fz_bond + if GRAVITY_INTERACTION: + n = len(m) + for i in range(n): + for j in range(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 <= 1e-12: + continue + mi = m[i]; mj = m[j] + f_mag = GRAVITY_STRENGTH * mi * mj / r2 + r = np.sqrt(r2) + fx_i = f_mag * dx / r + fy_i = f_mag * dy / r + fz_i = f_mag * dz / r + fx[i] += fx_i; fy[i] += fy_i; fz[i] += fz_i + fx[j] -= fx_i; fy[j] -= fy_i; fz[j] -= fz_i + return fx, fy, fz @@ -857,14 +980,13 @@ def Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b): y = y + vy_half * dt z = z + vz_half * dt - gamma_x = b[0] / m - gamma_y = b[1] / m - gamma_z = b[2] / m - vx_next = (vx_half + 0.5 * g[0] * dt) / (1.0 + 0.5 * gamma_x * dt) - vy_next = (vy_half + 0.5 * g[1] * dt) / (1.0 + 0.5 * gamma_y * dt) - vz_next = (vz_half + 0.5 * g[2] * dt) / (1.0 + 0.5 * gamma_z * dt) + # 显式预测器：用原始加速度外推半步（标准 Velocity-Verlet 预测步） + # v_pred = v_half + 0.5 * a(t)*dt，包含重力+阻尼+弹簧的所有贡献 + vx_pred = vx_half + 0.5 * ax * dt + vy_pred = vy_half + 0.5 * ay * dt + vz_pred = vz_half + 0.5 * az * dt ax_next, ay_next, az_next = compute_acceleration( - x, y, z, vx_next, vy_next, vz_next, m, g, b) + x, y, z, vx_pred, vy_pred, vz_pred, m, g, b) vx = vx_half + 0.5 * ax_next * dt vy = vy_half + 0.5 * ay_next * dt vz = vz_half + 0.5 * az_next * dt diff --git a/doc/index.html b/doc/index.html new file mode 100644 index 0000000..f27c577 --- /dev/null +++ b/doc/index.html @@ -0,0 +1,772 @@ + + + + + +Dynamics — 多语言蛙跳法实现对比分析 + + + + +
+

Dynamics — 多语言蛙跳法实现对比分析

+

Python · C · C++ · Fortran 四个引擎的蛙跳 (Velocity-Verlet) 算法与边界条件实现对比

+

最后更新：2026-05-20（力开关系统、万有引力、逐自由度约束）

+
+ + + +
+ + +
+

1. 概览

+

本框架实现了四套语言的计算引擎。各引擎共享相同的物理模型：

+
    +
  • 均匀重力场：Fg = m · g（常数矢量）—— 独立开关 gravity_field
    • +
  • 线性阻尼：Fd = −B · v —— 独立开关 damping_force
    • +
  • 弹簧键：胡克定律 F = k(r − r0) —— 独立开关 elastic_force
    • +
  • 万有引力（原子间）：F = G · mi · mj / r² —— 独立开关 gravity_interaction
    • +
  • 边界：硬壁反射，粒子限制在 [−a, a]3 内
    • +
  • 逐自由度固定约束：coord.txt 中 fix_x fix_y fix_z 列控制各自由度是否参与迭代
    • +
+

Python 引擎支持 4 种积分算法，C/C++/Fortran 引擎于 2026-05-20 重写，现支持全部 4 种算法，且边界条件与 Python 完全一致。经验证，四种引擎对于同一输入（case01, NT=1000, leapfrog, B=0）输出完全一致（至双精度浮点精度）。

+ +

分析范围

+
    +
  • 蛙跳法（Velocity-Verlet）积分器的四语言实现代码
    • +
  • 边界条件（硬壁反射：clamp 位置 + 反转速度）的实现
    • +
  • 2026-05-20 重写修复的 Bug 及版本差异
    • +
+ +

引擎简介

+ + + + + + +
引擎源文件语言标准算法支持
Pythoncompute.pyPython 3.8+显式欧拉、隐式欧拉、中点法、蛙跳法
Cengines/c/main.cC11同上（2026-05-20 新增 3 种）
C++engines/cpp/main.cppC++17同上（2026-05-20 新增 3 种）
Fortranengines/fortran/main.f90Fortran 90同上（2026-05-20 新增 3 种）
+
+ + +
+

2. Python 引擎（标准参考实现）

+

compute.py

+ +

2.1 蛙跳法积分步骤

+
+

Velocity-Verlet 三步法

+

Python 的 Leapfrog_Method 实现了标准 velocity-Verlet 格式。三步法对应于：

+
    +
  1. 速度半步推（第919-921行）：v(t+Δt/2) = v(t) + ½ · a(t) · Δt
    2. +
  2. 全推位置（第923-925行）：x(t+Δt) = x(t) + v(t+Δt/2) · Δt
    3. +
  3. 速度后半步（第930-937行）：v(t+Δt) = v(t+Δt/2) + ½ · a(t+Δt) · Δt
    4. +
+
+ +
Python917def Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b):
+918    ax, ay, az = compute_acceleration(x, y, z, vx, vy, vz, m, g, b)
+919    vx_half = vx + 0.5 * ax * dt
+920    vy_half = vy + 0.5 * ay * dt
+921    vz_half = vz + 0.5 * az * dt
+922
+923    x = x + vx_half * dt
+924    y = y + vy_half * dt
+925    z = z + vz_half * dt
+926
+927    gamma_x = b[0] / m
+928    gamma_y = b[1] / m
+929    gamma_z = b[2] / m
+930    vx_next = (vx_half + 0.5 * g[0] * dt) / (1.0 + 0.5 * gamma_x * dt)
+931    vy_next = (vy_half + 0.5 * g[1] * dt) / (1.0 + 0.5 * gamma_y * dt)
+932    vz_next = (vz_half + 0.5 * g[2] * dt) / (1.0 + 0.5 * gamma_z * dt)
+933    ax_next, ay_next, az_next = compute_acceleration(
+934        x, y, z, vx_next, vy_next, vz_next, m, g, b)
+935    vx = vx_half + 0.5 * ax_next * dt
+936    vy = vy_half + 0.5 * ay_next * dt
+937    vz = vz_half + 0.5 * az_next * dt
+938    return x, y, z, vx, vy, vz
+
+ +
+

⚠️ 阻尼的隐式处理（第927-932行）

+

Python 的蛙跳法中，阻尼项用隐式格式处理。gamma = B / m，然后 v_next = (v_half + ½·g·Δt) / (1 + ½·γ·Δt)。这个 v_next 是阻尼修正后的速度估计值，被传给第二次 compute_acceleration 作为入参。当 B = 0 时退化为显式 v_next = v_half + ½·g·Δt。

+

注意：vx_next 是独立变量，不覆盖 vx_half。最终速度 vx = vx_half + ½·a_next·dt 使用的是原始的半步速度。这是正确的做法。

+
+ +

2.2 边界条件

+
Python940def Limit_in_box(a, amin, amax, va):
+941    """限制物体在边界内，发生碰撞时反弹。"""
+942    over = a > amax
+943    under = a < amin
+944    a = np.where(over, amax, a)
+945    a = np.where(under, amin, a)
+946    va = np.where(over | under, -va, va)
+947    return a, va
+
+ +
Python949def apply_motion_update(x, y, z, vx, vy, vz, dt, m, g, b):
+957    elif METHOD == "leapfrog":
+958        x, y, z, vx, vy, vz = Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b)
+962    x, vx = Limit_in_box(x, X_MIN, X_MAX, vx)
+963    y, vy = Limit_in_box(y, Y_MIN, Y_MAX, vy)
+964    z, vz = Limit_in_box(z, Z_MIN, Z_MAX, vz)
+965    return x, y, z, vx, vy, vz
+
+ +
+

边界处理逻辑

+
    +
  • 积分方法调用（第957-958行）：只做速度/位置更新，不含边界
    • +
  • 边界后处理（第962-964行）：三个方向各调用一次 Limit_in_box
    • +
  • Limit_in_box：clamp 位置到边界 + 越界方向速度取反（弹性碰撞模型）
    • +
+
+
+ + +
+

3. 力开关与能量分析系统

+

2026-05-20 添加了独立的力开关系统，每种力可通过 input.txt 中的 0/1 开关独立控制。各引擎的受力计算均基于这些开关条件累加。

+ +

3.1 力开关配置

+ + + + + + +
开关名默认值控制的力物理公式
gravity_field1均匀重力场（常数矢量 g）F = m · g
gravity_interaction0原子间牛顿万有引力F = G · mi · mj / r²
elastic_force1弹簧键胡克力F = k(r − r0)
damping_force0线性阻尼F = −B · v
+ +

3.2 代码实现模式

+

Python compute_force() / C/C++ compute_acceleration() / Fortran accel() 统一采用以下模式：

+ +
伪代码// 1. 清零
+ax = 0; ay = 0; az = 0
+
+// 2. 按开关累加
+if (gravity_field):       ax += G[0]; ay += G[1]; az += G[2]
+if (damping_force):       ax -= B[0] * vx / m; ...
+if (elastic_force):       ax += bond_loop(...)
+if (gravity_interaction): ax += pair_gravity_loop(...)
+
+ +
+

✅ 开关关闭时零开销

+

每个力都包裹在 if (开关) 条件中。开关关闭时，对应力的计算代码完全跳过，不影响运行效率。

+
+ +

3.3 能量图联动

+

绘图时 dynamics.py 的 run_case() 根据力开关状态动态计算对应的势能分量：

+ + + + + + +
力开关势能公式关=0 时
gravity_field=1Ug = −m · G · rUg 归零
gravity_interaction=1Ugg = −∑i<j G · mi · mj / rijUgg 归零
elastic_force=1Us = ½ ∑ k(r − r0)²Us 归零
+ +
+

能量图线颜色

+

动能=蓝色、均匀重力势能=绿色、弹性势能=橙色、万有引力势能=紫色、总能量=红色虚线。

+
+
+ + +
+

4. C 引擎

+

engines/c/main.c

+

注意：2026-05-20 重写。原实现为 Symplectic Euler（非 Leapfrog），现已修正为 Velocity-Verlet。

+ +

3.1 蛙跳法积分步骤

+
C420/* ── 蛙跳法（Velocity-Verlet）── */
+421static void leapfrog_step(
+422    int n, double *x, double *y, double *z,
+423    double *vx, double *vy, double *vz,
+424    const double *m, const double G[3], const double B[3],
+425    const BondData *bonds, const int *fixed, double dt)
+426{
+427    double *ax = (double*)alloca(n * sizeof(double));
+428    double *ay = (double*)alloca(n * sizeof(double));
+429    double *az = (double*)alloca(n * sizeof(double));
+430    compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az);
+431
+432    /* 半推速度：v_half = v + 0.5*a*dt */
+433    for (int i = 0; i < n; i++) {
+434        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+435        vx[i] += ax[i] * dt * 0.5;
+436        vy[i] += ay[i] * dt * 0.5;
+437        vz[i] += az[i] * dt * 0.5;
+438    }
+439
+440    /* 全推位置（不含边界）*/
+441    for (int i = 0; i < n; i++) {
+442        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+443        x[i] += vx[i] * dt;   /* vx 此时是 v_half */
+444        y[i] += vy[i] * dt;
+445        z[i] += vz[i] * dt;
+446    }
+447
+448    /* 隐式阻尼处理：v_next = (v_half + 0.5*g*dt) / (1 + 0.5*gamma*dt)
+449       不覆盖 vx/vy/vz，用临时数组存入 */
+450    double *dmp_vx = (double*)alloca(n * sizeof(double));
+451    double *dmp_vy = (double*)alloca(n * sizeof(double));
+452    double *dmp_vz = (double*)alloca(n * sizeof(double));
+453    for (int i = 0; i < n; i++) {
+454        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+455        double gx = B[0] / m[i], gy = B[1] / m[i], gz = B[2] / m[i];
+456        dmp_vx[i] = (vx[i] + 0.5 * G[0] * dt) / (1.0 + 0.5 * gx * dt);
+457        dmp_vy[i] = (vy[i] + 0.5 * G[1] * dt) / (1.0 + 0.5 * gy * dt);
+458        dmp_vz[i] = (vz[i] + 0.5 * G[2] * dt) / (1.0 + 0.5 * gz * dt);
+459    }
+460
+461    /* 用新位置 + 阻尼处理后的速度重算加速度 */
+462    compute_acceleration(n, x, y, z, dmp_vx, dmp_vy, dmp_vz, m, G, B, bonds, ax, ay, az);
+463
+464    /* 速度后半步：v = v_half + 0.5*a_next*dt
+465       vx 仍为 v_half（未被覆盖）*/
+466    for (int i = 0; i < n; i++) {
+467        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+468        vx[i] += ax[i] * dt * 0.5;
+469        vy[i] += ay[i] * dt * 0.5;
+470        vz[i] += az[i] * dt * 0.5;
+471    }
+472}
+
+ +
+

✅ Velocity-Verlet 标准实现（与 Python 一致）

+

C 引擎的 leapfrog_step 与 Python 的 Leapfrog_Method 步骤完全相同：

+
    +
  1. 计算加速度 a(t)（第430行）
    2. +
  2. 速度半步推 v½（第432-438行）
    3. +
  3. 全推位置 x(t+Δt)（第440-446行）
    4. +
  4. 隐式阻尼 → 临时数组 dmp_vx（第448-459行）
    5. +
  5. 用新位置 + 阻尼速度重算加速度（第461-462行）
    6. +
  6. 速度后半步（第464-471行）
    7. +
+

关键：阻尼值存入 dmp_vx 而不覆盖 vx，与 Python 的 vx_next 独立变量做法一致。

+
+ +

3.2 边界条件

+
C309/* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
+310static void limit_in_box(double *pos, double *vel, double lo, double hi) {
+311    if (*pos > hi)  { *pos = hi;  *vel = -*vel; }
+312    if (*pos < lo)  { *pos = lo;  *vel = -*vel; }
+313}
+
+ +
C473/* ── 分发器：调用对应积分方法 + 边界条件 ── */
+488    } else if (strcmp(method, "leapfrog") == 0) {
+489        leapfrog_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
+490    }
+495    /* 边界条件（与 Python Limit_in_box 一致） */
+496    for (int i = 0; i < n; i++) {
+497        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+498        limit_in_box(&x[i],  &vx[i], -box_a, box_a);
+499        limit_in_box(&y[i],  &vy[i], -box_a, box_a);
+500        limit_in_box(&z[i],  &vz[i], -box_a, box_a);
+501    }
+502}
+
+ +
+

✅ 边界条件已修复

+

原 C 引擎仅有 clamp（只钳位置不反转速度），现已改为 limit_in_box（钳位置 + 反转速度），与 Python 一致。

+
+
+ + +
+

5. C++ 引擎

+

engines/cpp/main.cpp

+ +

4.1 蛙跳法积分步骤

+
C++360/* ── 蛙跳法（Velocity-Verlet）——与 Python Leapfrog_Method 一致 ── */
+361static void leapfrog_full_step(
+362    int n, double *x, double *y, double *z,
+363    double *vx, double *vy, double *vz,
+364    const double *m, const double G[3], const double B[3],
+365    const BondData &bonds, const int *fixed, double dt)
+366{
+367    // 第一次加速度
+368    std::vector ax(n), ay(n), az(n);
+369    compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax.data(), ay.data(), az.data());
+370
+371    // 半推速度：v_half = v + 0.5*a*dt  (存入 vx, vy, vz)
+372    for (int i = 0; i < n; i++) {
+373        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+374        vx[i] += ax[i] * dt * 0.5;
+375        vy[i] += ay[i] * dt * 0.5;
+376        vz[i] += az[i] * dt * 0.5;
+377    }
+378
+379    // 全推位置（不含边界，边界在外层统一处理）
+380    for (int i = 0; i < n; i++) {
+381        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+382        x[i] += vx[i] * dt;   // vx 此时是 v_half
+383        y[i] += vy[i] * dt;
+384        z[i] += vz[i] * dt;
+385    }
+386
+387    // 隐式阻尼处理得到 v_next（不覆盖 vx/vy/vz，Python 第929-931行 vx_next）
+388    std::vector dmp_vx(n), dmp_vy(n), dmp_vz(n);
+389    for (int i = 0; i < n; i++) {
+390        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+391        double gamma_x = B[0] / m[i];
+392        double gamma_y = B[1] / m[i];
+393        double gamma_z = B[2] / m[i];
+394        dmp_vx[i] = (vx[i] + 0.5 * G[0] * dt) / (1.0 + 0.5 * gamma_x * dt);
+395        dmp_vy[i] = (vy[i] + 0.5 * G[1] * dt) / (1.0 + 0.5 * gamma_y * dt);
+396        dmp_vz[i] = (vz[i] + 0.5 * G[2] * dt) / (1.0 + 0.5 * gamma_z * dt);
+397    }
+398
+399    // 用新位置 + 阻尼处理后的速度重算加速度
+400    compute_acceleration(n, x, y, z, dmp_vx.data(), dmp_vy.data(), dmp_vz.data(),
+401                         m, G, B, bonds, ax.data(), ay.data(), az.data());
+402
+403    // 速度后半步：v = v_half + 0.5*a_next*dt
+404    // vx 仍为 v_half（未被覆盖），直接加上 0.5*a_next*dt
+405    for (int i = 0; i < n; i++) {
+406        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+407        vx[i] += ax[i] * dt * 0.5;
+408        vy[i] += ay[i] * dt * 0.5;
+409        vz[i] += az[i] * dt * 0.5;
+410    }
+411}
+
+ +

4.2 边界条件

+
C++265/* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
+266static void limit_in_box(double &pos, double &vel, double lo, double hi) {
+267    if (pos > hi)  { pos = hi;  vel = -vel; }
+268    if (pos < lo)  { pos = lo;  vel = -vel; }
+269}
+
+ +
C++413/* ── 分发器：调用对应积分方法 + 边界条件 ── */
+429    } else if (method == "leapfrog") {
+430        leapfrog_full_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
+431    }
+435    // 边界条件（与 Python Limit_in_box 一致）
+436    for (int i = 0; i < n; i++) {
+437        if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
+438        limit_in_box(x[i],  vx[i], -box_a, box_a);
+439        limit_in_box(y[i],  vy[i], -box_a, box_a);
+440        limit_in_box(z[i],  vz[i], -box_a, box_a);
+441    }
+442}
+
+ +
+

✅ 与 Python 算法一致

+

C++ 引擎的 leapfrog_full_step 在重写后已与 Python 完全对齐，阻尼使用 std::vector<double> dmp_vx 临时数组而非原处覆盖。边界条件新增了速度反转。

+
+
+ + +
+

6. Fortran 引擎

+

engines/fortran/main.f90

+ +

5.1 蛙跳法积分步骤

+
Fortran498! ── 蛙跳法（Velocity-Verlet）──
+499subroutine leapfrog_full(n, x, y, z, vx, vy, vz, m, G, B, &
+500                         nb, bp, bk, br, fixed, dt)
+501    integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3)
+502    double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
+503    double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), dt
+504    double precision :: ax(n), ay(n), az(n)
+505    double precision :: dmp_vx(n), dmp_vy(n), dmp_vz(n)
+506    double precision :: gx, gy, gz
+507    integer :: i
+508
+509    call accel(n, x, y, z, vx, vy, vz, m, G, B, nb, bp, bk, br, ax, ay, az)
+510
+511    ! 速度半步推
+512    do i = 1, n
+513        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
+514        vx(i) = vx(i) + ax(i) * dt * 0.5d0
+515        vy(i) = vy(i) + ay(i) * dt * 0.5d0
+516        vz(i) = vz(i) + az(i) * dt * 0.5d0
+517    end do
+518
+519    ! 全推位置（不含边界）
+520    do i = 1, n
+521        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
+522        x(i) = x(i) + vx(i) * dt
+523        y(i) = y(i) + vy(i) * dt
+524        z(i) = z(i) + vz(i) * dt
+525    end do
+526
+527    ! 隐式阻尼处理（用临时数组 dmp_v，不覆盖 vx/vy/vz）
+528    do i = 1, n
+529        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then
+530            dmp_vx(i) = 0; dmp_vy(i) = 0; dmp_vz(i) = 0; cycle
+531        end if
+532        gx = B(1) / m(i); gy = B(2) / m(i); gz = B(3) / m(i)
+533        dmp_vx(i) = (vx(i) + 0.5d0 * G(1) * dt) / (1.0d0 + 0.5d0 * gx * dt)
+534        dmp_vy(i) = (vy(i) + 0.5d0 * G(2) * dt) / (1.0d0 + 0.5d0 * gy * dt)
+535        dmp_vz(i) = (vz(i) + 0.5d0 * G(3) * dt) / (1.0d0 + 0.5d0 * gz * dt)
+536    end do
+537
+538    ! 用新位置 + 阻尼速度重算加速度
+539    call accel(n, x, y, z, dmp_vx, dmp_vy, dmp_vz, m, G, B, nb, bp, bk, br, ax, ay, az)
+540
+541    ! 速度后半步：v = v_half + 0.5*a_next*dt（vx 仍为 v_half）
+542    do i = 1, n
+543        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
+544        vx(i) = vx(i) + ax(i) * dt * 0.5d0
+545        vy(i) = vy(i) + ay(i) * dt * 0.5d0
+546        vz(i) = vz(i) + az(i) * dt * 0.5d0
+547    end do
+548end subroutine leapfrog_full
+
+ +

5.2 边界条件

+
Fortran395! 边界条件：clamp 位置 + 速度反转
+396subroutine limit_in_box(pos, vel, lo, hi)
+397    double precision, intent(inout) :: pos, vel
+398    double precision, intent(in) :: lo, hi
+399    if (pos > hi) then
+400        pos = hi; vel = -vel
+401    else if (pos < lo) then
+402        pos = lo; vel = -vel
+403    end if
+404end subroutine limit_in_box
+
+ +
Fortran571    else if (trim(method) == 'leapfrog') then
+572        call leapfrog_full(n, x, y, z, vx, vy, vz, mm, G, B, &
+573                           nb, bp, bk, br, fixed, dt)
+574    end if
+579    ! 边界条件
+580    do i = 1, n
+581        if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
+582        call limit_in_box(x(i),  vx(i), -box_a, box_a)
+583        call limit_in_box(y(i),  vy(i), -box_a, box_a)
+584        call limit_in_box(z(i),  vz(i), -box_a, box_a)
+585    end do
+586end subroutine apply_step
+
+ +
+

✅ 与 Python 算法一致

+

Fortran 引擎的 leapfrog_full 在重写中修复了阻尼覆盖问题（使用 dmp_vx 临时数组），边界也改为速度反转。

+
+
+ + +
+

7. 版本变更记录

+

以下列出了 2026-05-20 三次迭代中所有修复和改进：

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
#问题描述影响引擎修复方式
1C 引擎蛙跳算法错误
原 leapfrog_step 实现的是 Symplectic Euler：
v += a·Δt; x += v·Δt（一次加速度，使用新速度推位置）。
这比 Velocity-Verlet 的截断误差大一阶（O(Δt) vs O(Δt²)），且只用了一次加速度计算。		C重写为 velocity-Verlet 三步法：
v½ → x → anew → v
2阻尼值覆盖速度变量
蛙跳法的隐式阻尼处理直接写回 vx[i]，导致后面的速度后半步 vx[i] += ½·a·dt 在覆盖后的值上叠加。等效于重力加速度被多算了一次，引入 ½·g·Δt 的系统偏差（Fortran）或 g·Δt 的偏差（C）。		C, Fortran用临时数组 dmp_vx 存储阻尼估计值，不覆盖原始速度变量
3bond.txt 表头未跳过
C/C++ 引擎的 read_bonds 用 while (fb >> bn >> bk >> br) 读取 bond.txt，但第一行是表头 bond_name k rest_length。解析 k 时尝试将字符串 "k" 转 double 失败，使用默认值 1.0（应为 100.0），弹簧劲度系数差 100 倍，完全改变了轨迹。Python 和 Fortran 正确跳过了表头。		C, C++在数据读取循环前用 std::getline(fb, header)（C++）或 fgets(header, ...)（C）跳过表头行
4C++ JSON 输出精度不足
C++ 的 JSON 输出用默认 << 运算符，只有 6 位有效数字，而 C 的 %.15g 和 Fortran 的 g0 都用全精度。导致 C++ 与其他引擎的输出到第 6 位之后截断。		C++添加 std::setprecision(15) 到输出流
5C 引擎 record_steps 忽略 warmup_steps
C 引擎的轨迹缓冲区大小写死为 params.NT 而非 params.NT - params.warmup_steps。Python/C++/Fortran 正确。		C改为 int record_steps = params.NT - params.warmup_steps;
6力开关系统
新增 4 个独立力开关 gravity_field/gravity_interaction/elastic_force/damping_force，通过 input.txt 的 0/1 控制各力是否参与计算。所有引擎统一实现：清零 → 按条件累加。		四引擎Python: 重写 compute_force()；C/C++: 重写 compute_acceleration()；Fortran: 重写 accel()
7万有引力引擎
原子间牛顿万有引力 F = G · mi · mj / r²，O(N²) 双循环对计算，gravity_interaction=0 时零开销。		四引擎所有引擎的受力计算函数中添加多原子对循环
8逐自由度固定约束
coord.txt 末尾增加 fix_x fix_y fix_z 三列，每列 0/1 控制该方向是否锁定。锁定后位置回复初始值、速度置零。		四引擎C/C++/Fortran: apply_step() 新增 pos_0 参数，边界条件后添加约束；Python: 已有 apply_fixed_constraints() 完美支持
9能量图随力开关联动
绘图时根据力开关状态动态计算势能分量。仅 gravity_interaction=1 时计算万有引力势能并以紫色曲线绘制。		Pythondynamics.py 能量计算添加 ug_grav 和开关判断
+
+ + +
+

8. 边界条件深度分析

+ +

8.1 硬壁反射实现

+

四引擎的边界条件在重写后完全一致：

+ + + + + + +
引擎函数名钳位置反转速度调用时机
PythonLimit_in_box✅✅蛙跳步之后（后处理）
Climit_in_box✅✅积分步之后（apply_step）
C++limit_in_box✅✅积分步之后（apply_step）
Fortranlimit_in_box✅✅积分步之后（apply_step）
+ +
+

算法合理性

+

速度反转是必须的。 硬壁反射的物理模型要求速度在碰撞时反向以保持运动方向的一致性。没有速度反转时，粒子到达边界后会被持续钳制在边界上（"黏壁"现象），这是功能性缺陷。

+
+ +

8.2 逐自由度固定约束

+

2026-05-20 新增了逐自由度固定约束。在 coord.txt 末尾添加 fix_x fix_y fix_z 三列：

+ +
coord.txt 格式n mass radius x y z vx vy vz fix_x fix_y fix_z
+1 1 0.28 -1 0 0 0 0 0 1 0 0   ← 仅固定 x 方向
+2 1 0.28 1 0 1 0 0 0 0 1 1   ← 固定 y 和 z 方向
+
+ +

Python 的 apply_fixed_constraints() 使用 NumPy 逐元素操作：

+
Pythonfixed = ATOM_FIXED != 0           # (n_atoms, 3) 布尔矩阵
+positions = np.where(fixed, ATOM_POSITIONS, positions)   # 逐度约束
+velocities = np.where(fixed, 0.0, velocities)
+
+ +

C/C++/Fortran 引擎在 apply_step() 的边界条件之后添加了相同的约束逻辑：

+
C/C++/Fortran 伪代码// 边界条件后，逐自由度固定约束
+for i = 1..n:
+    if fixed_x[i]:  x[i] = x0[i];  vx[i] = 0
+    if fixed_y[i]:  y[i] = y0[i];  vy[i] = 0
+    if fixed_z[i]:  z[i] = z0[i];  vz[i] = 0
+
+ +
+

✅ 约束在积分步与边界条件之后，保证正确性

+

约束在积分器完成位置/速度更新、边界钳位之后执行，确保固定自由度不参与任何物理演化。

+
+ +

8.3 能量局限性

+
+

⚠️ 硬壁反射与 Symplectic 积分器的不兼容

+

Velocity-Verlet 是 symplectic 积分器，在平滑势场中能量几乎精确守恒（前 20 步漂移仅 -0.000078%）。但硬壁反射是一个非光滑势（无限高势垒），每个碰撞步的工作流程如下：

+
    +
  1. 蛙跳步假定粒子在光滑势场中自由运动了 Δt 时间
    2. +
  2. 后处理发现粒子越界，将位置钳回边界 + 反转速度
    3. +
  3. 这个"钳回"动作改变了重力势能（因为 z 位置被改变）但没有相应调整动能
    4. +
  4. 每次碰撞产生约 0.03~0.5 J 的能量变化
    5. +
+

在典型 case01 参数下（100,000 步，83 次碰底），累积能量漂移约 36%。这是硬壁边界的固有限制，不是代码 Bug。

+
+ +
+

✅ 修复建议：软壁势

+

用光滑的排斥势替换硬壁 clamp 可保持 symplectic 守恒性：

+
    +
  • WCA 势（Weeks-Chandler-Andersen）：截断平移的 Lennard-Jones
    • +
  • 指数排斥：U(z) = A · exp(−|z|/z0)
    • +
  • 半壁谐振子：|z| > a 时受到线性回复力 F = −k(z − sign(z)·a)
    • +
+
+
+ + +
+

9. 总结

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
维度PythonCC++Fortran
积分器类型Velocity-VerletVelocity-Verlet
原为 Symplectic Euler		Velocity-VerletVelocity-Verlet
算法数量4 种4 种
新增 3 种		4 种
新增 3 种		4 种
新增 3 种
阻尼处理隐式 + 显式混合同 Python同 Python同 Python
边界：钳位置✅✅✅✅
边界：反转速度✅✅
原缺失		✅
原缺失		✅
原缺失
输出精度（JSON）双精度%.15gsetprecision(15)
原默认 6 位		g0
算法选择方式YAML configparam.jsonparam.jsonparam.json
引擎一致性四个引擎输出完全一致（双精度验证通过）
力开关系统✅✅✅✅
万有引力✅✅✅✅
逐自由度约束✅
np.where		✅
pos_0 参数		✅
pos_0 参数		✅
pos_0 参数
+
+ +
+ + + + + diff --git a/dynamics.py b/dynamics.py index 4c6a83e..0cd3504 100644 --- a/dynamics.py +++ b/dynamics.py @@ -4,9 +4,9 @@ dynamics.py 统一入口：读取 YAML 配置文件 → 运行模拟 → 抽帧 → 绘图（可选） 用法： - python dynamics.py # 使用 input/parameters.yaml - python dynamics.py input/parameters.yaml # 指定配置文件 - python dynamics.py --config input/parameters.yaml --no-plot + python dynamics.py # 使用 input/input.txt + python dynamics.py input/input.txt # 指定配置文件（YAML 格式） + python dynamics.py --config input/input.txt --no-plot """ import os @@ -115,6 +115,14 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", config = load_yaml_config(config_path) print(f"[run] 已加载配置: {config_path}") + # ── T_total → NT 自动转换 ────────────────────── + if "T_total" in config and "NT" not in config: + dt = float(config["DT"]) + config["NT"] = int(float(config["T_total"]) / dt) + print(f"[run] T_total={config['T_total']} → NT={config['NT']} (DT={dt})") + elif "T_total" in config and "NT" in config: + print(f"[run] 同时指定了 T_total 和 NT，使用 NT={config['NT']}") + # 显示步骤控制信息 steps_info = {k: config.get(k, 1) for k in ["step_simulate", "step_sample", "step_plot", "step_animation"]} step_flags = ", ".join(f"{k}={v}" for k, v in steps_info.items()) @@ -136,27 +144,42 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", disp_path = os.path.join(output_dir_abs, "display.txt") # ── 自动缓存检测 ─────────────────────────────────────── - # 若 output/ 中 trajectory.txt 和 display.txt 均已存在， - # 自动跳过模拟和抽帧，直接使用已有结果。 - # 如需强制重新计算，删除 output/ 目录或设 step_simulate: 1 即可。 - output_exists = ( - os.path.isdir(output_dir_abs) - and os.path.exists(traj_path) - and os.path.exists(disp_path) - ) - if output_exists: - if config.get("step_simulate", 1): - print(f"[run] 检测到已有输出（{traj_path}），自动跳过模拟与抽帧，直接进入后续步骤") - config["step_simulate"] = 0 - config["step_sample"] = 0 - else: - print(f"[run] 已有输出，步骤已被跳过") + # force_calc=1: 强制重新计算，忽略缓存 + # force_calc=0: 尊重 step_simulate 设置，不自动覆盖 + force_calc = int(config.get("force_calc", 0)) + if force_calc: + print(f"[run] force_calc=1，跳过缓存，强制重新计算") + config["step_simulate"] = 1 + config["step_sample"] = 1 + elif config.get("step_simulate", 1): + # step_simulate=1 且 force_calc=0 → 按用户要求执行计算 + # 但检测一下参数是否已变更（NT），如果变了则自动更新 step_sample + if os.path.isdir(output_dir_abs) and os.path.exists(traj_path) and os.path.exists(disp_path): + cached_nt = None + try: + with open(traj_path, 'rb') as _f: + _f.seek(-4096, 2) + _tail = _f.read().decode('utf-8', errors='replace') + import re as _re + _m = _re.search(r'"NT":\s*(\d+)', _tail) + if _m: + cached_nt = int(_m.group(1)) + except Exception: + pass + config_nt = int(config.get("NT", 0)) + if cached_nt is not None and cached_nt != config_nt: + print(f"[run] 参数已变更（缓存 NT={cached_nt}，配置 NT={config_nt}），" + f"将重新计算") + config["step_sample"] = 1 + else: + # 参数一致，按 step_simulate=1 执行，step_sample 由用户设置决定 + pass else: - # 目录存在但文件不全 → 强制重新计算 - if os.path.isdir(output_dir_abs): - print(f"[run] output/ 目录存在但文件不完整，将重新计算") + # step_simulate=0 → 检测缓存是否存在并提示 + if os.path.isdir(output_dir_abs) and os.path.exists(traj_path) and os.path.exists(disp_path): + print(f"[run] 已有输出，步骤已被跳过") else: - print(f"[run] output/ 目录不存在，将执行完整流程") + print(f"[run] 没有可用的缓存输出，但 step_simulate=0，将跳过模拟") # 2. 运行物理模拟 → output/trajectory.txt if config.get("step_simulate", 1): @@ -164,9 +187,12 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", total_steps = config["NT"] record_steps = total_steps - (config.get("warmup_steps") or 0) print(f"[run] 开始计算 总步数={total_steps} 记录步数={record_steps} DT={config['DT']}") + + import time as _time + _t0 = _time.time() + if engine == "python": traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = compute.run_from_config(config, str(runtime_base)) - print(f"[run] 计算完成，记录 {record_steps} 步") compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base)) else: # 外部引擎：先加载配置到全局变量，再运行引擎，再用 save_trajectory_txt 补全 metadata @@ -178,6 +204,9 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = compute.run_engine( engine, input_dir_abs, output_dir_abs, config) compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base)) + + _elapsed = _time.time() - _t0 + print(f"[run] 引擎: {engine} 计算完成: {record_steps} 步 {_elapsed:.3f} s") else: print("[run] 步骤 [模拟] 已跳过，直接加载已有轨迹") if not os.path.exists(traj_path): @@ -282,6 +311,11 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", "box_color_r": float(data["box_color_r"]), "box_color_g": float(data["box_color_g"]), "box_color_b": float(data["box_color_b"]), + "gravity_field": int(data.get("gravity_field", 1)), + "gravity_interaction": int(data.get("gravity_interaction", 0)), + "elastic_force": int(data.get("elastic_force", 1)), + "damping_force": int(data.get("damping_force", 0)), + "gravity_strength": float(data.get("gravity_strength", 1.0)), } save_display_txt(disp_data, str(runtime_base)) print(f"[run] 抽帧完成: {sample_end - sample_start} 步 -> {n_frames} 帧") @@ -328,21 +362,28 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", # ── 能量计算 ───────────────────────────────────── masses = np.array(data["atom_masses"]) # (n_atoms,) G_vec = np.array(data.get("G", [0.0, 0.0, -9.8])) # [gx, gy, gz] + gravity_field_enabled = int(data.get("gravity_field", 1)) + gravity_interaction_enabled = int(data.get("gravity_interaction", 0)) + gravity_strength = float(data.get("gravity_strength", 1.0)) + elastic_force_enabled = int(data.get("elastic_force", 1)) + damping_force_enabled = int(data.get("damping_force", 0)) # 动能 Ek = ½ m v² ek = 0.5 * masses[np.newaxis, :] * (all_vx**2 + all_vy**2 + all_vz**2) - # 重力势能 Ug = -m G·r - ug = -masses[np.newaxis, :] * ( - G_vec[0] * all_x + G_vec[1] * all_y + G_vec[2] * all_z - ) + # 均匀重力场势能 Ug = -m G·r + ug = np.zeros_like(ek) + if gravity_field_enabled: + ug = -masses[np.newaxis, :] * ( + G_vec[0] * all_x + G_vec[1] * all_y + G_vec[2] * all_z + ) # 弹性势能 Us = ½ k (d - d₀)² us = np.zeros_like(ek) bond_pairs = data.get("bond_pairs") bond_stiffness = data.get("bond_stiffness") bond_rest_lengths = data.get("bond_rest_lengths") - if bond_pairs is not None and len(bond_pairs) > 0: + if elastic_force_enabled and bond_pairs is not None and len(bond_pairs) > 0: for b_idx in range(len(bond_pairs)): i, j = bond_pairs[b_idx] dx = all_x[:, j] - all_x[:, i] @@ -350,17 +391,33 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", dz = all_z[:, j] - all_z[:, i] dist = np.sqrt(dx**2 + dy**2 + dz**2) stretch = dist - bond_rest_lengths[b_idx] - us[:, i] += 0.5 * bond_stiffness[b_idx] * stretch**2 - us[:, j] += 0.5 * bond_stiffness[b_idx] * stretch**2 + us_each = 0.5 * bond_stiffness[b_idx] * stretch**2 + us[:, i] += us_each # 将整根键的势能记给 i + + # 万有引力势能 Ug_grav = -G_grav * m_i * m_j / r + ug_grav = np.zeros_like(ek) + if gravity_interaction_enabled: + n_atoms_en = len(masses) + for i in range(n_atoms_en): + for j in range(i + 1, n_atoms_en): + dx = all_x[:, j] - all_x[:, i] + dy = all_y[:, j] - all_y[:, i] + dz = all_z[:, j] - all_z[:, i] + dist = np.sqrt(dx**2 + dy**2 + dz**2) + dist = np.maximum(dist, 1e-12) + pair_pe = -gravity_strength * masses[i] * masses[j] / dist + ug_grav[:, i] += 0.5 * pair_pe + ug_grav[:, j] += 0.5 * pair_pe # 各原子总能量 - e_total = ek + ug + us # (NT, n_atoms) + e_total = ek + ug + us + ug_grav # (NT, n_atoms) # 系统能量分量 ek_sys = np.sum(ek, axis=1) ug_sys = np.sum(ug, axis=1) us_sys = np.sum(us, axis=1) - e_sys = ek_sys + ug_sys + us_sys + ug_grav_sys = np.sum(ug_grav, axis=1) + e_sys = ek_sys + ug_sys + us_sys + ug_grav_sys # ── 第 3 行左：各原子总能量 ── ax_e = axes[2, 0] @@ -376,9 +433,11 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", # ── 第 3 行右：系统总能量 ── ax_sys = axes[2, 1] ax_sys.plot(time, ek_sys, 'b-', linewidth=1.5, label="系统动能") - ax_sys.plot(time, ug_sys, 'g-', linewidth=1.5, label="系统重力势能") - if bond_pairs is not None and len(bond_pairs) > 0: + ax_sys.plot(time, ug_sys, 'g-', linewidth=1.5, label="均匀重力势能") + if elastic_force_enabled and bond_pairs is not None and len(bond_pairs) > 0: ax_sys.plot(time, us_sys, color='orange', linewidth=1.5, label="系统弹性势能") + if gravity_interaction_enabled: + ax_sys.plot(time, ug_grav_sys, color='purple', linewidth=1.5, label="万有引力势能") ax_sys.plot(time, e_sys, 'r--', linewidth=1.5, label="系统总能量") ax_sys.set_title("系统总能量") ax_sys.set_xlabel("时间 (s)") @@ -421,8 +480,8 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", def main(): parser = argparse.ArgumentParser(description="物理模拟统一入口") - parser.add_argument("config_file", nargs="?", default=os.path.join("input", "parameters.yaml"), - help="YAML 配置文件路径（默认: input/parameters.yaml）") + parser.add_argument("config_file", nargs="?", default=os.path.join("input", "input.txt"), + help="YAML 配置文件路径（默认: input/input.txt，虽然是 .txt 后缀但使用 YAML 格式）") parser.add_argument("--config", dest="config_override", help="YAML 配置文件路径（可选，优先于位置参数）") parser.add_argument("--input-dir", default="input", diff --git a/engines/c/main.c b/engines/c/main.c index 695927b..7c2bc19 100644 --- a/engines/c/main.c +++ b/engines/c/main.c @@ -2,13 +2,15 @@ * engines/c/main.c * ----------------- * C 语言动力学模拟引擎。 + * 与 Python 版 (compute.py) 算法保持一致。 + * * 输入: param.json 数值参数（Python 从 YAML 转换得来） * /coord.txt * /connection.txt * /bond.txt * 输出: /trajectory.txt (JSON 格式，与 Python 版兼容) * - * 编译: make + * 编译: cmake --build build --target dynamics_c * 用法: ./build/dynamics_c */ @@ -27,8 +29,14 @@ typedef struct { 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; /* 万有引力强度 */ } SimParams; /* ======================================================================== @@ -79,7 +87,7 @@ static void *xmalloc(size_t sz) { return p; } -/* 从 JSON 中读取 double 值。简单实现：搜索 key 并解析 */ +/* 从 JSON 中读取 double 值 */ static double json_read_double(const char *json, const char *key) { char search[256]; snprintf(search, sizeof(search), "\"%s\"", key); @@ -96,6 +104,23 @@ 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]; @@ -112,6 +137,12 @@ static void json_read_double3(const char *json, const char *key, double out[3]) } /* 读取 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"); @@ -129,8 +160,20 @@ static SimParams read_params(const char *path) { 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"); + 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; @@ -147,7 +190,6 @@ static AtomData read_coord(const char *input_dir) { char line[1024]; if (!fgets(line, sizeof(line), f)) die("coord.txt 为空"); - /* 先读入所有行到内存 */ int capacity = 16; AtomData a; a.n_atoms = 0; @@ -173,10 +215,9 @@ static AtomData read_coord(const char *input_dir) { 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; /* 跳过空行/格式不符的行 */ + continue; } int i = a.n_atoms; a.atom_ids[i] = id; @@ -204,13 +245,11 @@ static BondData read_bonds(const char *input_dir, const AtomData *atoms) { snprintf(path, sizeof(path), "%s/connection.txt", input_dir); FILE *f = fopen(path, "r"); - if (!f) return b; /* 无成键 */ + 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++; @@ -223,24 +262,21 @@ static BondData read_bonds(const char *input_dir, const AtomData *atoms) { b.stiffness = (double*)xmalloc(n_lines * sizeof(double)); b.rest_lengths = (double*)xmalloc(n_lines * sizeof(double)); - /* 读取 bond.txt 获取劲度系数和平衡长度 */ char bond_path[512]; snprintf(bond_path, sizeof(bond_path), "%s/bond.txt", input_dir); FILE *fb = fopen(bond_path, "r"); - /* bond.txt 格式: name stiffness rest_length */ for (int i = 0; i < n_lines; i++) { fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name); - /* 转换原子编号为 0-based 索引 */ b.pairs[i*2+0] = tmp_a - 1; b.pairs[i*2+1] = tmp_b - 1; - /* 查找 bond.txt 中对应的参数 */ b.stiffness[i] = 1.0; - b.rest_lengths[i] = 2.0; /* 默认值 */ + b.rest_lengths[i] = 2.0; if (fb) { - char name[256]; + 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; @@ -256,9 +292,10 @@ static BondData read_bonds(const char *input_dir, const AtomData *atoms) { } /* ======================================================================== - * 物理核心 + * 物理核心（与 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, @@ -266,50 +303,124 @@ static void compute_acceleration( const BondData *bonds, double *ax, double *ay, double *az) { + /* 先清零 */ for (int i = 0; i < n; i++) { - double inv_m = 1.0 / m[i]; - ax[i] = (m[i] * G[0] - B[0] * vx[i]) * inv_m; - ay[i] = (m[i] * G[1] - B[1] * vy[i]) * inv_m; - az[i] = (m[i] * G[2] - B[2] * vz[i]) * inv_m; + ax[i] = 0.0; ay[i] = 0.0; az[i] = 0.0; } - /* 弹簧力 */ - 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_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]; + } + } } } -/* clamp */ -static double clamp(double v, double lo, double hi) { - return fmin(fmax(v, lo), hi); +/* 边界条件：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; } } -/* 蛙跳法一步 */ -static void leapfrog_step( +/* ── 显式欧拉法 ──────────── */ +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 box_a, double dt) + 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; @@ -318,16 +429,148 @@ static void leapfrog_step( x[i] += vx[i] * dt; y[i] += vy[i] * dt; z[i] += vz[i] * dt; - /* wrap */ - x[i] = clamp(x[i], -box_a, box_a); - y[i] = clamp(y[i], -box_a, box_a); - z[i] = clamp(z[i], -box_a, box_a); } } -/* ======================================================================== - * JSON 输出 - * ======================================================================== */ +/* ── 中点法 ──────────── */ +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）── */ +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)); + compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az); + + /* 半推速度：v_half = v + 0.5*a*dt */ + 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 * 0.5; + vy[i] += ay[i] * dt * 0.5; + vz[i] += az[i] * dt * 0.5; + } + + /* 全推位置（不含边界）*/ + 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; /* vx 此时是 v_half */ + y[i] += vy[i] * dt; + z[i] += vz[i] * dt; + } + + /* 显式预测器：v_pred = v_half + 0.5*a_old*dt，用第一次加速度外推半步 + 包含重力+阻尼+弹簧的所有贡献（标准 Velocity-Verlet 预测步）*/ + double *pred_vx = (double*)alloca(n * sizeof(double)); + double *pred_vy = (double*)alloca(n * sizeof(double)); + double *pred_vz = (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]) continue; + pred_vx[i] = vx[i] + 0.5 * ax[i] * dt; + pred_vy[i] = vy[i] + 0.5 * ay[i] * dt; + pred_vz[i] = vz[i] + 0.5 * az[i] * dt; + } + + /* 用新位置 + 预测速度重算加速度 */ + compute_acceleration(n, x, y, z, pred_vx, pred_vy, pred_vz, m, G, B, bonds, ax, ay, az); + + /* 速度后半步：v = v_half + 0.5*a_next*dt + vx 仍为 v_half（未被覆盖）*/ + 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 * 0.5; + vy[i] += ay[i] * dt * 0.5; + vz[i] += az[i] * dt * 0.5; + } +} + +/* ── 分发器：调用对应积分方法 + 边界条件 + 自由度约束（与 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 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, @@ -338,7 +581,6 @@ static void write_trajectory_json(const char *path, const Trajectory *traj, fprintf(f, "{\n"); - /* traj_x, traj_y, ... */ 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}; @@ -355,7 +597,6 @@ static void write_trajectory_json(const char *path, const Trajectory *traj, fputc('\n', f); } fprintf(f, " ]"); - // 每个轨迹数组后加逗号（后面还有标量参数） fputc(',', f); fputc('\n', f); } @@ -364,45 +605,43 @@ static void write_trajectory_json(const char *path, const Trajectory *traj, fprintf(f, " \"NT\": %d,\n", params->NT); fprintf(f, " \"DT\": %.15g,\n", params->DT); fprintf(f, " \"NSTEP\": %d,\n", params->NSTEP); - fprintf(f, " \"method\": \"leapfrog\",\n"); + fprintf(f, " \"method\": \"%s\",\n", params->method); fprintf(f, " \"warmup_steps\": %d,\n", params->warmup_steps); fprintf(f, " \"G\": [%.15g, %.15g, %.15g],\n", params->G[0], params->G[1], params->G[2]); fprintf(f, " \"B\": [%.15g, %.15g, %.15g],\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]); - if (i < atoms->n_atoms - 1) fputc(',', f); } fprintf(f, "],\n"); fprintf(f, " \"atom_masses\": ["); for (int i = 0; i < atoms->n_atoms; i++) { + if (i > 0) fputc(',', f); fprintf(f, "%.15g", atoms->masses[i]); - if (i < atoms->n_atoms - 1) fputc(',', f); } 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]); - if (b < bonds->n_bonds - 1) fputc(',', f); } fprintf(f, "],\n"); fprintf(f, " \"bond_stiffness\": ["); for (int b = 0; b < bonds->n_bonds; b++) { + if (b > 0) fputc(',', f); fprintf(f, "%.15g", bonds->stiffness[b]); - if (b < bonds->n_bonds - 1) fputc(',', f); } fprintf(f, "],\n"); fprintf(f, " \"bond_rest_lengths\": ["); for (int b = 0; b < bonds->n_bonds; b++) { + if (b > 0) fputc(',', f); fprintf(f, "%.15g", bonds->rest_lengths[b]); - if (b < bonds->n_bonds - 1) fputc(',', f); } fprintf(f, "]\n"); @@ -410,9 +649,9 @@ static void write_trajectory_json(const char *path, const Trajectory *traj, fclose(f); } -/* ======================================================================== - * 主函数 - * ======================================================================== */ +// ======================================================================== +// 主函数 +// ======================================================================== int main(int argc, char **argv) { if (argc < 4) { @@ -425,15 +664,13 @@ int main(int argc, char **argv) { clock_t t0 = clock(); - /* 读取参数和输入 */ SimParams params = read_params(param_path); AtomData atoms = read_coord(input_dir); BondData bonds = read_bonds(input_dir, &atoms); - printf("[C-engine] 原子数=%d, 键数=%d, NT=%d, DT=%.6g\n", - atoms.n_atoms, bonds.n_bonds, params.NT, params.DT); + printf("[C-engine] 原子数=%d, 键数=%d, NT=%d, DT=%.6g, method=%s\n", + atoms.n_atoms, bonds.n_bonds, 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)); @@ -450,30 +687,28 @@ int main(int argc, char **argv) { vz[i] = atoms.vel_0[i*3+2]; } - /* 分配轨迹缓冲区 */ + /* 分配轨迹缓冲区：改用 record_steps */ + int record_steps = params.NT - params.warmup_steps; Trajectory traj; - traj.n_steps = params.NT; + traj.n_steps = record_steps; traj.n_atoms = n; - traj.x = (double*)xmalloc(params.NT * n * sizeof(double) * 6); - traj.y = traj.x + params.NT * n; - traj.z = traj.y + params.NT * n; - traj.vx = traj.z + params.NT * n; - traj.vy = traj.vx + params.NT * n; - traj.vz = traj.vy + params.NT * n; + 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; /* 预热 */ for (int s = 0; s < params.warmup_steps; s++) { - leapfrog_step(n, x, y, z, vx, vy, vz, - atoms.masses, params.G, params.B, &bonds, atoms.fixed, - params.box_a, params.DT); + 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); } /* 记录 */ - for (int s = 0; s < params.NT; s++) { - leapfrog_step(n, x, y, z, vx, vy, vz, - atoms.masses, params.G, params.B, &bonds, atoms.fixed, - params.box_a, params.DT); - /* 保存这一帧 */ + for (int s = 0; s < record_steps; s++) { for (int i = 0; i < n; i++) { traj.x[ s * n + i] = x[i]; traj.y[ s * n + i] = y[i]; @@ -482,16 +717,19 @@ int main(int argc, char **argv) { traj.vy[s * n + i] = vy[i]; traj.vz[s * n + i] = vz[i]; } + 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]; snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir); write_trajectory_json(out_path, &traj, ¶ms, &atoms, &bonds); clock_t t1 = clock(); double elapsed = (double)(t1 - t0) / CLOCKS_PER_SEC; - printf("[C-engine] 计算完成: %d 步, %.3f s\n", params.NT, elapsed); + printf("[C-engine] 计算完成: %d 步, %.3f s\n", record_steps, elapsed); /* 清理 */ free(traj.x); diff --git a/engines/cpp/main.cpp b/engines/cpp/main.cpp index 42be97f..65eeaf3 100644 --- a/engines/cpp/main.cpp +++ b/engines/cpp/main.cpp @@ -1,54 +1,633 @@ /** * engines/cpp/main.cpp * -------------------- - * C++ 动力学模拟引擎（模板/参考实现）。 + * C++ 动力学模拟引擎。 + * 与 Python 版 (compute.py) 算法保持一致。 * - * 接口与 C 引擎一致： - * 输入: param.json, /coord.txt, connection.txt, bond.txt - * 输出: /trajectory.txt (JSON) + * 输入: param.json, /coord.txt, connection.txt, bond.txt + * 输出: /trajectory.txt (JSON, 与 Python 版兼容) * * 编译: - * g++ -O3 -march=native -o build/dynamics_cpp main.cpp + * 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 #include -#include -#include -#include +// ======================================================================== +// 配置参数（从 param.json 读取） +// ======================================================================== struct SimParams { - double box_a; - int NT; - double DT; - int NSTEP; - int warmup_steps; - double G[3]; - double B[3]; -}; - -struct Atom { - int id; - double mass, radius; - double x, y, z; - double vx, vy, vz; - bool fx, fy, fz; // fixed constraints -}; - -struct Bond { - int i, j; - double stiffness, rest_length; + 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; }; // ======================================================================== -// 请在下方实现 Leapfrog 算法 -// 参考 engines/c/main.c 中的物理核心实现 +// 原子数据 +// ======================================================================== +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; +}; + +// ======================================================================== +// 辅助函数 +// ======================================================================== + +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 数组 (如 "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(); + } +} + +/* 解析 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"); + 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; +} + +// ======================================================================== +// 物理核心 +// ======================================================================== + +/* 加速度计算（各力独立开关控制）——与 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]; + } + } + } +} + +/* 边界条件：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: 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 一致 ── */ +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); + 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()); + + // 半推速度：v_half = v + 0.5*a*dt (存入 vx, vy, vz) + 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 * 0.5; + vy[i] += ay[i] * dt * 0.5; + vz[i] += az[i] * dt * 0.5; + } + + // 全推位置（不含边界，边界在外层统一处理） + 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; // vx 此时是 v_half + y[i] += vy[i] * dt; + z[i] += vz[i] * dt; + } + + // 显式预测器：v_pred = v_half + 0.5*a_old*dt，用第一次加速度外推半步 + // 包含所有力的贡献（标准 Velocity-Verlet 预测步） + std::vector pred_vx(n), pred_vy(n), pred_vz(n); + for (int i = 0; i < n; i++) { + if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; + pred_vx[i] = vx[i] + 0.5 * ax[i] * dt; + pred_vy[i] = vy[i] + 0.5 * ay[i] * dt; + pred_vz[i] = vz[i] + 0.5 * az[i] * dt; + } + + // 用新位置 + 预测速度重算加速度 + compute_acceleration(n, x, y, z, pred_vx.data(), pred_vy.data(), pred_vz.data(), + m, G, B, bonds, + gravity_field, gravity_interaction, + elastic_force, damping_force, gravity_strength, + ax.data(), ay.data(), az.data()); + + // 速度后半步：v = v_half + 0.5*a_next*dt + // vx 仍为 v_half（未被覆盖），直接加上 0.5*a_next*dt + 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 * 0.5; + vy[i] += ay[i] * dt * 0.5; + vz[i] += az[i] * dt * 0.5; + } +} + +/* ── 分发器：调用对应积分方法 + 边界条件（与 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 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; } + } +} + +// ======================================================================== +// 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); + f << std::setprecision(15); + + 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 << "}\n"; +} + +// ======================================================================== +// 主函数 // ======================================================================== int main(int argc, char **argv) { @@ -63,15 +642,80 @@ int main(int argc, char **argv) { auto t0 = std::chrono::high_resolution_clock::now(); - // TODO: 读取 param.json、coord.txt、connection.txt、bond.txt - // TODO: 执行 Leapfrog 模拟 - // TODO: 输出 trajectory.txt (JSON 格式) - // - // 提示：输出格式与 Python 版兼容，参考 engines/c/main.c 中的 write_trajectory_json + // 读取参数和输入 + SimParams params = read_params(param_path); + AtomData atoms = read_coord(input_dir); + BondData bonds = read_bonds(input_dir); + + 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; + std::vector traj_x(record_steps * n); + std::vector traj_y(record_steps * n); + std::vector traj_z(record_steps * n); + std::vector traj_vx(record_steps * n); + std::vector traj_vy(record_steps * n); + std::vector traj_vz(record_steps * n); + + // 预热 + for (int s = 0; s < params.warmup_steps; s++) { + 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); + } + + // 记录 + for (int s = 0; s < record_steps; s++) { + // 保存当前帧 + for (int i = 0; i < n; i++) { + traj_x[s * n + i] = x[i]; + traj_y[s * n + i] = y[i]; + traj_z[s * n + i] = z[i]; + traj_vx[s * n + i] = vx[i]; + traj_vy[s * n + i] = vy[i]; + traj_vz[s * n + i] = vz[i]; + } + + 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); + } + + // 输出轨迹 + 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; - std::cout << "[C++-engine] 占位: NT=" << 0 << ", " << elapsed << " s (待实现)" << std::endl; return 0; } diff --git a/engines/fortran/main.f90 b/engines/fortran/main.f90 index ada2964..6edfdb9 100644 --- a/engines/fortran/main.f90 +++ b/engines/fortran/main.f90 @@ -1,22 +1,49 @@ ! engines/fortran/main.f90 ! ------------------------ -! Fortran 动力学模拟引擎（模板/参考实现）。 +! Fortran 动力学模拟引擎。 +! 与 Python 版 (compute.py) 算法保持一致。 ! -! 接口与 C 引擎一致： -! 输入: param.json, /coord.txt, connection.txt, bond.txt -! 输出: /trajectory.txt (JSON) +! 输入: param.json, /coord.txt, connection.txt, bond.txt +! 输出: /trajectory.txt (JSON) ! -! 编译: -! gfortran -O3 -march=native -o build/dynamics_f90 main.f90 -! -! 用法: -! ./build/dynamics_f90 +! 编译: 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 + 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 + double precision :: gravity_strength + character(len=32) :: method + double precision :: t0, t1, elapsed + + ! 原子数据 + 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(:) + + ! 运行时位置/速度 + double precision, allocatable :: x(:), y(:), z(:) + double precision, allocatable :: vx(:), vy(:), vz(:) + + ! 轨迹缓冲区 + integer :: record_steps + 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 " @@ -27,14 +54,808 @@ program dynamics_f90 call get_command_argument(2, output_dir) call get_command_argument(3, param_path) - ! TODO: 实现与 C 引擎相同的物理模拟 - ! - ! 参考实现: - ! 1. 读取 param.json → NT, DT, G, B, ... - ! 2. 读取 coord.txt → 原子数, 质量, 位置, 速度 - ! 3. 读取 connection.txt + bond.txt → 成键参数 - ! 4. Leapfrog 循环 NT 步 - ! 5. 输出 trajectory.txt (JSON 格式) + 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) + + ! 读取 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) + + 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 + allocate(traj_x(record_steps, n), traj_y(record_steps, n), traj_z(record_steps, n)) + allocate(traj_vx(record_steps, n), traj_vy(record_steps, n), traj_vz(record_steps, n)) + + ! 预热 + do s = 1, warmup_steps + 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 + + ! 记录 + do s = 1, record_steps + traj_x(s, :) = x; traj_y(s, :) = y; traj_z(s, :) = z + traj_vx(s, :) = vx; traj_vy(s, :) = vy; traj_vz(s, :) = vz + call apply_step(method, n, x, y, z, vx, vy, vz, masses, G, B, & + 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 + + ! 输出轨迹 + call write_json(output_dir, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, & + record_steps, n_atoms, atom_ids, masses, & + NT, DT, NSTEP, warmup_steps, method, G, B, & + n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths) + + 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) + +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) + 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 + 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 + + 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) +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) /= ' ' .and. (i == 1 .or. line(i-1:i-1) == ' ')) ncols = ncols + 1 + 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 + +! 边界条件：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 + +! ── 显式欧拉法 ──────────── +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 + +! ── 蛙跳法（Velocity-Verlet）── +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 :: dmp_vx(n), dmp_vy(n), dmp_vz(n) + double precision :: gx, gy, gz + 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 + vx(i) = vx(i) + ax(i) * dt * 0.5d0 + vy(i) = vy(i) + ay(i) * dt * 0.5d0 + vz(i) = vz(i) + az(i) * dt * 0.5d0 + end do + + ! 全推位置（不含边界） + 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 + end do + + ! 隐式阻尼处理（用临时数组 dmp_v，不覆盖 vx/vy/vz） + do i = 1, n + if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then + dmp_vx(i) = 0; dmp_vy(i) = 0; dmp_vz(i) = 0; cycle + end if + gx = B(1) / m(i); gy = B(2) / m(i); gz = B(3) / m(i) + dmp_vx(i) = (vx(i) + 0.5d0 * G(1) * dt) / (1.0d0 + 0.5d0 * gx * dt) + dmp_vy(i) = (vy(i) + 0.5d0 * G(2) * dt) / (1.0d0 + 0.5d0 * gy * dt) + dmp_vz(i) = (vz(i) + 0.5d0 * G(3) * dt) / (1.0d0 + 0.5d0 * gz * dt) + end do + + ! 用新位置 + 阻尼速度重算加速度 + call accel(n, x, y, z, dmp_vx, dmp_vy, dmp_vz, m, G, B, nb, bp, bk, br, & + gravity_field, gravity_interaction, & + elastic_force, damping_force, gravity_strength, ax, ay, az) + + ! 速度后半步：v = v_half + 0.5*a_next*dt（vx 仍为 v_half） + 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 * 0.5d0 + vy(i) = vy(i) + ay(i) * dt * 0.5d0 + vz(i) = vz(i) + az(i) * dt * 0.5d0 + 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 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 + +! ======================================================================== +! JSON 输出 +! ======================================================================== + +subroutine write_json(outdir, tx, ty, tz, tvx, tvy, tvz, & + nsteps, nat, aid, amass, & + NT, DT, NSTEP, warmup, method, G, B, & + nb, bp, bk, br) + character(len=*), intent(in) :: outdir, method + integer, intent(in) :: nsteps, nat, NT, NSTEP, warmup, nb, bp(nb, 2), aid(nat) + double precision, intent(in) :: tx(nsteps, nat), ty(nsteps, nat), tz(nsteps, nat) + double precision, intent(in) :: tvx(nsteps, nat), tvy(nsteps, nat), tvz(nsteps, nat) + double precision, intent(in) :: DT, G(3), B(3), bk(nb), br(nb), amass(nat) + + character(len=512) :: path, buf + integer :: u, s, i, ib, ios + + path = trim(outdir) // '/trajectory.txt' + open(newunit=u, file=trim(path), status='replace', action='write', iostat=ios) + if (ios /= 0) then + write(*, '("[Fortran-engine] 错误: 无法写入 ", a)') trim(path) + stop + end if + + write(u, '(a)') '{' + + ! traj_x + write(u, '(a)') ' "traj_x": [' + do s = 1, nsteps + call json_arr(u, tx(s, :), nat, s < nsteps, ' ') + end do + write(u, '(a)') ' ],' + + ! traj_y + write(u, '(a)') ' "traj_y": [' + do s = 1, nsteps + call json_arr(u, ty(s, :), nat, s < nsteps, ' ') + end do + write(u, '(a)') ' ],' + + ! traj_z + write(u, '(a)') ' "traj_z": [' + do s = 1, nsteps + call json_arr(u, tz(s, :), nat, s < nsteps, ' ') + end do + write(u, '(a)') ' ],' + + ! traj_vx + write(u, '(a)') ' "traj_vx": [' + do s = 1, nsteps + call json_arr(u, tvx(s, :), nat, s < nsteps, ' ') + end do + write(u, '(a)') ' ],' + + ! traj_vy + write(u, '(a)') ' "traj_vy": [' + do s = 1, nsteps + call json_arr(u, tvy(s, :), nat, s < nsteps, ' ') + end do + write(u, '(a)') ' ],' + + ! traj_vz + write(u, '(a)') ' "traj_vz": [' + do s = 1, nsteps + call json_arr(u, tvz(s, :), nat, s < nsteps, ' ') + end do + write(u, '(a)') ' ],' + + ! 标量参数 + write(buf, '(a, i0, a)') ' "NT": ', NT, ',' + write(u, '(a)') trim(buf) + write(buf, '(a, g0, a)') ' "DT": ', DT, ',' + write(u, '(a)') trim(buf) + write(buf, '(a, i0, a)') ' "NSTEP": ', NSTEP, ',' + write(u, '(a)') trim(buf) + write(buf, '(a, a, a)') ' "method": "', trim(method), '",' + write(u, '(a)') trim(buf) + write(buf, '(a, i0, a)') ' "warmup_steps": ', warmup, ',' + write(u, '(a)') trim(buf) + + write(buf, '(a, g0, a, g0, a, g0, a)') & + ' "G": [', G(1), ', ', G(2), ', ', G(3), '],' + write(u, '(a)') trim(buf) + write(buf, '(a, g0, a, g0, a, g0, a)') & + ' "B": [', B(1), ', ', B(2), ', ', B(3), '],' + write(u, '(a)') trim(buf) + + ! 原子信息 + write(u, '(a)', advance='no') ' "atom_ids": [' + do i = 1, nat + if (i > 1) write(u, '(a)', advance='no') ',' + write(u, '(i0)', advance='no') aid(i) + end do + write(u, '(a)') '],' + + write(u, '(a)', advance='no') ' "atom_masses": [' + do i = 1, nat + if (i > 1) write(u, '(a)', advance='no') ',' + write(u, '(g0)', advance='no') amass(i) + end do + write(u, '(a)') '],' + + ! 成键 + if (nb > 0) then + call write_int2_arr(u, 'bond_pairs', bp, nb, .true.) + call write_dbl_arr(u, 'bond_stiffness', bk, nb, .true.) + call write_dbl_arr(u, 'bond_rest_lengths', br, nb, .false.) + else + write(u, '(a)') ' "bond_pairs": [],' + write(u, '(a)') ' "bond_stiffness": [],' + write(u, '(a)') ' "bond_rest_lengths": []' + end if + + write(u, '(a)') '}' + close(u) +end subroutine write_json + +! 写出单行 JSON 数组 [v1, v2, ...] +subroutine json_arr(u, vals, n, has_next, indent) + integer, intent(in) :: u, n + double precision, intent(in) :: vals(n) + logical, intent(in) :: has_next + character(len=*), intent(in) :: indent + integer :: i + write(u, '(a)', advance='no') indent // '[' + do i = 1, n + if (i > 1) write(u, '(a)', advance='no') ',' + write(u, '(g0)', advance='no') vals(i) + end do + if (has_next) then + write(u, '(a)') '],' + else + write(u, '(a)') ']' + end if +end subroutine json_arr + +subroutine write_int2_arr(u, name, arr, n, has_next) + integer, intent(in) :: u, n, arr(n, 2) + character(len=*), intent(in) :: name + logical, intent(in) :: has_next + character(len=65536) :: buf + integer :: i, pos + write(u, '(a)', advance='no') ' "' // trim(name) // '": [' + do i = 1, n + if (i > 1) write(u, '(a)', advance='no') ',' + write(buf, '(a, i0, a, i0, a)') '[', arr(i, 1), ',', arr(i, 2), ']' + write(u, '(a)', advance='no') trim(buf) + end do + if (has_next) then + write(u, '(a)') '],' + else + write(u, '(a)') ']' + end if +end subroutine write_int2_arr + +subroutine write_dbl_arr(u, name, arr, n, has_next) + integer, intent(in) :: u, n + double precision, intent(in) :: arr(n) + character(len=*), intent(in) :: name + logical, intent(in) :: has_next + integer :: i + write(u, '(a)', advance='no') ' "' // trim(name) // '": [' + do i = 1, n + if (i > 1) write(u, '(a)', advance='no') ',' + write(u, '(g0)', advance='no') arr(i) + end do + if (has_next) then + write(u, '(a)') '],' + else + write(u, '(a)') ']' + end if +end subroutine write_dbl_arr - write(*,*) "[Fortran-engine] 待实现" end program dynamics_f90 diff --git a/examples/case01/input/coord.txt b/examples/case01/input/coord.txt index 9a2e805..23b9cff 100644 --- a/examples/case01/input/coord.txt +++ b/examples/case01/input/coord.txt @@ -1,3 +1,3 @@ -n mass radius x y z vx vy vz -1 1 0.28 -1 0 0 0 0 0 -2 1 0.28 1 0 1 0 0 0 +n mass radius x y z vx vy vz fix_x fix_y fix_z +1 1 0.28 -1 0 0 0 0 0 0 0 0 +2 1 0.28 1 0 1 0 0 0 0 0 0 diff --git a/examples/case01/input/parameters.yaml b/examples/case01/input/input.txt similarity index 79% rename from examples/case01/input/parameters.yaml rename to examples/case01/input/input.txt index 0ca9378..8be2c8d 100644 --- a/examples/case01/input/parameters.yaml +++ b/examples/case01/input/input.txt @@ -7,18 +7,20 @@ # 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧 step_simulate: 1 # 运行物理模拟 → output/trajectory.txt step_sample: 1 # 抽帧 → output/display.txt -step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_plot: 1 # 绘制轨迹/能量图 → output/trajectory_plots.png step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy） +force_calc: 0 # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出 # ── 计算引擎 ────────────────────────────────── # 可选: python, c, cpp, fortran, java # python = Python 参考实现(compute.py) # c = C 引擎 (engines/c/build/dynamics_c) # cpp = C++ 引擎 (engines/cpp/build/dynamics_cpp) -engine: python # 默认使用 Python 引擎 +# fortran = Fortran 引擎 (engines/fortran/build/dynamics_f90) +engine: python # 默认使用 Python 引擎 # ── 盒子 ────────────────────────────────────── -box_a: 10.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内 +box_a: 20.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内 # ── 初始构型 ────────────────────────────────── # 坐标文件格式： @@ -37,6 +39,14 @@ G: [0.0, 0.0, -9.8] # 重力场分量 (m/s²) # B: [0.5, 0.5, 0.5] # 阻尼分量 B: [0.0, 0.0, 0.0] # 阻尼分量 +# ── 力开关（0=关闭, 1=开启）────────────────── +gravity_field: 1 # 均匀重力场 (G) +gravity_interaction: 0 # 原子间万有引力 +elastic_force: 1 # 弹簧键力 +damping_force: 0 # 阻尼 (B) +# +gravity_strength: 1.0 # 万有引力强度（仅 gravity_interaction=1 时有效） + # ── 数值算法 ────────────────────────────────── # 可选： # explicit_euler 显式欧拉法 @@ -51,18 +61,21 @@ method: leapfrog # 预热步数：模拟开始时跳过不保存的步数（用于稳定初始状态） warmup_steps: 0 # 默认 0（立即开始记录） -# 总计算步数 -NT: 10000 +# 总模拟时间（秒），程序自动计算 NT = T_total / DT +# 如果同时指定了 NT，以 NT 为准 +T_total: 10.0 # 抽帧间隔（每 NSTEP 步取一帧用于动画） NSTEP: 100 +# ── 时间步长 ────────────────────────────────── +DT: 0.001 # 时间步长 (s) + # 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧 sample_start: null # null 表示从头开始（帧索引从 0 起） sample_end: null # null 表示到末尾 -# ── 时间步长 ────────────────────────────────── -DT: 0.001 # 时间步长 (s) + # ── 显示参数 ────────────────────────────────── # 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序 diff --git a/examples/case01/run_dynamics.py b/examples/case01/run_dynamics.py index 6e1c34a..d092d49 100644 --- a/examples/case01/run_dynamics.py +++ b/examples/case01/run_dynamics.py @@ -18,7 +18,7 @@ CASE_DIR = Path(__file__).resolve().parent DYNAMICS_PATH = Path("..") / ".." / "dynamics.py" INPUT_DIR = Path("input") OUTPUT_DIR = Path("output") -CONFIG_FILE = INPUT_DIR / "parameters.yaml" +CONFIG_FILE = INPUT_DIR / "input.txt" def load_dynamics_module(module_path: Path): diff --git a/examples/case02/input/bond.txt b/examples/case02/input/bond.txt new file mode 100644 index 0000000..2b183d5 --- /dev/null +++ b/examples/case02/input/bond.txt @@ -0,0 +1,2 @@ +bond_name k rest_length +k1 100.0 2.0 diff --git a/examples/case02/input/connection.txt b/examples/case02/input/connection.txt new file mode 100644 index 0000000..c845a64 --- /dev/null +++ b/examples/case02/input/connection.txt @@ -0,0 +1,2 @@ +n1 n2 bond_name + diff --git a/examples/case02/input/coord.txt b/examples/case02/input/coord.txt new file mode 100644 index 0000000..417b237 --- /dev/null +++ b/examples/case02/input/coord.txt @@ -0,0 +1,3 @@ +n mass radius x y z vx vy vz fix_x fix_y fix_z +1 1 0.28 0 0 0 0 0 0 1 1 1 +2 1 0.28 0 0 4 0 4 0 0 0 0 diff --git a/examples/case02/input/input.txt b/examples/case02/input/input.txt new file mode 100644 index 0000000..222af60 --- /dev/null +++ b/examples/case02/input/input.txt @@ -0,0 +1,93 @@ +# 物理模拟参数配置 +# 格式：YAML +# 用法：python run_dynamics.py + +# ── 流程控制 ────────────────────────────────── +# 每步用 0/1 单独开关，1=执行，0=跳过 +# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧 +step_simulate: 1 # 运行物理模拟 → output/trajectory.txt +step_sample: 1 # 抽帧 → output/display.txt +step_plot: 1 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy） +force_calc: 0 # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出 + +# ── 计算引擎 ────────────────────────────────── +# 可选: python, c, cpp, fortran, java +# python = Python 参考实现(compute.py) +# c = C 引擎 (engines/c/build/dynamics_c) +# cpp = C++ 引擎 (engines/cpp/build/dynamics_cpp) +# fortran = Fortran 引擎 (engines/fortran/build/dynamics_f90) +engine: python # 默认使用 Python 引擎 + +# ── 盒子 ────────────────────────────────────── +box_a: 20.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内 + +# ── 初始构型 ────────────────────────────────── +# 坐标文件格式： +# 第一行：n mass radius x y z vx vy vz +# 后续行：原子序号 质量 半径 x y z vx vy vz +coord_file: input/coord.txt +connection_file: input/connection.txt +bond_file: input/bond.txt + +# 绘图/动画展示的原子序号（对应 coord_file 第一列 n） +plot_atom: 1 + +# ── 物理参数 ────────────────────────────────── +# 三个方向分量分别对应 x, y, z +G: [0.0, 0.0, -9.8] # 重力场分量 (m/s²) +# B: [0.5, 0.5, 0.5] # 阻尼分量 +B: [0.0, 0.0, 0.0] # 阻尼分量 + +# ── 力开关（0=关闭, 1=开启）────────────────── +gravity_field: 0 # 均匀重力场 (G) +gravity_interaction: 1 # 原子间万有引力 +elastic_force: 0 # 弹簧键力 +damping_force: 0 # 阻尼 (B) +# +gravity_strength: 100.0 # 万有引力强度（仅 gravity_interaction=1 时有效） + +# ── 数值算法 ────────────────────────────────── +# 可选： +# explicit_euler 显式欧拉法 +# implicit_euler 隐式欧拉法 +# midpoint 中点法 +# leapfrog 蛙跳法 +method: leapfrog + +# ── 步骤控制 ────────────────────────────────── +# 以下参数控制哪些步骤被执行和保存 + +# 预热步数：模拟开始时跳过不保存的步数（用于稳定初始状态） +warmup_steps: 0 # 默认 0（立即开始记录） + +# 总模拟时间（秒），程序自动计算 NT = T_total / DT +# 如果同时指定了 NT，以 NT 为准 +T_total: 10.0 + +# 抽帧间隔（每 NSTEP 步取一帧用于动画） +NSTEP: 100 + +# ── 时间步长 ────────────────────────────────── +DT: 0.001 # 时间步长 (s) + +# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧 +sample_start: null # null 表示从头开始（帧索引从 0 起） +sample_end: null # null 表示到末尾 + + + +# ── 显示参数 ────────────────────────────────── +# 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序 +alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5] + +# 小球颜色 +# 小球半径从 coord_file 的 radius 列读取 +ball_color_r: 0.90 # R 分量 (0~1) +ball_color_g: 0.20 # G 分量 +ball_color_b: 0.20 # B 分量 + +# 盒子面颜色 +box_color_r: 0.80 +box_color_g: 0.80 +box_color_b: 0.85 diff --git a/examples/case02/run_dynamics.py b/examples/case02/run_dynamics.py new file mode 100644 index 0000000..d092d49 --- /dev/null +++ b/examples/case02/run_dynamics.py @@ -0,0 +1,54 @@ +""" +Case runner for Dynamics case01. + +This script keeps program and data separated: + - program: ../../dynamics.py + - input: ./input + - output: ./output +""" + +from __future__ import annotations + +import argparse +import importlib.util +from pathlib import Path + + +CASE_DIR = Path(__file__).resolve().parent +DYNAMICS_PATH = Path("..") / ".." / "dynamics.py" +INPUT_DIR = Path("input") +OUTPUT_DIR = Path("output") +CONFIG_FILE = INPUT_DIR / "input.txt" + + +def load_dynamics_module(module_path: Path): + spec = importlib.util.spec_from_file_location("dynamics_module", module_path) + if spec is None or spec.loader is None: + raise ImportError(f"无法加载 dynamics.py: {module_path}") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def main(): + parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01") + parser.add_argument("--no-plot", action="store_true", help="跳过 matplotlib 绘图") + args = parser.parse_args() + + dynamics_path = (CASE_DIR / DYNAMICS_PATH).resolve() + input_dir = (CASE_DIR / INPUT_DIR).resolve() + output_dir = (CASE_DIR / OUTPUT_DIR).resolve() + config_path = (CASE_DIR / CONFIG_FILE).resolve() + + module = load_dynamics_module(dynamics_path) + module.run_case( + config_path=config_path, + runtime_base=CASE_DIR, + input_dir=input_dir, + output_dir=output_dir, + no_plot=args.no_plot, + ) + + +if __name__ == "__main__": + main()