From 45513fe33458b8a3f9c12ca32b1eb1cb32efe461 Mon Sep 17 00:00:00 2001 From: Ying-Li Niu <64801511@qq.com> Date: Sun, 17 May 2026 08:47:25 +0800 Subject: [PATCH] init --- .gitattributes | 7 + .gitignore | 101 +++ CMakeLists.txt | 113 +++ INSTALL.md | 568 ++++++++++++++ README.md | 218 +++++- build_engines.py | 93 +++ build_release_zip.py | 36 + cmake/toolchain-linux.cmake | 16 + cmake/toolchain-macos.cmake | 15 + cmake/toolchain-mingw.cmake | 17 + compute.py | 1032 +++++++++++++++++++++++++ draw.py | 492 ++++++++++++ dynamics.py | 454 +++++++++++ engines/c/Makefile | 63 ++ engines/c/main.c | 507 ++++++++++++ engines/cpp/main.cpp | 77 ++ engines/fortran/main.f90 | 40 + examples/case01/input/bond.txt | 2 + examples/case01/input/connection.txt | 2 + examples/case01/input/coord.txt | 3 + examples/case01/input/parameters.yaml | 80 ++ examples/case01/run_dynamics.py | 54 ++ export_web_data.py | 125 +++ migrate_npz_outputs.py | 160 ++++ plot_trajectory.py | 265 +++++++ requirements.txt | 12 + sample.py | 184 +++++ 27 files changed, 4734 insertions(+), 2 deletions(-) create mode 100644 .gitattributes create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 INSTALL.md create mode 100644 build_engines.py create mode 100644 build_release_zip.py create mode 100644 cmake/toolchain-linux.cmake create mode 100644 cmake/toolchain-macos.cmake create mode 100644 cmake/toolchain-mingw.cmake create mode 100644 compute.py create mode 100644 draw.py create mode 100644 dynamics.py create mode 100644 engines/c/Makefile create mode 100644 engines/c/main.c create mode 100644 engines/cpp/main.cpp create mode 100644 engines/fortran/main.f90 create mode 100644 examples/case01/input/bond.txt create mode 100644 examples/case01/input/connection.txt create mode 100644 examples/case01/input/coord.txt create mode 100644 examples/case01/input/parameters.yaml create mode 100644 examples/case01/run_dynamics.py create mode 100644 export_web_data.py create mode 100644 migrate_npz_outputs.py create mode 100644 plot_trajectory.py create mode 100644 requirements.txt create mode 100644 sample.py diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..0b45a66 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,7 @@ +# 全部使用 LF 换行，仓库内外一致，不随系统自动转换 +* text eol=lf + +# 二进制文件不转换 +*.png binary +*.jpg binary +*.ico binary diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..370a6d8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,101 @@ +# ============================================================ +# dynamics 项目 .gitignore +# ============================================================ + +# ── Python ────────────────────────────────────────────────── +__pycache__/ +*.py[cod] +*.pyo +*.pyd +.Python +*.egg-info/ +dist/ +build/ +*.egg +pip-wheel-metadata/ +.env +.venv +venv/ +ENV/ + +# ── C / C++ 编译产物 ───────────────────────────────────────── +# Makefile 构建输出（engines/c/build/） +engines/c/build/ +engines/cpp/build/ + +# CMake 构建目录（根目录或自定义 build 目录） +CMakeCache.txt +CMakeFiles/ +cmake_install.cmake +Makefile.cmake +CTestTestfile.cmake +_CPack_Packages/ +*.cmake.bak +build/ +build_*/ + +# 目标文件 / 静态库 / 共享库 +*.o +*.obj +*.a +*.lib +*.so +*.so.* +*.dylib +*.dll + +# 可执行文件（保留源码，排除编译出的二进制） +# 注意：Windows 下 .exe 后缀的可执行文件 +*.exe +# 但 engines/c/Makefile 里指定了 build/ 目录，已由上面覆盖 + +# 运行时生成的引擎参数文件（每次运行都会覆盖） +engines/*/param.json + +# Fortran 模块文件 +*.mod +*.smod + +# ── examples 输出目录 ──────────────────────────────────────── +# 所有案例的 output/ 目录下均为运行时产物，不提交 +examples/*/output/ + +# 如果希望保留 output/ 目录结构（占位用）， +# 可在各 output/ 目录中放一个 .gitkeep 文件， +# 并在此处改为只忽略具体文件类型： +# examples/*/output/*.txt +# examples/*/output/*.log +# examples/*/output/*.png +# examples/*/output/*.js +# examples/*/output/*.json + +# ── 运行时日志与调试文件 ───────────────────────────────────── +*.log +draw_debug.log + +# ── 根目录运行时输出 ───────────────────────────────────────── +# output/ 目录（根目录下的全局输出） +output/ + +# ── 临时 / 系统文件 ────────────────────────────────────────── +.DS_Store +Thumbs.db +desktop.ini +*.tmp +*.bak +*.swp +*~ + +# ── IDE / 编辑器 ───────────────────────────────────────────── +.vscode/ +.idea/ +*.sublime-project +*.sublime-workspace + +# ── WorkBuddy 工作记忆（本地使用，不提交）─────────────────── +.workbuddy/ + +# ── 发布压缩包 ─────────────────────────────────────────────── +*.zip +*.tar.gz +*.tar.bz2 diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..a1354ad --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,113 @@ +# ============================================================================ +# Dynamics — 多语言引擎 CMake 构建系统 +# ============================================================================ +# 用法: +# cmake -B build # 配置（使用默认编译器） +# cmake --build build # 编译 +# cmake --build build --target dynamics_c # 只编译 C 引擎 +# +# 交叉编译: +# cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-linux.cmake +# cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-mingw.cmake +# cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-macos.cmake +# ============================================================================ + +cmake_minimum_required(VERSION 3.15) +project(DynamicsEngines + VERSION 1.0.0 + DESCRIPTION "Multi-language molecular dynamics simulation engines" + LANGUAGES C CXX Fortran +) + +# ── 全局编译选项 ────────────────────────────── +set(CMAKE_C_STANDARD 11) +set(CMAKE_C_STANDARD_REQUIRED ON) +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +# 优化等级（可被 -DCMAKE_BUILD_TYPE=Debug 覆盖） +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() + +# 统一编译选项 +set(COMMON_C_FLAGS -O3 -Wall -Wextra) +set(COMMON_CXX_FLAGS -O3 -Wall -Wextra) + +# ── 统一输出目录：所有引擎可执行文件输出到 engines//build/ ──── +# 这样 Python 的 run_engine() 可以在固定路径找到它们 +function(add_engine_target name lang src) + if(lang STREQUAL "C") + add_executable(${name} ${src}) + target_compile_options(${name} PRIVATE ${COMMON_C_FLAGS}) + target_link_libraries(${name} PRIVATE m) + elseif(lang STREQUAL "CXX") + add_executable(${name} ${src}) + target_compile_options(${name} PRIVATE ${COMMON_CXX_FLAGS}) + elseif(lang STREQUAL "Fortran") + enable_language(Fortran OPTIONAL) + if(NOT CMAKE_Fortran_COMPILER) + message(STATUS "Fortran compiler not found — skipping ${name}") + return() + endif() + add_executable(${name} ${src}) + endif() + + # 输出到 engines//build/（Python 调度器期望的位置） + set_target_properties(${name} PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/engines/${name}/build" + ) + # macOS 上禁止生成 .app 包（命令行工具不需要） + set_target_properties(${name} PROPERTIES + MACOSX_BUNDLE FALSE + ) +endfunction() + +# ============================================================================ +# 引擎目标定义 +# ============================================================================ + +# ── C 引擎 ──────────────────────────────────── +add_engine_target(dynamics_c C "engines/c/main.c") + +# ── C++ 引擎 ────────────────────────────────── +add_engine_target(dynamics_cpp CXX "engines/cpp/main.cpp") + +# ── Fortran 引擎 ────────────────────────────── +add_engine_target(dynamics_f90 Fortran "engines/fortran/main.f90") + +# ============================================================================ +# 自定义目标 +# ============================================================================ + +# 编译所有引擎 +add_custom_target(all-engines + DEPENDS dynamics_c dynamics_cpp dynamics_f90 +) + +# 显示引擎信息 +add_custom_target(info + COMMAND ${CMAKE_COMMAND} -E echo "=== Dynamics Engines ===" + COMMAND ${CMAKE_COMMAND} -E echo "C engine: dynamics_c" + COMMAND ${CMAKE_COMMAND} -E echo "C++ engine: dynamics_cpp" + COMMAND ${CMAKE_COMMAND} -E echo "Fortran engine: dynamics_f90" + COMMAND ${CMAKE_COMMAND} -E echo "" + COMMAND ${CMAKE_COMMAND} -E echo "Build: cmake --build ." + COMMAND ${CMAKE_COMMAND} -E echo "Single: cmake --build . --target dynamics_c" +) + +message(STATUS "DynamicsEngines ${PROJECT_VERSION}") +message(STATUS " Build type: ${CMAKE_BUILD_TYPE}") +message(STATUS " C compiler: ${CMAKE_C_COMPILER}") +message(STATUS " C++ compiler: ${CMAKE_CXX_COMPILER}") +if(CMAKE_Fortran_COMPILER) + message(STATUS " Fortran compiler: ${CMAKE_Fortran_COMPILER}") +else() + message(STATUS " Fortran compiler: (not found, skipped)") +endif() +message(STATUS "") +message(STATUS " Targets:") +message(STATUS " cmake --build . --target dynamics_c") +message(STATUS " cmake --build . --target dynamics_cpp") +message(STATUS " cmake --build . --target dynamics_f90") +message(STATUS " cmake --build . --target all-engines") diff --git a/INSTALL.md b/INSTALL.md new file mode 100644 index 0000000..4d237f5 --- /dev/null +++ b/INSTALL.md @@ -0,0 +1,568 @@ +# 安装指南 + +## 目录 + +- [快速安装（5 分钟）](#快速安装5-分钟) +- [详细步骤](#详细步骤) + - [1. 安装 Python](#1-安装-python) + - [2. 安装 Python 依赖](#2-安装-python-依赖) + - [3. 安装 C 编译器（可选）](#3-安装-c-编译器可选) + - [4. 安装 CMake（可选）](#4-安装-cmake可选) + - [5. 编译 C 引擎](#5-编译-c-引擎) + - [6. 安装 VisPy（可选）](#6-安装-vispy可选) +- [验证安装](#验证安装) +- [平台对照表](#平台对照表) +- [常见问题](#常见问题) + +--- + +## 快速安装（5 分钟） + +```bash +# 1. 安装 Python 依赖 +pip install -r requirements.txt + +# 2. 运行案例验证 +cd Dynamics +py -3 examples/case01/run_dynamics.py --no-plot +``` + +或手动逐个安装： + +```bash +pip install numpy pyyaml tqdm matplotlib +``` + +看到 `[run] 完成！输出目录: ...` 即安装成功。 + +--- + +## 详细步骤 + +### 1. 安装 Python + +**要求：Python 3.9 或更高版本。** + +| 系统 | 安装方式 | +|------|---------| +| **Windows** | 从 [python.org](https://www.python.org/downloads/) 下载安装包，安装时勾选 **"Add Python to PATH"** | +| **macOS** | `brew install python` 或从 python.org 下载 | +| **Linux (Ubuntu/Debian)** | `sudo apt install python3 python3-pip` | +| **Linux (Fedora)** | `sudo dnf install python3 python3-pip` | +| **Linux (Arch)** | `sudo pacman -S python python-pip` | + +验证： + +```bash +py -3 --version # Windows +# 或 +python3 --version # Linux/macOS +``` + +### 2. 安装 Python 依赖 + +```bash +# 一键安装所有必需 + 推荐依赖 +pip install -r requirements.txt + +# 或手动逐个安装 +pip install numpy pyyaml tqdm + +# 绘图（可选，用于生成轨迹图） +pip install matplotlib + +# 根据你的 pip 命令，可能是 pip3： +pip3 install numpy pyyaml tqdm matplotlib +``` + +> 如果遇到网络慢，使用国内镜像： +> ```bash +> pip install -i https://pypi.tuna.tsinghua.edu.cn/simple numpy pyyaml tqdm matplotlib +> ``` + +### 3. 安装 C 编译器（可选） + +> 仅当你想**编译 C/C++ 引擎**时才需要。如果只用默认的 Python 引擎可以跳过。 + +#### Windows — 下载安装 MSYS2 + MinGW64 + +MSYS2 是一个在 Windows 上提供 Linux 风格开发环境的软件包集合， +内含 MinGW64 编译器（gcc、g++、gfortran 等）。 + +##### 步骤 1：下载 MSYS2 + +访问 [msys2.org](https://www.msys2.org/)，点击页面上的 **"msys2-x86_64-…….exe"** 下载安装包。 + +> 国内用户如果下载慢，可尝试镜像： +> - 清华大学: https://mirrors.tuna.tsinghua.edu.cn/msys2/distrib/msys2-x86_64-latest.exe +> - 中科大: https://mirrors.ustc.edu.cn/msys2/distrib/msys2-x86_64-latest.exe + +##### 步骤 2：安装 MSYS2 + +运行下载的安装包： + +1. 一路点 **"Next"**（下一步） +2. 安装路径保持默认：`C:\msys64` +3. 勾选 **"Run MSYS2 now"**（立即运行），点 Finish + +##### 步骤 3：安装 MinGW64 编译器 + +安装完成后会弹出 **MSYS2 UCRT64 终端**（黑窗口）。在终端中执行： + +```bash +pacman -Syu +``` + +> 这步会更新 MSYS2 自身。更新完毕后终端可能自动关闭，请从开始菜单重新打开 **"MSYS2 UCRT64"**。 + +然后安装编译器： + +```bash +pacman -S mingw-w64-ucrt-x86_64-gcc +``` + +安装 g++（可选，编译 C++ 引擎时需要）： + +```bash +pacman -S mingw-w64-ucrt-x86_64-gcc +``` + +安装 gfortran（可选，编译 Fortran 引擎时需要）： + +```bash +pacman -S mingw-w64-ucrt-x86_64-gcc-fortran +``` + +安装 make（可选，直接用 Makefile 编译时需要）： + +```bash +pacman -S make +``` + +> 如果安装时提示选择软件包版本，直接回车选默认即可。 + +##### 步骤 4：添加环境变量 + +MSYS2 安装完成后，需要把编译器所在目录添加到系统的 **PATH** 环境变量中，这样在 CMD 或 PowerShell 中直接输入 `gcc` 就能调用。 + +**操作步骤：** + +1. 打开 **系统属性** → **高级** → **环境变量** + - 快捷键：Win + R，输入 `sysdm.cpl`，回车 → 高级 → 环境变量 + - 或：右键"此电脑" → 属性 → 高级系统设置 → 环境变量 + +2. 在 **"系统变量"** 列表中找到 `Path`，选中后点 **"编辑"** + +3. 点 **"新建"**，添加以下路径： + + ``` + C:\msys64\ucrt64\bin + ``` + + > 如果安装时修改了路径，请替换为实际安装路径。 + +4. 点 **"确定"** 关闭所有对话框 + +##### 步骤 5：验证安装 + +打开一个新的 **CMD** 或 **PowerShell** 窗口（不是 MSYS2 终端），运行： + +```bash +gcc --version +``` + +如果显示类似以下信息，说明安装成功： + +``` +gcc.exe (Rev13, Built by MSYS2 project) 15.2.0 +Copyright (C) 2025 Free Software Foundation, Inc. +``` + +如果提示 `'gcc' 不是内部或外部命令`，说明环境变量未生效： +- 检查是否添加了正确的路径（`C:\msys64\ucrt64\bin`） +- 检查是否打开了**新的**命令行窗口（旧的窗口不识别新的环境变量） +- 重启电脑使环境变量生效 + +##### 步骤 6：编译 C 引擎 + +```bash +cd Dynamics\engines\c +gcc -O3 -march=native -o build\dynamics_c.exe main.c -lm +``` + +> 如果 MSYS2 安装在非默认路径，编译时可能需要指定完整路径： +> ```bash +> C:\msys64\ucrt64\bin\gcc -O3 -o build\dynamics_c.exe main.c -lm +> ``` + +##### 完整安装 MSYS2 + 编译器一键脚本 + +将以下内容保存为 `install_msys2.bat`，以管理员身份运行： + +```batch +@echo off +echo 本脚本将安装 MSYS2 和 MinGW64 编译器 +echo 请确保已下载 MSYS2 安装包: https://www.msys2.org/ +echo. + +REM 检查是否已安装 +if exist "C:\msys64\ucrt64\bin\gcc.exe" ( + echo 检测到 gcc 已安装，跳过编译安装步骤 + goto check_path +) + +echo 请手动运行 MSYS2 UCRT64 终端，执行以下命令： +echo. +echo pacman -Syu +echo pacman -S mingw-w64-ucrt-x86_64-gcc +echo. +echo 安装完成后重新运行此脚本。 +pause +exit /b + +:check_path +echo 检查环境变量... +echo %PATH% | findstr /C:"C:\msys64\ucrt64\bin" >nul +if %errorlevel% equ 0 ( + echo 环境变量已设置 +) else ( + echo 警告: C:\msys64\ucrt64\bin 不在 PATH 中 + echo 请手动添加：系统属性 → 高级 → 环境变量 → Path → 新建 +) +echo. +gcc --version +echo. +echo 安装验证完成! +pause +``` + +#### macOS + +```bash +# 安装 Xcode Command Line Tools（包含 clang） +xcode-select --install + +# 或通过 Homebrew 安装 gcc +brew install gcc +``` + +#### Linux (Ubuntu/Debian) + +```bash +sudo apt install gcc g++ gfortran +``` + +#### Linux (Fedora) + +```bash +sudo dnf install gcc gcc-c++ gcc-gfortran +``` + +### 4. 安装 CMake（可选） + +> 仅当你想用 **CMake 方式编译引擎** 时才需要。如果只用 gcc 直接编译可跳过。 + +#### 下载安装 + +| 系统 | 安装方式 | +|------|---------| +| **Windows** | ① 从 [cmake.org/download](https://cmake.org/download/) 下载 `.msi` 安装包 → 安装时勾选 **"Add CMake to system PATH"**
② 或使用 winget：`winget install Kitware.CMake` | +| **macOS** | `brew install cmake`
或从 cmake.org 下载 `.dmg` 安装包 | +| **Linux (Ubuntu/Debian)** | `sudo apt install cmake` | +| **Linux (Fedora)** | `sudo dnf install cmake` | +| **Linux (Arch)** | `sudo pacman -S cmake` | + +#### 验证安装 + +打开新的终端，运行： + +```bash +cmake --version +``` + +预期输出： + +``` +cmake version 3.29.0 +CMake suite maintained and supported by Kitware (kitware.com/cmake). +``` + +> 如果提示 `'cmake' 不是内部或外部命令`，说明未添加到 PATH： +> - **Windows**：重新运行安装包，选择 **"Add CMake to the system PATH for all users"** +> - 或手动将 CMake 安装目录（如 `C:\Program Files\CMake\bin`）添加到环境变量 + +#### CMake 基本用法 + +CMake 采用 **源码外构建**（out-of-source build），所有编译中间文件生成在 `build/` 目录中， +不会污染源码目录。 + +```bash +# 1. 在项目根目录下配置（只需执行一次） +cd Dynamics +cmake -B build + +# 2. 编译（--build 会自动调用 make / nmake / MSBuild） +cmake --build build + +# 3. 只编译特定引擎 +cmake --build build --target dynamics_c # C 引擎 +cmake --build build --target dynamics_cpp # C++ 引擎 +cmake --build build --target dynamics_f90 # Fortran 引擎 +cmake --build build --target all-engines # 全部（默认） + +# 4. 指定编译配置 +cmake -B build -DCMAKE_BUILD_TYPE=Debug # Debug 模式（含符号表） +cmake -B build -DCMAKE_BUILD_TYPE=Release # Release 模式（-O3 优化） + +# 5. 清理后重新编译 +rm -rf build # Windows: rmdir /s /q build +cmake -B build +cmake --build build +``` + +#### CMake 在 Windows 上的注意事项 + +**情况 A：使用 MinGW（MSYS2）** + +```bash +# 确保 gcc 在 PATH 中 +gcc --version + +# 配置时指定 MinGW Makefiles 生成器 +cmake -B build -G "MinGW Makefiles" + +# 编译 +cmake --build build + +# 或者用 ninja（如果安装了 ninja） +cmake -B build -G "Ninja" +cmake --build build +``` + +**情况 B：使用 Visual Studio** + +如果安装了 Visual Studio，CMake 会自动检测并使用 MSVC 编译器： + +```bash +# 在 "Developer Command Prompt for VS" 中运行 +cmake -B build +cmake --build build --config Release +``` + +> 也可在 CMake GUI 中操作： +> 1. 打开 `cmake-gui` +> 2. 源码目录：`D:/.../Dynamics` +> 3. 构建目录：`D:/.../Dynamics/build` +> 4. 点 **Configure** → 选择 VC 或 MinGW → 点 **Generate** +> 5. 点 **Open Project** 在 VS 中打开，或直接点 **Build** + +#### CMake 交叉编译 + +从任意平台编译到其他平台（需安装对应交叉编译器）： + +```bash +# 编译到 Windows（需 mingw-w64） +cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-mingw.cmake +cmake --build build + +# 编译到 Linux（需 Linux 交叉编译器） +cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-linux.cmake +cmake --build build + +# 编译到 macOS（需 osxcross） +cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-macos.cmake +cmake --build build +``` + +### 5. 编译 C 引擎 + +#### 方式一：用 gcc 直接编译（最快） + +```bash +cd engines/c +gcc -O3 -march=native -o build/dynamics_c.exe main.c -lm +``` + +#### 方式二：用 CMake（跨平台推荐） + +```bash +# 配置 +cmake -B build + +# 只编译 C 引擎 +cmake --build build --target dynamics_c + +# 或编译全部引擎 +cmake --build build +``` + +#### 方式三：用 Python 脚本 + +```bash +py -3 build_engines.py --target dynamics_c +``` + +#### 验证 C 引擎 + +```bash +# 确认可执行文件已生成 +ls engines/c/build/dynamics_c.exe + +# 直接运行测试 +engines/c/build/dynamics_c.exe examples/case01/input examples/case01/output engines/c/param.json +``` + +> `ls` 在 Windows 上可用 `dir` 替代。 + +### 6. 安装 VisPy（可选） + +> 仅当你想观看 **3D 动画** 时才需要。 + +```bash +pip install vispy PyQt5 +``` + +验证： + +```bash +py -3 -c "from vispy import app; print('VisPy 安装成功')" +``` + +--- + +## 验证安装 + +### 基础验证（Python 引擎） + +```bash +cd Dynamics +py -3 examples/case01/run_dynamics.py --no-plot +``` + +预期输出末尾： + +``` +[run] 完成！输出目录: .../examples/case01/output +[run] 运行 python draw.py 查看动画。 +``` + +### 完整验证（包括绘图） + +```bash +py -3 examples/case01/run_dynamics.py +``` + +会在 `examples/case01/output/` 下生成： +- `trajectory.txt` — 完整轨迹 +- `display.txt` — 动画数据 +- `trajectory_plots.png` — 轨迹图 +- `dynamics.log` — 计算日志 + +### C 引擎验证 + +```bash +# 修改 parameters.yaml 中 engine: python → engine: c +# 或用命令行指定 +py -3 examples/case01/run_dynamics.py +``` + +预期输出中出现 `[C-engine]` 字样。 + +### 动画验证 + +```bash +py -3 draw.py examples/case01/output +``` + +会弹出一个 3D 窗口，显示原子运动。 + +--- + +## 平台对照表 + +| 组件 | Windows | macOS | Linux | +|------|---------|-------|-------| +| **Python** (必需) | python.org 安装包 | `brew install python` | `apt install python3` | +| **numpy** (必需) | `pip install numpy` | 同上 | 同上 | +| **PyYAML** (必需) | `pip install pyyaml` | 同上 | 同上 | +| **tqdm** (推荐) | `pip install tqdm` | 同上 | 同上 | +| **matplotlib** (可选) | `pip install matplotlib` | 同上 | 同上 | +| **C 编译器** (可选) | MSYS2 MinGW | Xcode CLT / brew gcc | `apt install gcc` | +| **CMake** (可选) | cmake.org / winget | `brew install cmake` | `apt install cmake` | +| **VisPy** (可选) | `pip install vispy PyQt5` | `pip install vispy PyQt5` | `pip install vispy` | +| **交叉编译→Windows** | — | `brew install mingw-w64` | `apt install mingw-w64` | +| **交叉编译→Linux** | MSYS2 中 `pacman -S mingw-w64-x86_64-linux-gnu` | `brew install x86_64-elf-gcc` | — | +| **交叉编译→macOS** | 需 osxcross | — | 需 osxcross | + +--- + +## 常见问题 + +### Q: `pip` 命令找不到？ + +Python 3 的包管理命令： + +```bash +# Windows +py -3 -m pip install numpy + +# Linux/macOS +python3 -m pip install numpy +``` + +### Q: `gcc` 命令找不到？ + +| 系统 | 解决 | +|------|------| +| **Windows (MSYS2)** | 将 `C:\msys64\ucrt64\bin` 添加到 PATH | +| **macOS** | 运行 `xcode-select --install` | +| **Linux** | `sudo apt install gcc` | + +### Q: 运行 `draw.py` 闪退？ + +确保已安装 VisPy 和 PyQt5： + +```bash +pip install vispy PyQt5 +``` + +如果仍闪退，直接运行查看错误信息： + +```bash +py -3 draw.py examples/case01/output +``` + +### Q: 网络慢，pip 安装失败？ + +使用国内镜像： + +```bash +pip install -i https://pypi.tuna.tsinghua.edu.cn/simple numpy pyyaml tqdm matplotlib vispy PyQt5 +``` + +### Q: 如何使用 C 引擎？ + +```yaml +# 1. 修改 examples/case01/input/parameters.yaml +engine: c + +# 2. 先编译 +cd engines/c +gcc -O3 -march=native -o build/dynamics_c.exe main.c -lm + +# 3. 运行 +cd ../.. +py -3 examples/case01/run_dynamics.py +``` + +### Q: 如何安装旧版 Python？ + +项目要求 Python 3.9+。如果系统自带 Python 2（某些旧版 Linux），请单独安装 Python 3： + +```bash +# Ubuntu +sudo apt install python3 python3-pip + +# 然后用 python3 和 pip3 代替 python 和 pip +python3 examples/case01/run_dynamics.py +``` diff --git a/README.md b/README.md index 9ec0103..db022f5 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,217 @@ -# dynamics +# Dynamics -动力学程序 \ No newline at end of file +多语言动力学数值积分与轨迹可视化框架。 + +支持 Python、C、C++、Fortran 多种语言实现同一套物理引擎，共享同一套可视化管线 +（抽帧、绘图、3D 动画），直接在 `parameters.yaml` 中切换引擎即可对比性能。 + +## 快速开始 + +```bash +# 1. 安装依赖 +pip install numpy pyyaml matplotlib tqdm + +# 2. 运行案例（计算 + 抽帧 + 绘图 + 动画） +py -3 examples/case01/run_dynamics.py + +# 3. 仅计算（跳过绘图，最小依赖） +py -3 examples/case01/run_dynamics.py --no-plot + +# 4. 仅播放动画（需先安装 VisPy） +pip install vispy PyQt5 +py -3 draw.py examples/case01/output +``` + +## 目录结构 + +``` +Dynamics/ +├── dynamics.py ← 统一入口：调度引擎、抽帧、绘图、动画 +├── compute.py ← Python 参考引擎 + 外部引擎调度 +├── draw.py ← VisPy 3D 动画（所有引擎共用） +├── plot_trajectory.py ← Matplotlib 静态绘图 +├── sample.py ← 抽帧生成显示数据 +├── export_web_data.py ← 导出 Web 可视化数据 +│ +├── engines/ ← 多语言计算引擎 +│ ├── c/ ← C 语言 (完整实现，已编译) +│ │ ├── main.c +│ │ ├── Makefile +│ │ └── build/dynamics_c.exe +│ ├── cpp/ ← C++ (模板) +│ │ └── main.cpp +│ └── fortran/ ← Fortran (模板) +│ └── main.f90 +│ +├── examples/ +│ └── case01/ +│ ├── input/ +│ │ ├── parameters.yaml ← 唯一配置文件 +│ │ ├── coord.txt +│ │ ├── connection.txt +│ │ └── bond.txt +│ └── output/ +│ ├── trajectory.txt ← 完整轨迹 +│ ├── display.txt ← 抽帧数据（动画用） +│ ├── trajectory_plots.png +│ └── dynamics.log ← 计算日志 +│ +└── .workbuddy/memory/ ← 工作记忆（AI 辅助开发用） +``` + +## 配置文件 (`parameters.yaml`) + +所有控制集中在一个文件里： + +```yaml +# ── 计算引擎 ────────────────────────────────── +engine: python # python | c | cpp | fortran + +# ── 流程控制（每步独立开关） ───────────────── +step_simulate: 1 # 运行物理模拟 +step_sample: 1 # 抽帧 +step_plot: 1 # 绘图 +step_animation: 0 # VisPy 3D 动画 + +# ── 物理参数 ────────────────────────────────── +G: [0.0, 0.0, -9.8] # 重力 +B: [0.0, 0.0, 0.0] # 阻尼 +method: leapfrog # explicit_euler | implicit_euler | midpoint | leapfrog +NT: 100000 +DT: 0.0001 +NSTEP: 1000 + +# ── 显示参数 ────────────────────────────────── +# 六个面的透明度，按 [-x,+x,-y,+y,-z,+z] 顺序 +alpha: [0.0, 0.0, 0.2, 0.2, 0.0, 0.0] +``` + +## 多语言引擎 + +### 一键切换 + +```yaml +# parameters.yaml 中改一行 +engine: c # 切换到 C 引擎 +engine: python # 切回 Python +``` + +### 编译 C 引擎 + +**方式一：CMake（推荐，跨平台）** + +```bash +# 安装 CMake (https://cmake.org/download/) +# 然后： +cmake -B build # 配置 +cmake --build build # 编译全部 +cmake --build build --target dynamics_c # 只编译 C 引擎 + +# 编译完成后： +py -3 examples/case01/run_dynamics.py # 配合 Python 调度器使用 +``` + +**方式二：直接 gcc（最快）** + +```bash +cd engines/c +gcc -O3 -march=native -o build/dynamics_c.exe main.c -lm +``` + +### 交叉编译 + +从任意平台编译到其他平台： + +```bash +# 编译到 Windows（需 mingw-w64） +cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-mingw.cmake +cmake --build build + +# 编译到 Linux（需 Linux 交叉编译器） +cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-linux.cmake +cmake --build build + +# 编译到 macOS（需 osxcross） +cmake -B build -DCMAKE_TOOLCHAIN_FILE=cmake/toolchain-macos.cmake +cmake --build build +``` + +或用 Python 一键脚本（自动调用 CMake）： + +```bash +py -3 build_engines.py # 本地编译全部 +py -3 build_engines.py --target dynamics_c # 只编译 C 引擎 +py -3 build_engines.py --clean # 清理后重新编译 +``` + +所有平台统一输出到 `engines/*/build/dynamics_c.exe`，Python 自动识别调用。 + +### 性能对比 (case01, 10000步) + +| 引擎 | 耗时 | 相对速度 | +|------|------|----------| +| Python | ~1.3 s | 1× | +| C | ~0.05 s | 26× | + +## 3D 动画功能 + +运行 `py -3 draw.py examples/case01/output` 后： + +| 操作 | 说明 | +|------|------| +| 鼠标拖拽 | 旋转视角 | +| 滚轮 | 缩放 | +| 点击 **`reset`** | 复位相机 + 从头播放 | +| 点击 **`info`** | 显示/隐藏信息面板 | +| 点击 **`axes`** | 显示/隐藏坐标轴 | +| `Q` / `E` 键 | 绕屏幕法线旋转 90° | + +所有原子用 tab10 调色板着色，成键用灰色线段连接。 + +## 依赖 + +| 用途 | 包 | 是否必需 | +|------|----|---------| +| 数值计算 | `numpy` | 是 | +| 配置解析 | `PyYAML` | 是 | +| 进度条 | `tqdm` | 推荐 | +| 静态绘图 | `matplotlib` | 可选 | +| 3D 动画 | `vispy` + `PyQt5` | 可选 | +| C/C++/Fortran 编译 | `cmake` (>= 3.15) | 可选（编译引擎时） | +| C 快速编译 | `gcc` (MSYS2/MinGW) | 可选 | + +```bash +pip install numpy pyyaml tqdm matplotlib vispy PyQt5 +``` + +## 工作原理 + +``` +parameters.yaml + │ + ▼ + dynamics.py ──→ run_from_config() 加载全局参数 + │ + ├── engine=python? → compute.py → run_simulation() + │ (Leapfrog / Euler) + └── engine=c? → subprocess → engines/c/dynamics_c.exe + (Leapfrog in C) + │ + ▼ + trajectory.txt (JSON, 统一格式) + │ + ▼ + sample.py → display.txt (抽帧数据) + │ + ├── plot_trajectory.py → trajectory_plots.png + └── draw.py → VisPy 3D 窗口 +``` + +所有引擎的输出都是相同格式的 `trajectory.txt`（JSON 嵌套数组）， +后续的抽帧、绘图、动画完全复用，每个引擎只需实现纯计算逻辑。 + +## 项目初衷 + +- **教学**：对比不同语言的数值计算性能 +- **验证**：确保 C/C++/Fortran 实现与 Python 参考实现结果一致 +- **可视化**：一套代码搞定所有语言的轨迹展示 diff --git a/build_engines.py b/build_engines.py new file mode 100644 index 0000000..abd3731 --- /dev/null +++ b/build_engines.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +""" +build_engines.py — 跨平台一键编译所有引擎。 + +用法: + py -3 build_engines.py # 本地编译 + py -3 build_engines.py --target c # 只编译 C 引擎 + py -3 build_engines.py --release # Release 模式 + py -3 build_engines.py --clean # 清理构建目录 + +依赖: cmake (>= 3.15) +""" + +import argparse +import os +import subprocess +import sys +import platform + +SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) +BUILD_DIR = os.path.join(SCRIPT_DIR, "build") + +def run(cmd, **kwargs): + print(f"[build] {' '.join(cmd)}") + result = subprocess.run(cmd, cwd=SCRIPT_DIR, **kwargs) + if result.returncode != 0: + sys.exit(result.returncode) + return result + +def main(): + parser = argparse.ArgumentParser(description="一键编译 Dynamics 引擎") + parser.add_argument("--target", "-t", default="all-engines", + help="编译目标: dynamics_c, dynamics_cpp, dynamics_f90, all-engines (默认)") + parser.add_argument("--release", action="store_true", help="Release 模式") + parser.add_argument("--debug", action="store_true", help="Debug 模式") + parser.add_argument("--clean", action="store_true", help="清理后重新编译") + parser.add_argument("--toolchain", help="交叉编译工具链文件路径") + args = parser.parse_args() + + # 检查 cmake + try: + subprocess.run(["cmake", "--version"], capture_output=True, check=True) + except (FileNotFoundError, subprocess.CalledProcessError): + print("[build] 错误: 未找到 cmake，请先安装:") + system = platform.system().lower() + if system == "windows": + print(" 下载安装: https://cmake.org/download/") + print(" 或: winget install Kitware.CMake") + elif system == "darwin": + print(" brew install cmake") + else: + print(" apt install cmake # Ubuntu/Debian") + print(" dnf install cmake # Fedora") + sys.exit(1) + + # 确定构建类型 + if args.debug: + build_type = "Debug" + elif args.release: + build_type = "Release" + else: + build_type = "Release" + + # 清理 + if args.clean and os.path.exists(BUILD_DIR): + import shutil + shutil.rmtree(BUILD_DIR) + print(f"[build] 已清理: {BUILD_DIR}") + + # cmake 配置 + config_cmd = ["cmake", "-B", BUILD_DIR, + f"-DCMAKE_BUILD_TYPE={build_type}"] + if args.toolchain: + config_cmd.append(f"-DCMAKE_TOOLCHAIN_FILE={args.toolchain}") + + # 并行编译 + nproc = os.cpu_count() or 4 + build_cmd = ["cmake", "--build", BUILD_DIR, + "--target", args.target, + "--parallel", str(nproc)] + + # 执行 + print(f"[build] 目标: {args.target}") + print(f"[build] 模式: {build_type}") + print(f"[build] 系统: {platform.system()} {platform.machine()}") + run(config_cmd) + run(build_cmd) + + print(f"\n[build] 完成!") + print(f"[build] 可执行文件位置: engines/*/build/") + +if __name__ == "__main__": + main() diff --git a/build_release_zip.py b/build_release_zip.py new file mode 100644 index 0000000..f7775b7 --- /dev/null +++ b/build_release_zip.py @@ -0,0 +1,36 @@ +""" +Build a lean Dynamics.zip package for the website download button. +""" + +from __future__ import annotations + +import zipfile +from pathlib import Path + + +HERE = Path(__file__).resolve().parent +ZIP_PATH = HERE.parent / "Dynamics-demo.zip" +ARC_ROOT = "Dynamics" + +INCLUDE_FILES = [ + "README.md", + "dynamics.py", + "compute.py", + "sample.py", + "draw.py", + "plot_trajectory.py", + "examples\\case01\\input\\parameters.yaml", + "examples\\case01\\input\\coord.txt", +] + + +def main() -> None: + with zipfile.ZipFile(ZIP_PATH, "w", compression=zipfile.ZIP_DEFLATED) as archive: + for relative_name in INCLUDE_FILES: + path = HERE / relative_name + archive.write(path, arcname=f"{ARC_ROOT}/{relative_name}") + print(f"[release] wrote {ZIP_PATH}") + + +if __name__ == "__main__": + main() diff --git a/cmake/toolchain-linux.cmake b/cmake/toolchain-linux.cmake new file mode 100644 index 0000000..5a078df --- /dev/null +++ b/cmake/toolchain-linux.cmake @@ -0,0 +1,16 @@ +# cmake/toolchain-linux.cmake +# 交叉编译 → Linux (x86_64) +# 需要安装 Linux 交叉编译器： +# Ubuntu/Debian: apt install gcc-x86-64-linux-gnu g++-x86-64-linux-gnu +# macOS: brew install x86_64-elf-gcc + +set(CMAKE_SYSTEM_NAME Linux) +set(CMAKE_SYSTEM_PROCESSOR x86_64) + +set(CMAKE_C_COMPILER x86_64-linux-gnu-gcc) +set(CMAKE_CXX_COMPILER x86_64-linux-gnu-g++) +set(CMAKE_Fortran_COMPILER x86_64-linux-gnu-gfortran) + +set(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM NEVER) +set(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY) +set(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY) diff --git a/cmake/toolchain-macos.cmake b/cmake/toolchain-macos.cmake new file mode 100644 index 0000000..fc273f3 --- /dev/null +++ b/cmake/toolchain-macos.cmake @@ -0,0 +1,15 @@ +# cmake/toolchain-macos.cmake +# 交叉编译 → macOS (x86_64) +# 需要安装 osxcross： +# https://github.com/tpoechtrager/osxcross +# 或使用 Apple 官方工具链（仅在 macOS 上原生编译） + +set(CMAKE_SYSTEM_NAME Darwin) +set(CMAKE_SYSTEM_PROCESSOR x86_64) + +set(CMAKE_C_COMPILER x86_64-apple-darwin-clang) +set(CMAKE_CXX_COMPILER x86_64-apple-darwin-clang++) + +set(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM NEVER) +set(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY) +set(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY) diff --git a/cmake/toolchain-mingw.cmake b/cmake/toolchain-mingw.cmake new file mode 100644 index 0000000..f41498c --- /dev/null +++ b/cmake/toolchain-mingw.cmake @@ -0,0 +1,17 @@ +# cmake/toolchain-mingw.cmake +# 交叉编译 → Windows (x86_64) MinGW +# 需要安装 mingw-w64： +# Ubuntu/Debian: apt install mingw-w64 +# macOS: brew install mingw-w64 +# Fedora: dnf install mingw64-gcc + +set(CMAKE_SYSTEM_NAME Windows) +set(CMAKE_SYSTEM_PROCESSOR x86_64) + +set(CMAKE_C_COMPILER x86_64-w64-mingw32-gcc) +set(CMAKE_CXX_COMPILER x86_64-w64-mingw32-g++) +set(CMAKE_Fortran_COMPILER x86_64-w64-mingw32-gfortran) + +set(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM NEVER) +set(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY) +set(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY) diff --git a/compute.py b/compute.py new file mode 100644 index 0000000..f297418 --- /dev/null +++ b/compute.py @@ -0,0 +1,1032 @@ +""" +compute.py +---------- +纯 Python 计算脚本，不依赖 VisPy。 + +功能： + 1. 运行 NT 步物理模拟（ kinematics / dynamics 等运动模式） + 2. 将每一步的 (x, y, z, vx, vy, vz) 保存到 output/trajectory.txt + 3. 同时保存所有模拟参数元数据 +""" + +import json +import numpy as np +import os +import platform +import subprocess +import sys +import time +import datetime +from tqdm import trange + +# =========================================================================== +# 物理参数（必须与 draw.py 加载的边界常量保持一致） +# =========================================================================== +# 全局参数变量，将在 load_parameters 中填充 +box_a = None +alpha = None +COORD_FILE = None +X0 = None +Y0 = None +Z0 = None +VX0 = None +VY0 = None +VZ0 = None +M = None +ATOM_IDS = None +ATOM_MASSES = None +ATOM_RADII = None +ATOM_POSITIONS = None +ATOM_VELOCITIES = None +ATOM_FIXED = None +BOND_CONNECTION_FILE = None +BOND_PARAMETER_FILE = None +BOND_PAIRS = None +BOND_NAMES = None +BOND_STIFFNESS = None +BOND_REST_LENGTHS = None +PLOT_ATOM_ID = None +PLOT_ATOM_ROW = None +G = None +B = None +METHOD = None +NT = None +DT = None +NSTEP = None +warmup_steps = None # 预热步数（跳过不保存） +sample_start = None # 抽帧起始索引 +sample_end = None # 抽帧结束索引 +ball_radius = None +ball_color_r = None +ball_color_g = None +ball_color_b = None +box_color_r = None +box_color_g = None +box_color_b = None + +# 派生边界（根据 box_a 计算） +X_MIN = None +X_MAX = None +Y_MIN = None +Y_MAX = None +Z_MIN = None +Z_MAX = None + + +def _to_text_value(value): + """Convert numpy-heavy objects into JSON-friendly plain Python values.""" + if isinstance(value, np.ndarray): + return value.tolist() + if isinstance(value, (np.integer,)): + return int(value) + if isinstance(value, (np.floating,)): + return float(value) + if isinstance(value, dict): + return {key: _to_text_value(val) for key, val in value.items()} + if isinstance(value, (list, tuple)): + return [_to_text_value(item) for item in value] + return value + + +def _from_text_value(value): + """Convert nested numeric lists back into numpy arrays when appropriate.""" + if isinstance(value, list): + if not value: + return np.array([], dtype=np.float64) + if all(not isinstance(item, (list, dict)) for item in value): + if all(isinstance(item, (int, float, bool)) for item in value): + return np.array(value) + return value + if all(isinstance(item, list) for item in value): + return np.array(value) + return [_from_text_value(item) for item in value] + if isinstance(value, dict): + return {key: _from_text_value(val) for key, val in value.items()} + return value + + +def save_text_data(path, data): + """Save structured simulation data as formatted JSON text.""" + os.makedirs(os.path.dirname(path), exist_ok=True) + with open(path, "w", encoding="utf-8") as f: + json.dump(_to_text_value(data), f, ensure_ascii=False, indent=2) + return path + + +def load_text_data(path): + """Load structured simulation data from JSON text.""" + with open(path, "r", encoding="utf-8") as f: + data = json.load(f) + return _from_text_value(data) + + +def get_output_dir(base_dir=None): + """Return the output directory used for generated artifacts.""" + override = os.environ.get("DYNAMICS_OUTPUT_DIR") + if override: + output_dir = os.path.abspath(override) + else: + if base_dir is None: + base_dir = os.path.dirname(os.path.abspath(__file__)) + output_dir = os.path.join(base_dir, "output") + os.makedirs(output_dir, exist_ok=True) + return output_dir + + +def get_input_dir(base_dir=None): + """Return the input directory used for configuration and source data.""" + override = os.environ.get("DYNAMICS_INPUT_DIR") + if override: + input_dir = os.path.abspath(override) + else: + if base_dir is None: + base_dir = os.path.dirname(os.path.abspath(__file__)) + input_dir = os.path.join(base_dir, "input") + os.makedirs(input_dir, exist_ok=True) + return input_dir + + +def load_coord_file(coord_path): + """Load atom ids, masses, radii, coordinates, and velocities.""" + if not os.path.exists(coord_path): + raise FileNotFoundError(f"坐标文件不存在: {coord_path}") + + rows = [] + with open(coord_path, "r", encoding="utf-8") as f: + 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): + raise ValueError( + f"坐标文件表头应为: {' '.join(expected)}，实际为: {' '.join(header)}") + + for line_no, line in enumerate(f, start=2): + stripped = line.strip() + if not stripped or stripped.startswith("#"): + continue + 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: + raise ValueError(f"{coord_path}:{line_no} 应有 12 列，实际为 {len(parts)} 列") + rows.append(parts) + + if not rows: + raise ValueError(f"坐标文件没有原子数据: {coord_path}") + + atom_ids = np.array([int(row[0]) for row in rows], dtype=np.int64) + masses = np.array([float(row[1]) for row in rows], dtype=np.float64) + radii = np.array([float(row[2]) for row in rows], dtype=np.float64) + positions = np.array([[float(row[3]), float(row[4]), float(row[5])] for row in rows], + dtype=np.float64) + 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: + fixed = np.array([[int(row[9]), int(row[10]), int(row[11])] for row in rows], + dtype=np.int64) + + if np.any(masses <= 0): + raise ValueError("坐标文件中的质量必须为正数") + if np.any(radii <= 0): + raise ValueError("坐标文件中的半径必须为正数") + + print(f"[compute] 已加载坐标文件: {coord_path}，原子数={len(atom_ids)}") + return atom_ids, masses, radii, positions, velocities, fixed + + +def load_bond_parameters(bond_path): + """Load bond stiffness and optional rest length by bond name. + + File format (表头): + bond_name k # 2 列：从初始坐标计算键长 + bond_name k rest_length # 3 列：显式指定键长 + """ + if not os.path.exists(bond_path): + raise FileNotFoundError(f"键参数文件不存在: {bond_path}") + + bond_map = {} + with open(bond_path, "r", encoding="utf-8") as f: + header = f.readline().strip().split() + + if header == ["bond_name", "k"]: + has_rest_length = False + elif header == ["bond_name", "k", "rest_length"]: + has_rest_length = True + else: + expected = ["bond_name", "k"] + (["rest_length"] if len(header) == 3 else []) + raise ValueError( + f"键参数文件表头应为: {' '.join(expected)}，实际为: {' '.join(header)}") + + ncols = 3 if has_rest_length else 2 + + for line_no, line in enumerate(f, start=2): + stripped = line.strip() + if not stripped or stripped.startswith("#"): + continue + parts = stripped.split() + if len(parts) != ncols: + raise ValueError( + f"{bond_path}:{line_no} 应有 {ncols} 列，实际为 {len(parts)} 列") + bond_name = parts[0] + stiffness = float(parts[1]) + if stiffness < 0: + raise ValueError(f"{bond_path}:{line_no} 劲度系数必须非负") + + rest_length = None + if has_rest_length: + rest_length = float(parts[2]) + if rest_length <= 0: + raise ValueError(f"{bond_path}:{line_no} 键长必须为正数") + + bond_map[bond_name] = {"stiffness": stiffness, "rest_length": rest_length} + + if not bond_map: + raise ValueError(f"键参数文件没有有效数据: {bond_path}") + return bond_map + + +def load_bond_connections(connection_path, atom_ids, atom_positions, bond_map): + """Load bonded atom pairs and derive rest lengths. + + 键长优先级： + 1. bond_map 中显式指定的 rest_length（bond.txt 第三列） + 2. 否则从初始坐标计算 + """ + if not os.path.exists(connection_path): + raise FileNotFoundError(f"成键连接文件不存在: {connection_path}") + + atom_index = {int(atom_id): idx for idx, atom_id in enumerate(atom_ids)} + pairs = [] + bond_names = [] + stiffness = [] + rest_lengths = [] + + with open(connection_path, "r", encoding="utf-8") as f: + header = f.readline().strip().split() + expected = ["n1", "n2", "bond_name"] + if header != expected: + raise ValueError( + f"成键连接文件表头应为: {' '.join(expected)}，实际为: {' '.join(header)}") + + for line_no, line in enumerate(f, start=2): + stripped = line.strip() + if not stripped or stripped.startswith("#"): + continue + parts = stripped.split() + if len(parts) != 3: + raise ValueError(f"{connection_path}:{line_no} 应有 3 列，实际为 {len(parts)} 列") + + atom_1 = int(parts[0]) + atom_2 = int(parts[1]) + bond_name = parts[2] + + if atom_1 not in atom_index or atom_2 not in atom_index: + raise ValueError( + f"{connection_path}:{line_no} 中的原子序号 {atom_1}, {atom_2} 不在坐标文件里") + if bond_name not in bond_map: + raise ValueError( + f"{connection_path}:{line_no} 使用了未定义的键类型 {bond_name}") + + idx_1 = atom_index[atom_1] + idx_2 = atom_index[atom_2] + + bond_entry = bond_map[bond_name] + rest_length = bond_entry["rest_length"] + if rest_length is None: + # 未指定键长时从初始坐标计算 + delta = atom_positions[idx_2] - atom_positions[idx_1] + rest_length = float(np.linalg.norm(delta)) + + pairs.append((idx_1, idx_2)) + bond_names.append(bond_name) + stiffness.append(bond_entry["stiffness"]) + rest_lengths.append(rest_length) + + if not pairs: + return ( + np.zeros((0, 2), dtype=np.int64), + [], + np.zeros(0, dtype=np.float64), + np.zeros(0, dtype=np.float64), + ) + + return ( + np.array(pairs, dtype=np.int64), + bond_names, + np.array(stiffness, dtype=np.float64), + np.array(rest_lengths, dtype=np.float64), + ) + + +def parse_gravity_vector(value): + """Parse gravity into a 3D acceleration vector. + + Backward compatibility: + - scalar `G: 9.8` -> [0, 0, -9.8] + - vector `G: [gx, gy, gz]` -> [gx, gy, gz] + """ + if isinstance(value, (int, float, np.integer, np.floating)): + return np.array([0.0, 0.0, -float(value)], dtype=np.float64) + vector = np.asarray(value, dtype=np.float64) + if vector.shape != (3,): + raise ValueError(f"G 必须是长度为 3 的分量数组，实际为 {value}") + return vector + + +def parse_damping_vector(value): + """Parse damping into a 3D component vector. + + Backward compatibility: + - scalar `B: 0.5` -> [0.5, 0.5, 0.5] + - vector `B: [bx, by, bz]` -> [bx, by, bz] + """ + if isinstance(value, (int, float, np.integer, np.floating)): + scalar = float(value) + return np.array([scalar, scalar, scalar], dtype=np.float64) + vector = np.asarray(value, dtype=np.float64) + if vector.shape != (3,): + raise ValueError(f"B 必须是长度为 3 的分量数组，实际为 {value}") + return vector + + +def run_from_config(config, out_dir=None): + """直接接收配置字典，运行模拟并保存结果。 + + config 字典结构： + box_a, alpha, coord_file, connection_file, bond_file, plot_atom, G, B, method, NT, DT, NSTEP, + warmup_steps, sample_start, sample_end（可选）, + ball_radius, ball_color_r/g/b, box_color_r/g/b（可选） + + out_dir: 输出目录，默认同脚本所在目录 + 返回：(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz) + """ + global box_a, alpha, COORD_FILE, X0, Y0, Z0, VX0, VY0, VZ0, M, G, B, METHOD, NT, DT, NSTEP + global PLOT_ATOM_ID, PLOT_ATOM_ROW + global ATOM_IDS, ATOM_MASSES, ATOM_RADII, ATOM_POSITIONS, ATOM_VELOCITIES, ATOM_FIXED + global BOND_CONNECTION_FILE, BOND_PARAMETER_FILE, BOND_PAIRS, BOND_NAMES, BOND_STIFFNESS, BOND_REST_LENGTHS + global X_MIN, X_MAX, Y_MIN, Y_MAX, Z_MIN, Z_MAX + 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 + + box_a = float(config.get("box_a", 10.0)) + raw_alpha = config.get("alpha", 0.2) + if isinstance(raw_alpha, (list, tuple)): + alpha = [float(a) for a in raw_alpha] + else: + alpha = float(raw_alpha) + COORD_FILE = str(config.get("coord_file", os.path.join("input", "coord.txt"))) + coord_path = COORD_FILE + if out_dir is not None and not os.path.isabs(coord_path): + coord_path = os.path.join(out_dir, coord_path) + (ATOM_IDS, ATOM_MASSES, ATOM_RADII, ATOM_POSITIONS, + ATOM_VELOCITIES, ATOM_FIXED) = load_coord_file(coord_path) + BOND_CONNECTION_FILE = str(config.get("connection_file", os.path.join("input", "connection.txt"))) + BOND_PARAMETER_FILE = str(config.get("bond_file", os.path.join("input", "bond.txt"))) + connection_path = BOND_CONNECTION_FILE + bond_path = BOND_PARAMETER_FILE + if out_dir is not None and not os.path.isabs(connection_path): + connection_path = os.path.join(out_dir, connection_path) + if out_dir is not None and not os.path.isabs(bond_path): + bond_path = os.path.join(out_dir, bond_path) + bond_map = load_bond_parameters(bond_path) + (BOND_PAIRS, BOND_NAMES, BOND_STIFFNESS, + BOND_REST_LENGTHS) = load_bond_connections( + connection_path, ATOM_IDS, ATOM_POSITIONS, bond_map) + PLOT_ATOM_ID = int(config.get("plot_atom", int(ATOM_IDS[0]))) + matches = np.where(ATOM_IDS == PLOT_ATOM_ID)[0] + if len(matches) == 0: + raise ValueError(f"plot_atom={PLOT_ATOM_ID} 不在坐标文件原子序号中") + PLOT_ATOM_ROW = int(matches[0]) + + M = float(ATOM_MASSES[PLOT_ATOM_ROW]) + X0, Y0, Z0 = [float(v) for v in ATOM_POSITIONS[PLOT_ATOM_ROW]] + VX0, VY0, VZ0 = [float(v) for v in ATOM_VELOCITIES[PLOT_ATOM_ROW]] + G = parse_gravity_vector(config["G"]) + B = parse_damping_vector(config["B"]) + METHOD = normalize_method_name(config.get("method", "explicit_euler")) + NT = int(config["NT"]) + DT = float(config["DT"]) + NSTEP = int(config["NSTEP"]) + + # 步骤控制参数（可选，有默认值） + warmup_steps = int(config.get("warmup_steps", 0)) + sample_start = config.get("sample_start") # None 表示从头开始 + sample_end = config.get("sample_end") # None 表示到末尾 + + X_MIN = -box_a; X_MAX = box_a + Y_MIN = -box_a; Y_MAX = box_a + Z_MIN = -box_a; Z_MAX = box_a + + ball_radius = float(config.get("ball_radius", ATOM_RADII[PLOT_ATOM_ROW])) + ball_color_r = float(config.get("ball_color_r", 0.9)) + ball_color_g = float(config.get("ball_color_g", 0.2)) + ball_color_b = float(config.get("ball_color_b", 0.2)) + box_color_r = float(config.get("box_color_r", 0.8)) + box_color_g = float(config.get("box_color_g", 0.8)) + box_color_b = float(config.get("box_color_b", 0.85)) + + print(f"[compute] 使用算法: {METHOD}") + print(f"[compute] 已加载成键信息: {len(BOND_PAIRS)} 条键") + if config.get("_skip_run", False): + return None, None, None, None, None, None + t_start = time.time() + t_start_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = run_simulation() + t_end = time.time() + t_end_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + elapsed = t_end - t_start + + # 写入 dynamics.log + log_dir = get_output_dir(out_dir) + log_path = os.path.join(log_dir, "dynamics.log") + with open(log_path, "w", encoding="utf-8") as f: + f.write("=" * 50 + "\n") + f.write("Dynamics 计算日志\n") + f.write("=" * 50 + "\n") + f.write(f"计算开始时间: {t_start_str}\n") + f.write(f"计算结束时间: {t_end_str}\n") + f.write(f"计算耗时: {elapsed:.3f} s\n") + f.write("-" * 50 + "\n") + f.write("计算参数:\n") + f.write(f" box_a: {box_a}\n") + _alpha_str = alpha if isinstance(alpha, str) else (str(alpha) if isinstance(alpha, list) else str(alpha)) + f.write(f" alpha: {_alpha_str}\n") + f.write(f" coord_file: {COORD_FILE}\n") + f.write(f" method: {METHOD}\n") + f.write(f" NT: {NT}\n") + f.write(f" DT: {DT}\n") + f.write(f" NSTEP: {NSTEP}\n") + f.write(f" G: ({G[0]:.2f}, {G[1]:.2f}, {G[2]:.2f})\n") + f.write(f" B: ({B[0]:.2f}, {B[1]:.2f}, {B[2]:.2f})\n") + f.write(f" warmup: {warmup_steps}\n") + f.write(f" 原子数: {len(ATOM_IDS)}\n") + f.write(f" 键数: {len(BOND_PAIRS)}\n") + f.write("-" * 50 + "\n") + f.write("初始状态 (plot_atom):\n") + f.write(f" position: ({X0:.6f}, {Y0:.6f}, {Z0:.6f})\n") + f.write(f" velocity: ({VX0:.6f}, {VY0:.6f}, {VZ0:.6f})\n") + f.write(f" mass: {M}\n") + f.write("=" * 50 + "\n") + print(f"[compute] 日志已保存至: {log_path}") + print(f"[compute] 计算耗时: {elapsed:.3f} s") + + return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz + + +def run_engine(engine, input_dir, output_dir, config): + """调用外部计算引擎（C/C++/Fortran），生成 trajectory.txt。 + + Args: + engine: 引擎名称 ("c", "cpp", "fortran") + input_dir: 输入目录（含 coord.txt, connection.txt, bond.txt） + output_dir: 输出目录（轨迹文件写入位置） + config: YAML 配置字典 + Returns: + (traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz) + """ + script_dir = os.path.dirname(os.path.abspath(__file__)) + system = platform.system().lower() + engine_map = { + "c": "engines/c/build/dynamics_c", + "cpp": "engines/cpp/build/dynamics_cpp", + "fortran": "engines/fortran/build/dynamics_f90", + } + if engine not in engine_map: + raise ValueError(f"不支持的引擎: {engine}，可选: {list(engine_map.keys())}") + + engine_rel = engine_map[engine] + engine_path = os.path.join(script_dir, engine_rel) + + # 自动检测可执行文件后缀和平台专用版本 + candidates = [ + engine_path, # 无后缀 + engine_path + ".exe", # Windows .exe + engine_path + f"_{system}.exe", # 平台专用 (c_linux.exe, c_darwin.exe) + ] + found = None + for p in candidates: + if os.path.exists(p): + found = p + break + if found is None: + raise FileNotFoundError( + f"引擎可执行文件不存在: 尝试了 {candidates}\n" + f"请先编译: cd engines/{engine} && make\n" + f"或安装交叉编译器后: cd engines/{engine} && make {system}") + + # 构造 param.json（数值参数） + G = parse_gravity_vector(config.get("G", [0, 0, -9.8])) + B = parse_damping_vector(config.get("B", [0, 0, 0])) + param_json = { + "box_a": float(config.get("box_a", 10.0)), + "NT": int(config["NT"]), + "DT": float(config["DT"]), + "NSTEP": int(config.get("NSTEP", 1)), + "warmup_steps": int(config.get("warmup_steps", 0)), + "G": [float(v) for v in G], + "B": [float(v) for v in B], + } + param_path = os.path.join(script_dir, "engines", engine, "param.json") + os.makedirs(os.path.dirname(param_path), exist_ok=True) + with open(param_path, "w", encoding="utf-8") as f: + json.dump(param_json, f, indent=2) + + # 确保输出目录存在 + os.makedirs(output_dir, exist_ok=True) + + print(f"[compute] 引擎: {engine} → {os.path.basename(engine_path)}") + print(f"[compute] input: {input_dir}") + print(f"[compute] output: {output_dir}") + + 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}") + + # 加载输出的 trajectory.txt + traj_path = os.path.join(os.path.abspath(output_dir), "trajectory.txt") + if not os.path.exists(traj_path): + raise FileNotFoundError(f"引擎未生成 trajectory.txt: {traj_path}") + data = load_text_data(traj_path) + traj_x = data["traj_x"] + traj_y = data["traj_y"] + traj_z = data["traj_z"] + traj_vx = data["traj_vx"] + traj_vy = data["traj_vy"] + traj_vz = data["traj_vz"] + + # 写入日志文件 + log_path = os.path.join(output_dir, "dynamics.log") + with open(log_path, "w", encoding="utf-8") as f: + f.write("=" * 50 + "\n") + f.write("Dynamics 计算日志\n") + f.write("=" * 50 + "\n") + f.write(f"引擎: {engine}\n") + f.write(f"计算开始时间: {t_start_str}\n") + f.write(f"计算结束时间: {t_end_str}\n") + f.write(f"计算耗时: {elapsed:.3f} s\n") + + print(f"[compute] 引擎完成: {len(traj_x)} 步, {traj_x.shape[1]} 原子, {elapsed:.3f} s") + return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz + + +def save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, out_dir=None): + """将轨迹数据保存到结构化 txt 文件（含所有参数元数据）。""" + if out_dir is None: + out_dir = os.path.dirname(os.path.abspath(__file__)) + out_path = os.path.join(get_output_dir(out_dir), "trajectory.txt") + + # 记录的实际步数（扣除预热后） + record_steps = len(traj_x) + sample_start_to_save = -1 if sample_start is None else int(sample_start) + sample_end_to_save = -1 if sample_end is None else int(sample_end) + + payload = { + "traj_x": traj_x, "traj_y": traj_y, "traj_z": traj_z, + "traj_vx": traj_vx, "traj_vy": traj_vy, "traj_vz": traj_vz, + "NT": record_steps, "DT": DT, "NSTEP": NSTEP, + "method": METHOD, + "coord_file": COORD_FILE, + "connection_file": BOND_CONNECTION_FILE, + "bond_file": BOND_PARAMETER_FILE, + "atom_ids": ATOM_IDS, + "atom_masses": ATOM_MASSES, + "atom_radii": ATOM_RADII, + "atom_positions": ATOM_POSITIONS, + "atom_velocities": ATOM_VELOCITIES, + "atom_fixed": ATOM_FIXED, + "bond_pairs": BOND_PAIRS, + "bond_names": BOND_NAMES, + "bond_stiffness": BOND_STIFFNESS, + "bond_rest_lengths": BOND_REST_LENGTHS, + "G": G.tolist() if G is not None else None, + "B": B.tolist() if B is not None else None, + "plot_atom_id": PLOT_ATOM_ID, + "plot_atom_row": PLOT_ATOM_ROW, + "warmup_steps": warmup_steps or 0, + "sample_start": sample_start_to_save, + "sample_end": sample_end_to_save, + "X_MIN": X_MIN, "X_MAX": X_MAX, + "Y_MIN": Y_MIN, "Y_MAX": Y_MAX, + "Z_MIN": Z_MIN, "Z_MAX": Z_MAX, + "X0": X0, "Y0": Y0, "Z0": Z0, + "VX0": VX0, "VY0": VY0, "VZ0": VZ0, + "M": M, + "alpha": alpha, + "ball_radius": ball_radius, + "ball_color_r": ball_color_r, + "ball_color_g": ball_color_g, + "ball_color_b": ball_color_b, + "box_color_r": box_color_r, + "box_color_g": box_color_g, + "box_color_b": box_color_b, + } + save_text_data(out_path, payload) + print(f"[compute] 轨迹数据已保存至: {out_path}") + return out_path + + +def load_parameters(param_file): + """从文本文件加载参数，文件格式：参数 = 值""" + global box_a, alpha, X0, Y0, Z0, VX0, VY0, VZ0, M, G, B, METHOD, NT, DT, NSTEP + global X_MIN, X_MAX, Y_MIN, Y_MAX, Z_MIN, Z_MAX + global ball_radius, ball_color_r, ball_color_g, ball_color_b + global box_color_r, box_color_g, box_color_b + + print(f"[compute] 正在加载参数文件: {param_file}") + with open(param_file, 'r', encoding='utf-8') as f: + lines = f.readlines() + + for line in lines: + line = line.strip() + if not line or line.startswith('#'): + continue + if '=' not in line: + continue + # 分割键值，并去除行内注释 + key, value = line.split('=', 1) + key = key.strip() + value = value.strip() + # 去除值中可能存在的注释（# 之后的部分） + if '#' in value: + value = value.split('#', 1)[0].strip() + # 尝试转换为 int 或 float + try: + if '.' in value or 'e' in value.lower(): + value = float(value) + else: + value = int(value) + except ValueError: + pass # 保持字符串 + + # 设置全局变量 + if key in globals(): + globals()[key] = value + else: + print(f"警告: 未知参数 '{key}'，忽略") + + # 计算派生边界 + if box_a is not None: + X_MIN = -box_a + X_MAX = box_a + Y_MIN = -box_a + Y_MAX = box_a + Z_MIN = -box_a + Z_MAX = box_a + else: + print("错误: 未找到 box_a 参数") + sys.exit(1) + + # 检查必需参数 + required = ['alpha', 'X0', 'Y0', 'Z0', 'VX0', 'VY0', 'VZ0', 'M', 'G', 'B', 'NT', 'DT', 'NSTEP'] + for param in required: + if globals()[param] is None: + print(f"错误: 未找到必需参数 '{param}'") + sys.exit(1) + + METHOD = normalize_method_name(METHOD or "explicit_euler") + + # 设置绘图参数的默认值（如果未提供） + if ball_radius is None: + ball_radius = 0.28 + if ball_color_r is None: + ball_color_r = 0.9 + if ball_color_g is None: + ball_color_g = 0.2 + if ball_color_b is None: + ball_color_b = 0.2 + if box_color_r is None: + box_color_r = 0.8 + if box_color_g is None: + box_color_g = 0.8 + if box_color_b is None: + box_color_b = 0.85 + + print(f"[compute] 参数加载完成: box_a={box_a}, alpha={alpha}, NT={NT}, DT={DT}, method={METHOD}") + +def normalize_method_name(method): + """Normalize method aliases from YAML/text config.""" + aliases = { + "explicit_euler": "explicit_euler", + "explicit": "explicit_euler", + "euler_explicit": "explicit_euler", + "显式欧拉": "explicit_euler", + "显式欧拉法": "explicit_euler", + "implicit_euler": "implicit_euler", + "implicit": "implicit_euler", + "euler_implicit": "implicit_euler", + "隐式欧拉": "implicit_euler", + "隐式欧拉法": "implicit_euler", + "midpoint": "midpoint", + "mid_point": "midpoint", + "midpoint_method": "midpoint", + "中点法": "midpoint", + "中点差分法": "midpoint", + "leapfrog": "leapfrog", + "leap_frog": "leapfrog", + "蛙跳法": "leapfrog", + } + key = str(method).strip().lower() + if key not in aliases: + valid = ", ".join(["explicit_euler", "implicit_euler", "midpoint", "leapfrog"]) + raise ValueError(f"未知算法 '{method}'。可选值: {valid}") + return aliases[key] + + +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 + """ + fx = m * g[0] - b[0] * vx + fy = m * g[1] - b[1] * vy + fz = m * g[2] - b[2] * vz + + if 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] + dz = z[idx_2] - z[idx_1] + dist = np.sqrt(dx * dx + dy * dy + dz * dz) + if dist <= 1e-12: + continue + + stretch = dist - BOND_REST_LENGTHS[bond_idx] + force_mag = BOND_STIFFNESS[bond_idx] * stretch + ux = dx / dist + uy = dy / dist + uz = dz / dist + + fx_bond = force_mag * ux + fy_bond = force_mag * uy + fz_bond = force_mag * uz + + fx[idx_1] += fx_bond + fy[idx_1] += fy_bond + fz[idx_1] += fz_bond + fx[idx_2] -= fx_bond + fy[idx_2] -= fy_bond + fz[idx_2] -= fz_bond + + return fx, fy, fz + + +def compute_acceleration(x, y, z, vx, vy, vz, m, g, b): + """Compute acceleration from the shared force model.""" + fx, fy, fz = compute_force(x, y, z, vx, vy, vz, m, g, b) + return fx / m, fy / m, fz / m + + +def Explicit_Euler_Method(x, y, z, vx, vy, vz, dt, m, g, b): + ax, ay, az = compute_acceleration(x, y, z, vx, vy, vz, m, g, b) + x = x + vx * dt + y = y + vy * dt + z = z + vz * dt + vx = vx + ax * dt + vy = vy + ay * dt + vz = vz + az * dt + return x, y, z, vx, vy, vz + + +def Implicit_Euler_Method(x, y, z, vx, vy, vz, dt, m, g, b): + gamma_x = b[0] / m + gamma_y = b[1] / m + gamma_z = b[2] / m + vx_next = (vx + g[0] * dt) / (1.0 + gamma_x * dt) + vy_next = (vy + g[1] * dt) / (1.0 + gamma_y * dt) + vz_next = (vz + g[2] * dt) / (1.0 + gamma_z * dt) + ax_next, ay_next, az_next = compute_acceleration( + x, y, z, vx_next, vy_next, vz_next, m, g, b) + + vx = vx + ax_next * dt + vy = vy + ay_next * dt + vz = vz + az_next * dt + x = x + vx * dt + y = y + vy * dt + z = z + vz * dt + return x, y, z, vx, vy, vz + + +def Midpoint_Method(x, y, z, vx, vy, vz, dt, m, g, b): + ax, ay, az = compute_acceleration(x, y, z, vx, vy, vz, m, g, b) + x_mid = x + 0.5 * vx * dt + y_mid = y + 0.5 * vy * dt + z_mid = z + 0.5 * vz * dt + vx_mid = vx + 0.5 * ax * dt + vy_mid = vy + 0.5 * ay * dt + vz_mid = vz + 0.5 * az * dt + + x = x + vx_mid * dt + y = y + vy_mid * dt + z = z + vz_mid * dt + + ax_mid, ay_mid, az_mid = compute_acceleration( + x_mid, y_mid, z_mid, vx_mid, vy_mid, vz_mid, m, g, b) + vx = vx + ax_mid * dt + vy = vy + ay_mid * dt + vz = vz + az_mid * dt + return x, y, z, vx, vy, vz + + +def Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b): + ax, ay, az = compute_acceleration(x, y, z, vx, vy, vz, m, g, b) + vx_half = vx + 0.5 * ax * dt + vy_half = vy + 0.5 * ay * dt + vz_half = vz + 0.5 * az * dt + + x = x + vx_half * dt + 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) + ax_next, ay_next, az_next = compute_acceleration( + x, y, z, vx_next, vy_next, vz_next, 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 + return x, y, z, vx, vy, vz + +def Limit_in_box(a, amin, amax, va): + """限制物体在边界内，发生碰撞时反弹。""" + over = a > amax + under = a < amin + a = np.where(over, amax, a) + a = np.where(under, amin, a) + va = np.where(over | under, -va, va) + return a, va + +def apply_motion_update(x, y, z, vx, vy, vz, dt, m, g, b): + """按配置选择的位置更新算法推进一步。""" + if METHOD == "explicit_euler": + x, y, z, vx, vy, vz = Explicit_Euler_Method(x, y, z, vx, vy, vz, dt, m, g, b) + elif METHOD == "implicit_euler": + x, y, z, vx, vy, vz = Implicit_Euler_Method(x, y, z, vx, vy, vz, dt, m, g, b) + elif METHOD == "midpoint": + x, y, z, vx, vy, vz = Midpoint_Method(x, y, z, vx, vy, vz, dt, m, g, b) + elif METHOD == "leapfrog": + x, y, z, vx, vy, vz = Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b) + else: + raise ValueError(f"未知算法: {METHOD}") + + x, vx = Limit_in_box(x, X_MIN, X_MAX, vx) + y, vy = Limit_in_box(y, Y_MIN, Y_MAX, vy) + z, vz = Limit_in_box(z, Z_MIN, Z_MAX, vz) + return x, y, z, vx, vy, vz + + +def apply_fixed_constraints(x, y, z, vx, vy, vz): + """Keep fixed degrees of freedom at their initial coordinate with zero speed.""" + fixed = ATOM_FIXED != 0 + positions = np.column_stack((x, y, z)) + velocities = np.column_stack((vx, vy, vz)) + positions = np.where(fixed, ATOM_POSITIONS, positions) + velocities = np.where(fixed, 0.0, velocities) + return ( + positions[:, 0], positions[:, 1], positions[:, 2], + velocities[:, 0], velocities[:, 1], velocities[:, 2], + ) + +def wrap_position(x, y, z): + """边界回绕（ dynamics 模式）。""" + x = np.where(x > X_MAX, X_MIN, x) + x = np.where(x < X_MIN, X_MAX, x) + y = np.where(y > Y_MAX, Y_MIN, y) + y = np.where(y < Y_MIN, Y_MAX, y) + z = np.where(z > Z_MAX, Z_MIN, z) + z = np.where(z < Z_MIN, Z_MAX, z) + return x, y, z + + +# =========================================================================== +# 主计算流程 +# =========================================================================== + +def run_simulation(): + """计算 NT 步轨迹，返回位置/速度数组。 + + 步骤控制： + - warmup_steps: 预热步数，跳过不记录（用于稳定初始状态） + - 实际记录步数 = NT - warmup_steps + """ + # 预热阶段 + x = ATOM_POSITIONS[:, 0].copy() + y = ATOM_POSITIONS[:, 1].copy() + z = ATOM_POSITIONS[:, 2].copy() + vx = ATOM_VELOCITIES[:, 0].copy() + vy = ATOM_VELOCITIES[:, 1].copy() + vz = ATOM_VELOCITIES[:, 2].copy() + x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz) + + if warmup_steps is not None and warmup_steps > 0: + print(f"[compute] 预热阶段: 前 {warmup_steps} 步不记录") + for step in trange(warmup_steps, desc="[compute] 预热"): + t = (step + 1) * DT + x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B) + x, y, z = wrap_position(x, y, z) + x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz) + print( + f"[compute] 预热完成，展示原子位置: " + f"({x[PLOT_ATOM_ROW]:.4f}, {y[PLOT_ATOM_ROW]:.4f}, {z[PLOT_ATOM_ROW]:.4f})" + ) + + # 记录阶段 + record_steps = NT - (warmup_steps or 0) + n_atoms = len(ATOM_IDS) + traj_x = np.zeros((record_steps, n_atoms), dtype=np.float64) + traj_y = np.zeros((record_steps, n_atoms), dtype=np.float64) + traj_z = np.zeros((record_steps, n_atoms), dtype=np.float64) + traj_vx = np.zeros((record_steps, n_atoms), dtype=np.float64) + traj_vy = np.zeros((record_steps, n_atoms), dtype=np.float64) + traj_vz = np.zeros((record_steps, n_atoms), dtype=np.float64) + + for step in trange(record_steps, desc="[compute] 计算中"): + traj_x[step] = x + traj_y[step] = y + traj_z[step] = z + traj_vx[step] = vx + traj_vy[step] = vy + traj_vz[step] = vz + + x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B) + x, y, z = wrap_position(x, y, z) + x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz) + + return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz + + +def save_trajectory_table_txt(txt_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT): + """将轨迹数据保存为逐行表格文本，便于教学查看。""" + if np.ndim(traj_x) > 1: + row = PLOT_ATOM_ROW if PLOT_ATOM_ROW is not None else 0 + traj_x = traj_x[:, row] + traj_y = traj_y[:, row] + traj_z = traj_z[:, row] + traj_vx = traj_vx[:, row] + traj_vy = traj_vy[:, row] + traj_vz = traj_vz[:, row] + + with open(txt_path, 'w', encoding='utf-8') as f: + # 写入表头 + f.write("# 步数, 时间(s), x, y, z, vx, vy, vz\n") + for step in range(NT): + t = step * DT + x = traj_x[step] + y = traj_y[step] + z = traj_z[step] + vx = traj_vx[step] + vy = traj_vy[step] + vz = traj_vz[step] + f.write(f"{step+1:6d}, {t:.6f}, {x:.6f}, {y:.6f}, {z:.6f}, {vx:.6f}, {vy:.6f}, {vz:.6f}\n") + print(f"[compute] 轨迹文本文件已保存至: {txt_path}") + + +def main(): + # 默认参数文件 + script_dir = os.path.dirname(os.path.abspath(__file__)) + param_file = os.path.join(get_input_dir(script_dir), "input.txt") + if len(sys.argv) > 1: + param_file = sys.argv[1] + + # 加载参数 + load_parameters(param_file) + + script_dir = os.path.dirname(os.path.abspath(__file__)) + output_dir = get_output_dir(script_dir) + + print(f"[compute] 开始计算 NT={NT} DT={DT} ") + traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = run_simulation() + print(f"[compute] 计算完成，共 {NT} 步") + + save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, script_dir) + + # 同时保存为逐行表格，便于直接查看 + txt_path = os.path.join(output_dir, "trajectory_table.txt") + save_trajectory_table_txt(txt_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT) + + +if __name__ == "__main__": + main() diff --git a/draw.py b/draw.py new file mode 100644 index 0000000..530d221 --- /dev/null +++ b/draw.py @@ -0,0 +1,492 @@ +"""VisPy 演示：加载预计算轨迹数据，驱动小球运动动画。 + +计算与显示完全分离： + 1. 先运行 compute.py → 生成 output/trajectory.txt（全量 NT 步轨迹） + 2. 再运行 sample.py → 从 output/trajectory.txt 抽帧生成 output/display.txt + 3. 本文件加载 output/display.txt，按帧播放动画 + +用法： + python draw.py # 使用 dynamics 根目录下的 output/ + python draw.py examples/case01/output # 指定案例输出目录 +""" + +import numpy as np +import os +import sys +from vispy import app, scene +from vispy.visuals.transforms import STTransform + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import compute + +# =========================================================================== +# 加载预计算轨迹 +# =========================================================================== +script_dir = os.path.dirname(os.path.abspath(__file__)) +if len(sys.argv) > 1: + # 用户指定了输出目录 + output_dir = os.path.abspath(sys.argv[1]) +else: + output_dir = compute.get_output_dir(script_dir) +os.environ["DYNAMICS_OUTPUT_DIR"] = output_dir +disp_path = os.path.join(output_dir, "display.txt") + +if not os.path.exists(disp_path): + raise FileNotFoundError( + f"找不到 display.txt！\n" + f"期望路径: {disp_path}\n" + f"请先运行 compute.py 计算轨迹，再运行 sample.py 生成显示数组。\n" + f"用法: python draw.py [案例输出目录]" + ) + +disp_data = compute.load_text_data(disp_path) + +# 单原子数据（plot_atom：用于信息显示） +DISP_X = disp_data["disp_x"] +DISP_Y = disp_data["disp_y"] +DISP_Z = disp_data["disp_z"] +DISP_VX = disp_data["disp_vx"] +DISP_VY = disp_data["disp_vy"] +DISP_VZ = disp_data["disp_vz"] + +# 全原子数据（用于多球绘制） +DISP_ALL_X = disp_data["disp_all_x"] # (n_frames, n_atoms) +DISP_ALL_Y = disp_data["disp_all_y"] +DISP_ALL_Z = disp_data["disp_all_z"] +DISP_ALL_VX = disp_data["disp_all_vx"] +DISP_ALL_VY = disp_data["disp_all_vy"] +DISP_ALL_VZ = disp_data["disp_all_vz"] + +DISP_T = disp_data["disp_t"] +DISP_STEP = disp_data["disp_step"] +N_FRAMES = int(disp_data["n_frames"]) +NT = int(disp_data["NT"]) +DT = float(disp_data["DT"]) +NSTEP = int(disp_data["NSTEP"]) + +# 原子信息 +ATOM_IDS = disp_data.get("atom_ids", np.array([1])) +ATOM_RADII = disp_data.get("atom_radii", np.array([float(disp_data["ball_radius"])])) +N_ATOMS = len(ATOM_IDS) +PLOT_ATOM_ROW = int(disp_data.get("plot_atom_row", 0)) +PLOT_ATOM_ID = int(disp_data.get("plot_atom_id", ATOM_IDS[0])) +BOND_PAIRS = disp_data.get("bond_pairs", []) + +if N_FRAMES <= 0: + raise ValueError( + "output/display.txt 中没有可播放的帧，请检查 sample_start/sample_end/NSTEP 配置。") + +# 保留模拟边界常量（用于场景缩放、相机等），从 output/display.txt 中读取 +X_MIN = float(disp_data["X_MIN"]); X_MAX = float(disp_data["X_MAX"]) +Y_MIN = float(disp_data["Y_MIN"]); Y_MAX = float(disp_data["Y_MAX"]) +Z_MIN = float(disp_data["Z_MIN"]); Z_MAX = float(disp_data["Z_MAX"]) +X0 = float(disp_data["X0"]); Y0 = float(disp_data["Y0"]); Z0 = float(disp_data["Z0"]) +raw_alpha = disp_data["alpha"] +if isinstance(raw_alpha, (list, tuple, np.ndarray)): + alpha_list = [float(a) for a in raw_alpha] + if len(alpha_list) != 6: + raise ValueError(f"alpha 数组长度须为 6，实际为 {len(alpha_list)}") +else: + alpha_list = [float(raw_alpha)] * 6 +# 绘图参数 +ball_radius = float(disp_data["ball_radius"]) +ball_color_r = float(disp_data["ball_color_r"]) +ball_color_g = float(disp_data["ball_color_g"]) +ball_color_b = float(disp_data["ball_color_b"]) +box_color_r = float(disp_data["box_color_r"]) +box_color_g = float(disp_data["box_color_g"]) +box_color_b = float(disp_data["box_color_b"]) + + +# =========================================================================== +# 图形界面无关的几何参数（不参与物理计算，仅用于场景外观） +# =========================================================================== + +info_margin = 36 +axis_length = 10.0 + +initial_camera = { + "distance": 40.0, + "elevation": 0, + "azimuth": 0, + "center": (0, 0, 0), +} + + +# =========================================================================== +# 创建画布与相机 +# =========================================================================== +canvas = scene.SceneCanvas( + keys="interactive", + size=(1000, 700), + bgcolor=(0.08, 0.08, 0.10, 1.0), + show=True, +) + +view = canvas.central_widget.add_view() +view.camera = "turntable" +view.camera.distance = initial_camera["distance"] +view.camera.elevation = initial_camera["elevation"] +view.camera.azimuth = initial_camera["azimuth"] +view.camera.center = initial_camera["center"] + + +# =========================================================================== +# 场景对象 +# =========================================================================== +axis_width = 3 +arrow_size = 14 +axes_visible = True +axes_group = [] + +axes_group.append(scene.visuals.Arrow( + pos=np.array([[0, 0, 0], [axis_length, 0, 0]], dtype=np.float32), + color=(1.0, 0.2, 0.2, 1.0), + width=axis_width, + arrows=np.array([[0, 0, 0, axis_length, 0, 0]], dtype=np.float32), + arrow_size=arrow_size, + parent=view.scene, +)) +axes_group.append(scene.visuals.Arrow( + pos=np.array([[0, 0, 0], [0, axis_length, 0]], dtype=np.float32), + color=(0.2, 1.0, 0.2, 1.0), + width=axis_width, + arrows=np.array([[0, 0, 0, 0, axis_length, 0]], dtype=np.float32), + arrow_size=arrow_size, + parent=view.scene, +)) +axes_group.append(scene.visuals.Arrow( + pos=np.array([[0, 0, 0], [0, 0, axis_length]], dtype=np.float32), + color=(0.3, 0.6, 1.0, 1.0), + width=axis_width, + arrows=np.array([[0, 0, 0, 0, 0, axis_length]], dtype=np.float32), + arrow_size=arrow_size, + parent=view.scene, +)) + +axes_group.append(scene.visuals.Text(text="x", color=(1.0, 0.2, 0.2, 1.0), font_size=18, + pos=(axis_length + 0.2, 0, 0), anchor_x="left", anchor_y="center", parent=view.scene)) +axes_group.append(scene.visuals.Text(text="y", color=(0.2, 1.0, 0.2, 1.0), font_size=18, + pos=(0, axis_length + 0.2, 0), anchor_x="left", anchor_y="bottom", parent=view.scene)) +axes_group.append(scene.visuals.Text(text="z", color=(0.3, 0.6, 1.0, 1.0), font_size=18, + pos=(0, 0, axis_length + 0.2), anchor_x="left", anchor_y="bottom", parent=view.scene)) + +# 所有小球（每个原子一个球，不同颜色） +TAB10_RGB = np.array([ + [0.1216, 0.4667, 0.7059], # 蓝 + [0.8902, 0.4667, 0.1137], # 橙 + [0.1725, 0.6275, 0.1725], # 绿 + [0.8392, 0.1529, 0.1569], # 红 + [0.5804, 0.4039, 0.7412], # 紫 + [0.5490, 0.3373, 0.2941], # 棕 + [0.8902, 0.4667, 0.7608], # 粉 + [0.4980, 0.4980, 0.4980], # 灰 + [0.7373, 0.7412, 0.1333], # 黄绿 + [0.0902, 0.7451, 0.8118], # 青 +]) +balls = [] +for i in range(N_ATOMS): + r, g, b = TAB10_RGB[i % len(TAB10_RGB)] + s = scene.visuals.Sphere( + radius=float(ATOM_RADII[i]), method="latitude", + color=(r, g, b, 1.0), edge_color=None, parent=view.scene) + s.mesh.shading = "smooth" + balls.append(s) + +# 成键线（原子之间的连接） +if len(BOND_PAIRS) > 0: + n_bonds = len(BOND_PAIRS) + bond_pos = np.zeros((n_bonds * 2, 3), dtype=np.float32) + bond_lines = scene.visuals.Line( + pos=bond_pos, color=(0.7, 0.7, 0.7, 0.8), width=3, + connect="segments", parent=view.scene) +else: + bond_lines = None + +# 六个面形成立方体边界（每个面独立透明度，alpha<=0 时隐藏该面） +box_size = X_MAX - X_MIN +faces = [ + ((X_MAX, 0, 0), "-x"), + ((X_MIN, 0, 0), "+x"), + ((0, Y_MAX, 0), "-y"), + ((0, Y_MIN, 0), "+y"), + ((0, 0, Z_MAX), "-z"), + ((0, 0, Z_MIN), "+z"), +] +for f_idx, (pos, direction) in enumerate(faces): + a = alpha_list[f_idx] + if a <= 0: + continue + face_color = (box_color_r, box_color_g, box_color_b, a) + plane = scene.visuals.Plane( + width=box_size, height=box_size, width_segments=1, height_segments=1, + direction=direction, color=face_color, parent=view.scene) + plane.set_gl_state(depth_test=False, blend=True) + plane.transform = STTransform(translate=pos) + +# 右上角：相机信息 +camera_info = scene.visuals.Text( + text="", color="white", font_size=14, + pos=(0, 0), anchor_x="right", anchor_y="top", parent=canvas.scene) + +# 左上角：小球信息 +ball_info = scene.visuals.Text( + text="", color=(0.2, 1.0, 0.2, 1.0), font_size=28, + pos=(0, 0), anchor_x="left", anchor_y="top", + face="黑体", bold=True, parent=canvas.scene) + +# 左上角：reset 按钮（在 info 上方） +reset_btn_size = (60, 30) +reset_button = scene.visuals.Rectangle( + center=(reset_btn_size[0] / 2 + 8, reset_btn_size[1] / 2 + 8), + width=reset_btn_size[0], height=reset_btn_size[1], + radius=6, color=(0.18, 0.35, 0.65, 0.85), + border_color="white", parent=canvas.scene) +reset_button_label = scene.visuals.Text( + text="reset", color="white", font_size=16, + pos=(reset_btn_size[0] / 2 + 8, reset_btn_size[1] / 2 + 8), + anchor_x="center", anchor_y="center", + bold=True, parent=canvas.scene) + +# 信息显示/隐藏 切换按钮（左上角小方块，在 reset 下方） +info_btn_size = (60, 30) +info_toggle_visible = True +info_button = scene.visuals.Rectangle( + center=(info_btn_size[0] / 2 + 8, info_btn_size[1] / 2 + 8), + width=info_btn_size[0], height=info_btn_size[1], + radius=6, color=(0.9, 0.3, 0.3, 0.9), + border_color="white", parent=canvas.scene) +info_button_label = scene.visuals.Text( + text="info", color="white", font_size=16, + pos=(info_btn_size[0] / 2 + 8, info_btn_size[1] / 2 + 8), + anchor_x="center", anchor_y="center", + bold=True, parent=canvas.scene) + +# 强制在初始化时定位到左上角（canvas.scene坐标系：左下角为原点） +_cw, _ch = canvas.size +reset_button.center = (reset_btn_size[0] / 2 + 8, _ch - reset_btn_size[1] / 2 - 8 - info_btn_size[1] - 4) +reset_button_label.pos = reset_button.center +info_button.center = (info_btn_size[0] / 2 + 8, _ch - info_btn_size[1] / 2 - 8) +info_button_label.pos = info_button.center + +# axes 显示/隐藏 按钮（在 info 下方） +axes_btn_size = (60, 30) +axes_button = scene.visuals.Rectangle( + center=(axes_btn_size[0] / 2 + 8, axes_btn_size[1] / 2 + 8), + width=axes_btn_size[0], height=axes_btn_size[1], + radius=6, color=(0.3, 0.7, 0.3, 0.9), + border_color="white", parent=canvas.scene) +axes_button_label = scene.visuals.Text( + text="axes", color="white", font_size=16, + pos=(axes_btn_size[0] / 2 + 8, axes_btn_size[1] / 2 + 8), + anchor_x="center", anchor_y="center", + bold=True, parent=canvas.scene) +axes_button.center = (axes_btn_size[0] / 2 + 8, _ch - axes_btn_size[1] / 2 - 8 - info_btn_size[1] - 4 - axes_btn_size[1] - 4) +axes_button_label.pos = axes_button.center + + +# =========================================================================== +# 回调函数 +# =========================================================================== + +def update_camera_info(event=None): + c = view.camera + camera_info.text = ( + "Camera\n" + f"center = ({c.center[0]:.2f}, {c.center[1]:.2f}, {c.center[2]:.2f})\n" + f"distance = {c.distance:.2f}\n" + f"elevation = {c.elevation:.2f}\n" + f"azimuth = {c.azimuth:.2f}" + ) + + +def update_ball_info(frame_idx, x, y, z, vx, vy, vz): + step = int(DISP_STEP[frame_idx]) + t = float(DISP_T[frame_idx]) + ball_info.text = ( + f"原子 {PLOT_ATOM_ID} (共 {N_ATOMS} 个原子)\n" + f"frame {frame_idx+1}/{N_FRAMES} | saved step {step}/{NT-1}\n" + f"t = {t:.2f} s | dt = {DT:.3f} s | nstep = {NSTEP}\n" + f"Position: ({x:.2f}, {y:.2f}, {z:.2f})\n" + f"Velocity: ({vx:.2f}, {vy:.2f}, {vz:.2f})\n" + ) + + +def reposition_camera_info(event=None): + width, height = canvas.size + camera_info.pos = (width - 20, height - 20) + # reset → info → axes 三按钮纵向排列 + gap = 4 + info_y = info_btn_size[1] / 2 + 8 + reset_y = info_y + info_btn_size[1] + gap + axes_y = info_y + info_btn_size[1] + gap + axes_btn_size[1] + gap + reset_button.center = (reset_btn_size[0] / 2 + 8, height - reset_y) + reset_button_label.pos = reset_button.center + info_button.center = (info_btn_size[0] / 2 + 8, height - info_y) + info_button_label.pos = info_button.center + axes_button.center = (axes_btn_size[0] / 2 + 8, height - axes_y) + axes_button_label.pos = axes_button.center + # info 文字放在所有按钮下方 + buttons_bottom = height - (axes_y + axes_btn_size[1] / 2) + ball_info.pos = (info_margin, buttons_bottom - 10) + update_ball_info(frame_idx, DISP_X[frame_idx], DISP_Y[frame_idx], DISP_Z[frame_idx], + DISP_VX[frame_idx], DISP_VY[frame_idx], DISP_VZ[frame_idx]) + update_camera_info() + + +def handle_view_interaction(event): + update_camera_info() + + +def rotate_about_screen_normal(angle): + if hasattr(view.camera, "roll"): + view.camera.roll = (view.camera.roll + angle) % 360 + else: + view.camera.azimuth = (view.camera.azimuth + angle) % 360 + update_camera_info() + + +def handle_key_press(event): + key_name = "" + if getattr(event, "text", None): + key_name = event.text.lower() + elif getattr(event, "key", None) is not None: + key_name = str(event.key).lower() + if key_name == "q": + rotate_about_screen_normal(-90) + elif key_name == "e": + rotate_about_screen_normal(90) + + +def reset_camera_view(): + global frame_idx + frame_idx = 0 + # 立即复位所有小球到第 0 帧 + for i in range(N_ATOMS): + balls[i].transform = STTransform(translate=( + float(DISP_ALL_X[frame_idx, i]), + float(DISP_ALL_Y[frame_idx, i]), + float(DISP_ALL_Z[frame_idx, i]), + )) + if bond_lines is not None and len(BOND_PAIRS) > 0: + pos = np.zeros((len(BOND_PAIRS) * 2, 3), dtype=np.float32) + for b_idx, (i, j) in enumerate(BOND_PAIRS): + pos[b_idx * 2] = [DISP_ALL_X[frame_idx, i], DISP_ALL_Y[frame_idx, i], DISP_ALL_Z[frame_idx, i]] + pos[b_idx * 2 + 1] = [DISP_ALL_X[frame_idx, j], DISP_ALL_Y[frame_idx, j], DISP_ALL_Z[frame_idx, j]] + bond_lines.set_data(pos=pos) + # 复位相机 + view.camera.distance = initial_camera["distance"] + view.camera.elevation = initial_camera["elevation"] + view.camera.azimuth = initial_camera["azimuth"] + view.camera.center = initial_camera["center"] + if hasattr(view.camera, "roll"): + view.camera.roll = 0 + update_camera_info() + + +def handle_mouse_press(event): + global info_toggle_visible, axes_visible + if event.pos is None: + return + ex, ey = event.pos + + # reset 按钮 + bx, by = reset_button.center + lw = reset_btn_size[0] / 2; lh = reset_btn_size[1] / 2 + if (bx - lw) <= ex <= (bx + lw) and (by - lh) <= ey <= (by + lh): + reset_camera_view() + return + + # info 按钮 + bx, by = info_button.center + if (bx - info_btn_size[0]/2) <= ex <= (bx + info_btn_size[0]/2) and (by - info_btn_size[1]/2) <= ey <= (by + info_btn_size[1]/2): + info_toggle_visible = not info_toggle_visible + ball_info.visible = info_toggle_visible + camera_info.visible = info_toggle_visible + return + + # axes 按钮 + bx, by = axes_button.center + if (bx - axes_btn_size[0]/2) <= ex <= (bx + axes_btn_size[0]/2) and (by - axes_btn_size[1]/2) <= ey <= (by + axes_btn_size[1]/2): + axes_visible = not axes_visible + for axis in axes_group: + axis.visible = axes_visible + + +# =========================================================================== +# 动画初始化 +# =========================================================================== +frame_idx = 0 + +# 初始帧：摆放所有小球并刷新 UI +for i in range(N_ATOMS): + balls[i].transform = STTransform(translate=( + float(DISP_ALL_X[frame_idx, i]), + float(DISP_ALL_Y[frame_idx, i]), + float(DISP_ALL_Z[frame_idx, i]), + )) +# 初始帧：更新成键线 +if bond_lines is not None and len(BOND_PAIRS) > 0: + pos = np.zeros((len(BOND_PAIRS) * 2, 3), dtype=np.float32) + for b_idx, (i, j) in enumerate(BOND_PAIRS): + pos[b_idx * 2] = [DISP_ALL_X[frame_idx, i], DISP_ALL_Y[frame_idx, i], DISP_ALL_Z[frame_idx, i]] + pos[b_idx * 2 + 1] = [DISP_ALL_X[frame_idx, j], DISP_ALL_Y[frame_idx, j], DISP_ALL_Z[frame_idx, j]] + bond_lines.set_data(pos=pos) +reposition_camera_info() +update_ball_info(frame_idx, + float(DISP_X[frame_idx]), float(DISP_Y[frame_idx]), float(DISP_Z[frame_idx]), + float(DISP_VX[frame_idx]), float(DISP_VY[frame_idx]), float(DISP_VZ[frame_idx])) + +print(f"[draw] 加载 output/display.txt: {N_FRAMES} 帧, {N_ATOMS} 个原子, NT={NT}, DT={DT}, NSTEP={NSTEP}") +print(f"[draw] 绘图参数: ball_radius={ball_radius}, box_color=({box_color_r:.2f},{box_color_g:.2f},{box_color_b:.2f}), alpha={alpha_list}") + + +# =========================================================================== +# 每帧回调：仅推进帧索引，从预存数组读取位置，零物理计算 +# =========================================================================== +def update(event): + global frame_idx + frame_idx = (frame_idx + 1) % N_FRAMES # 循环播放 + + # 更新所有小球位置 + for i in range(N_ATOMS): + x = float(DISP_ALL_X[frame_idx, i]) + y = float(DISP_ALL_Y[frame_idx, i]) + z = float(DISP_ALL_Z[frame_idx, i]) + balls[i].transform = STTransform(translate=(x, y, z)) + + # 更新成键线 + if bond_lines is not None and len(BOND_PAIRS) > 0: + pos = np.zeros((len(BOND_PAIRS) * 2, 3), dtype=np.float32) + for b_idx, (i, j) in enumerate(BOND_PAIRS): + pos[b_idx * 2] = [DISP_ALL_X[frame_idx, i], DISP_ALL_Y[frame_idx, i], DISP_ALL_Z[frame_idx, i]] + pos[b_idx * 2 + 1] = [DISP_ALL_X[frame_idx, j], DISP_ALL_Y[frame_idx, j], DISP_ALL_Z[frame_idx, j]] + bond_lines.set_data(pos=pos) + + # 信息面板显示 plot_atom 的数据 + x = float(DISP_X[frame_idx]) + y = float(DISP_Y[frame_idx]) + z = float(DISP_Z[frame_idx]) + vx = float(DISP_VX[frame_idx]) + vy = float(DISP_VY[frame_idx]) + vz = float(DISP_VZ[frame_idx]) + update_ball_info(frame_idx, x, y, z, vx, vy, vz) + update_camera_info() + + +# =========================================================================== +# 事件绑定与启动 +# =========================================================================== +timer = app.Timer(interval=0.02, connect=update, start=True) +canvas.events.mouse_move.connect(handle_view_interaction) +canvas.events.mouse_press.connect(handle_mouse_press) +canvas.events.mouse_wheel.connect(handle_view_interaction) +canvas.events.resize.connect(reposition_camera_info) +canvas.events.key_press.connect(handle_key_press) + +if hasattr(canvas, "native") and hasattr(canvas.native, "setFocus"): + canvas.native.setFocus() + + +if __name__ == "__main__": + app.run() diff --git a/dynamics.py b/dynamics.py new file mode 100644 index 0000000..4c6a83e --- /dev/null +++ b/dynamics.py @@ -0,0 +1,454 @@ +""" +dynamics.py +------ +统一入口：读取 YAML 配置文件 → 运行模拟 → 抽帧 → 绘图（可选） + +用法： + python dynamics.py # 使用 input/parameters.yaml + python dynamics.py input/parameters.yaml # 指定配置文件 + python dynamics.py --config input/parameters.yaml --no-plot +""" + +import os +import sys +import subprocess +import argparse +from contextlib import contextmanager +from pathlib import Path + +import yaml +import numpy as np + +# 导入同目录下的模块 +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import compute + + +def read_optional_index(data, key, default_value): + """Read an optional integer index from structured txt metadata.""" + if key not in data: + return default_value + value = data[key] + if value is None or int(value) < 0: + return default_value + return int(value) + + +def load_yaml_config(config_path): + """从 YAML 文件加载配置字典。""" + with open(config_path, 'r', encoding='utf-8') as f: + config = yaml.safe_load(f) + return config + + +def resolve_path(base_dir, path_value): + """Resolve a path relative to the runtime base directory.""" + path = Path(path_value) + if path.is_absolute(): + return path + return (Path(base_dir) / path).resolve() + + +@contextmanager +def runtime_dir_overrides(input_dir=None, output_dir=None): + """Temporarily override runtime input/output directories for compute helpers.""" + old_input = os.environ.get("DYNAMICS_INPUT_DIR") + old_output = os.environ.get("DYNAMICS_OUTPUT_DIR") + try: + if input_dir is not None: + os.environ["DYNAMICS_INPUT_DIR"] = str(Path(input_dir).resolve()) + if output_dir is not None: + os.environ["DYNAMICS_OUTPUT_DIR"] = str(Path(output_dir).resolve()) + yield + finally: + if old_input is None: + os.environ.pop("DYNAMICS_INPUT_DIR", None) + else: + os.environ["DYNAMICS_INPUT_DIR"] = old_input + if old_output is None: + os.environ.pop("DYNAMICS_OUTPUT_DIR", None) + else: + os.environ["DYNAMICS_OUTPUT_DIR"] = old_output + + +def build_sample_indices(total_steps, sample_step, sample_start, sample_end): + """Validate sampling settings and build frame indices.""" + if sample_step <= 0: + raise ValueError(f"NSTEP 必须为正整数，实际为 {sample_step}") + if sample_start < 0: + raise ValueError(f"sample_start 不能小于 0，实际为 {sample_start}") + if sample_end > total_steps: + raise ValueError( + f"sample_end 不能大于记录步数 {total_steps}，实际为 {sample_end}") + if sample_start >= sample_end: + raise ValueError( + f"sample_start 必须小于 sample_end，实际为 [{sample_start}, {sample_end})") + + n_frames = (sample_end - sample_start) // sample_step + if n_frames <= 0: + raise ValueError( + f"抽帧范围 [{sample_start}, {sample_end}) 过短，按 NSTEP={sample_step} 无法抽出任何帧") + + indices = np.arange(n_frames, dtype=np.int64) * sample_step + sample_start + return indices + + +def save_display_txt(data, out_dir=None): + """将抽帧数据保存到 output/display.txt（含所有参数元数据）。""" + if out_dir is None: + out_dir = os.path.dirname(os.path.abspath(__file__)) + disp_path = os.path.join(compute.get_output_dir(out_dir), "display.txt") + compute.save_text_data(disp_path, data) + print(f"[sample] 显示数组已保存至: {disp_path}") + return disp_path + + +def run_case(config_path, runtime_base, input_dir="input", output_dir="output", no_plot=False): + """Run one case with explicit program path, input path, and output path.""" + runtime_base = Path(runtime_base).resolve() + input_dir_path = resolve_path(runtime_base, input_dir) + output_dir_path = resolve_path(runtime_base, output_dir) + config_path = resolve_path(runtime_base, config_path) + + with runtime_dir_overrides(input_dir_path, output_dir_path): + # 1. 加载 YAML 配置 + config = load_yaml_config(config_path) + print(f"[run] 已加载配置: {config_path}") + + # 显示步骤控制信息 + 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()) + print(f"[run] 步骤开关: {step_flags}") + + warmup = config.get("warmup_steps", 0) + ss = config.get("sample_start", "从头") + se = config.get("sample_end", "到尾") + method = config.get("method", "explicit_euler") + coord_file = config.get("coord_file", os.path.join("input", "coord.txt")) + plot_atom = config.get("plot_atom", "第一个原子") + print(f"[run] 算法: {method}") + print(f"[run] 坐标文件: {coord_file}") + print(f"[run] 绘图/动画原子序号: {plot_atom}") + print(f"[run] 步骤控制: 预热={warmup}步, 抽帧范围=[{ss}, {se})") + + output_dir_abs = compute.get_output_dir(str(runtime_base)) + traj_path = os.path.join(output_dir_abs, "trajectory.txt") + 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] 已有输出，步骤已被跳过") + else: + # 目录存在但文件不全 → 强制重新计算 + if os.path.isdir(output_dir_abs): + print(f"[run] output/ 目录存在但文件不完整，将重新计算") + else: + print(f"[run] output/ 目录不存在，将执行完整流程") + + # 2. 运行物理模拟 → output/trajectory.txt + if config.get("step_simulate", 1): + engine = config.get("engine", "python") + total_steps = config["NT"] + record_steps = total_steps - (config.get("warmup_steps") or 0) + print(f"[run] 开始计算 总步数={total_steps} 记录步数={record_steps} DT={config['DT']}") + 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 + config["_skip_run"] = True + compute.run_from_config(config, str(runtime_base)) + config.pop("_skip_run", None) + input_dir_abs = str(input_dir_path.resolve()) + output_dir_abs = str(output_dir_path.resolve()) + 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)) + else: + print("[run] 步骤 [模拟] 已跳过，直接加载已有轨迹") + if not os.path.exists(traj_path): + print(f"[run] 错误: trajectory.txt 不存在，无法跳过模拟") + sys.exit(1) + + # 3. 抽帧 → output/display.txt + traj_path = os.path.join(output_dir_abs, "trajectory.txt") + data = compute.load_text_data(traj_path) + + NT = int(data["NT"]); DT = float(data["DT"]); NSTEP = int(data["NSTEP"]) + warmup_steps = int(data.get("warmup_steps", 0)) + plot_atom_row = int(data["plot_atom_row"]) if "plot_atom_row" in data else 0 + plot_atom_id = int(data["plot_atom_id"]) if "plot_atom_id" in data else int(data["atom_ids"][plot_atom_row]) + + # 抽帧范围控制 + sample_start = read_optional_index(data, "sample_start", 0) + sample_end = read_optional_index(data, "sample_end", NT) + + indices = build_sample_indices(NT, NSTEP, sample_start, sample_end) + n_frames = len(indices) + + print(f"[run] 抽帧范围: [{sample_start}, {sample_end}), 共 {n_frames} 帧") + + traj_x = data["traj_x"] + traj_y = data["traj_y"] + traj_z = data["traj_z"] + traj_vx = data["traj_vx"] + traj_vy = data["traj_vy"] + traj_vz = data["traj_vz"] + + if traj_x.ndim == 1: + selected_x = traj_x + selected_y = traj_y + selected_z = traj_z + selected_vx = traj_vx + selected_vy = traj_vy + selected_vz = traj_vz + all_x = traj_x[:, None] + all_y = traj_y[:, None] + all_z = traj_z[:, None] + all_vx = traj_vx[:, None] + all_vy = traj_vy[:, None] + all_vz = traj_vz[:, None] + else: + selected_x = traj_x[:, plot_atom_row] + selected_y = traj_y[:, plot_atom_row] + selected_z = traj_z[:, plot_atom_row] + selected_vx = traj_vx[:, plot_atom_row] + selected_vy = traj_vy[:, plot_atom_row] + selected_vz = traj_vz[:, plot_atom_row] + all_x = traj_x + all_y = traj_y + all_z = traj_z + all_vx = traj_vx + all_vy = traj_vy + all_vz = traj_vz + + if config.get("step_sample", 1): + disp_data = { + "disp_x": selected_x[indices], + "disp_y": selected_y[indices], + "disp_z": selected_z[indices], + "disp_vx": selected_vx[indices], + "disp_vy": selected_vy[indices], + "disp_vz": selected_vz[indices], + "disp_all_x": all_x[indices], + "disp_all_y": all_y[indices], + "disp_all_z": all_z[indices], + "disp_all_vx": all_vx[indices], + "disp_all_vy": all_vy[indices], + "disp_all_vz": all_vz[indices], + "disp_t": indices * DT, + "disp_step": indices, + "n_frames": n_frames, + "NT": NT, "DT": DT, "NSTEP": NSTEP, + "plot_atom_id": plot_atom_id, + "plot_atom_row": plot_atom_row, + "method": str(data["method"]) if "method" in data else "explicit_euler", + "coord_file": str(data["coord_file"]) if "coord_file" in data else os.path.join("input", "coord.txt"), + "atom_ids": data["atom_ids"] if "atom_ids" in data else np.array([1]), + "atom_masses": data["atom_masses"] if "atom_masses" in data else np.array([float(data["M"])]), + "atom_radii": data["atom_radii"] if "atom_radii" in data else np.array([float(data["ball_radius"])]), + "atom_positions": data["atom_positions"] if "atom_positions" in data else np.array([[float(data["X0"]), float(data["Y0"]), float(data["Z0"])]]), + "atom_velocities": data["atom_velocities"] if "atom_velocities" in data else np.array([[float(data["VX0"]), float(data["VY0"]), float(data["VZ0"])]]), + "atom_fixed": data["atom_fixed"] if "atom_fixed" in data else np.array([[0, 0, 0]]), + "bond_pairs": data.get("bond_pairs", np.zeros((0, 2), dtype=np.int64)).tolist(), + "warmup_steps": warmup_steps, + "sample_start": sample_start, + "sample_end": sample_end, + "X_MIN": float(data["X_MIN"]), "X_MAX": float(data["X_MAX"]), + "Y_MIN": float(data["Y_MIN"]), "Y_MAX": float(data["Y_MAX"]), + "Z_MIN": float(data["Z_MIN"]), "Z_MAX": float(data["Z_MAX"]), + "X0": float(data["X0"]), "Y0": float(data["Y0"]), "Z0": float(data["Z0"]), + "VX0": float(data["VX0"]), "VY0": float(data["VY0"]), "VZ0": float(data["VZ0"]), + "M": float(data["M"]) if "M" in data else 1.0, + "alpha": data["alpha"], + "ball_radius": float(data["ball_radius"]), + "ball_color_r": float(data["ball_color_r"]), + "ball_color_g": float(data["ball_color_g"]), + "ball_color_b": float(data["ball_color_b"]), + "box_color_r": float(data["box_color_r"]), + "box_color_g": float(data["box_color_g"]), + "box_color_b": float(data["box_color_b"]), + } + save_display_txt(disp_data, str(runtime_base)) + print(f"[run] 抽帧完成: {sample_end - sample_start} 步 -> {n_frames} 帧") + else: + print("[run] 步骤 [抽帧] 已跳过") + + # 4. 绘图（可选） + if not no_plot and config.get("step_plot", 1): + try: + import matplotlib.pyplot as plt + + # 配置中文字体支持 + plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'WenQuanYi Micro Hei', 'DejaVu Sans'] + plt.rcParams['axes.unicode_minus'] = False + + time = np.arange(NT) * DT + n_atoms = all_x.shape[1] + atom_ids_list = data.get("atom_ids", np.arange(n_atoms) + 1) + + fig, axes = plt.subplots(3, 3, figsize=(15, 13)) + fig.suptitle("轨迹与能量分析", fontsize=16) + + # ── 位置 / 速度 6 子图（前 2 行，每行 3 列） ── + plot_configs = [ + (axes[0, 0], all_x, "x - 时间"), + (axes[0, 1], all_y, "y - 时间"), + (axes[0, 2], all_z, "z - 时间"), + (axes[1, 0], all_vx, "vx - 时间"), + (axes[1, 1], all_vy, "vy - 时间"), + (axes[1, 2], all_vz, "vz - 时间"), + ] + + colors = plt.cm.tab10(np.linspace(0, 1, n_atoms)) + + for ax, data_arr, title in plot_configs: + for i in range(n_atoms): + atom_id = int(atom_ids_list[i]) + ax.plot(time, data_arr[:, i], color=colors[i], linewidth=1.5, label=f"原子 {atom_id}") + ax.set_title(title) + ax.set_xlabel("时间 (s)") + ax.grid(True, alpha=0.3) + ax.legend() + + # ── 能量计算 ───────────────────────────────────── + masses = np.array(data["atom_masses"]) # (n_atoms,) + G_vec = np.array(data.get("G", [0.0, 0.0, -9.8])) # [gx, gy, gz] + + # 动能 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 + ) + + # 弹性势能 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: + for b_idx in range(len(bond_pairs)): + i, j = bond_pairs[b_idx] + 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) + 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 + + # 各原子总能量 + e_total = ek + ug + us # (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 + + # ── 第 3 行左：各原子总能量 ── + ax_e = axes[2, 0] + for i in range(n_atoms): + aid = int(atom_ids_list[i]) + ax_e.plot(time, e_total[:, i], color=colors[i], linewidth=1.5, label=f"原子 {aid}") + ax_e.set_title("各原子总能量") + ax_e.set_xlabel("时间 (s)") + ax_e.set_ylabel("能量") + ax_e.grid(True, alpha=0.3) + ax_e.legend(loc="upper right") + + # ── 第 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, us_sys, color='orange', linewidth=1.5, label="系统弹性势能") + ax_sys.plot(time, e_sys, 'r--', linewidth=1.5, label="系统总能量") + ax_sys.set_title("系统总能量") + ax_sys.set_xlabel("时间 (s)") + ax_sys.set_ylabel("能量") + ax_sys.grid(True, alpha=0.3) + ax_sys.legend(loc="upper right") + + # 隐藏第 3 行第 3 列空子图 + axes[2, 2].set_visible(False) + + plt.tight_layout(rect=[0, 0.03, 1, 0.95]) + plot_path = os.path.join(output_dir_abs, "trajectory_plots.png") + plt.savefig(plot_path, dpi=300, bbox_inches="tight") + print(f"[run] 轨迹图表已保存至: {plot_path}") + plt.show() + except ImportError: + print("[run] 警告: 未安装 matplotlib，跳过绘图") + + print(f"[run] 完成！输出目录: {output_dir_abs}") + + # 5. 自动播放动画（可选） + if config.get("step_animation", 0): + draw_script = os.path.join(os.path.dirname(os.path.abspath(__file__)), "draw.py") + if os.path.exists(draw_script): + try: + print("[run] 正在启动 VisPy 3D 动画窗口…") + subprocess.Popen( + [sys.executable, draw_script], + cwd=runtime_base, + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + ) + except Exception as e: + print(f"[run] 启动动画失败: {e}") + else: + print(f"[run] 未找到动画脚本: {draw_script}") + else: + print("[run] 运行 python draw.py 查看动画。") + + +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", dest="config_override", + help="YAML 配置文件路径（可选，优先于位置参数）") + parser.add_argument("--input-dir", default="input", + help="输入目录路径（默认: input）") + parser.add_argument("--output-dir", default="output", + help="输出目录路径（默认: output）") + parser.add_argument("--runtime-base", default=".", + help="案例运行根目录（默认: 当前目录）") + parser.add_argument("--no-plot", action="store_true", + help="跳过 matplotlib 绘图") + args = parser.parse_args() + + config_path = args.config_override or args.config_file + runtime_base = resolve_path(os.getcwd(), args.runtime_base) + config_path_abs = resolve_path(runtime_base, config_path) + if not os.path.exists(config_path_abs): + print(f"错误: 配置文件不存在: {config_path_abs}") + sys.exit(1) + run_case( + config_path=config_path_abs, + runtime_base=runtime_base, + input_dir=args.input_dir, + output_dir=args.output_dir, + no_plot=args.no_plot, + ) + + +if __name__ == "__main__": + main() diff --git a/engines/c/Makefile b/engines/c/Makefile new file mode 100644 index 0000000..79cfbe4 --- /dev/null +++ b/engines/c/Makefile @@ -0,0 +1,63 @@ +# engines/c/Makefile +# 跨平台编译：make → 本地系统编译 +# make linux → Linux 交叉编译（需 x86_64-linux-gnu-gcc） +# make windows → Windows 交叉编译（需 x86_64-w64-mingw32-gcc） +# make macos → macOS 交叉编译（需 osxcross 工具链） + +CC = gcc +CFLAGS = -O3 -march=native -Wall -Wextra +LDFLAGS = -lm +SRCS = main.c + +# 自动检测系统 +UNAME_S := $(shell uname -s 2>/dev/null || echo Windows) + +# 目标文件名：统一使用 .exe 后缀（方便 Python 跨平台调用） +TARGET = build/dynamics_c.exe + +# ── 本地编译 ───────────────────────────────── +.PHONY: all clean linux windows macos + +all: $(TARGET) + +$(TARGET): $(SRCS) | build + $(CC) $(CFLAGS) -o $@ $(SRCS) $(LDFLAGS) + @echo " === C engine built: $@ ===" + +build: + mkdir -p build + +# ── 交叉编译 ───────────────────────────────── +# Linux → Linux (x86_64) +linux: CROSS_PREFIX = x86_64-linux-gnu- +linux: CC = $(CROSS_PREFIX)gcc +linux: CFLAGS = -O3 -march=x86-64 -Wall -Wextra +linux: $(SRCS) | build + $(CC) $(CFLAGS) -o build/dynamics_c_linux.exe $(SRCS) $(LDFLAGS) + @echo " === Linux binary: build/dynamics_c_linux.exe ===" + +# 任意平台 → Windows (x86_64) +# 需要安装 MinGW 交叉编译器： +# apt install mingw-w64 (Debian/Ubuntu) +# brew install mingw-w64 (macOS) +windows: CROSS_PREFIX = x86_64-w64-mingw32- +windows: CC = $(CROSS_PREFIX)gcc +windows: CFLAGS = -O3 -march=x86-64 -Wall -Wextra +windows: $(SRCS) | build + $(CC) $(CFLAGS) -o build/dynamics_c_win.exe $(SRCS) $(LDFLAGS) + @echo " === Windows binary: build/dynamics_c_win.exe ===" + +# 任意平台 → macOS (x86_64) +# 需要安装 osxcross 工具链 +macos: CROSS_PREFIX = x86_64-apple-darwin- +macos: CC = $(CROSS_PREFIX)gcc +macos: CFLAGS = -O3 -march=x86-64 -Wall -Wextra +macos: $(SRCS) | build + $(CC) $(CFLAGS) -o build/dynamics_c_mac.exe $(SRCS) $(LDFLAGS) + @echo " === macOS binary: build/dynamics_c_mac.exe ===" + +# ── 编译所有平台 ────────────────────────────── +all-platforms: linux windows macos + +clean: + rm -rf build *.o diff --git a/engines/c/main.c b/engines/c/main.c new file mode 100644 index 0000000..695927b --- /dev/null +++ b/engines/c/main.c @@ -0,0 +1,507 @@ +/** + * engines/c/main.c + * ----------------- + * C 语言动力学模拟引擎。 + * 输入: param.json 数值参数（Python 从 YAML 转换得来） + * /coord.txt + * /connection.txt + * /bond.txt + * 输出: /trajectory.txt (JSON 格式，与 Python 版兼容) + * + * 编译: make + * 用法: ./build/dynamics_c + */ + +#include +#include +#include +#include +#include + +/* ======================================================================== + * 配置参数（从 param.json 读取） + * ======================================================================== */ +typedef struct { + double box_a; /* 盒子半边长 */ + int NT; /* 总步数 */ + double DT; /* 时间步长 */ + int NSTEP; /* 抽帧间隔 */ + int warmup_steps; /* 预热步数 */ + double G[3]; /* 重力分量 */ + double B[3]; /* 阻尼分量 */ +} SimParams; + +/* ======================================================================== + * 原子数据 + * ======================================================================== */ +typedef struct { + int n_atoms; + int *atom_ids; + double *masses; + double *radii; + double *pos_0; /* 初始位置 (n_atoms*3) */ + double *vel_0; /* 初始速度 (n_atoms*3) */ + int *fixed; /* 固定约束 (n_atoms*3) */ +} AtomData; + +/* ======================================================================== + * 成键数据 + * ======================================================================== */ +typedef struct { + int n_bonds; + int *pairs; /* (n_bonds*2) */ + double *stiffness; + double *rest_lengths; +} BondData; + +/* ======================================================================== + * 轨迹缓冲区 + * ======================================================================== */ +typedef struct { + int n_steps; + int n_atoms; + double *x, *y, *z; + double *vx, *vy, *vz; +} Trajectory; + +/* ======================================================================== + * 辅助函数 + * ======================================================================== */ + +static void die(const char *msg) { + fprintf(stderr, "[C-engine] 错误: %s\n", msg); + exit(1); +} + +static void *xmalloc(size_t sz) { + void *p = malloc(sz); + if (!p) die("内存分配失败"); + return p; +} + +/* 从 JSON 中读取 double 值。简单实现：搜索 key 并解析 */ +static double json_read_double(const char *json, const char *key) { + char search[256]; + snprintf(search, sizeof(search), "\"%s\"", key); + const char *p = strstr(json, search); + if (!p) return 0.0; + p = strchr(p, ':'); + if (!p) return 0.0; + p++; + while (*p == ' ' || *p == '\t' || *p == '\n') p++; + return strtod(p, NULL); +} + +static int json_read_int(const char *json, const char *key) { + return (int)json_read_double(json, key); +} + +/* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */ +static void json_read_double3(const char *json, const char *key, double out[3]) { + char search[256]; + snprintf(search, sizeof(search), "\"%s\"", key); + const char *p = strstr(json, search); + if (!p) { out[0]=out[1]=out[2]=0; return; } + p = strchr(p, '['); + if (!p) { out[0]=out[1]=out[2]=0; return; } + p++; + for (int i = 0; i < 3; i++) { + while (*p == ' ' || *p == '\t' || *p == '\n' || *p == ',') p++; + out[i] = strtod(p, (char**)&p); + } +} + +/* 读取 param.json */ +static SimParams read_params(const char *path) { + FILE *f = fopen(path, "rb"); + if (!f) die("无法打开 param.json"); + fseek(f, 0, SEEK_END); + long sz = ftell(f); + fseek(f, 0, SEEK_SET); + char *buf = (char*)xmalloc(sz + 1); + fread(buf, 1, sz, f); + buf[sz] = '\0'; + fclose(f); + + SimParams p; + p.box_a = json_read_double(buf, "box_a"); + p.NT = json_read_int(buf, "NT"); + p.DT = json_read_double(buf, "DT"); + p.NSTEP = json_read_int(buf, "NSTEP"); + p.warmup_steps = json_read_int(buf, "warmup_steps"); + json_read_double3(buf, "G", p.G); + json_read_double3(buf, "B", p.B); + + free(buf); + return p; +} + +/* 读取 coord.txt */ +static AtomData read_coord(const char *input_dir) { + char path[512]; + snprintf(path, sizeof(path), "%s/coord.txt", input_dir); + FILE *f = fopen(path, "r"); + if (!f) die("无法打开 coord.txt"); + + /* 跳过第一行表头 */ + char line[1024]; + if (!fgets(line, sizeof(line), f)) die("coord.txt 为空"); + + /* 先读入所有行到内存 */ + int capacity = 16; + AtomData a; + a.n_atoms = 0; + a.atom_ids = (int*)xmalloc(capacity * sizeof(int)); + a.masses = (double*)xmalloc(capacity * sizeof(double)); + a.radii = (double*)xmalloc(capacity * sizeof(double)); + a.pos_0 = (double*)xmalloc(capacity * 3 * sizeof(double)); + a.vel_0 = (double*)xmalloc(capacity * 3 * sizeof(double)); + a.fixed = (int*)xmalloc(capacity * 3 * sizeof(int)); + + while (fgets(line, sizeof(line), f)) { + if (a.n_atoms >= capacity) { + capacity *= 2; + a.atom_ids = realloc(a.atom_ids, capacity * sizeof(int)); + a.masses = realloc(a.masses, capacity * sizeof(double)); + a.radii = realloc(a.radii, capacity * sizeof(double)); + a.pos_0 = realloc(a.pos_0, capacity * 3 * sizeof(double)); + a.vel_0 = realloc(a.vel_0, capacity * 3 * sizeof(double)); + a.fixed = realloc(a.fixed, capacity * 3 * sizeof(int)); + } + int id, fx, fy, fz; + double mass, rad, px, py, pz, vx, vy, vz; + int n_parsed = sscanf(line, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %d", + &id, &mass, &rad, &px, &py, &pz, &vx, &vy, &vz, &fx, &fy, &fz); + if (n_parsed == 9) { + /* 无固定约束列 */ + fx = fy = fz = 0; + } else if (n_parsed != 12) { + continue; /* 跳过空行/格式不符的行 */ + } + int i = a.n_atoms; + a.atom_ids[i] = id; + a.masses[i] = mass; + a.radii[i] = rad; + a.pos_0[i*3+0] = px; a.pos_0[i*3+1] = py; a.pos_0[i*3+2] = pz; + a.vel_0[i*3+0] = vx; a.vel_0[i*3+1] = vy; a.vel_0[i*3+2] = vz; + a.fixed[i*3+0] = fx; a.fixed[i*3+1] = fy; a.fixed[i*3+2] = fz; + a.n_atoms++; + } + fclose(f); + + if (a.n_atoms <= 0) die("coord.txt 原子数无效"); + return a; +} + +/* 读取 connection.txt */ +static BondData read_bonds(const char *input_dir, const AtomData *atoms) { + char path[512]; + BondData b; + b.n_bonds = 0; + b.pairs = NULL; + b.stiffness = NULL; + b.rest_lengths = NULL; + + snprintf(path, sizeof(path), "%s/connection.txt", input_dir); + FILE *f = fopen(path, "r"); + if (!f) return b; /* 无成键 */ + + /* 跳过第一行表头 */ + char line[256]; + if (!fgets(line, sizeof(line), f)) { fclose(f); return b; } + + /* 先数行数 */ + int n_lines = 0, tmp_a, tmp_b; + char bond_name[256]; + while (fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name) == 3) n_lines++; + rewind(f); + + if (n_lines == 0) { fclose(f); return b; } + + b.n_bonds = n_lines; + b.pairs = (int*)xmalloc(n_lines * 2 * sizeof(int)); + b.stiffness = (double*)xmalloc(n_lines * sizeof(double)); + b.rest_lengths = (double*)xmalloc(n_lines * sizeof(double)); + + /* 读取 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; /* 默认值 */ + if (fb) { + char name[256]; + double k, r0; + rewind(fb); + while (fscanf(fb, "%s %lf %lf", name, &k, &r0) == 3) { + if (strcmp(name, bond_name) == 0) { + b.stiffness[i] = k; + b.rest_lengths[i] = r0; + break; + } + } + } + } + fclose(f); + if (fb) fclose(fb); + return b; +} + +/* ======================================================================== + * 物理核心 + * ======================================================================== */ + +static void compute_acceleration( + int n, const double *x, const double *y, const double *z, + const double *vx, const double *vy, const double *vz, + const double *m, const double G[3], const double B[3], + const BondData *bonds, + double *ax, double *ay, double *az) +{ + for (int i = 0; i < n; i++) { + 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; + } + + /* 弹簧力 */ + 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]; + } +} + +/* clamp */ +static double clamp(double v, double lo, double hi) { + return fmin(fmax(v, lo), hi); +} + +/* 蛙跳法一步 */ +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 box_a, 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; + 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; + /* 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 write_trajectory_json(const char *path, const Trajectory *traj, + const SimParams *params, const AtomData *atoms, + const BondData *bonds) +{ + FILE *f = fopen(path, "w"); + if (!f) die("无法写入 trajectory.txt"); + + fprintf(f, "{\n"); + + /* 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}; + + for (int a = 0; a < 6; a++) { + fprintf(f, " \"%s\": [\n", names[a]); + for (int t = 0; t < traj->n_steps; t++) { + fprintf(f, " ["); + for (int i = 0; i < traj->n_atoms; i++) { + fprintf(f, "%.15g", arrs[a][t * traj->n_atoms + i]); + if (i < traj->n_atoms - 1) fputc(',', f); + } + fprintf(f, "]"); + if (t < traj->n_steps - 1) fputc(',', f); + fputc('\n', f); + } + fprintf(f, " ]"); + // 每个轨迹数组后加逗号（后面还有标量参数） + fputc(',', f); + fputc('\n', f); + } + + /* 标量参数 */ + fprintf(f, " \"NT\": %d,\n", params->NT); + fprintf(f, " \"DT\": %.15g,\n", params->DT); + fprintf(f, " \"NSTEP\": %d,\n", params->NSTEP); + fprintf(f, " \"method\": \"leapfrog\",\n"); + 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++) { + 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++) { + 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++) { + 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++) { + 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++) { + fprintf(f, "%.15g", bonds->rest_lengths[b]); + if (b < bonds->n_bonds - 1) fputc(',', f); + } + fprintf(f, "]\n"); + + fprintf(f, "}\n"); + fclose(f); +} + +/* ======================================================================== + * 主函数 + * ======================================================================== */ + +int main(int argc, char **argv) { + if (argc < 4) { + fprintf(stderr, "用法: %s \n", argv[0]); + return 1; + } + const char *input_dir = argv[1]; + const char *output_dir = argv[2]; + const char *param_path = argv[3]; + + clock_t t0 = clock(); + + /* 读取参数和输入 */ + SimParams params = read_params(param_path); + AtomData atoms = read_coord(input_dir); + BondData bonds = read_bonds(input_dir, &atoms); + + printf("[C-engine] 原子数=%d, 键数=%d, NT=%d, DT=%.6g\n", + atoms.n_atoms, bonds.n_bonds, params.NT, params.DT); + + /* 初始化位置/速度 */ + int n = atoms.n_atoms; + double *x = (double*)xmalloc(n * sizeof(double)); + double *y = (double*)xmalloc(n * sizeof(double)); + double *z = (double*)xmalloc(n * sizeof(double)); + double *vx = (double*)xmalloc(n * sizeof(double)); + double *vy = (double*)xmalloc(n * sizeof(double)); + double *vz = (double*)xmalloc(n * sizeof(double)); + for (int i = 0; i < n; i++) { + x[i] = atoms.pos_0[i*3+0]; + y[i] = atoms.pos_0[i*3+1]; + z[i] = atoms.pos_0[i*3+2]; + vx[i] = atoms.vel_0[i*3+0]; + vy[i] = atoms.vel_0[i*3+1]; + vz[i] = atoms.vel_0[i*3+2]; + } + + /* 分配轨迹缓冲区 */ + Trajectory traj; + traj.n_steps = params.NT; + 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; + + /* 预热 */ + 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); + } + + /* 记录 */ + 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 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]; + } + } + + /* 输出轨迹 */ + 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); + + /* 清理 */ + free(traj.x); + free(x); free(y); free(z); + free(vx); free(vy); free(vz); + free(atoms.atom_ids); + free(atoms.masses); free(atoms.radii); + free(atoms.pos_0); free(atoms.vel_0); + free(atoms.fixed); + if (bonds.pairs) { free(bonds.pairs); free(bonds.stiffness); free(bonds.rest_lengths); } + + return 0; +} diff --git a/engines/cpp/main.cpp b/engines/cpp/main.cpp new file mode 100644 index 0000000..42be97f --- /dev/null +++ b/engines/cpp/main.cpp @@ -0,0 +1,77 @@ +/** + * engines/cpp/main.cpp + * -------------------- + * C++ 动力学模拟引擎（模板/参考实现）。 + * + * 接口与 C 引擎一致： + * 输入: param.json, /coord.txt, connection.txt, bond.txt + * 输出: /trajectory.txt (JSON) + * + * 编译: + * g++ -O3 -march=native -o build/dynamics_cpp main.cpp + * + * 用法: + * ./build/dynamics_cpp + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +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; +}; + +// ======================================================================== +// 请在下方实现 Leapfrog 算法 +// 参考 engines/c/main.c 中的物理核心实现 +// ======================================================================== + +int main(int argc, char **argv) { + if (argc < 4) { + std::cerr << "用法: " << argv[0] << " " << std::endl; + return 1; + } + + std::string input_dir = argv[1]; + std::string output_dir = argv[2]; + std::string param_path = argv[3]; + + auto t0 = std::chrono::high_resolution_clock::now(); + + // TODO: 读取 param.json、coord.txt、connection.txt、bond.txt + // TODO: 执行 Leapfrog 模拟 + // TODO: 输出 trajectory.txt (JSON 格式) + // + // 提示：输出格式与 Python 版兼容，参考 engines/c/main.c 中的 write_trajectory_json + + auto t1 = std::chrono::high_resolution_clock::now(); + double elapsed = std::chrono::duration(t1 - t0).count(); + + std::cout << "[C++-engine] 占位: NT=" << 0 << ", " << elapsed << " s (待实现)" << std::endl; + return 0; +} diff --git a/engines/fortran/main.f90 b/engines/fortran/main.f90 new file mode 100644 index 0000000..ada2964 --- /dev/null +++ b/engines/fortran/main.f90 @@ -0,0 +1,40 @@ +! engines/fortran/main.f90 +! ------------------------ +! Fortran 动力学模拟引擎（模板/参考实现）。 +! +! 接口与 C 引擎一致： +! 输入: param.json, /coord.txt, connection.txt, bond.txt +! 输出: /trajectory.txt (JSON) +! +! 编译: +! gfortran -O3 -march=native -o build/dynamics_f90 main.f90 +! +! 用法: +! ./build/dynamics_f90 + +program dynamics_f90 + implicit none + character(len=256) :: input_dir, output_dir, param_path + integer :: narg + + narg = command_argument_count() + if (narg < 3) then + write(*,*) "用法: dynamics_f90 " + stop + end if + + call get_command_argument(1, input_dir) + call get_command_argument(2, output_dir) + call get_command_argument(3, param_path) + + ! TODO: 实现与 C 引擎相同的物理模拟 + ! + ! 参考实现: + ! 1. 读取 param.json → NT, DT, G, B, ... + ! 2. 读取 coord.txt → 原子数, 质量, 位置, 速度 + ! 3. 读取 connection.txt + bond.txt → 成键参数 + ! 4. Leapfrog 循环 NT 步 + ! 5. 输出 trajectory.txt (JSON 格式) + + write(*,*) "[Fortran-engine] 待实现" +end program dynamics_f90 diff --git a/examples/case01/input/bond.txt b/examples/case01/input/bond.txt new file mode 100644 index 0000000..2b183d5 --- /dev/null +++ b/examples/case01/input/bond.txt @@ -0,0 +1,2 @@ +bond_name k rest_length +k1 100.0 2.0 diff --git a/examples/case01/input/connection.txt b/examples/case01/input/connection.txt new file mode 100644 index 0000000..689f729 --- /dev/null +++ b/examples/case01/input/connection.txt @@ -0,0 +1,2 @@ +n1 n2 bond_name +1 2 k1 diff --git a/examples/case01/input/coord.txt b/examples/case01/input/coord.txt new file mode 100644 index 0000000..9a2e805 --- /dev/null +++ b/examples/case01/input/coord.txt @@ -0,0 +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 diff --git a/examples/case01/input/parameters.yaml b/examples/case01/input/parameters.yaml new file mode 100644 index 0000000..0ca9378 --- /dev/null +++ b/examples/case01/input/parameters.yaml @@ -0,0 +1,80 @@ +# 物理模拟参数配置 +# 格式：YAML +# 用法：python run_dynamics.py + +# ── 流程控制 ────────────────────────────────── +# 每步用 0/1 单独开关，1=执行，0=跳过 +# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧 +step_simulate: 1 # 运行物理模拟 → output/trajectory.txt +step_sample: 1 # 抽帧 → output/display.txt +step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy） + +# ── 计算引擎 ────────────────────────────────── +# 可选: 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 引擎 + +# ── 盒子 ────────────────────────────────────── +box_a: 10.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] # 阻尼分量 + +# ── 数值算法 ────────────────────────────────── +# 可选： +# explicit_euler 显式欧拉法 +# implicit_euler 隐式欧拉法 +# midpoint 中点法 +# leapfrog 蛙跳法 +method: leapfrog + +# ── 步骤控制 ────────────────────────────────── +# 以下参数控制哪些步骤被执行和保存 + +# 预热步数：模拟开始时跳过不保存的步数（用于稳定初始状态） +warmup_steps: 0 # 默认 0（立即开始记录） + +# 总计算步数 +NT: 10000 + +# 抽帧间隔（每 NSTEP 步取一帧用于动画） +NSTEP: 100 + +# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧 +sample_start: null # null 表示从头开始（帧索引从 0 起） +sample_end: null # null 表示到末尾 + +# ── 时间步长 ────────────────────────────────── +DT: 0.001 # 时间步长 (s) + +# ── 显示参数 ────────────────────────────────── +# 盒子透明度：单个数值（统一）或 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/case01/run_dynamics.py b/examples/case01/run_dynamics.py new file mode 100644 index 0000000..6e1c34a --- /dev/null +++ b/examples/case01/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 / "parameters.yaml" + + +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() diff --git a/export_web_data.py b/export_web_data.py new file mode 100644 index 0000000..b00e028 --- /dev/null +++ b/export_web_data.py @@ -0,0 +1,125 @@ +""" +Export selected output/display.txt fields into browser-friendly JSON files. +""" + +from __future__ import annotations + +import json +import math +import os +import sys +from pathlib import Path + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import compute + + +SCRIPT_DIR = Path(__file__).resolve().parent +OUTPUT_DIR = Path(compute.get_output_dir(SCRIPT_DIR)) +DISPLAY_TXT = OUTPUT_DIR / "display.txt" +DISPLAY_JSON = OUTPUT_DIR / "display.json" +DISPLAY_JS = OUTPUT_DIR / "display.js" + + +def rounded_list(values, digits=6): + rounded = [] + for value in values: + if isinstance(value, float): + rounded.append(round(value, digits)) + else: + rounded.append(value) + return rounded + + +def build_payload(arrays): + disp_t = rounded_list(arrays["disp_t"]) + disp_x = rounded_list(arrays["disp_x"]) + disp_y = rounded_list(arrays["disp_y"]) + disp_z = rounded_list(arrays["disp_z"]) + disp_vx = rounded_list(arrays["disp_vx"]) + disp_vy = rounded_list(arrays["disp_vy"]) + disp_vz = rounded_list(arrays["disp_vz"]) + + z_values = [float(value) for value in arrays["disp_z"]] + speed_values = [ + math.sqrt( + float(vx) ** 2 + float(vy) ** 2 + float(vz) ** 2 + ) + for vx, vy, vz in zip(arrays["disp_vx"], arrays["disp_vy"], arrays["disp_vz"]) + ] + + return { + "metadata": { + "method": arrays["method"], + "coord_file": arrays["coord_file"], + "plot_atom_id": int(arrays["plot_atom_id"]), + "plot_atom_row": int(arrays["plot_atom_row"]), + "n_frames": int(arrays["n_frames"]), + "NT": int(arrays["NT"]), + "DT": float(arrays["DT"]), + "NSTEP": int(arrays["NSTEP"]), + "warmup_steps": int(arrays["warmup_steps"]), + "sample_start": int(arrays["sample_start"]), + "sample_end": int(arrays["sample_end"]), + "bounds": { + "x": [float(arrays["X_MIN"]), float(arrays["X_MAX"])], + "y": [float(arrays["Y_MIN"]), float(arrays["Y_MAX"])], + "z": [float(arrays["Z_MIN"]), float(arrays["Z_MAX"])], + }, + "ball_radius": float(arrays["ball_radius"]), + "ball_color": [ + float(arrays["ball_color_r"]), + float(arrays["ball_color_g"]), + float(arrays["ball_color_b"]), + ], + "box_color": [ + float(arrays["box_color_r"]), + float(arrays["box_color_g"]), + float(arrays["box_color_b"]), + ], + "alpha": float(arrays["alpha"]), + }, + "series": { + "t": disp_t, + "x": disp_x, + "y": disp_y, + "z": disp_z, + "vx": disp_vx, + "vy": disp_vy, + "vz": disp_vz, + "step": [int(value) for value in arrays["disp_step"]], + }, + "summary": { + "z_min": round(min(z_values), 6), + "z_max": round(max(z_values), 6), + "max_speed": round(max(speed_values), 6), + "final_position": [ + disp_x[-1], + disp_y[-1], + disp_z[-1], + ], + "final_velocity": [ + disp_vx[-1], + disp_vy[-1], + disp_vz[-1], + ], + }, + } + + +def main(): + OUTPUT_DIR.mkdir(exist_ok=True) + arrays = compute.load_text_data(DISPLAY_TXT) + payload = build_payload(arrays) + json_text = json.dumps(payload, ensure_ascii=False, separators=(",", ":")) + DISPLAY_JSON.write_text(json_text, encoding="utf-8") + DISPLAY_JS.write_text( + f"window.__DYNAMICS_DISPLAY__ = {json_text};\n", + encoding="utf-8", + ) + print(f"[export] Wrote {DISPLAY_JSON}") + print(f"[export] Wrote {DISPLAY_JS}") + + +if __name__ == "__main__": + main() diff --git a/migrate_npz_outputs.py b/migrate_npz_outputs.py new file mode 100644 index 0000000..b82835f --- /dev/null +++ b/migrate_npz_outputs.py @@ -0,0 +1,160 @@ +""" +migrate_npz_outputs.py +---------------------- +将旧版 `.npz` 运行产物迁移为当前使用的 `.txt` 文本格式。 + +不依赖 numpy，便于在教学机器上直接运行： + + py -3 migrate_npz_outputs.py +""" + +from __future__ import annotations + +import ast +import json +import os +import struct +import zipfile + + +BASE_DIR = os.path.dirname(os.path.abspath(__file__)) +OUTPUT_DIR = os.path.join(BASE_DIR, "output") + + +def parse_npy(raw): + """Parse a small subset of the NPY format used by this project.""" + if raw[:6] != b"\x93NUMPY": + raise ValueError("无效的 NPY 文件头") + + major = raw[6] + minor = raw[7] + if major == 1: + header_len = struct.unpack("" + return list(struct.unpack(prefix + fmt * count, payload)) + + if dtype == "i": + fmt = {"1": "b", "2": "h", "4": "i", "8": "q"}[item_chars] + size = int(item_chars) + count = len(payload) // size + prefix = "<" if endian in ("<", "|") else ">" + return list(struct.unpack(prefix + fmt * count, payload)) + + if dtype == "u": + fmt = {"1": "B", "2": "H", "4": "I", "8": "Q"}[item_chars] + size = int(item_chars) + count = len(payload) // size + prefix = "<" if endian in ("<", "|") else ">" + return list(struct.unpack(prefix + fmt * count, payload)) + + if dtype == "b": + return [bool(x) for x in payload] + + if dtype == "U": + char_count = int(item_chars) + item_size = char_count * 4 + out = [] + for offset in range(0, len(payload), item_size): + chunk = payload[offset:offset + item_size] + out.append(chunk.decode("utf-32le").rstrip("\x00")) + return out + + raise ValueError(f"暂不支持的数据类型: {descr}") + + +def reshape_values(values, shape): + if shape == (): + return values[0] + if len(shape) == 1: + return values[:shape[0]] + + step = 1 + for dim in shape[1:]: + step *= dim + + return [ + reshape_values(values[index:index + step], shape[1:]) + for index in range(0, len(values), step) + ] + + +def load_npz(path): + data = {} + with zipfile.ZipFile(path, "r") as zf: + for name in zf.namelist(): + if not name.endswith(".npy"): + continue + key = os.path.splitext(os.path.basename(name))[0] + data[key] = parse_npy(zf.read(name)) + return data + + +def dump_json(path, payload): + os.makedirs(os.path.dirname(path), exist_ok=True) + with open(path, "w", encoding="utf-8") as f: + json.dump(payload, f, ensure_ascii=False, indent=2) + + +def maybe_move_legacy_table(): + table_like = os.path.join(BASE_DIR, "trajectory.txt") + table_target = os.path.join(OUTPUT_DIR, "trajectory_table.txt") + if not os.path.exists(table_like) or os.path.exists(table_target): + return + + with open(table_like, "r", encoding="utf-8") as f: + first_line = f.readline() + if first_line.startswith("#"): + os.replace(table_like, table_target) + + +def main(): + os.makedirs(OUTPUT_DIR, exist_ok=True) + maybe_move_legacy_table() + + mappings = [ + ("trajectory.npz", os.path.join("output", "trajectory.txt")), + ("display.npz", os.path.join("output", "display.txt")), + ] + + for src_name, dst_name in mappings: + src = os.path.join(BASE_DIR, src_name) + dst = os.path.join(BASE_DIR, dst_name) + if not os.path.exists(src): + continue + payload = load_npz(src) + dump_json(dst, payload) + print(f"[migrate] 已生成: {dst_name}") + + +if __name__ == "__main__": + main() diff --git a/plot_trajectory.py b/plot_trajectory.py new file mode 100644 index 0000000..125dfba --- /dev/null +++ b/plot_trajectory.py @@ -0,0 +1,265 @@ +""" +plot_trajectory.py +------------------ +绘制轨迹数据的x-t, y-t, z-t, vx-t, vy-t, vz-t函数图像，支持多原子同时绘制。 + +用法： + python plot_trajectory.py # 绘制所有原子 + python plot_trajectory.py output/trajectory.txt # 指定文件 + python plot_trajectory.py output/trajectory.txt --atom-id 1 # 仅绘制原子1 + +参数： + trajectory_file: 轨迹数据文件，支持结构化 .txt 格式（默认：output/trajectory.txt） + --atom-id: 可选，指定只绘制单个原子的轨迹 +""" + +import numpy as np +import matplotlib.pyplot as plt +import sys +import os + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import compute + +def select_atom_series(data, atom_id=None): + """Select one atom from trajectory arrays when the file stores many atoms.""" + traj_x = data['traj_x'] + traj_y = data['traj_y'] + traj_z = data['traj_z'] + traj_vx = data['traj_vx'] + traj_vy = data['traj_vy'] + traj_vz = data['traj_vz'] + + if traj_x.ndim == 1: + selected_id = int(data["plot_atom_id"]) if "plot_atom_id" in data else 1 + return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, selected_id + + atom_ids = data["atom_ids"] if "atom_ids" in data else np.arange(traj_x.shape[1]) + 1 + if atom_id is None: + row = int(data["plot_atom_row"]) if "plot_atom_row" in data else 0 + else: + matches = np.where(atom_ids == int(atom_id))[0] + if len(matches) == 0: + raise ValueError(f"原子序号 {atom_id} 不在轨迹文件中") + row = int(matches[0]) + selected_id = int(atom_ids[row]) + return ( + traj_x[:, row], traj_y[:, row], traj_z[:, row], + traj_vx[:, row], traj_vy[:, row], traj_vz[:, row], + selected_id, + ) + + +def load_trajectory_txt(file_path, atom_id=None): + """加载结构化 .txt 格式的轨迹数据""" + data = compute.load_text_data(file_path) + traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, selected_id = select_atom_series(data, atom_id) + NT = int(data['NT']) + DT = float(data['DT']) + return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT, selected_id + + +def load_full_trajectory(file_path): + """加载完整轨迹数据（所有原子），返回 2D 数组及原始数据字典。""" + data = compute.load_text_data(file_path) + traj_x = data['traj_x'] + traj_y = data['traj_y'] + traj_z = data['traj_z'] + traj_vx = data['traj_vx'] + traj_vy = data['traj_vy'] + traj_vz = data['traj_vz'] + NT = int(data['NT']) + DT = float(data['DT']) + n_atoms = traj_x.shape[1] if traj_x.ndim > 1 else 1 + atom_ids = data.get("atom_ids", np.arange(n_atoms) + 1) + return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT, atom_ids, data + +def create_plots(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT, output_dir='.', atom_ids=None, extra_data=None): + """创建六个子图的图表，支持多原子绘制与能量曲线。""" + # 配置中文字体支持 + plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'WenQuanYi Micro Hei', 'DejaVu Sans'] + plt.rcParams['axes.unicode_minus'] = False + + # 生成时间数组 + time = np.arange(NT) * DT + + # 判断单原子（1D）或多原子（2D） + is_multi = traj_x.ndim == 2 + n_atoms = traj_x.shape[1] if is_multi else 1 + + # 是否绘制能量子图（需要 multi-atom + extra_data 中有 G） + has_energy = is_multi and extra_data is not None and "G" in extra_data + n_rows = 3 if has_energy else 2 + n_cols = 3 + figsize = (15, 13) if has_energy else (15, 9) + + # 创建图形和子图 + fig, axes = plt.subplots(n_rows, n_cols, figsize=figsize) + title = '轨迹与能量分析' if has_energy else '轨迹分析' + if atom_ids is not None and not is_multi: + title += f' - 原子 {int(atom_ids)}' + fig.suptitle(title, fontsize=16) + + plot_configs = [ + (axes[0, 0], traj_x, 'x'), + (axes[0, 1], traj_y, 'y'), + (axes[0, 2], traj_z, 'z'), + (axes[1, 0], traj_vx, 'vx'), + (axes[1, 1], traj_vy, 'vy'), + (axes[1, 2], traj_vz, 'vz'), + ] + + if is_multi: + colors = plt.cm.tab10(np.linspace(0, 1, n_atoms)) + if atom_ids is None: + atom_ids = np.arange(n_atoms) + 1 + + for ax, data_arr, label in plot_configs: + if is_multi: + for i in range(n_atoms): + aid = int(atom_ids[i]) + ax.plot(time, data_arr[:, i], color=colors[i], linewidth=1.5, label=f"原子 {aid}") + ax.legend() + else: + ax.plot(time, data_arr, 'b-', linewidth=1.5) + ax.set_title(f'{label} - 时间') + ax.set_xlabel('时间 (s)') + ax.set_ylabel(label) + ax.grid(True, alpha=0.3) + + # ── 能量子图 ───────────────────────────────────── + if has_energy: + masses = np.array(extra_data["atom_masses"]) # (n_atoms,) + G_vec = np.array(extra_data["G"]) # [gx, gy, gz] + + # 动能 Ek = ½ m v² + ek = 0.5 * masses[np.newaxis, :] * (traj_vx**2 + traj_vy**2 + traj_vz**2) + # 重力势能 Ug = -m G·r + ug = -masses[np.newaxis, :] * (G_vec[0] * traj_x + G_vec[1] * traj_y + G_vec[2] * traj_z) + # 弹性势能 Us = ½ k (d - d₀)² + us = np.zeros_like(ek) + bond_pairs = extra_data.get("bond_pairs") + bond_stiffness = extra_data.get("bond_stiffness") + bond_rest_lengths = extra_data.get("bond_rest_lengths") + if 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 = traj_x[:, j] - traj_x[:, i] + dy = traj_y[:, j] - traj_y[:, i] + dz = traj_z[:, j] - traj_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 + + e_total = ek + ug + us # (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 + + # 第 3 行左：各原子总能量 + ax_e = axes[2, 0] + for i in range(n_atoms): + aid = int(atom_ids[i]) + ax_e.plot(time, e_total[:, i], color=colors[i], linewidth=1.5, label=f"原子 {aid}") + ax_e.set_title("各原子总能量") + ax_e.set_xlabel("时间 (s)") + ax_e.set_ylabel("能量") + ax_e.grid(True, alpha=0.3) + ax_e.legend(loc="upper right") + + # 第 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, us_sys, color='orange', linewidth=1.5, label="系统弹性势能") + ax_sys.plot(time, e_sys, 'r--', linewidth=1.5, label="系统总能量") + ax_sys.set_title("系统总能量") + ax_sys.set_xlabel("时间 (s)") + ax_sys.set_ylabel("能量") + ax_sys.grid(True, alpha=0.3) + ax_sys.legend(loc="upper right") + + # 隐藏第 3 行第 3 列空子图 + axes[2, 2].set_visible(False) + + # 调整布局 + plt.tight_layout(rect=[0, 0.03, 1, 0.95]) + + # 保存图像 + output_path = os.path.join(output_dir, 'trajectory_plots.png') + plt.savefig(output_path, dpi=300, bbox_inches='tight') + print(f"图表已保存至: {output_path}") + + # 显示图像（如果环境支持） + plt.show() + + return output_path + +def main(): + import argparse + + parser = argparse.ArgumentParser(description="绘制轨迹数据图表") + parser.add_argument("trajectory_file", nargs="?", + help="轨迹数据文件路径（默认: output/trajectory.txt）") + parser.add_argument("--atom-id", type=int, default=None, + help="指定只绘制某个原子的轨迹（默认绘制所有原子）") + args = parser.parse_args() + + # 默认文件 + script_dir = os.path.dirname(os.path.abspath(__file__)) + default_file = os.path.join(compute.get_output_dir(script_dir), "trajectory.txt") + input_file = args.trajectory_file or default_file + + # 检查文件是否存在 + if not os.path.exists(input_file): + print(f"错误: 文件 '{input_file}' 不存在") + print(f"请先运行 compute.py 生成轨迹数据") + sys.exit(1) + + file_ext = os.path.splitext(input_file)[1].lower() + + try: + if args.atom_id is not None: + # 单原子模式（旧版兼容） + traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT, selected_id = \ + load_trajectory_txt(input_file, args.atom_id) + ids_for_plot = selected_id + extra_data = None + n_atoms_str = f" 绘图原子序号: {selected_id}" + else: + # 全原子模式 + traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, NT, DT, atom_ids, raw_data = \ + load_full_trajectory(input_file) + ids_for_plot = atom_ids + extra_data = raw_data + n_atoms = traj_x.shape[1] if traj_x.ndim > 1 else 1 + n_atoms_str = f" 绘图原子数: {n_atoms}" + except Exception as e: + print(f"加载文件时出错: {e}") + sys.exit(1) + + print(f"加载轨迹数据: NT={NT}, DT={DT}") + if args.atom_id is not None: + print(n_atoms_str) + else: + print(n_atoms_str) + print(f" 位置范围: x [{traj_x.min():.3f}, {traj_x.max():.3f}], " + f"y [{traj_y.min():.3f}, {traj_y.max():.3f}], " + f"z [{traj_z.min():.3f}, {traj_z.max():.3f}]") + print(f" 速度范围: vx [{traj_vx.min():.3f}, {traj_vx.max():.3f}], " + f"vy [{traj_vy.min():.3f}, {traj_vy.max():.3f}], " + f"vz [{traj_vz.min():.3f}, {traj_vz.max():.3f}]") + + # 创建图表 + output_path = create_plots( + traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, + NT, DT, compute.get_output_dir(script_dir), ids_for_plot, extra_data, + ) + + print("绘图完成！") + +if __name__ == "__main__": + main() diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..8aae06f --- /dev/null +++ b/requirements.txt @@ -0,0 +1,12 @@ +numpy>=1.20 +PyYAML>=5.4 + +# 进度条（推荐） +tqdm>=4.60 + +# 静态绘图（可选——用于生成 trajectory_plots.png） +# matplotlib>=3.4 + +# 3D 动画（可选——用于播放 VisPy 动画） +# vispy>=0.14 +# PyQt5>=5.15 diff --git a/sample.py b/sample.py new file mode 100644 index 0000000..d57339a --- /dev/null +++ b/sample.py @@ -0,0 +1,184 @@ +""" +sample.py +--------- +纯 Python 抽帧脚本，不依赖 VisPy。 + +功能： + 1. 从 output/trajectory.txt 读取完整轨迹 + 2. 每隔 NSTEP 帧抽取一帧，生成 output/display.txt + 3. output/display.txt 可直接被 draw.py 加载驱动动画 +""" + +import numpy as np +import os +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import compute + + +def build_sample_indices(total_steps, sample_step, sample_start, sample_end): + """Validate sampling settings and build frame indices.""" + if sample_step <= 0: + raise ValueError(f"NSTEP 必须为正整数，实际为 {sample_step}") + if sample_start < 0: + raise ValueError(f"sample_start 不能小于 0，实际为 {sample_start}") + if sample_end > total_steps: + raise ValueError( + f"sample_end 不能大于记录步数 {total_steps}，实际为 {sample_end}") + if sample_start >= sample_end: + raise ValueError( + f"sample_start 必须小于 sample_end，实际为 [{sample_start}, {sample_end})") + + n_frames = (sample_end - sample_start) // sample_step + if n_frames <= 0: + raise ValueError( + f"抽帧范围 [{sample_start}, {sample_end}) 过短，按 NSTEP={sample_step} 无法抽出任何帧") + + return np.arange(n_frames, dtype=np.int64) * sample_step + sample_start + + +def read_optional_index(data, key, default_value): + """Read an optional integer index from txt metadata.""" + if key not in data: + return default_value + value = data[key] + if value is None or int(value) < 0: + return default_value + return int(value) + + +def main(): + script_dir = os.path.dirname(os.path.abspath(__file__)) + output_dir = compute.get_output_dir(script_dir) + traj_path = os.path.join(output_dir, "trajectory.txt") + disp_path = os.path.join(output_dir, "display.txt") + + # ------------------------------------------------------------------------- + # 1. 加载完整轨迹 + # ------------------------------------------------------------------------- + data = compute.load_text_data(traj_path) + NT = int(data["NT"]) + DT = float(data["DT"]) + NSTEP = int(data["NSTEP"]) + traj_x = data["traj_x"] + traj_y = data["traj_y"] + traj_z = data["traj_z"] + traj_vx = data["traj_vx"] + traj_vy = data["traj_vy"] + traj_vz = data["traj_vz"] + plot_atom_row = int(data["plot_atom_row"]) if "plot_atom_row" in data else 0 + plot_atom_id = int(data["plot_atom_id"]) if "plot_atom_id" in data else int(data["atom_ids"][plot_atom_row]) + + print(f"[sample] 加载轨迹数据: NT={NT}, DT={DT}, NSTEP={NSTEP}") + + # ------------------------------------------------------------------------- + # 2. 抽帧：支持配置文件中保存的 [sample_start, sample_end) 区间 + # ------------------------------------------------------------------------- + sample_start = read_optional_index(data, "sample_start", 0) + sample_end = read_optional_index(data, "sample_end", NT) + indices = build_sample_indices(NT, NSTEP, sample_start, sample_end) + n_frames = len(indices) + + if traj_x.ndim == 1: + selected_x = traj_x + selected_y = traj_y + selected_z = traj_z + selected_vx = traj_vx + selected_vy = traj_vy + selected_vz = traj_vz + all_x = traj_x[:, None] + all_y = traj_y[:, None] + all_z = traj_z[:, None] + all_vx = traj_vx[:, None] + all_vy = traj_vy[:, None] + all_vz = traj_vz[:, None] + else: + selected_x = traj_x[:, plot_atom_row] + selected_y = traj_y[:, plot_atom_row] + selected_z = traj_z[:, plot_atom_row] + selected_vx = traj_vx[:, plot_atom_row] + selected_vy = traj_vy[:, plot_atom_row] + selected_vz = traj_vz[:, plot_atom_row] + all_x = traj_x + all_y = traj_y + all_z = traj_z + all_vx = traj_vx + all_vy = traj_vy + all_vz = traj_vz + + disp_x = selected_x[indices] + disp_y = selected_y[indices] + disp_z = selected_z[indices] + disp_vx = selected_vx[indices] + disp_vy = selected_vy[indices] + disp_vz = selected_vz[indices] + disp_t = indices * DT # 物理时间 = step * DT + disp_step = indices # 对应保存轨迹的步编号（0-based） + + print(f"[sample] 抽帧完成: [{sample_start}, {sample_end}) -> {n_frames} 帧") + print(f"[sample] 每帧时间跨度: {NSTEP*DT:.3f} s (即每隔 {NSTEP} 步取一帧)") + + # ------------------------------------------------------------------------- + # 3. 保存显示数组 + # ------------------------------------------------------------------------- + payload = { + "disp_x": disp_x, + "disp_y": disp_y, + "disp_z": disp_z, + "disp_vx": disp_vx, + "disp_vy": disp_vy, + "disp_vz": disp_vz, + "disp_all_x": all_x[indices], + "disp_all_y": all_y[indices], + "disp_all_z": all_z[indices], + "disp_all_vx": all_vx[indices], + "disp_all_vy": all_vy[indices], + "disp_all_vz": all_vz[indices], + "disp_t": disp_t, + "disp_step": disp_step, + "n_frames": n_frames, + "NT": NT, + "DT": DT, + "NSTEP": NSTEP, + "plot_atom_id": plot_atom_id, + "plot_atom_row": plot_atom_row, + "method": str(data["method"]) if "method" in data else "explicit_euler", + "coord_file": str(data["coord_file"]) if "coord_file" in data else os.path.join("input", "coord.txt"), + "atom_ids": data["atom_ids"] if "atom_ids" in data else np.array([1]), + "atom_masses": data["atom_masses"] if "atom_masses" in data else np.array([float(data["M"])]), + "atom_radii": data["atom_radii"] if "atom_radii" in data else np.array([float(data["ball_radius"])]), + "atom_positions": data["atom_positions"] if "atom_positions" in data else np.array([[float(data["X0"]), float(data["Y0"]), float(data["Z0"])]]), + "atom_velocities": data["atom_velocities"] if "atom_velocities" in data else np.array([[float(data["VX0"]), float(data["VY0"]), float(data["VZ0"])]]), + "atom_fixed": data["atom_fixed"] if "atom_fixed" in data else np.array([[0, 0, 0]]), + "warmup_steps": int(data["warmup_steps"]) if "warmup_steps" in data else 0, + "sample_start": sample_start, + "sample_end": sample_end, + "X_MIN": float(data["X_MIN"]), + "X_MAX": float(data["X_MAX"]), + "Y_MIN": float(data["Y_MIN"]), + "Y_MAX": float(data["Y_MAX"]), + "Z_MIN": float(data["Z_MIN"]), + "Z_MAX": float(data["Z_MAX"]), + "X0": float(data["X0"]), + "Y0": float(data["Y0"]), + "Z0": float(data["Z0"]), + "VX0": float(data["VX0"]), + "VY0": float(data["VY0"]), + "VZ0": float(data["VZ0"]), + "M": float(data["M"]) if "M" in data else 1.0, + "alpha": float(data["alpha"]), + "ball_radius": float(data["ball_radius"]), + "ball_color_r": float(data["ball_color_r"]), + "ball_color_g": float(data["ball_color_g"]), + "ball_color_b": float(data["ball_color_b"]), + "box_color_r": float(data["box_color_r"]), + "box_color_g": float(data["box_color_g"]), + "box_color_b": float(data["box_color_b"]), + } + compute.save_text_data(disp_path, payload) + print(f"[sample] 显示数组已保存至: {disp_path}") + + +if __name__ == "__main__": + main()