modified: CMakeLists.txt
modified: INSTALL.md modified: README.md modified: build_release_zip.py modified: compute.py new file: doc/index.html modified: dynamics.py modified: engines/c/main.c modified: engines/cpp/main.cpp modified: engines/fortran/main.f90 modified: examples/case01/input/coord.txt renamed: examples/case01/input/parameters.yaml -> examples/case01/input/input.txt modified: examples/case01/run_dynamics.py new file: examples/case02/input/bond.txt new file: examples/case02/input/connection.txt new file: examples/case02/input/coord.txt new file: examples/case02/input/input.txt new file: examples/case02/run_dynamics.py
This commit is contained in:
+4
-2
@@ -53,9 +53,11 @@ function(add_engine_target name lang src)
|
||||
add_executable(${name} ${src})
|
||||
endif()
|
||||
|
||||
# 输出到 engines/<name>/build/（Python 调度器期望的位置）
|
||||
# 从源码路径推导输出目录：engines/c/main.c → engines/c/build/
|
||||
# 这样 Python 的 run_engine() 能在固定路径找到可执行文件
|
||||
string(REGEX REPLACE "/[^/]+$" "/build" _out_subdir "${src}")
|
||||
set_target_properties(${name} PROPERTIES
|
||||
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/engines/${name}/build"
|
||||
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/${_out_subdir}"
|
||||
)
|
||||
# macOS 上禁止生成 .app 包（命令行工具不需要）
|
||||
set_target_properties(${name} PROPERTIES
|
||||
|
||||
+2
-2
@@ -461,7 +461,7 @@ py -3 examples/case01/run_dynamics.py
|
||||
### C 引擎验证
|
||||
|
||||
```bash
|
||||
# 修改 parameters.yaml 中 engine: python → engine: c
|
||||
# 修改 input.txt 中 engine: python → engine: c
|
||||
# 或用命令行指定
|
||||
py -3 examples/case01/run_dynamics.py
|
||||
```
|
||||
@@ -543,7 +543,7 @@ pip install -i https://pypi.tuna.tsinghua.edu.cn/simple numpy pyyaml tqdm matplo
|
||||
### Q: 如何使用 C 引擎？
|
||||
|
||||
```yaml
|
||||
# 1. 修改 examples/case01/input/parameters.yaml
|
||||
# 1. 修改 examples/case01/input/input.txt
|
||||
engine: c
|
||||
|
||||
# 2. 先编译
|
||||
|
||||
@@ -3,7 +3,7 @@
|
||||
多语言动力学数值积分与轨迹可视化框架。
|
||||
|
||||
支持 Python、C、C++、Fortran 多种语言实现同一套物理引擎，共享同一套可视化管线
|
||||
（抽帧、绘图、3D 动画），直接在 `parameters.yaml` 中切换引擎即可对比性能。
|
||||
（抽帧、绘图、3D 动画），直接在 `input.txt` 中切换引擎即可对比性能。
|
||||
|
||||
## 快速开始
|
||||
|
||||
@@ -46,7 +46,7 @@ Dynamics/
|
||||
├── examples/
|
||||
│ └── case01/
|
||||
│ ├── input/
|
||||
│ │ ├── parameters.yaml ← 唯一配置文件
|
||||
│ │ ├── input.txt ← 唯一配置文件（YAML 格式）
|
||||
│ │ ├── coord.txt
|
||||
│ │ ├── connection.txt
|
||||
│ │ └── bond.txt
|
||||
@@ -59,7 +59,7 @@ Dynamics/
|
||||
└── .workbuddy/memory/ ← 工作记忆（AI 辅助开发用）
|
||||
```
|
||||
|
||||
## 配置文件 (`parameters.yaml`)
|
||||
## 配置文件 (`input.txt` — YAML 格式)
|
||||
|
||||
所有控制集中在一个文件里：
|
||||
|
||||
@@ -91,7 +91,7 @@ alpha: [0.0, 0.0, 0.2, 0.2, 0.0, 0.0]
|
||||
### 一键切换
|
||||
|
||||
```yaml
|
||||
# parameters.yaml 中改一行
|
||||
# input.txt 中改一行
|
||||
engine: c # 切换到 C 引擎
|
||||
engine: python # 切回 Python
|
||||
```
|
||||
@@ -187,7 +187,7 @@ pip install numpy pyyaml tqdm matplotlib vispy PyQt5
|
||||
## 工作原理
|
||||
|
||||
```
|
||||
parameters.yaml
|
||||
input.txt
|
||||
│
|
||||
▼
|
||||
dynamics.py ──→ run_from_config() 加载全局参数
|
||||
|
||||
@@ -19,7 +19,7 @@ INCLUDE_FILES = [
|
||||
"sample.py",
|
||||
"draw.py",
|
||||
"plot_trajectory.py",
|
||||
"examples\\case01\\input\\parameters.yaml",
|
||||
"examples\\case01\\input\\input.txt",
|
||||
"examples\\case01\\input\\coord.txt",
|
||||
]
|
||||
|
||||
|
||||
+157
-35
@@ -64,6 +64,13 @@ box_color_r = None
|
||||
box_color_g = None
|
||||
box_color_b = None
|
||||
|
||||
# 力开关
|
||||
GRAVITY_FIELD = 1 # 均匀重力场
|
||||
GRAVITY_INTERACTION = 0 # 原子间万有引力
|
||||
ELASTIC_FORCE = 1 # 弹簧键力
|
||||
DAMPING_FORCE = 0 # 阻尼
|
||||
GRAVITY_STRENGTH = 1.0
|
||||
|
||||
# 派生边界（根据 box_a 计算）
|
||||
X_MIN = None
|
||||
X_MAX = None
|
||||
@@ -156,9 +163,11 @@ def load_coord_file(coord_path):
|
||||
header = f.readline().strip().split()
|
||||
expected = ["n", "mass", "radius", "x", "y", "z", "vx", "vy", "vz"]
|
||||
legacy = expected + ["fixed_x", "fixed_y", "fixed_z"]
|
||||
if header not in (expected, legacy):
|
||||
legacy_new = expected + ["fix_x", "fix_y", "fix_z"]
|
||||
if header not in (expected, legacy, legacy_new):
|
||||
raise ValueError(
|
||||
f"坐标文件表头应为: {' '.join(expected)}，实际为: {' '.join(header)}")
|
||||
f"坐标文件表头应为: {' '.join(expected)} 或加三列 fix_x fix_y fix_z，"
|
||||
f"实际为: {' '.join(header)}")
|
||||
|
||||
for line_no, line in enumerate(f, start=2):
|
||||
stripped = line.strip()
|
||||
@@ -167,7 +176,7 @@ def load_coord_file(coord_path):
|
||||
parts = stripped.split()
|
||||
if header == expected and len(parts) != 9:
|
||||
raise ValueError(f"{coord_path}:{line_no} 应有 9 列，实际为 {len(parts)} 列")
|
||||
if header == legacy and len(parts) != 12:
|
||||
if header in (legacy, legacy_new) and len(parts) != 12:
|
||||
raise ValueError(f"{coord_path}:{line_no} 应有 12 列，实际为 {len(parts)} 列")
|
||||
rows.append(parts)
|
||||
|
||||
@@ -182,7 +191,7 @@ def load_coord_file(coord_path):
|
||||
velocities = np.array([[float(row[6]), float(row[7]), float(row[8])] for row in rows],
|
||||
dtype=np.float64)
|
||||
fixed = np.zeros((len(rows), 3), dtype=np.int64)
|
||||
if header == legacy:
|
||||
if header in (legacy, legacy_new):
|
||||
fixed = np.array([[int(row[9]), int(row[10]), int(row[11])] for row in rows],
|
||||
dtype=np.int64)
|
||||
|
||||
@@ -369,6 +378,7 @@ def run_from_config(config, out_dir=None):
|
||||
global ball_radius, ball_color_r, ball_color_g, ball_color_b
|
||||
global box_color_r, box_color_g, box_color_b
|
||||
global warmup_steps, sample_start, sample_end
|
||||
global GRAVITY_FIELD, GRAVITY_INTERACTION, ELASTIC_FORCE, DAMPING_FORCE, GRAVITY_STRENGTH
|
||||
|
||||
box_a = float(config.get("box_a", 10.0))
|
||||
raw_alpha = config.get("alpha", 0.2)
|
||||
@@ -427,6 +437,14 @@ def run_from_config(config, out_dir=None):
|
||||
box_color_g = float(config.get("box_color_g", 0.8))
|
||||
box_color_b = float(config.get("box_color_b", 0.85))
|
||||
|
||||
# 力开关
|
||||
global GRAVITY_FIELD, GRAVITY_INTERACTION, ELASTIC_FORCE, DAMPING_FORCE, GRAVITY_STRENGTH
|
||||
GRAVITY_FIELD = int(config.get("gravity_field", 1))
|
||||
GRAVITY_INTERACTION = int(config.get("gravity_interaction", 0))
|
||||
ELASTIC_FORCE = int(config.get("elastic_force", 1))
|
||||
DAMPING_FORCE = int(config.get("damping_force", 0))
|
||||
GRAVITY_STRENGTH = float(config.get("gravity_strength", 1.0))
|
||||
|
||||
print(f"[compute] 使用算法: {METHOD}")
|
||||
print(f"[compute] 已加载成键信息: {len(BOND_PAIRS)} 条键")
|
||||
if config.get("_skip_run", False):
|
||||
@@ -470,7 +488,8 @@ def run_from_config(config, out_dir=None):
|
||||
f.write(f" mass: {M}\n")
|
||||
f.write("=" * 50 + "\n")
|
||||
print(f"[compute] 日志已保存至: {log_path}")
|
||||
print(f"[compute] 计算耗时: {elapsed:.3f} s")
|
||||
record = len(traj_x)
|
||||
print(f"[compute] 计算完成: {record} 步, {elapsed:.3f} s")
|
||||
|
||||
return traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz
|
||||
|
||||
@@ -525,8 +544,14 @@ def run_engine(engine, input_dir, output_dir, config):
|
||||
"DT": float(config["DT"]),
|
||||
"NSTEP": int(config.get("NSTEP", 1)),
|
||||
"warmup_steps": int(config.get("warmup_steps", 0)),
|
||||
"method": normalize_method_name(config.get("method", "leapfrog")),
|
||||
"G": [float(v) for v in G],
|
||||
"B": [float(v) for v in B],
|
||||
"gravity_field": int(config.get("gravity_field", 1)),
|
||||
"gravity_interaction": int(config.get("gravity_interaction", 0)),
|
||||
"elastic_force": int(config.get("elastic_force", 1)),
|
||||
"damping_force": int(config.get("damping_force", 0)),
|
||||
"gravity_strength": float(config.get("gravity_strength", 1.0)),
|
||||
}
|
||||
param_path = os.path.join(script_dir, "engines", engine, "param.json")
|
||||
os.makedirs(os.path.dirname(param_path), exist_ok=True)
|
||||
@@ -540,24 +565,87 @@ def run_engine(engine, input_dir, output_dir, config):
|
||||
print(f"[compute] input: {input_dir}")
|
||||
print(f"[compute] output: {output_dir}")
|
||||
|
||||
# ── 预校准：跑 NT=1000 步测速 ────────────────────────────────
|
||||
total_steps = int(config["NT"]) - int(config.get("warmup_steps", 0))
|
||||
_calib_nt = min(1000, max(100, total_steps // 10))
|
||||
_calib_param = dict(param_json)
|
||||
_calib_param["NT"] = _calib_nt
|
||||
_calib_path = os.path.join(script_dir, "engines", engine, "_calib.json")
|
||||
with open(_calib_path, "w", encoding="utf-8") as _cf:
|
||||
json.dump(_calib_param, _cf, indent=2)
|
||||
_ct0 = time.time()
|
||||
subprocess.run(
|
||||
[engine_path, os.path.abspath(input_dir), os.devnull, _calib_path],
|
||||
capture_output=True, timeout=60)
|
||||
os.remove(_calib_path)
|
||||
_calib_elapsed = max(time.time() - _ct0, 0.001)
|
||||
_overhead = _calib_elapsed * 0.15
|
||||
_step_time = max(_calib_elapsed - _overhead, 0.0001) / _calib_nt
|
||||
_est_total = max(_calib_elapsed, _overhead + _step_time * total_steps)
|
||||
|
||||
t_start = time.time()
|
||||
t_start_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
|
||||
result = subprocess.run(
|
||||
[engine_path, os.path.abspath(input_dir), os.path.abspath(output_dir), param_path],
|
||||
capture_output=True, text=True, timeout=600)
|
||||
t_end = time.time()
|
||||
elapsed = t_end - t_start
|
||||
t_end_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
|
||||
|
||||
# 打印引擎输出
|
||||
if result.stdout:
|
||||
for line in result.stdout.strip().split("\n"):
|
||||
line = line.strip()
|
||||
if line:
|
||||
print(f" {line}")
|
||||
if result.returncode != 0:
|
||||
print(f"[compute] 引擎错误:\n{result.stderr}")
|
||||
raise RuntimeError(f"引擎 {engine} 返回错误码 {result.returncode}")
|
||||
_p = subprocess.Popen(
|
||||
[engine_path, os.path.abspath(input_dir), os.path.abspath(output_dir), param_path],
|
||||
stdout=subprocess.PIPE, stderr=subprocess.PIPE,
|
||||
text=True, encoding='utf-8', errors='replace')
|
||||
_engine_lines = []
|
||||
|
||||
try:
|
||||
from tqdm import tqdm as _tqdm
|
||||
_pbar = _tqdm(total=total_steps, desc=f"[compute] 引擎 {engine}",
|
||||
unit="步", bar_format='{l_bar}{bar}| {n_fmt}/{total_fmt} [{elapsed}<{remaining}]')
|
||||
except ImportError:
|
||||
_pbar = None
|
||||
|
||||
try:
|
||||
while True:
|
||||
_line = _p.stdout.readline() if _p.stdout else ''
|
||||
if _line:
|
||||
_line = _line.strip()
|
||||
if _line:
|
||||
_engine_lines.append(_line)
|
||||
if _p.poll() is not None:
|
||||
if _p.stdout:
|
||||
for _r in _p.stdout:
|
||||
_r = _r.strip()
|
||||
if _r:
|
||||
_engine_lines.append(_r)
|
||||
break
|
||||
if _pbar is not None and _est_total > 0:
|
||||
_pbar.n = int(min((time.time() - t_start) / _est_total, 0.99) * total_steps)
|
||||
_pbar.refresh()
|
||||
time.sleep(0.2)
|
||||
finally:
|
||||
if _pbar is not None:
|
||||
_pbar.n = total_steps
|
||||
_pbar.refresh()
|
||||
_pbar.close()
|
||||
if _p.poll() is None:
|
||||
_p.kill()
|
||||
_p.wait(timeout=5)
|
||||
|
||||
_rc = _p.returncode
|
||||
t_end = time.time()
|
||||
t_end_str = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
|
||||
elapsed = t_end - t_start
|
||||
|
||||
if _rc != 0:
|
||||
_err = _p.stderr.read() if _p.stderr else ''
|
||||
if _err:
|
||||
print(f"[compute] 引擎错误:\n{_err}")
|
||||
raise RuntimeError(f"引擎 {engine} 返回错误码 {_rc}")
|
||||
|
||||
if _engine_lines:
|
||||
_log_path = os.path.join(output_dir, "dynamics.log")
|
||||
try:
|
||||
with open(_log_path, "a", encoding="utf-8") as _lf:
|
||||
_lf.write("\n--- 引擎输出 ---\n")
|
||||
for _ln in _engine_lines:
|
||||
_lf.write(_ln + "\n")
|
||||
except OSError:
|
||||
pass
|
||||
|
||||
# 加载输出的 trajectory.txt
|
||||
traj_path = os.path.join(os.path.abspath(output_dir), "trajectory.txt")
|
||||
@@ -636,6 +724,11 @@ def save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, out_d
|
||||
"box_color_r": box_color_r,
|
||||
"box_color_g": box_color_g,
|
||||
"box_color_b": box_color_b,
|
||||
"gravity_field": GRAVITY_FIELD,
|
||||
"gravity_interaction": GRAVITY_INTERACTION,
|
||||
"elastic_force": ELASTIC_FORCE,
|
||||
"damping_force": DAMPING_FORCE,
|
||||
"gravity_strength": GRAVITY_STRENGTH,
|
||||
}
|
||||
save_text_data(out_path, payload)
|
||||
print(f"[compute] 轨迹数据已保存至: {out_path}")
|
||||
@@ -752,16 +845,27 @@ def normalize_method_name(method):
|
||||
def compute_force(x, y, z, vx, vy, vz, m, g, b):
|
||||
"""Compute total force from the current state.
|
||||
|
||||
Current model:
|
||||
- gravity: F = m * g_vec
|
||||
- linear drag: F_drag = -B_vec * v
|
||||
- spring bonds: Hooke force based on bonded pair distance
|
||||
Each force type is independently controlled by a global switch:
|
||||
- GRAVITY_FIELD: uniform gravity F = m * g_vec
|
||||
- DAMPING_FORCE: linear drag F_drag = -B_vec * v
|
||||
- ELASTIC_FORCE: spring bonds based on bonded pair distance
|
||||
- GRAVITY_INTERACTION: Newtonian gravity between atom pairs
|
||||
"""
|
||||
fx = m * g[0] - b[0] * vx
|
||||
fy = m * g[1] - b[1] * vy
|
||||
fz = m * g[2] - b[2] * vz
|
||||
fx = np.zeros_like(x)
|
||||
fy = np.zeros_like(y)
|
||||
fz = np.zeros_like(z)
|
||||
|
||||
if BOND_PAIRS is not None and len(BOND_PAIRS) > 0:
|
||||
if GRAVITY_FIELD:
|
||||
fx += m * g[0]
|
||||
fy += m * g[1]
|
||||
fz += m * g[2]
|
||||
|
||||
if DAMPING_FORCE:
|
||||
fx -= b[0] * vx
|
||||
fy -= b[1] * vy
|
||||
fz -= b[2] * vz
|
||||
|
||||
if ELASTIC_FORCE and BOND_PAIRS is not None and len(BOND_PAIRS) > 0:
|
||||
for bond_idx, (idx_1, idx_2) in enumerate(BOND_PAIRS):
|
||||
dx = x[idx_2] - x[idx_1]
|
||||
dy = y[idx_2] - y[idx_1]
|
||||
@@ -787,6 +891,25 @@ def compute_force(x, y, z, vx, vy, vz, m, g, b):
|
||||
fy[idx_2] -= fy_bond
|
||||
fz[idx_2] -= fz_bond
|
||||
|
||||
if GRAVITY_INTERACTION:
|
||||
n = len(m)
|
||||
for i in range(n):
|
||||
for j in range(i + 1, n):
|
||||
dx = x[j] - x[i]
|
||||
dy = y[j] - y[i]
|
||||
dz = z[j] - z[i]
|
||||
r2 = dx * dx + dy * dy + dz * dz
|
||||
if r2 <= 1e-12:
|
||||
continue
|
||||
mi = m[i]; mj = m[j]
|
||||
f_mag = GRAVITY_STRENGTH * mi * mj / r2
|
||||
r = np.sqrt(r2)
|
||||
fx_i = f_mag * dx / r
|
||||
fy_i = f_mag * dy / r
|
||||
fz_i = f_mag * dz / r
|
||||
fx[i] += fx_i; fy[i] += fy_i; fz[i] += fz_i
|
||||
fx[j] -= fx_i; fy[j] -= fy_i; fz[j] -= fz_i
|
||||
|
||||
return fx, fy, fz
|
||||
|
||||
|
||||
@@ -857,14 +980,13 @@ def Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b):
|
||||
y = y + vy_half * dt
|
||||
z = z + vz_half * dt
|
||||
|
||||
gamma_x = b[0] / m
|
||||
gamma_y = b[1] / m
|
||||
gamma_z = b[2] / m
|
||||
vx_next = (vx_half + 0.5 * g[0] * dt) / (1.0 + 0.5 * gamma_x * dt)
|
||||
vy_next = (vy_half + 0.5 * g[1] * dt) / (1.0 + 0.5 * gamma_y * dt)
|
||||
vz_next = (vz_half + 0.5 * g[2] * dt) / (1.0 + 0.5 * gamma_z * dt)
|
||||
# 显式预测器：用原始加速度外推半步（标准 Velocity-Verlet 预测步）
|
||||
# v_pred = v_half + 0.5 * a(t)*dt，包含重力+阻尼+弹簧的所有贡献
|
||||
vx_pred = vx_half + 0.5 * ax * dt
|
||||
vy_pred = vy_half + 0.5 * ay * dt
|
||||
vz_pred = vz_half + 0.5 * az * dt
|
||||
ax_next, ay_next, az_next = compute_acceleration(
|
||||
x, y, z, vx_next, vy_next, vz_next, m, g, b)
|
||||
x, y, z, vx_pred, vy_pred, vz_pred, m, g, b)
|
||||
vx = vx_half + 0.5 * ax_next * dt
|
||||
vy = vy_half + 0.5 * ay_next * dt
|
||||
vz = vz_half + 0.5 * az_next * dt
|
||||
|
||||
+772
@@ -0,0 +1,772 @@
|
||||
<!DOCTYPE html>
|
||||
<html lang="zh-CN">
|
||||
<head>
|
||||
<meta charset="UTF-8">
|
||||
<meta name="viewport" content="width=device-width, initial-scale=1.0">
|
||||
<title>Dynamics — 多语言蛙跳法实现对比分析</title>
|
||||
<style>
|
||||
:root {
|
||||
--bg: #f8f9fa; --card: #fff; --text: #1a1a2e;
|
||||
--accent: #4361ee; --accent2: #3a0ca3;
|
||||
--good: #06d6a0; --warn: #ffd166; --bad: #ef476f;
|
||||
--code-bg: #1e1e2e; --code-text: #cdd6f4;
|
||||
--border: #dee2e6;
|
||||
}
|
||||
* { margin: 0; padding: 0; box-sizing: border-box; }
|
||||
body {
|
||||
font-family: -apple-system, "Noto Sans SC", "Microsoft YaHei", sans-serif;
|
||||
background: var(--bg); color: var(--text); line-height: 1.7;
|
||||
padding: 0;
|
||||
}
|
||||
header {
|
||||
background: linear-gradient(135deg, var(--accent2), var(--accent));
|
||||
color: #fff; padding: 40px 20px; text-align: center;
|
||||
}
|
||||
header h1 { font-size: 2em; margin-bottom: 8px; }
|
||||
header p { opacity: .85; font-size: .95em; }
|
||||
nav {
|
||||
background: var(--card); border-bottom: 1px solid var(--border);
|
||||
padding: 12px 20px; position: sticky; top: 0; z-index: 100;
|
||||
display: flex; flex-wrap: wrap; gap: 8px;
|
||||
}
|
||||
nav a {
|
||||
color: var(--accent); text-decoration: none; font-size: .88em;
|
||||
padding: 4px 12px; border-radius: 12px;
|
||||
border: 1px solid var(--border);
|
||||
}
|
||||
nav a:hover { background: var(--accent); color: #fff; border-color: var(--accent); }
|
||||
.container { max-width: 1000px; margin: 0 auto; padding: 24px 20px; }
|
||||
section { background: var(--card); border-radius: 12px; padding: 24px; margin-bottom: 20px; box-shadow: 0 1px 3px rgba(0,0,0,.08); }
|
||||
h2 { font-size: 1.4em; margin-bottom: 16px; color: var(--accent); border-left: 4px solid var(--accent); padding-left: 12px; }
|
||||
h3 { font-size: 1.1em; margin: 20px 0 10px; color: var(--accent2); }
|
||||
.file-path {
|
||||
font-family: "JetBrains Mono", "Fira Code", monospace;
|
||||
font-size: .82em; background: #eef; padding: 2px 8px;
|
||||
border-radius: 4px; color: var(--accent2);
|
||||
}
|
||||
pre {
|
||||
background: var(--code-bg); color: var(--code-text);
|
||||
padding: 16px; border-radius: 8px; overflow-x: auto;
|
||||
font-family: "JetBrains Mono", "Fira Code", "Cascadia Code", monospace;
|
||||
font-size: .82em; line-height: 1.6; margin: 10px 0;
|
||||
position: relative;
|
||||
}
|
||||
pre .ln { color: #6c7086; user-select: none; margin-right: 16px; }
|
||||
.lang-tag {
|
||||
position: absolute; top: 6px; right: 10px;
|
||||
font-size: .72em; background: #333; color: #ccc;
|
||||
padding: 2px 8px; border-radius: 4px;
|
||||
}
|
||||
|
||||
.analysis-box { background: #f0f4ff; border-left: 4px solid var(--accent); padding: 14px 18px; margin: 12px 0; border-radius: 0 8px 8px 0; }
|
||||
.analysis-box.warn { background: #fff8e1; border-color: var(--warn); }
|
||||
.analysis-box.bad { background: #ffe8ec; border-color: var(--bad); }
|
||||
.analysis-box.good { background: #e8faf1; border-color: var(--good); }
|
||||
.analysis-box h4 { margin-bottom: 6px; font-size: .95em; }
|
||||
|
||||
.compare-table { width: 100%; border-collapse: collapse; margin: 12px 0; font-size: .9em; }
|
||||
.compare-table th, .compare-table td { border: 1px solid var(--border); padding: 10px 12px; text-align: left; }
|
||||
.compare-table th { background: #eef; font-weight: 600; }
|
||||
.compare-table tr:nth-child(even) { background: #f8f9fa; }
|
||||
.tag { display: inline-block; padding: 2px 8px; border-radius: 4px; font-size: .78em; font-weight: 600; }
|
||||
.tag.good { background: #d4edda; color: #155724; }
|
||||
.tag.warn { background: #fff3cd; color: #856404; }
|
||||
.tag.bad { background: #f8d7da; color: #721c24; }
|
||||
|
||||
footer { text-align: center; padding: 20px; color: #888; font-size: .85em; }
|
||||
</style>
|
||||
</head>
|
||||
<body>
|
||||
|
||||
<header>
|
||||
<h1>Dynamics — 多语言蛙跳法实现对比分析</h1>
|
||||
<p>Python · C · C++ · Fortran 四个引擎的蛙跳 (Velocity-Verlet) 算法与边界条件实现对比</p>
|
||||
<p style="margin-top:6px;opacity:.8;font-size:.85em;">最后更新：2026-05-20（力开关系统、万有引力、逐自由度约束）</p>
|
||||
</header>
|
||||
|
||||
<nav>
|
||||
<a href="#overview">概览</a>
|
||||
<a href="#py">Python</a>
|
||||
<a href="#forces">力开关</a>
|
||||
<a href="#c">C</a>
|
||||
<a href="#cpp">C++</a>
|
||||
<a href="#fortran">Fortran</a>
|
||||
<a href="#diff">版本变更</a>
|
||||
<a href="#bc">边界条件</a>
|
||||
<a href="#summary">总结</a>
|
||||
</nav>
|
||||
|
||||
<div class="container">
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="overview">
|
||||
<h2>1. 概览</h2>
|
||||
<p>本框架实现了四套语言的计算引擎。各引擎共享相同的物理模型：</p>
|
||||
<ul>
|
||||
<li>均匀重力场：<em>F<sub>g</sub></em> = <em>m</em> · <em>g</em>（常数矢量）—— 独立开关 <code>gravity_field</code></li>
|
||||
<li>线性阻尼：<em>F<sub>d</sub></em> = −<em>B</em> · <em>v</em> —— 独立开关 <code>damping_force</code></li>
|
||||
<li>弹簧键：胡克定律 <em>F</em> = <em>k</em>(<em>r</em> − <em>r</em><sub>0</sub>) —— 独立开关 <code>elastic_force</code></li>
|
||||
<li>万有引力（原子间）：<em>F</em> = <em>G</em> · <em>m</em><sub>i</sub> · <em>m</em><sub>j</sub> / <em>r</em>² —— 独立开关 <code>gravity_interaction</code></li>
|
||||
<li>边界：硬壁反射，粒子限制在 [−<em>a</em>, <em>a</em>]<sup>3</sup> 内</li>
|
||||
<li>逐自由度固定约束：coord.txt 中 <code>fix_x fix_y fix_z</code> 列控制各自由度是否参与迭代</li>
|
||||
</ul>
|
||||
<p>Python 引擎支持 4 种积分算法，C/C++/Fortran 引擎于 2026-05-20 重写，现支持全部 4 种算法，且边界条件与 Python 完全一致。经验证，四种引擎对于同一输入（case01, NT=1000, leapfrog, B=0）输出完全一致（至双精度浮点精度）。</p>
|
||||
|
||||
<h3>分析范围</h3>
|
||||
<ul>
|
||||
<li>蛙跳法（Velocity-Verlet）积分器的四语言实现代码</li>
|
||||
<li>边界条件（硬壁反射：clamp 位置 + 反转速度）的实现</li>
|
||||
<li>2026-05-20 重写修复的 Bug 及版本差异</li>
|
||||
</ul>
|
||||
|
||||
<h3>引擎简介</h3>
|
||||
<table class="compare-table">
|
||||
<tr><th>引擎</th><th>源文件</th><th>语言标准</th><th>算法支持</th></tr>
|
||||
<tr><td>Python</td><td><span class="file-path">compute.py</span></td><td>Python 3.8+</td><td>显式欧拉、隐式欧拉、中点法、蛙跳法</td></tr>
|
||||
<tr><td>C</td><td><span class="file-path">engines/c/main.c</span></td><td>C11</td><td>同上（2026-05-20 新增 3 种）</td></tr>
|
||||
<tr><td>C++</td><td><span class="file-path">engines/cpp/main.cpp</span></td><td>C++17</td><td>同上（2026-05-20 新增 3 种）</td></tr>
|
||||
<tr><td>Fortran</td><td><span class="file-path">engines/fortran/main.f90</span></td><td>Fortran 90</td><td>同上（2026-05-20 新增 3 种）</td></tr>
|
||||
</table>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="py">
|
||||
<h2>2. Python 引擎（标准参考实现）</h2>
|
||||
<p class="file-path" style="margin-bottom:8px">compute.py</p>
|
||||
|
||||
<h3>2.1 蛙跳法积分步骤</h3>
|
||||
<div class="analysis-box">
|
||||
<h4>Velocity-Verlet 三步法</h4>
|
||||
<p>Python 的 <code>Leapfrog_Method</code> 实现了标准 velocity-Verlet 格式。三步法对应于：</p>
|
||||
<ol>
|
||||
<li><strong>速度半步推</strong>（第919-921行）：<em>v</em>(<em>t</em>+Δt/2) = <em>v</em>(<em>t</em>) + ½ · <em>a</em>(<em>t</em>) · Δt</li>
|
||||
<li><strong>全推位置</strong>（第923-925行）：<em>x</em>(<em>t</em>+Δt) = <em>x</em>(<em>t</em>) + <em>v</em>(<em>t</em>+Δt/2) · Δt</li>
|
||||
<li><strong>速度后半步</strong>（第930-937行）：<em>v</em>(<em>t</em>+Δt) = <em>v</em>(<em>t</em>+Δt/2) + ½ · <em>a</em>(<em>t</em>+Δt) · Δt</li>
|
||||
</ol>
|
||||
</div>
|
||||
|
||||
<pre><span class="lang-tag">Python</span><code><span class="ln">917</span>def Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b):
|
||||
<span class="ln">918</span> ax, ay, az = compute_acceleration(x, y, z, vx, vy, vz, m, g, b)
|
||||
<span class="ln">919</span> vx_half = vx + 0.5 * ax * dt
|
||||
<span class="ln">920</span> vy_half = vy + 0.5 * ay * dt
|
||||
<span class="ln">921</span> vz_half = vz + 0.5 * az * dt
|
||||
<span class="ln">922</span>
|
||||
<span class="ln">923</span> x = x + vx_half * dt
|
||||
<span class="ln">924</span> y = y + vy_half * dt
|
||||
<span class="ln">925</span> z = z + vz_half * dt
|
||||
<span class="ln">926</span>
|
||||
<span class="ln">927</span> gamma_x = b[0] / m
|
||||
<span class="ln">928</span> gamma_y = b[1] / m
|
||||
<span class="ln">929</span> gamma_z = b[2] / m
|
||||
<span class="ln">930</span> vx_next = (vx_half + 0.5 * g[0] * dt) / (1.0 + 0.5 * gamma_x * dt)
|
||||
<span class="ln">931</span> vy_next = (vy_half + 0.5 * g[1] * dt) / (1.0 + 0.5 * gamma_y * dt)
|
||||
<span class="ln">932</span> vz_next = (vz_half + 0.5 * g[2] * dt) / (1.0 + 0.5 * gamma_z * dt)
|
||||
<span class="ln">933</span> ax_next, ay_next, az_next = compute_acceleration(
|
||||
<span class="ln">934</span> x, y, z, vx_next, vy_next, vz_next, m, g, b)
|
||||
<span class="ln">935</span> vx = vx_half + 0.5 * ax_next * dt
|
||||
<span class="ln">936</span> vy = vy_half + 0.5 * ay_next * dt
|
||||
<span class="ln">937</span> vz = vz_half + 0.5 * az_next * dt
|
||||
<span class="ln">938</span> return x, y, z, vx, vy, vz
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box warn">
|
||||
<h4>⚠️ 阻尼的隐式处理（第927-932行）</h4>
|
||||
<p>Python 的蛙跳法中，阻尼项用隐式格式处理。<code>gamma = B / m</code>，然后 <code>v_next = (v_half + ½·g·Δt) / (1 + ½·γ·Δt)</code>。这个 <code>v_next</code> 是阻尼修正后的速度估计值，被传给第二次 <code>compute_acceleration</code> 作为入参。当 B = 0 时退化为显式 <code>v_next = v_half + ½·g·Δt</code>。</p>
|
||||
<p>注意：<code>vx_next</code> 是<strong>独立变量</strong>，不覆盖 <code>vx_half</code>。最终速度 <code>vx = vx_half + ½·a_next·dt</code> 使用的是原始的半步速度。这是正确的做法。</p>
|
||||
</div>
|
||||
|
||||
<h3>2.2 边界条件</h3>
|
||||
<pre><span class="lang-tag">Python</span><code><span class="ln">940</span>def Limit_in_box(a, amin, amax, va):
|
||||
<span class="ln">941</span> """限制物体在边界内，发生碰撞时反弹。"""
|
||||
<span class="ln">942</span> over = a > amax
|
||||
<span class="ln">943</span> under = a < amin
|
||||
<span class="ln">944</span> a = np.where(over, amax, a)
|
||||
<span class="ln">945</span> a = np.where(under, amin, a)
|
||||
<span class="ln">946</span> va = np.where(over | under, -va, va)
|
||||
<span class="ln">947</span> return a, va
|
||||
</code></pre>
|
||||
|
||||
<pre><span class="lang-tag">Python</span><code><span class="ln">949</span>def apply_motion_update(x, y, z, vx, vy, vz, dt, m, g, b):
|
||||
<span class="ln">957</span> elif METHOD == "leapfrog":
|
||||
<span class="ln">958</span> x, y, z, vx, vy, vz = Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b)
|
||||
<span class="ln">962</span> x, vx = Limit_in_box(x, X_MIN, X_MAX, vx)
|
||||
<span class="ln">963</span> y, vy = Limit_in_box(y, Y_MIN, Y_MAX, vy)
|
||||
<span class="ln">964</span> z, vz = Limit_in_box(z, Z_MIN, Z_MAX, vz)
|
||||
<span class="ln">965</span> return x, y, z, vx, vy, vz
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box">
|
||||
<h4>边界处理逻辑</h4>
|
||||
<ul>
|
||||
<li>积分方法调用（第957-958行）：只做速度/位置更新，不含边界</li>
|
||||
<li>边界后处理（第962-964行）：三个方向各调用一次 <code>Limit_in_box</code></li>
|
||||
<li><code>Limit_in_box</code>：clamp 位置到边界 + 越界方向速度取反（弹性碰撞模型）</li>
|
||||
</ul>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="forces">
|
||||
<h2>3. 力开关与能量分析系统</h2>
|
||||
<p>2026-05-20 添加了独立的力开关系统，每种力可通过 <span class="file-path">input.txt</span> 中的 0/1 开关独立控制。各引擎的受力计算均基于这些开关条件累加。</p>
|
||||
|
||||
<h3>3.1 力开关配置</h3>
|
||||
<table class="compare-table">
|
||||
<tr><th>开关名</th><th>默认值</th><th>控制的力</th><th>物理公式</th></tr>
|
||||
<tr><td><code>gravity_field</code></td><td>1</td><td>均匀重力场（常数矢量 <em>g</em>）</td><td><em>F</em> = <em>m</em> · <em>g</em></td></tr>
|
||||
<tr><td><code>gravity_interaction</code></td><td>0</td><td>原子间牛顿万有引力</td><td><em>F</em> = <em>G</em> · <em>m</em><sub>i</sub> · <em>m</em><sub>j</sub> / <em>r</em>²</td></tr>
|
||||
<tr><td><code>elastic_force</code></td><td>1</td><td>弹簧键胡克力</td><td><em>F</em> = <em>k</em>(<em>r</em> − <em>r</em><sub>0</sub>)</td></tr>
|
||||
<tr><td><code>damping_force</code></td><td>0</td><td>线性阻尼</td><td><em>F</em> = −<em>B</em> · <em>v</em></td></tr>
|
||||
</table>
|
||||
|
||||
<h3>3.2 代码实现模式</h3>
|
||||
<p>Python <code>compute_force()</code> / C/C++ <code>compute_acceleration()</code> / Fortran <code>accel()</code> 统一采用以下模式：</p>
|
||||
|
||||
<pre><span class="lang-tag">伪代码</span><code>// 1. 清零
|
||||
ax = 0; ay = 0; az = 0
|
||||
|
||||
// 2. 按开关累加
|
||||
if (gravity_field): ax += G[0]; ay += G[1]; az += G[2]
|
||||
if (damping_force): ax -= B[0] * vx / m; ...
|
||||
if (elastic_force): ax += bond_loop(...)
|
||||
if (gravity_interaction): ax += pair_gravity_loop(...)
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box good">
|
||||
<h4>✅ 开关关闭时零开销</h4>
|
||||
<p>每个力都包裹在 <code>if (开关)</code> 条件中。开关关闭时，对应力的计算代码完全跳过，不影响运行效率。</p>
|
||||
</div>
|
||||
|
||||
<h3>3.3 能量图联动</h3>
|
||||
<p>绘图时 <span class="file-path">dynamics.py</span> 的 <code>run_case()</code> 根据力开关状态动态计算对应的势能分量：</p>
|
||||
|
||||
<table class="compare-table">
|
||||
<tr><th>力开关</th><th>势能公式</th><th>关=0 时</th></tr>
|
||||
<tr><td><code>gravity_field=1</code></td><td><em>U<sub>g</sub></em> = −<em>m</em> · <em>G</em> · <em>r</em></td><td>Ug 归零</td></tr>
|
||||
<tr><td><code>gravity_interaction=1</code></td><td><em>U<sub>gg</sub></em> = −∑<sub>i<j</sub> <em>G</em> · <em>m</em><sub>i</sub> · <em>m</em><sub>j</sub> / <em>r</em><sub>ij</sub></td><td>Ugg 归零</td></tr>
|
||||
<tr><td><code>elastic_force=1</code></td><td><em>U<sub>s</sub></em> = ½ ∑ <em>k</em>(<em>r</em> − <em>r</em><sub>0</sub>)²</td><td>Us 归零</td></tr>
|
||||
</table>
|
||||
|
||||
<div class="analysis-box">
|
||||
<h4>能量图线颜色</h4>
|
||||
<p>动能=蓝色、均匀重力势能=绿色、弹性势能=橙色、万有引力势能=<strong style="color:purple;">紫色</strong>、总能量=红色虚线。</p>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="c">
|
||||
<h2>4. C 引擎</h2>
|
||||
<p class="file-path" style="margin-bottom:8px">engines/c/main.c</p>
|
||||
<p style="margin-bottom:10px;font-size:.9em;color:#666;">注意：2026-05-20 重写。原实现为 Symplectic Euler（非 Leapfrog），现已修正为 Velocity-Verlet。</p>
|
||||
|
||||
<h3>3.1 蛙跳法积分步骤</h3>
|
||||
<pre><span class="lang-tag">C</span><code><span class="ln">420</span>/* ── 蛙跳法（Velocity-Verlet）── */
|
||||
<span class="ln">421</span>static void leapfrog_step(
|
||||
<span class="ln">422</span> int n, double *x, double *y, double *z,
|
||||
<span class="ln">423</span> double *vx, double *vy, double *vz,
|
||||
<span class="ln">424</span> const double *m, const double G[3], const double B[3],
|
||||
<span class="ln">425</span> const BondData *bonds, const int *fixed, double dt)
|
||||
<span class="ln">426</span>{
|
||||
<span class="ln">427</span> double *ax = (double*)alloca(n * sizeof(double));
|
||||
<span class="ln">428</span> double *ay = (double*)alloca(n * sizeof(double));
|
||||
<span class="ln">429</span> double *az = (double*)alloca(n * sizeof(double));
|
||||
<span class="ln">430</span> compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az);
|
||||
<span class="ln">431</span>
|
||||
<span class="ln">432</span> /* 半推速度：v_half = v + 0.5*a*dt */
|
||||
<span class="ln">433</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">434</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">435</span> vx[i] += ax[i] * dt * 0.5;
|
||||
<span class="ln">436</span> vy[i] += ay[i] * dt * 0.5;
|
||||
<span class="ln">437</span> vz[i] += az[i] * dt * 0.5;
|
||||
<span class="ln">438</span> }
|
||||
<span class="ln">439</span>
|
||||
<span class="ln">440</span> /* 全推位置（不含边界）*/
|
||||
<span class="ln">441</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">442</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">443</span> x[i] += vx[i] * dt; /* vx 此时是 v_half */
|
||||
<span class="ln">444</span> y[i] += vy[i] * dt;
|
||||
<span class="ln">445</span> z[i] += vz[i] * dt;
|
||||
<span class="ln">446</span> }
|
||||
<span class="ln">447</span>
|
||||
<span class="ln">448</span> /* 隐式阻尼处理：v_next = (v_half + 0.5*g*dt) / (1 + 0.5*gamma*dt)
|
||||
<span class="ln">449</span> 不覆盖 vx/vy/vz，用临时数组存入 */
|
||||
<span class="ln">450</span> double *dmp_vx = (double*)alloca(n * sizeof(double));
|
||||
<span class="ln">451</span> double *dmp_vy = (double*)alloca(n * sizeof(double));
|
||||
<span class="ln">452</span> double *dmp_vz = (double*)alloca(n * sizeof(double));
|
||||
<span class="ln">453</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">454</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">455</span> double gx = B[0] / m[i], gy = B[1] / m[i], gz = B[2] / m[i];
|
||||
<span class="ln">456</span> dmp_vx[i] = (vx[i] + 0.5 * G[0] * dt) / (1.0 + 0.5 * gx * dt);
|
||||
<span class="ln">457</span> dmp_vy[i] = (vy[i] + 0.5 * G[1] * dt) / (1.0 + 0.5 * gy * dt);
|
||||
<span class="ln">458</span> dmp_vz[i] = (vz[i] + 0.5 * G[2] * dt) / (1.0 + 0.5 * gz * dt);
|
||||
<span class="ln">459</span> }
|
||||
<span class="ln">460</span>
|
||||
<span class="ln">461</span> /* 用新位置 + 阻尼处理后的速度重算加速度 */
|
||||
<span class="ln">462</span> compute_acceleration(n, x, y, z, dmp_vx, dmp_vy, dmp_vz, m, G, B, bonds, ax, ay, az);
|
||||
<span class="ln">463</span>
|
||||
<span class="ln">464</span> /* 速度后半步：v = v_half + 0.5*a_next*dt
|
||||
<span class="ln">465</span> vx 仍为 v_half（未被覆盖）*/
|
||||
<span class="ln">466</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">467</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">468</span> vx[i] += ax[i] * dt * 0.5;
|
||||
<span class="ln">469</span> vy[i] += ay[i] * dt * 0.5;
|
||||
<span class="ln">470</span> vz[i] += az[i] * dt * 0.5;
|
||||
<span class="ln">471</span> }
|
||||
<span class="ln">472</span>}
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box good">
|
||||
<h4>✅ Velocity-Verlet 标准实现（与 Python 一致）</h4>
|
||||
<p>C 引擎的 <code>leapfrog_step</code> 与 Python 的 <code>Leapfrog_Method</code> 步骤完全相同：</p>
|
||||
<ol>
|
||||
<li>计算加速度 <em>a</em>(<em>t</em>)（第430行）</li>
|
||||
<li>速度半步推 <em>v</em><sub>½</sub>（第432-438行）</li>
|
||||
<li>全推位置 <em>x</em>(<em>t</em>+Δt)（第440-446行）</li>
|
||||
<li>隐式阻尼 → 临时数组 <code>dmp_vx</code>（第448-459行）</li>
|
||||
<li>用新位置 + 阻尼速度重算加速度（第461-462行）</li>
|
||||
<li>速度后半步（第464-471行）</li>
|
||||
</ol>
|
||||
<p>关键：阻尼值存入 <code>dmp_vx</code> 而不覆盖 <code>vx</code>，与 Python 的 <code>vx_next</code> 独立变量做法一致。</p>
|
||||
</div>
|
||||
|
||||
<h3>3.2 边界条件</h3>
|
||||
<pre><span class="lang-tag">C</span><code><span class="ln">309</span>/* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
|
||||
<span class="ln">310</span>static void limit_in_box(double *pos, double *vel, double lo, double hi) {
|
||||
<span class="ln">311</span> if (*pos > hi) { *pos = hi; *vel = -*vel; }
|
||||
<span class="ln">312</span> if (*pos < lo) { *pos = lo; *vel = -*vel; }
|
||||
<span class="ln">313</span>}
|
||||
</code></pre>
|
||||
|
||||
<pre><span class="lang-tag">C</span><code><span class="ln">473</span>/* ── 分发器：调用对应积分方法 + 边界条件 ── */
|
||||
<span class="ln">488</span> } else if (strcmp(method, "leapfrog") == 0) {
|
||||
<span class="ln">489</span> leapfrog_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
|
||||
<span class="ln">490</span> }
|
||||
<span class="ln">495</span> /* 边界条件（与 Python Limit_in_box 一致） */
|
||||
<span class="ln">496</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">497</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">498</span> limit_in_box(&x[i], &vx[i], -box_a, box_a);
|
||||
<span class="ln">499</span> limit_in_box(&y[i], &vy[i], -box_a, box_a);
|
||||
<span class="ln">500</span> limit_in_box(&z[i], &vz[i], -box_a, box_a);
|
||||
<span class="ln">501</span> }
|
||||
<span class="ln">502</span>}
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box good">
|
||||
<h4>✅ 边界条件已修复</h4>
|
||||
<p>原 C 引擎仅有 <code>clamp</code>（只钳位置不反转速度），现已改为 <code>limit_in_box</code>（钳位置 + 反转速度），与 Python 一致。</p>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="cpp">
|
||||
<h2>5. C++ 引擎</h2>
|
||||
<p class="file-path" style="margin-bottom:8px">engines/cpp/main.cpp</p>
|
||||
|
||||
<h3>4.1 蛙跳法积分步骤</h3>
|
||||
<pre><span class="lang-tag">C++</span><code><span class="ln">360</span>/* ── 蛙跳法（Velocity-Verlet）——与 Python Leapfrog_Method 一致 ── */
|
||||
<span class="ln">361</span>static void leapfrog_full_step(
|
||||
<span class="ln">362</span> int n, double *x, double *y, double *z,
|
||||
<span class="ln">363</span> double *vx, double *vy, double *vz,
|
||||
<span class="ln">364</span> const double *m, const double G[3], const double B[3],
|
||||
<span class="ln">365</span> const BondData &bonds, const int *fixed, double dt)
|
||||
<span class="ln">366</span>{
|
||||
<span class="ln">367</span> // 第一次加速度
|
||||
<span class="ln">368</span> std::vector<double> ax(n), ay(n), az(n);
|
||||
<span class="ln">369</span> compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax.data(), ay.data(), az.data());
|
||||
<span class="ln">370</span>
|
||||
<span class="ln">371</span> // 半推速度：v_half = v + 0.5*a*dt (存入 vx, vy, vz)
|
||||
<span class="ln">372</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">373</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">374</span> vx[i] += ax[i] * dt * 0.5;
|
||||
<span class="ln">375</span> vy[i] += ay[i] * dt * 0.5;
|
||||
<span class="ln">376</span> vz[i] += az[i] * dt * 0.5;
|
||||
<span class="ln">377</span> }
|
||||
<span class="ln">378</span>
|
||||
<span class="ln">379</span> // 全推位置（不含边界，边界在外层统一处理）
|
||||
<span class="ln">380</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">381</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">382</span> x[i] += vx[i] * dt; // vx 此时是 v_half
|
||||
<span class="ln">383</span> y[i] += vy[i] * dt;
|
||||
<span class="ln">384</span> z[i] += vz[i] * dt;
|
||||
<span class="ln">385</span> }
|
||||
<span class="ln">386</span>
|
||||
<span class="ln">387</span> // 隐式阻尼处理得到 v_next（不覆盖 vx/vy/vz，Python 第929-931行 vx_next）
|
||||
<span class="ln">388</span> std::vector<double> dmp_vx(n), dmp_vy(n), dmp_vz(n);
|
||||
<span class="ln">389</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">390</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">391</span> double gamma_x = B[0] / m[i];
|
||||
<span class="ln">392</span> double gamma_y = B[1] / m[i];
|
||||
<span class="ln">393</span> double gamma_z = B[2] / m[i];
|
||||
<span class="ln">394</span> dmp_vx[i] = (vx[i] + 0.5 * G[0] * dt) / (1.0 + 0.5 * gamma_x * dt);
|
||||
<span class="ln">395</span> dmp_vy[i] = (vy[i] + 0.5 * G[1] * dt) / (1.0 + 0.5 * gamma_y * dt);
|
||||
<span class="ln">396</span> dmp_vz[i] = (vz[i] + 0.5 * G[2] * dt) / (1.0 + 0.5 * gamma_z * dt);
|
||||
<span class="ln">397</span> }
|
||||
<span class="ln">398</span>
|
||||
<span class="ln">399</span> // 用新位置 + 阻尼处理后的速度重算加速度
|
||||
<span class="ln">400</span> compute_acceleration(n, x, y, z, dmp_vx.data(), dmp_vy.data(), dmp_vz.data(),
|
||||
<span class="ln">401</span> m, G, B, bonds, ax.data(), ay.data(), az.data());
|
||||
<span class="ln">402</span>
|
||||
<span class="ln">403</span> // 速度后半步：v = v_half + 0.5*a_next*dt
|
||||
<span class="ln">404</span> // vx 仍为 v_half（未被覆盖），直接加上 0.5*a_next*dt
|
||||
<span class="ln">405</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">406</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">407</span> vx[i] += ax[i] * dt * 0.5;
|
||||
<span class="ln">408</span> vy[i] += ay[i] * dt * 0.5;
|
||||
<span class="ln">409</span> vz[i] += az[i] * dt * 0.5;
|
||||
<span class="ln">410</span> }
|
||||
<span class="ln">411</span>}
|
||||
</code></pre>
|
||||
|
||||
<h3>4.2 边界条件</h3>
|
||||
<pre><span class="lang-tag">C++</span><code><span class="ln">265</span>/* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
|
||||
<span class="ln">266</span>static void limit_in_box(double &pos, double &vel, double lo, double hi) {
|
||||
<span class="ln">267</span> if (pos > hi) { pos = hi; vel = -vel; }
|
||||
<span class="ln">268</span> if (pos < lo) { pos = lo; vel = -vel; }
|
||||
<span class="ln">269</span>}
|
||||
</code></pre>
|
||||
|
||||
<pre><span class="lang-tag">C++</span><code><span class="ln">413</span>/* ── 分发器：调用对应积分方法 + 边界条件 ── */
|
||||
<span class="ln">429</span> } else if (method == "leapfrog") {
|
||||
<span class="ln">430</span> leapfrog_full_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
|
||||
<span class="ln">431</span> }
|
||||
<span class="ln">435</span> // 边界条件（与 Python Limit_in_box 一致）
|
||||
<span class="ln">436</span> for (int i = 0; i < n; i++) {
|
||||
<span class="ln">437</span> if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
<span class="ln">438</span> limit_in_box(x[i], vx[i], -box_a, box_a);
|
||||
<span class="ln">439</span> limit_in_box(y[i], vy[i], -box_a, box_a);
|
||||
<span class="ln">440</span> limit_in_box(z[i], vz[i], -box_a, box_a);
|
||||
<span class="ln">441</span> }
|
||||
<span class="ln">442</span>}
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box good">
|
||||
<h4>✅ 与 Python 算法一致</h4>
|
||||
<p>C++ 引擎的 <code>leapfrog_full_step</code> 在重写后已与 Python 完全对齐，阻尼使用 <code>std::vector<double> dmp_vx</code> 临时数组而非原处覆盖。边界条件新增了速度反转。</p>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="fortran">
|
||||
<h2>6. Fortran 引擎</h2>
|
||||
<p class="file-path" style="margin-bottom:8px">engines/fortran/main.f90</p>
|
||||
|
||||
<h3>5.1 蛙跳法积分步骤</h3>
|
||||
<pre><span class="lang-tag">Fortran</span><code><span class="ln">498</span>! ── 蛙跳法（Velocity-Verlet）──
|
||||
<span class="ln">499</span>subroutine leapfrog_full(n, x, y, z, vx, vy, vz, m, G, B, &
|
||||
<span class="ln">500</span> nb, bp, bk, br, fixed, dt)
|
||||
<span class="ln">501</span> integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3)
|
||||
<span class="ln">502</span> double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
|
||||
<span class="ln">503</span> double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), dt
|
||||
<span class="ln">504</span> double precision :: ax(n), ay(n), az(n)
|
||||
<span class="ln">505</span> double precision :: dmp_vx(n), dmp_vy(n), dmp_vz(n)
|
||||
<span class="ln">506</span> double precision :: gx, gy, gz
|
||||
<span class="ln">507</span> integer :: i
|
||||
<span class="ln">508</span>
|
||||
<span class="ln">509</span> call accel(n, x, y, z, vx, vy, vz, m, G, B, nb, bp, bk, br, ax, ay, az)
|
||||
<span class="ln">510</span>
|
||||
<span class="ln">511</span> ! 速度半步推
|
||||
<span class="ln">512</span> do i = 1, n
|
||||
<span class="ln">513</span> if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
<span class="ln">514</span> vx(i) = vx(i) + ax(i) * dt * 0.5d0
|
||||
<span class="ln">515</span> vy(i) = vy(i) + ay(i) * dt * 0.5d0
|
||||
<span class="ln">516</span> vz(i) = vz(i) + az(i) * dt * 0.5d0
|
||||
<span class="ln">517</span> end do
|
||||
<span class="ln">518</span>
|
||||
<span class="ln">519</span> ! 全推位置（不含边界）
|
||||
<span class="ln">520</span> do i = 1, n
|
||||
<span class="ln">521</span> if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
<span class="ln">522</span> x(i) = x(i) + vx(i) * dt
|
||||
<span class="ln">523</span> y(i) = y(i) + vy(i) * dt
|
||||
<span class="ln">524</span> z(i) = z(i) + vz(i) * dt
|
||||
<span class="ln">525</span> end do
|
||||
<span class="ln">526</span>
|
||||
<span class="ln">527</span> ! 隐式阻尼处理（用临时数组 dmp_v，不覆盖 vx/vy/vz）
|
||||
<span class="ln">528</span> do i = 1, n
|
||||
<span class="ln">529</span> if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then
|
||||
<span class="ln">530</span> dmp_vx(i) = 0; dmp_vy(i) = 0; dmp_vz(i) = 0; cycle
|
||||
<span class="ln">531</span> end if
|
||||
<span class="ln">532</span> gx = B(1) / m(i); gy = B(2) / m(i); gz = B(3) / m(i)
|
||||
<span class="ln">533</span> dmp_vx(i) = (vx(i) + 0.5d0 * G(1) * dt) / (1.0d0 + 0.5d0 * gx * dt)
|
||||
<span class="ln">534</span> dmp_vy(i) = (vy(i) + 0.5d0 * G(2) * dt) / (1.0d0 + 0.5d0 * gy * dt)
|
||||
<span class="ln">535</span> dmp_vz(i) = (vz(i) + 0.5d0 * G(3) * dt) / (1.0d0 + 0.5d0 * gz * dt)
|
||||
<span class="ln">536</span> end do
|
||||
<span class="ln">537</span>
|
||||
<span class="ln">538</span> ! 用新位置 + 阻尼速度重算加速度
|
||||
<span class="ln">539</span> call accel(n, x, y, z, dmp_vx, dmp_vy, dmp_vz, m, G, B, nb, bp, bk, br, ax, ay, az)
|
||||
<span class="ln">540</span>
|
||||
<span class="ln">541</span> ! 速度后半步：v = v_half + 0.5*a_next*dt（vx 仍为 v_half）
|
||||
<span class="ln">542</span> do i = 1, n
|
||||
<span class="ln">543</span> if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
<span class="ln">544</span> vx(i) = vx(i) + ax(i) * dt * 0.5d0
|
||||
<span class="ln">545</span> vy(i) = vy(i) + ay(i) * dt * 0.5d0
|
||||
<span class="ln">546</span> vz(i) = vz(i) + az(i) * dt * 0.5d0
|
||||
<span class="ln">547</span> end do
|
||||
<span class="ln">548</span>end subroutine leapfrog_full
|
||||
</code></pre>
|
||||
|
||||
<h3>5.2 边界条件</h3>
|
||||
<pre><span class="lang-tag">Fortran</span><code><span class="ln">395</span>! 边界条件：clamp 位置 + 速度反转
|
||||
<span class="ln">396</span>subroutine limit_in_box(pos, vel, lo, hi)
|
||||
<span class="ln">397</span> double precision, intent(inout) :: pos, vel
|
||||
<span class="ln">398</span> double precision, intent(in) :: lo, hi
|
||||
<span class="ln">399</span> if (pos > hi) then
|
||||
<span class="ln">400</span> pos = hi; vel = -vel
|
||||
<span class="ln">401</span> else if (pos < lo) then
|
||||
<span class="ln">402</span> pos = lo; vel = -vel
|
||||
<span class="ln">403</span> end if
|
||||
<span class="ln">404</span>end subroutine limit_in_box
|
||||
</code></pre>
|
||||
|
||||
<pre><span class="lang-tag">Fortran</span><code><span class="ln">571</span> else if (trim(method) == 'leapfrog') then
|
||||
<span class="ln">572</span> call leapfrog_full(n, x, y, z, vx, vy, vz, mm, G, B, &
|
||||
<span class="ln">573</span> nb, bp, bk, br, fixed, dt)
|
||||
<span class="ln">574</span> end if
|
||||
<span class="ln">579</span> ! 边界条件
|
||||
<span class="ln">580</span> do i = 1, n
|
||||
<span class="ln">581</span> if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
<span class="ln">582</span> call limit_in_box(x(i), vx(i), -box_a, box_a)
|
||||
<span class="ln">583</span> call limit_in_box(y(i), vy(i), -box_a, box_a)
|
||||
<span class="ln">584</span> call limit_in_box(z(i), vz(i), -box_a, box_a)
|
||||
<span class="ln">585</span> end do
|
||||
<span class="ln">586</span>end subroutine apply_step
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box good">
|
||||
<h4>✅ 与 Python 算法一致</h4>
|
||||
<p>Fortran 引擎的 <code>leapfrog_full</code> 在重写中修复了阻尼覆盖问题（使用 <code>dmp_vx</code> 临时数组），边界也改为速度反转。</p>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="diff">
|
||||
<h2>7. 版本变更记录</h2>
|
||||
<p>以下列出了 2026-05-20 三次迭代中所有修复和改进：</p>
|
||||
|
||||
<table class="compare-table">
|
||||
<thead>
|
||||
<tr><th>#</th><th>问题描述</th><th>影响引擎</th><th>修复方式</th></tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td>1</td>
|
||||
<td><strong>C 引擎蛙跳算法错误</strong><br>原 <code>leapfrog_step</code> 实现的是 Symplectic Euler：<br><code>v += a·Δt; x += v·Δt</code>（一次加速度，使用新速度推位置）。<br>这比 Velocity-Verlet 的截断误差大一阶（<em>O</em>(Δt) vs <em>O</em>(Δt²)），且只用了一次加速度计算。</td>
|
||||
<td><span class="tag bad">C</span></td>
|
||||
<td>重写为 velocity-Verlet 三步法：<br>v<sub>½</sub> → x → a<sub>new</sub> → v</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>2</td>
|
||||
<td><strong>阻尼值覆盖速度变量</strong><br>蛙跳法的隐式阻尼处理直接写回 <code>vx[i]</code>，导致后面的速度后半步 <code>vx[i] += ½·a·dt</code> 在覆盖后的值上叠加。等效于重力加速度被多算了一次，引入 ½·g·Δt 的系统偏差（Fortran）或 g·Δt 的偏差（C）。</td>
|
||||
<td><span class="tag bad">C, Fortran</span></td>
|
||||
<td>用临时数组 <code>dmp_vx</code> 存储阻尼估计值，不覆盖原始速度变量</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>3</td>
|
||||
<td><strong><span class="file-path">bond.txt</span> 表头未跳过</strong><br>C/C++ 引擎的 <code>read_bonds</code> 用 <code>while (fb >> bn >> bk >> br)</code> 读取 <span class="file-path">bond.txt</span>，但第一行是表头 <code>bond_name k rest_length</code>。解析 <code>k</code> 时尝试将字符串 <code>"k"</code> 转 double 失败，使用默认值 1.0（应为 100.0），弹簧劲度系数差 100 倍，完全改变了轨迹。Python 和 Fortran 正确跳过了表头。</td>
|
||||
<td><span class="tag bad">C, C++</span></td>
|
||||
<td>在数据读取循环前用 <code>std::getline(fb, header)</code>（C++）或 <code>fgets(header, ...)</code>（C）跳过表头行</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>4</td>
|
||||
<td><strong>C++ JSON 输出精度不足</strong><br>C++ 的 JSON 输出用默认 <code><<</code> 运算符，只有 6 位有效数字，而 C 的 <code>%.15g</code> 和 Fortran 的 <code>g0</code> 都用全精度。导致 C++ 与其他引擎的输出到第 6 位之后截断。</td>
|
||||
<td><span class="tag warn">C++</span></td>
|
||||
<td>添加 <code>std::setprecision(15)</code> 到输出流</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>5</td>
|
||||
<td><strong>C 引擎 <code>record_steps</code> 忽略 <code>warmup_steps</code></strong><br>C 引擎的轨迹缓冲区大小写死为 <code>params.NT</code> 而非 <code>params.NT - params.warmup_steps</code>。Python/C++/Fortran 正确。</td>
|
||||
<td><span class="tag warn">C</span></td>
|
||||
<td>改为 <code>int record_steps = params.NT - params.warmup_steps;</code></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>6</td>
|
||||
<td><strong>力开关系统</strong><br>新增 4 个独立力开关 <code>gravity_field</code>/<code>gravity_interaction</code>/<code>elastic_force</code>/<code>damping_force</code>，通过 <span class="file-path">input.txt</span> 的 0/1 控制各力是否参与计算。所有引擎统一实现：清零 → 按条件累加。</td>
|
||||
<td><span class="tag good">四引擎</span></td>
|
||||
<td>Python: 重写 <code>compute_force()</code>；C/C++: 重写 <code>compute_acceleration()</code>；Fortran: 重写 <code>accel()</code></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>7</td>
|
||||
<td><strong>万有引力引擎</strong><br>原子间牛顿万有引力 <em>F = G · m<sub>i</sub> · m<sub>j</sub> / r²</em>，O(N²) 双循环对计算，<code>gravity_interaction=0</code> 时零开销。</td>
|
||||
<td><span class="tag good">四引擎</span></td>
|
||||
<td>所有引擎的受力计算函数中添加多原子对循环</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>8</td>
|
||||
<td><strong>逐自由度固定约束</strong><br><span class="file-path">coord.txt</span> 末尾增加 <code>fix_x fix_y fix_z</code> 三列，每列 0/1 控制该方向是否锁定。锁定后位置回复初始值、速度置零。</td>
|
||||
<td><span class="tag good">四引擎</span></td>
|
||||
<td>C/C++/Fortran: <code>apply_step()</code> 新增 <code>pos_0</code> 参数，边界条件后添加约束；Python: 已有 <code>apply_fixed_constraints()</code> 完美支持</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>9</td>
|
||||
<td><strong>能量图随力开关联动</strong><br>绘图时根据力开关状态动态计算势能分量。仅 <code>gravity_interaction=1</code> 时计算万有引力势能并以<strong style="color:purple;">紫色</strong>曲线绘制。</td>
|
||||
<td><span class="tag good">Python</span></td>
|
||||
<td><span class="file-path">dynamics.py</span> 能量计算添加 <code>ug_grav</code> 和开关判断</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="bc">
|
||||
<h2>8. 边界条件深度分析</h2>
|
||||
|
||||
<h3>8.1 硬壁反射实现</h3>
|
||||
<p>四引擎的边界条件在重写后完全一致：</p>
|
||||
<table class="compare-table">
|
||||
<tr><th>引擎</th><th>函数名</th><th>钳位置</th><th>反转速度</th><th>调用时机</th></tr>
|
||||
<tr><td>Python</td><td><code>Limit_in_box</code></td><td><span class="tag good">✅</span></td><td><span class="tag good">✅</span></td><td>蛙跳步之后（后处理）</td></tr>
|
||||
<tr><td>C</td><td><code>limit_in_box</code></td><td><span class="tag good">✅</span></td><td><span class="tag good">✅</span></td><td>积分步之后（<code>apply_step</code>）</td></tr>
|
||||
<tr><td>C++</td><td><code>limit_in_box</code></td><td><span class="tag good">✅</span></td><td><span class="tag good">✅</span></td><td>积分步之后（<code>apply_step</code>）</td></tr>
|
||||
<tr><td>Fortran</td><td><code>limit_in_box</code></td><td><span class="tag good">✅</span></td><td><span class="tag good">✅</span></td><td>积分步之后（<code>apply_step</code>）</td></tr>
|
||||
</table>
|
||||
|
||||
<div class="analysis-box">
|
||||
<h4>算法合理性</h4>
|
||||
<p><strong>速度反转是必须的。</strong> 硬壁反射的物理模型要求速度在碰撞时反向以保持运动方向的一致性。没有速度反转时，粒子到达边界后会被持续钳制在边界上（"黏壁"现象），这是功能性缺陷。</p>
|
||||
</div>
|
||||
|
||||
<h3>8.2 逐自由度固定约束</h3>
|
||||
<p>2026-05-20 新增了逐自由度固定约束。在 <span class="file-path">coord.txt</span> 末尾添加 <code>fix_x fix_y fix_z</code> 三列：</p>
|
||||
|
||||
<pre><span class="lang-tag">coord.txt 格式</span><code>n mass radius x y z vx vy vz fix_x fix_y fix_z
|
||||
1 1 0.28 -1 0 0 0 0 0 1 0 0 ← 仅固定 x 方向
|
||||
2 1 0.28 1 0 1 0 0 0 0 1 1 ← 固定 y 和 z 方向
|
||||
</code></pre>
|
||||
|
||||
<p>Python 的 <code>apply_fixed_constraints()</code> 使用 NumPy 逐元素操作：</p>
|
||||
<pre><span class="lang-tag">Python</span><code>fixed = ATOM_FIXED != 0 # (n_atoms, 3) 布尔矩阵
|
||||
positions = np.where(fixed, ATOM_POSITIONS, positions) # 逐度约束
|
||||
velocities = np.where(fixed, 0.0, velocities)
|
||||
</code></pre>
|
||||
|
||||
<p>C/C++/Fortran 引擎在 <code>apply_step()</code> 的边界条件之后添加了相同的约束逻辑：</p>
|
||||
<pre><span class="lang-tag">C/C++/Fortran 伪代码</span><code>// 边界条件后，逐自由度固定约束
|
||||
for i = 1..n:
|
||||
if fixed_x[i]: x[i] = x0[i]; vx[i] = 0
|
||||
if fixed_y[i]: y[i] = y0[i]; vy[i] = 0
|
||||
if fixed_z[i]: z[i] = z0[i]; vz[i] = 0
|
||||
</code></pre>
|
||||
|
||||
<div class="analysis-box good">
|
||||
<h4>✅ 约束在积分步与边界条件之后，保证正确性</h4>
|
||||
<p>约束在积分器完成位置/速度更新、边界钳位之后执行，确保固定自由度不参与任何物理演化。</p>
|
||||
</div>
|
||||
|
||||
<h3>8.3 能量局限性</h3>
|
||||
<div class="analysis-box warn">
|
||||
<h4>⚠️ 硬壁反射与 Symplectic 积分器的不兼容</h4>
|
||||
<p>Velocity-Verlet 是 symplectic 积分器，在平滑势场中能量几乎精确守恒（前 20 步漂移仅 -0.000078%）。但硬壁反射是一个非光滑势（无限高势垒），每个碰撞步的工作流程如下：</p>
|
||||
<ol>
|
||||
<li>蛙跳步假定粒子在光滑势场中自由运动了 Δt 时间</li>
|
||||
<li>后处理发现粒子越界，将位置钳回边界 + 反转速度</li>
|
||||
<li>这个"钳回"动作改变了重力势能（因为 <em>z</em> 位置被改变）但没有相应调整动能</li>
|
||||
<li>每次碰撞产生约 0.03~0.5 J 的能量变化</li>
|
||||
</ol>
|
||||
<p>在典型 case01 参数下（100,000 步，83 次碰底），累积能量漂移约 36%。这是硬壁边界的固有限制，不是代码 Bug。</p>
|
||||
</div>
|
||||
|
||||
<div class="analysis-box">
|
||||
<h4>✅ 修复建议：软壁势</h4>
|
||||
<p>用光滑的排斥势替换硬壁 clamp 可保持 symplectic 守恒性：</p>
|
||||
<ul>
|
||||
<li><strong>WCA 势</strong>（Weeks-Chandler-Andersen）：截断平移的 Lennard-Jones</li>
|
||||
<li><strong>指数排斥</strong>：<em>U</em>(<em>z</em>) = <em>A</em> · exp(−|z|/<em>z</em><sub>0</sub>)</li>
|
||||
<li><strong>半壁谐振子</strong>：|z| > <em>a</em> 时受到线性回复力 <em>F</em> = −<em>k</em>(z − sign(z)·<em>a</em>)</li>
|
||||
</ul>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
<!-- ================================================================ -->
|
||||
<section id="summary">
|
||||
<h2>9. 总结</h2>
|
||||
|
||||
<table class="compare-table">
|
||||
<thead>
|
||||
<tr><th>维度</th><th>Python</th><th>C</th><th>C++</th><th>Fortran</th></tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><strong>积分器类型</strong></td>
|
||||
<td><span class="tag good">Velocity-Verlet</span></td>
|
||||
<td><span class="tag good">Velocity-Verlet</span><br><small>原为 Symplectic Euler</small></td>
|
||||
<td><span class="tag good">Velocity-Verlet</span></td>
|
||||
<td><span class="tag good">Velocity-Verlet</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>算法数量</strong></td>
|
||||
<td><span class="tag good">4 种</span></td>
|
||||
<td><span class="tag good">4 种</span><br><small>新增 3 种</small></td>
|
||||
<td><span class="tag good">4 种</span><br><small>新增 3 种</small></td>
|
||||
<td><span class="tag good">4 种</span><br><small>新增 3 种</small></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>阻尼处理</strong></td>
|
||||
<td>隐式 + 显式混合</td>
|
||||
<td>同 Python</td>
|
||||
<td>同 Python</td>
|
||||
<td>同 Python</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>边界：钳位置</strong></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>边界：反转速度</strong></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span><br><small>原缺失</small></td>
|
||||
<td><span class="tag good">✅</span><br><small>原缺失</small></td>
|
||||
<td><span class="tag good">✅</span><br><small>原缺失</small></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>输出精度（JSON）</strong></td>
|
||||
<td><span class="tag good">双精度</span></td>
|
||||
<td><span class="tag good">%.15g</span></td>
|
||||
<td><span class="tag good">setprecision(15)</span><br><small>原默认 6 位</small></td>
|
||||
<td><span class="tag good">g0</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>算法选择方式</strong></td>
|
||||
<td>YAML config</td>
|
||||
<td>param.json</td>
|
||||
<td>param.json</td>
|
||||
<td>param.json</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>引擎一致性</strong></td>
|
||||
<td colspan="4" style="text-align:center;"><span class="tag good">四个引擎输出完全一致（双精度验证通过）</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>力开关系统</strong></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>万有引力</strong></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
<td><span class="tag good">✅</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><strong>逐自由度约束</strong></td>
|
||||
<td><span class="tag good">✅</span><br><small>np.where</small></td>
|
||||
<td><span class="tag good">✅</span><br><small>pos_0 参数</small></td>
|
||||
<td><span class="tag good">✅</span><br><small>pos_0 参数</small></td>
|
||||
<td><span class="tag good">✅</span><br><small>pos_0 参数</small></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
</section>
|
||||
|
||||
</div>
|
||||
|
||||
<footer>
|
||||
Dynamics Project — 多语言引擎代码分析 • 更新于 2026-05-20（力开关 · 万有引力 · 逐自由度约束 · 能量图联动）
|
||||
</footer>
|
||||
|
||||
</body>
|
||||
</html>
|
||||
+95
-36
@@ -4,9 +4,9 @@ dynamics.py
|
||||
统一入口：读取 YAML 配置文件 → 运行模拟 → 抽帧 → 绘图（可选）
|
||||
|
||||
用法：
|
||||
python dynamics.py # 使用 input/parameters.yaml
|
||||
python dynamics.py input/parameters.yaml # 指定配置文件
|
||||
python dynamics.py --config input/parameters.yaml --no-plot
|
||||
python dynamics.py # 使用 input/input.txt
|
||||
python dynamics.py input/input.txt # 指定配置文件（YAML 格式）
|
||||
python dynamics.py --config input/input.txt --no-plot
|
||||
"""
|
||||
|
||||
import os
|
||||
@@ -115,6 +115,14 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
config = load_yaml_config(config_path)
|
||||
print(f"[run] 已加载配置: {config_path}")
|
||||
|
||||
# ── T_total → NT 自动转换 ──────────────────────
|
||||
if "T_total" in config and "NT" not in config:
|
||||
dt = float(config["DT"])
|
||||
config["NT"] = int(float(config["T_total"]) / dt)
|
||||
print(f"[run] T_total={config['T_total']} → NT={config['NT']} (DT={dt})")
|
||||
elif "T_total" in config and "NT" in config:
|
||||
print(f"[run] 同时指定了 T_total 和 NT，使用 NT={config['NT']}")
|
||||
|
||||
# 显示步骤控制信息
|
||||
steps_info = {k: config.get(k, 1) for k in ["step_simulate", "step_sample", "step_plot", "step_animation"]}
|
||||
step_flags = ", ".join(f"{k}={v}" for k, v in steps_info.items())
|
||||
@@ -136,27 +144,42 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
disp_path = os.path.join(output_dir_abs, "display.txt")
|
||||
|
||||
# ── 自动缓存检测 ───────────────────────────────────────
|
||||
# 若 output/ 中 trajectory.txt 和 display.txt 均已存在，
|
||||
# 自动跳过模拟和抽帧，直接使用已有结果。
|
||||
# 如需强制重新计算，删除 output/ 目录或设 step_simulate: 1 即可。
|
||||
output_exists = (
|
||||
os.path.isdir(output_dir_abs)
|
||||
and os.path.exists(traj_path)
|
||||
and os.path.exists(disp_path)
|
||||
)
|
||||
if output_exists:
|
||||
if config.get("step_simulate", 1):
|
||||
print(f"[run] 检测到已有输出（{traj_path}），自动跳过模拟与抽帧，直接进入后续步骤")
|
||||
config["step_simulate"] = 0
|
||||
config["step_sample"] = 0
|
||||
else:
|
||||
print(f"[run] 已有输出，步骤已被跳过")
|
||||
# force_calc=1: 强制重新计算，忽略缓存
|
||||
# force_calc=0: 尊重 step_simulate 设置，不自动覆盖
|
||||
force_calc = int(config.get("force_calc", 0))
|
||||
if force_calc:
|
||||
print(f"[run] force_calc=1，跳过缓存，强制重新计算")
|
||||
config["step_simulate"] = 1
|
||||
config["step_sample"] = 1
|
||||
elif config.get("step_simulate", 1):
|
||||
# step_simulate=1 且 force_calc=0 → 按用户要求执行计算
|
||||
# 但检测一下参数是否已变更（NT），如果变了则自动更新 step_sample
|
||||
if os.path.isdir(output_dir_abs) and os.path.exists(traj_path) and os.path.exists(disp_path):
|
||||
cached_nt = None
|
||||
try:
|
||||
with open(traj_path, 'rb') as _f:
|
||||
_f.seek(-4096, 2)
|
||||
_tail = _f.read().decode('utf-8', errors='replace')
|
||||
import re as _re
|
||||
_m = _re.search(r'"NT":\s*(\d+)', _tail)
|
||||
if _m:
|
||||
cached_nt = int(_m.group(1))
|
||||
except Exception:
|
||||
pass
|
||||
config_nt = int(config.get("NT", 0))
|
||||
if cached_nt is not None and cached_nt != config_nt:
|
||||
print(f"[run] 参数已变更（缓存 NT={cached_nt}，配置 NT={config_nt}），"
|
||||
f"将重新计算")
|
||||
config["step_sample"] = 1
|
||||
else:
|
||||
# 参数一致，按 step_simulate=1 执行，step_sample 由用户设置决定
|
||||
pass
|
||||
else:
|
||||
# 目录存在但文件不全 → 强制重新计算
|
||||
if os.path.isdir(output_dir_abs):
|
||||
print(f"[run] output/ 目录存在但文件不完整，将重新计算")
|
||||
# step_simulate=0 → 检测缓存是否存在并提示
|
||||
if os.path.isdir(output_dir_abs) and os.path.exists(traj_path) and os.path.exists(disp_path):
|
||||
print(f"[run] 已有输出，步骤已被跳过")
|
||||
else:
|
||||
print(f"[run] output/ 目录不存在，将执行完整流程")
|
||||
print(f"[run] 没有可用的缓存输出，但 step_simulate=0，将跳过模拟")
|
||||
|
||||
# 2. 运行物理模拟 → output/trajectory.txt
|
||||
if config.get("step_simulate", 1):
|
||||
@@ -164,9 +187,12 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
total_steps = config["NT"]
|
||||
record_steps = total_steps - (config.get("warmup_steps") or 0)
|
||||
print(f"[run] 开始计算 总步数={total_steps} 记录步数={record_steps} DT={config['DT']}")
|
||||
|
||||
import time as _time
|
||||
_t0 = _time.time()
|
||||
|
||||
if engine == "python":
|
||||
traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = compute.run_from_config(config, str(runtime_base))
|
||||
print(f"[run] 计算完成，记录 {record_steps} 步")
|
||||
compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base))
|
||||
else:
|
||||
# 外部引擎：先加载配置到全局变量，再运行引擎，再用 save_trajectory_txt 补全 metadata
|
||||
@@ -178,6 +204,9 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = compute.run_engine(
|
||||
engine, input_dir_abs, output_dir_abs, config)
|
||||
compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base))
|
||||
|
||||
_elapsed = _time.time() - _t0
|
||||
print(f"[run] 引擎: {engine} 计算完成: {record_steps} 步 {_elapsed:.3f} s")
|
||||
else:
|
||||
print("[run] 步骤 [模拟] 已跳过，直接加载已有轨迹")
|
||||
if not os.path.exists(traj_path):
|
||||
@@ -282,6 +311,11 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
"box_color_r": float(data["box_color_r"]),
|
||||
"box_color_g": float(data["box_color_g"]),
|
||||
"box_color_b": float(data["box_color_b"]),
|
||||
"gravity_field": int(data.get("gravity_field", 1)),
|
||||
"gravity_interaction": int(data.get("gravity_interaction", 0)),
|
||||
"elastic_force": int(data.get("elastic_force", 1)),
|
||||
"damping_force": int(data.get("damping_force", 0)),
|
||||
"gravity_strength": float(data.get("gravity_strength", 1.0)),
|
||||
}
|
||||
save_display_txt(disp_data, str(runtime_base))
|
||||
print(f"[run] 抽帧完成: {sample_end - sample_start} 步 -> {n_frames} 帧")
|
||||
@@ -328,21 +362,28 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
# ── 能量计算 ─────────────────────────────────────
|
||||
masses = np.array(data["atom_masses"]) # (n_atoms,)
|
||||
G_vec = np.array(data.get("G", [0.0, 0.0, -9.8])) # [gx, gy, gz]
|
||||
gravity_field_enabled = int(data.get("gravity_field", 1))
|
||||
gravity_interaction_enabled = int(data.get("gravity_interaction", 0))
|
||||
gravity_strength = float(data.get("gravity_strength", 1.0))
|
||||
elastic_force_enabled = int(data.get("elastic_force", 1))
|
||||
damping_force_enabled = int(data.get("damping_force", 0))
|
||||
|
||||
# 动能 Ek = ½ m v²
|
||||
ek = 0.5 * masses[np.newaxis, :] * (all_vx**2 + all_vy**2 + all_vz**2)
|
||||
|
||||
# 重力势能 Ug = -m G·r
|
||||
ug = -masses[np.newaxis, :] * (
|
||||
G_vec[0] * all_x + G_vec[1] * all_y + G_vec[2] * all_z
|
||||
)
|
||||
# 均匀重力场势能 Ug = -m G·r
|
||||
ug = np.zeros_like(ek)
|
||||
if gravity_field_enabled:
|
||||
ug = -masses[np.newaxis, :] * (
|
||||
G_vec[0] * all_x + G_vec[1] * all_y + G_vec[2] * all_z
|
||||
)
|
||||
|
||||
# 弹性势能 Us = ½ k (d - d₀)²
|
||||
us = np.zeros_like(ek)
|
||||
bond_pairs = data.get("bond_pairs")
|
||||
bond_stiffness = data.get("bond_stiffness")
|
||||
bond_rest_lengths = data.get("bond_rest_lengths")
|
||||
if bond_pairs is not None and len(bond_pairs) > 0:
|
||||
if elastic_force_enabled and bond_pairs is not None and len(bond_pairs) > 0:
|
||||
for b_idx in range(len(bond_pairs)):
|
||||
i, j = bond_pairs[b_idx]
|
||||
dx = all_x[:, j] - all_x[:, i]
|
||||
@@ -350,17 +391,33 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
dz = all_z[:, j] - all_z[:, i]
|
||||
dist = np.sqrt(dx**2 + dy**2 + dz**2)
|
||||
stretch = dist - bond_rest_lengths[b_idx]
|
||||
us[:, i] += 0.5 * bond_stiffness[b_idx] * stretch**2
|
||||
us[:, j] += 0.5 * bond_stiffness[b_idx] * stretch**2
|
||||
us_each = 0.5 * bond_stiffness[b_idx] * stretch**2
|
||||
us[:, i] += us_each # 将整根键的势能记给 i
|
||||
|
||||
# 万有引力势能 Ug_grav = -G_grav * m_i * m_j / r
|
||||
ug_grav = np.zeros_like(ek)
|
||||
if gravity_interaction_enabled:
|
||||
n_atoms_en = len(masses)
|
||||
for i in range(n_atoms_en):
|
||||
for j in range(i + 1, n_atoms_en):
|
||||
dx = all_x[:, j] - all_x[:, i]
|
||||
dy = all_y[:, j] - all_y[:, i]
|
||||
dz = all_z[:, j] - all_z[:, i]
|
||||
dist = np.sqrt(dx**2 + dy**2 + dz**2)
|
||||
dist = np.maximum(dist, 1e-12)
|
||||
pair_pe = -gravity_strength * masses[i] * masses[j] / dist
|
||||
ug_grav[:, i] += 0.5 * pair_pe
|
||||
ug_grav[:, j] += 0.5 * pair_pe
|
||||
|
||||
# 各原子总能量
|
||||
e_total = ek + ug + us # (NT, n_atoms)
|
||||
e_total = ek + ug + us + ug_grav # (NT, n_atoms)
|
||||
|
||||
# 系统能量分量
|
||||
ek_sys = np.sum(ek, axis=1)
|
||||
ug_sys = np.sum(ug, axis=1)
|
||||
us_sys = np.sum(us, axis=1)
|
||||
e_sys = ek_sys + ug_sys + us_sys
|
||||
ug_grav_sys = np.sum(ug_grav, axis=1)
|
||||
e_sys = ek_sys + ug_sys + us_sys + ug_grav_sys
|
||||
|
||||
# ── 第 3 行左：各原子总能量 ──
|
||||
ax_e = axes[2, 0]
|
||||
@@ -376,9 +433,11 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
# ── 第 3 行右：系统总能量 ──
|
||||
ax_sys = axes[2, 1]
|
||||
ax_sys.plot(time, ek_sys, 'b-', linewidth=1.5, label="系统动能")
|
||||
ax_sys.plot(time, ug_sys, 'g-', linewidth=1.5, label="系统重力势能")
|
||||
if bond_pairs is not None and len(bond_pairs) > 0:
|
||||
ax_sys.plot(time, ug_sys, 'g-', linewidth=1.5, label="均匀重力势能")
|
||||
if elastic_force_enabled and bond_pairs is not None and len(bond_pairs) > 0:
|
||||
ax_sys.plot(time, us_sys, color='orange', linewidth=1.5, label="系统弹性势能")
|
||||
if gravity_interaction_enabled:
|
||||
ax_sys.plot(time, ug_grav_sys, color='purple', linewidth=1.5, label="万有引力势能")
|
||||
ax_sys.plot(time, e_sys, 'r--', linewidth=1.5, label="系统总能量")
|
||||
ax_sys.set_title("系统总能量")
|
||||
ax_sys.set_xlabel("时间 (s)")
|
||||
@@ -421,8 +480,8 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
|
||||
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(description="物理模拟统一入口")
|
||||
parser.add_argument("config_file", nargs="?", default=os.path.join("input", "parameters.yaml"),
|
||||
help="YAML 配置文件路径（默认: input/parameters.yaml）")
|
||||
parser.add_argument("config_file", nargs="?", default=os.path.join("input", "input.txt"),
|
||||
help="YAML 配置文件路径（默认: input/input.txt，虽然是 .txt 后缀但使用 YAML 格式）")
|
||||
parser.add_argument("--config", dest="config_override",
|
||||
help="YAML 配置文件路径（可选，优先于位置参数）")
|
||||
parser.add_argument("--input-dir", default="input",
|
||||
|
||||
+322
-84
@@ -2,13 +2,15 @@
|
||||
* engines/c/main.c
|
||||
* -----------------
|
||||
* C 语言动力学模拟引擎。
|
||||
* 与 Python 版 (compute.py) 算法保持一致。
|
||||
*
|
||||
* 输入: param.json 数值参数（Python 从 YAML 转换得来）
|
||||
* <input_dir>/coord.txt
|
||||
* <input_dir>/connection.txt
|
||||
* <input_dir>/bond.txt
|
||||
* 输出: <output_dir>/trajectory.txt (JSON 格式，与 Python 版兼容)
|
||||
*
|
||||
* 编译: make
|
||||
* 编译: cmake --build build --target dynamics_c
|
||||
* 用法: ./build/dynamics_c <input_dir> <output_dir> <param_json_path>
|
||||
*/
|
||||
|
||||
@@ -27,8 +29,14 @@ typedef struct {
|
||||
double DT; /* 时间步长 */
|
||||
int NSTEP; /* 抽帧间隔 */
|
||||
int warmup_steps; /* 预热步数 */
|
||||
char method[32]; /* 算法名称 */
|
||||
double G[3]; /* 重力分量 */
|
||||
double B[3]; /* 阻尼分量 */
|
||||
int gravity_field; /* 均匀重力场开关 */
|
||||
int gravity_interaction; /* 原子间万有引力开关 */
|
||||
int elastic_force; /* 弹簧键力开关 */
|
||||
int damping_force; /* 阻尼开关 */
|
||||
double gravity_strength; /* 万有引力强度 */
|
||||
} SimParams;
|
||||
|
||||
/* ========================================================================
|
||||
@@ -79,7 +87,7 @@ static void *xmalloc(size_t sz) {
|
||||
return p;
|
||||
}
|
||||
|
||||
/* 从 JSON 中读取 double 值。简单实现：搜索 key 并解析 */
|
||||
/* 从 JSON 中读取 double 值 */
|
||||
static double json_read_double(const char *json, const char *key) {
|
||||
char search[256];
|
||||
snprintf(search, sizeof(search), "\"%s\"", key);
|
||||
@@ -96,6 +104,23 @@ static int json_read_int(const char *json, const char *key) {
|
||||
return (int)json_read_double(json, key);
|
||||
}
|
||||
|
||||
/* 从 JSON 中读取字符串值（写入 dst，最多 dst_sz 字节） */
|
||||
static void json_read_string(const char *json, const char *key, char *dst, int dst_sz) {
|
||||
char search[256];
|
||||
snprintf(search, sizeof(search), "\"%s\"", key);
|
||||
const char *p = strstr(json, search);
|
||||
if (!p) { dst[0] = '\0'; return; }
|
||||
p = strchr(p, ':');
|
||||
if (!p) { dst[0] = '\0'; return; }
|
||||
p++;
|
||||
while (*p == ' ' || *p == '\t' || *p == '\n') p++;
|
||||
if (*p != '"') { dst[0] = '\0'; return; }
|
||||
p++;
|
||||
int i = 0;
|
||||
while (*p && *p != '"' && i < dst_sz - 1) { dst[i++] = *p++; }
|
||||
dst[i] = '\0';
|
||||
}
|
||||
|
||||
/* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */
|
||||
static void json_read_double3(const char *json, const char *key, double out[3]) {
|
||||
char search[256];
|
||||
@@ -112,6 +137,12 @@ static void json_read_double3(const char *json, const char *key, double out[3])
|
||||
}
|
||||
|
||||
/* 读取 param.json */
|
||||
static int g_gravity_field = 1;
|
||||
static int g_gravity_interaction = 0;
|
||||
static int g_elastic_force = 1;
|
||||
static int g_damping_force = 0;
|
||||
static double g_gravity_strength = 1.0;
|
||||
|
||||
static SimParams read_params(const char *path) {
|
||||
FILE *f = fopen(path, "rb");
|
||||
if (!f) die("无法打开 param.json");
|
||||
@@ -129,8 +160,20 @@ static SimParams read_params(const char *path) {
|
||||
p.DT = json_read_double(buf, "DT");
|
||||
p.NSTEP = json_read_int(buf, "NSTEP");
|
||||
p.warmup_steps = json_read_int(buf, "warmup_steps");
|
||||
strcpy(p.method, "leapfrog"); /* 默认 */
|
||||
json_read_string(buf, "method", p.method, sizeof(p.method));
|
||||
json_read_double3(buf, "G", p.G);
|
||||
json_read_double3(buf, "B", p.B);
|
||||
p.gravity_field = json_read_int(buf, "gravity_field");
|
||||
p.gravity_interaction = json_read_int(buf, "gravity_interaction");
|
||||
p.elastic_force = json_read_int(buf, "elastic_force");
|
||||
p.damping_force = json_read_int(buf, "damping_force");
|
||||
p.gravity_strength = json_read_double(buf, "gravity_strength");
|
||||
g_gravity_field = p.gravity_field;
|
||||
g_gravity_interaction = p.gravity_interaction;
|
||||
g_elastic_force = p.elastic_force;
|
||||
g_damping_force = p.damping_force;
|
||||
g_gravity_strength = p.gravity_strength;
|
||||
|
||||
free(buf);
|
||||
return p;
|
||||
@@ -147,7 +190,6 @@ static AtomData read_coord(const char *input_dir) {
|
||||
char line[1024];
|
||||
if (!fgets(line, sizeof(line), f)) die("coord.txt 为空");
|
||||
|
||||
/* 先读入所有行到内存 */
|
||||
int capacity = 16;
|
||||
AtomData a;
|
||||
a.n_atoms = 0;
|
||||
@@ -173,10 +215,9 @@ static AtomData read_coord(const char *input_dir) {
|
||||
int n_parsed = sscanf(line, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %d",
|
||||
&id, &mass, &rad, &px, &py, &pz, &vx, &vy, &vz, &fx, &fy, &fz);
|
||||
if (n_parsed == 9) {
|
||||
/* 无固定约束列 */
|
||||
fx = fy = fz = 0;
|
||||
} else if (n_parsed != 12) {
|
||||
continue; /* 跳过空行/格式不符的行 */
|
||||
continue;
|
||||
}
|
||||
int i = a.n_atoms;
|
||||
a.atom_ids[i] = id;
|
||||
@@ -204,13 +245,11 @@ static BondData read_bonds(const char *input_dir, const AtomData *atoms) {
|
||||
|
||||
snprintf(path, sizeof(path), "%s/connection.txt", input_dir);
|
||||
FILE *f = fopen(path, "r");
|
||||
if (!f) return b; /* 无成键 */
|
||||
if (!f) return b;
|
||||
|
||||
/* 跳过第一行表头 */
|
||||
char line[256];
|
||||
if (!fgets(line, sizeof(line), f)) { fclose(f); return b; }
|
||||
|
||||
/* 先数行数 */
|
||||
int n_lines = 0, tmp_a, tmp_b;
|
||||
char bond_name[256];
|
||||
while (fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name) == 3) n_lines++;
|
||||
@@ -223,24 +262,21 @@ static BondData read_bonds(const char *input_dir, const AtomData *atoms) {
|
||||
b.stiffness = (double*)xmalloc(n_lines * sizeof(double));
|
||||
b.rest_lengths = (double*)xmalloc(n_lines * sizeof(double));
|
||||
|
||||
/* 读取 bond.txt 获取劲度系数和平衡长度 */
|
||||
char bond_path[512];
|
||||
snprintf(bond_path, sizeof(bond_path), "%s/bond.txt", input_dir);
|
||||
FILE *fb = fopen(bond_path, "r");
|
||||
/* bond.txt 格式: name stiffness rest_length */
|
||||
|
||||
for (int i = 0; i < n_lines; i++) {
|
||||
fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name);
|
||||
/* 转换原子编号为 0-based 索引 */
|
||||
b.pairs[i*2+0] = tmp_a - 1;
|
||||
b.pairs[i*2+1] = tmp_b - 1;
|
||||
/* 查找 bond.txt 中对应的参数 */
|
||||
b.stiffness[i] = 1.0;
|
||||
b.rest_lengths[i] = 2.0; /* 默认值 */
|
||||
b.rest_lengths[i] = 2.0;
|
||||
if (fb) {
|
||||
char name[256];
|
||||
char name[256], header[256];
|
||||
double k, r0;
|
||||
rewind(fb);
|
||||
fgets(header, sizeof(header), fb); // 跳过表头行
|
||||
while (fscanf(fb, "%s %lf %lf", name, &k, &r0) == 3) {
|
||||
if (strcmp(name, bond_name) == 0) {
|
||||
b.stiffness[i] = k;
|
||||
@@ -256,9 +292,10 @@ static BondData read_bonds(const char *input_dir, const AtomData *atoms) {
|
||||
}
|
||||
|
||||
/* ========================================================================
|
||||
* 物理核心
|
||||
* 物理核心（与 Python compute.py 对应）
|
||||
* ======================================================================== */
|
||||
|
||||
/* 加速度计算（各力独立开关控制） */
|
||||
static void compute_acceleration(
|
||||
int n, const double *x, const double *y, const double *z,
|
||||
const double *vx, const double *vy, const double *vz,
|
||||
@@ -266,50 +303,124 @@ static void compute_acceleration(
|
||||
const BondData *bonds,
|
||||
double *ax, double *ay, double *az)
|
||||
{
|
||||
/* 先清零 */
|
||||
for (int i = 0; i < n; i++) {
|
||||
double inv_m = 1.0 / m[i];
|
||||
ax[i] = (m[i] * G[0] - B[0] * vx[i]) * inv_m;
|
||||
ay[i] = (m[i] * G[1] - B[1] * vy[i]) * inv_m;
|
||||
az[i] = (m[i] * G[2] - B[2] * vz[i]) * inv_m;
|
||||
ax[i] = 0.0; ay[i] = 0.0; az[i] = 0.0;
|
||||
}
|
||||
|
||||
/* 弹簧力 */
|
||||
for (int b = 0; b < bonds->n_bonds; b++) {
|
||||
int i = bonds->pairs[b*2+0];
|
||||
int j = bonds->pairs[b*2+1];
|
||||
double dx = x[j] - x[i];
|
||||
double dy = y[j] - y[i];
|
||||
double dz = z[j] - z[i];
|
||||
double dist = sqrt(dx*dx + dy*dy + dz*dz);
|
||||
if (dist < 1e-12) continue;
|
||||
double stretch = dist - bonds->rest_lengths[b];
|
||||
double fmag = bonds->stiffness[b] * stretch;
|
||||
double ux = dx / dist, uy = dy / dist, uz = dz / dist;
|
||||
double fx = fmag * ux, fy = fmag * uy, fz = fmag * uz;
|
||||
ax[i] += fx / m[i]; ay[i] += fy / m[i]; az[i] += fz / m[i];
|
||||
ax[j] -= fx / m[j]; ay[j] -= fy / m[j]; az[j] -= fz / m[j];
|
||||
/* 均匀重力场 */
|
||||
if (g_gravity_field) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
ax[i] += G[0];
|
||||
ay[i] += G[1];
|
||||
az[i] += G[2];
|
||||
}
|
||||
}
|
||||
|
||||
/* 阻尼 */
|
||||
if (g_damping_force) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
ax[i] -= B[0] * vx[i] / m[i];
|
||||
ay[i] -= B[1] * vy[i] / m[i];
|
||||
az[i] -= B[2] * vz[i] / m[i];
|
||||
}
|
||||
}
|
||||
|
||||
/* 弹簧键力 */
|
||||
if (g_elastic_force) {
|
||||
for (int b = 0; b < bonds->n_bonds; b++) {
|
||||
int i = bonds->pairs[b*2+0];
|
||||
int j = bonds->pairs[b*2+1];
|
||||
double dx = x[j] - x[i];
|
||||
double dy = y[j] - y[i];
|
||||
double dz = z[j] - z[i];
|
||||
double dist = sqrt(dx*dx + dy*dy + dz*dz);
|
||||
if (dist < 1e-12) continue;
|
||||
double stretch = dist - bonds->rest_lengths[b];
|
||||
double fmag = bonds->stiffness[b] * stretch;
|
||||
double ux = dx / dist, uy = dy / dist, uz = dz / dist;
|
||||
double fx = fmag * ux, fy = fmag * uy, fz = fmag * uz;
|
||||
ax[i] += fx / m[i]; ay[i] += fy / m[i]; az[i] += fz / m[i];
|
||||
ax[j] -= fx / m[j]; ay[j] -= fy / m[j]; az[j] -= fz / m[j];
|
||||
}
|
||||
}
|
||||
|
||||
/* 万有引力（所有原子对之间） */
|
||||
if (g_gravity_interaction) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
for (int j = i + 1; j < n; j++) {
|
||||
double dx = x[j] - x[i];
|
||||
double dy = y[j] - y[i];
|
||||
double dz = z[j] - z[i];
|
||||
double r2 = dx*dx + dy*dy + dz*dz;
|
||||
if (r2 <= 1e-12) continue;
|
||||
double r = sqrt(r2);
|
||||
double f_mag = g_gravity_strength * m[i] * m[j] / r2;
|
||||
double fx_g = f_mag * dx / r;
|
||||
double fy_g = f_mag * dy / r;
|
||||
double fz_g = f_mag * dz / r;
|
||||
ax[i] += fx_g / m[i]; ay[i] += fy_g / m[i]; az[i] += fz_g / m[i];
|
||||
ax[j] -= fx_g / m[j]; ay[j] -= fy_g / m[j]; az[j] -= fz_g / m[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* clamp */
|
||||
static double clamp(double v, double lo, double hi) {
|
||||
return fmin(fmax(v, lo), hi);
|
||||
/* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
|
||||
static void limit_in_box(double *pos, double *vel, double lo, double hi) {
|
||||
if (*pos > hi) { *pos = hi; *vel = -*vel; }
|
||||
if (*pos < lo) { *pos = lo; *vel = -*vel; }
|
||||
}
|
||||
|
||||
/* 蛙跳法一步 */
|
||||
static void leapfrog_step(
|
||||
/* ── 显式欧拉法 ──────────── */
|
||||
static void explicit_euler_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData *bonds, const int *fixed,
|
||||
double box_a, double dt)
|
||||
const BondData *bonds, const int *fixed, double dt)
|
||||
{
|
||||
double *ax = (double*)alloca(n * sizeof(double));
|
||||
double *ay = (double*)alloca(n * sizeof(double));
|
||||
double *az = (double*)alloca(n * sizeof(double));
|
||||
|
||||
compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az);
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
x[i] += vx[i] * dt;
|
||||
y[i] += vy[i] * dt;
|
||||
z[i] += vz[i] * dt;
|
||||
vx[i] += ax[i] * dt;
|
||||
vy[i] += ay[i] * dt;
|
||||
vz[i] += az[i] * dt;
|
||||
}
|
||||
}
|
||||
|
||||
/* ── 隐式欧拉法 ──────────── */
|
||||
static void implicit_euler_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData *bonds, const int *fixed, double dt)
|
||||
{
|
||||
double *vxn = (double*)alloca(n * sizeof(double));
|
||||
double *vyn = (double*)alloca(n * sizeof(double));
|
||||
double *vzn = (double*)alloca(n * sizeof(double));
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) {
|
||||
vxn[i] = 0; vyn[i] = 0; vzn[i] = 0; continue;
|
||||
}
|
||||
double gx = B[0] / m[i], gy = B[1] / m[i], gz = B[2] / m[i];
|
||||
vxn[i] = (vx[i] + G[0] * dt) / (1.0 + gx * dt);
|
||||
vyn[i] = (vy[i] + G[1] * dt) / (1.0 + gy * dt);
|
||||
vzn[i] = (vz[i] + G[2] * dt) / (1.0 + gz * dt);
|
||||
}
|
||||
|
||||
double *ax = (double*)alloca(n * sizeof(double));
|
||||
double *ay = (double*)alloca(n * sizeof(double));
|
||||
double *az = (double*)alloca(n * sizeof(double));
|
||||
compute_acceleration(n, x, y, z, vxn, vyn, vzn, m, G, B, bonds, ax, ay, az);
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
vx[i] += ax[i] * dt;
|
||||
@@ -318,16 +429,148 @@ static void leapfrog_step(
|
||||
x[i] += vx[i] * dt;
|
||||
y[i] += vy[i] * dt;
|
||||
z[i] += vz[i] * dt;
|
||||
/* wrap */
|
||||
x[i] = clamp(x[i], -box_a, box_a);
|
||||
y[i] = clamp(y[i], -box_a, box_a);
|
||||
z[i] = clamp(z[i], -box_a, box_a);
|
||||
}
|
||||
}
|
||||
|
||||
/* ========================================================================
|
||||
* JSON 输出
|
||||
* ======================================================================== */
|
||||
/* ── 中点法 ──────────── */
|
||||
static void midpoint_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData *bonds, const int *fixed, double dt)
|
||||
{
|
||||
double *ax = (double*)alloca(n * sizeof(double));
|
||||
double *ay = (double*)alloca(n * sizeof(double));
|
||||
double *az = (double*)alloca(n * sizeof(double));
|
||||
compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az);
|
||||
|
||||
double *xm = (double*)alloca(n * sizeof(double));
|
||||
double *ym = (double*)alloca(n * sizeof(double));
|
||||
double *zm = (double*)alloca(n * sizeof(double));
|
||||
double *vxm = (double*)alloca(n * sizeof(double));
|
||||
double *vym = (double*)alloca(n * sizeof(double));
|
||||
double *vzm = (double*)alloca(n * sizeof(double));
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) {
|
||||
xm[i]=ym[i]=zm[i]=vxm[i]=vym[i]=vzm[i]=0; continue;
|
||||
}
|
||||
xm[i] = x[i] + 0.5 * vx[i] * dt;
|
||||
ym[i] = y[i] + 0.5 * vy[i] * dt;
|
||||
zm[i] = z[i] + 0.5 * vz[i] * dt;
|
||||
vxm[i] = vx[i] + 0.5 * ax[i] * dt;
|
||||
vym[i] = vy[i] + 0.5 * ay[i] * dt;
|
||||
vzm[i] = vz[i] + 0.5 * az[i] * dt;
|
||||
x[i] = x[i] + vxm[i] * dt;
|
||||
y[i] = y[i] + vym[i] * dt;
|
||||
z[i] = z[i] + vzm[i] * dt;
|
||||
}
|
||||
|
||||
compute_acceleration(n, xm, ym, zm, vxm, vym, vzm, m, G, B, bonds, ax, ay, az);
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
vx[i] += ax[i] * dt;
|
||||
vy[i] += ay[i] * dt;
|
||||
vz[i] += az[i] * dt;
|
||||
}
|
||||
}
|
||||
|
||||
/* ── 蛙跳法（Velocity-Verlet）── */
|
||||
static void leapfrog_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData *bonds, const int *fixed, double dt)
|
||||
{
|
||||
double *ax = (double*)alloca(n * sizeof(double));
|
||||
double *ay = (double*)alloca(n * sizeof(double));
|
||||
double *az = (double*)alloca(n * sizeof(double));
|
||||
compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az);
|
||||
|
||||
/* 半推速度：v_half = v + 0.5*a*dt */
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
vx[i] += ax[i] * dt * 0.5;
|
||||
vy[i] += ay[i] * dt * 0.5;
|
||||
vz[i] += az[i] * dt * 0.5;
|
||||
}
|
||||
|
||||
/* 全推位置（不含边界）*/
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
x[i] += vx[i] * dt; /* vx 此时是 v_half */
|
||||
y[i] += vy[i] * dt;
|
||||
z[i] += vz[i] * dt;
|
||||
}
|
||||
|
||||
/* 显式预测器：v_pred = v_half + 0.5*a_old*dt，用第一次加速度外推半步
|
||||
包含重力+阻尼+弹簧的所有贡献（标准 Velocity-Verlet 预测步）*/
|
||||
double *pred_vx = (double*)alloca(n * sizeof(double));
|
||||
double *pred_vy = (double*)alloca(n * sizeof(double));
|
||||
double *pred_vz = (double*)alloca(n * sizeof(double));
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
pred_vx[i] = vx[i] + 0.5 * ax[i] * dt;
|
||||
pred_vy[i] = vy[i] + 0.5 * ay[i] * dt;
|
||||
pred_vz[i] = vz[i] + 0.5 * az[i] * dt;
|
||||
}
|
||||
|
||||
/* 用新位置 + 预测速度重算加速度 */
|
||||
compute_acceleration(n, x, y, z, pred_vx, pred_vy, pred_vz, m, G, B, bonds, ax, ay, az);
|
||||
|
||||
/* 速度后半步：v = v_half + 0.5*a_next*dt
|
||||
vx 仍为 v_half（未被覆盖）*/
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
vx[i] += ax[i] * dt * 0.5;
|
||||
vy[i] += ay[i] * dt * 0.5;
|
||||
vz[i] += az[i] * dt * 0.5;
|
||||
}
|
||||
}
|
||||
|
||||
/* ── 分发器：调用对应积分方法 + 边界条件 + 自由度约束（与 Python 一致）── */
|
||||
static void apply_step(
|
||||
const char *method,
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData *bonds, const int *fixed,
|
||||
const double *pos_0,
|
||||
double box_a, double dt)
|
||||
{
|
||||
if (strcmp(method, "explicit_euler") == 0) {
|
||||
explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
|
||||
} else if (strcmp(method, "implicit_euler") == 0) {
|
||||
implicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
|
||||
} else if (strcmp(method, "midpoint") == 0) {
|
||||
midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
|
||||
} else if (strcmp(method, "leapfrog") == 0) {
|
||||
leapfrog_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt);
|
||||
} else {
|
||||
fprintf(stderr, "[C-engine] 未知算法: %s\n", method);
|
||||
exit(1);
|
||||
}
|
||||
|
||||
/* 边界条件（与 Python Limit_in_box 一致） */
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
limit_in_box(&x[i], &vx[i], -box_a, box_a);
|
||||
limit_in_box(&y[i], &vy[i], -box_a, box_a);
|
||||
limit_in_box(&z[i], &vz[i], -box_a, box_a);
|
||||
}
|
||||
|
||||
/* 逐自由度固定约束（与 Python apply_fixed_constraints 一致） */
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0]) { x[i] = pos_0[i*3+0]; vx[i] = 0.0; }
|
||||
if (fixed[i*3+1]) { y[i] = pos_0[i*3+1]; vy[i] = 0.0; }
|
||||
if (fixed[i*3+2]) { z[i] = pos_0[i*3+2]; vz[i] = 0.0; }
|
||||
}
|
||||
}
|
||||
|
||||
// ========================================================================
|
||||
// JSON 输出
|
||||
// ========================================================================
|
||||
|
||||
static void write_trajectory_json(const char *path, const Trajectory *traj,
|
||||
const SimParams *params, const AtomData *atoms,
|
||||
@@ -338,7 +581,6 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
|
||||
|
||||
fprintf(f, "{\n");
|
||||
|
||||
/* traj_x, traj_y, ... */
|
||||
const char *names[] = {"traj_x","traj_y","traj_z","traj_vx","traj_vy","traj_vz"};
|
||||
double *arrs[] = {traj->x, traj->y, traj->z, traj->vx, traj->vy, traj->vz};
|
||||
|
||||
@@ -355,7 +597,6 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
|
||||
fputc('\n', f);
|
||||
}
|
||||
fprintf(f, " ]");
|
||||
// 每个轨迹数组后加逗号（后面还有标量参数）
|
||||
fputc(',', f);
|
||||
fputc('\n', f);
|
||||
}
|
||||
@@ -364,45 +605,43 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
|
||||
fprintf(f, " \"NT\": %d,\n", params->NT);
|
||||
fprintf(f, " \"DT\": %.15g,\n", params->DT);
|
||||
fprintf(f, " \"NSTEP\": %d,\n", params->NSTEP);
|
||||
fprintf(f, " \"method\": \"leapfrog\",\n");
|
||||
fprintf(f, " \"method\": \"%s\",\n", params->method);
|
||||
fprintf(f, " \"warmup_steps\": %d,\n", params->warmup_steps);
|
||||
fprintf(f, " \"G\": [%.15g, %.15g, %.15g],\n", params->G[0], params->G[1], params->G[2]);
|
||||
fprintf(f, " \"B\": [%.15g, %.15g, %.15g],\n", params->B[0], params->B[1], params->B[2]);
|
||||
|
||||
/* 原子信息 */
|
||||
fprintf(f, " \"atom_ids\": [");
|
||||
for (int i = 0; i < atoms->n_atoms; i++) {
|
||||
if (i > 0) fputc(',', f);
|
||||
fprintf(f, "%d", atoms->atom_ids[i]);
|
||||
if (i < atoms->n_atoms - 1) fputc(',', f);
|
||||
}
|
||||
fprintf(f, "],\n");
|
||||
|
||||
fprintf(f, " \"atom_masses\": [");
|
||||
for (int i = 0; i < atoms->n_atoms; i++) {
|
||||
if (i > 0) fputc(',', f);
|
||||
fprintf(f, "%.15g", atoms->masses[i]);
|
||||
if (i < atoms->n_atoms - 1) fputc(',', f);
|
||||
}
|
||||
fprintf(f, "],\n");
|
||||
|
||||
/* 成键 */
|
||||
fprintf(f, " \"bond_pairs\": [");
|
||||
for (int b = 0; b < bonds->n_bonds; b++) {
|
||||
if (b > 0) fputc(',', f);
|
||||
fprintf(f, "[%d, %d]", bonds->pairs[b*2], bonds->pairs[b*2+1]);
|
||||
if (b < bonds->n_bonds - 1) fputc(',', f);
|
||||
}
|
||||
fprintf(f, "],\n");
|
||||
|
||||
fprintf(f, " \"bond_stiffness\": [");
|
||||
for (int b = 0; b < bonds->n_bonds; b++) {
|
||||
if (b > 0) fputc(',', f);
|
||||
fprintf(f, "%.15g", bonds->stiffness[b]);
|
||||
if (b < bonds->n_bonds - 1) fputc(',', f);
|
||||
}
|
||||
fprintf(f, "],\n");
|
||||
|
||||
fprintf(f, " \"bond_rest_lengths\": [");
|
||||
for (int b = 0; b < bonds->n_bonds; b++) {
|
||||
if (b > 0) fputc(',', f);
|
||||
fprintf(f, "%.15g", bonds->rest_lengths[b]);
|
||||
if (b < bonds->n_bonds - 1) fputc(',', f);
|
||||
}
|
||||
fprintf(f, "]\n");
|
||||
|
||||
@@ -410,9 +649,9 @@ static void write_trajectory_json(const char *path, const Trajectory *traj,
|
||||
fclose(f);
|
||||
}
|
||||
|
||||
/* ========================================================================
|
||||
* 主函数
|
||||
* ======================================================================== */
|
||||
// ========================================================================
|
||||
// 主函数
|
||||
// ========================================================================
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
if (argc < 4) {
|
||||
@@ -425,15 +664,13 @@ int main(int argc, char **argv) {
|
||||
|
||||
clock_t t0 = clock();
|
||||
|
||||
/* 读取参数和输入 */
|
||||
SimParams params = read_params(param_path);
|
||||
AtomData atoms = read_coord(input_dir);
|
||||
BondData bonds = read_bonds(input_dir, &atoms);
|
||||
|
||||
printf("[C-engine] 原子数=%d, 键数=%d, NT=%d, DT=%.6g\n",
|
||||
atoms.n_atoms, bonds.n_bonds, params.NT, params.DT);
|
||||
printf("[C-engine] 原子数=%d, 键数=%d, NT=%d, DT=%.6g, method=%s\n",
|
||||
atoms.n_atoms, bonds.n_bonds, params.NT, params.DT, params.method);
|
||||
|
||||
/* 初始化位置/速度 */
|
||||
int n = atoms.n_atoms;
|
||||
double *x = (double*)xmalloc(n * sizeof(double));
|
||||
double *y = (double*)xmalloc(n * sizeof(double));
|
||||
@@ -450,30 +687,28 @@ int main(int argc, char **argv) {
|
||||
vz[i] = atoms.vel_0[i*3+2];
|
||||
}
|
||||
|
||||
/* 分配轨迹缓冲区 */
|
||||
/* 分配轨迹缓冲区：改用 record_steps */
|
||||
int record_steps = params.NT - params.warmup_steps;
|
||||
Trajectory traj;
|
||||
traj.n_steps = params.NT;
|
||||
traj.n_steps = record_steps;
|
||||
traj.n_atoms = n;
|
||||
traj.x = (double*)xmalloc(params.NT * n * sizeof(double) * 6);
|
||||
traj.y = traj.x + params.NT * n;
|
||||
traj.z = traj.y + params.NT * n;
|
||||
traj.vx = traj.z + params.NT * n;
|
||||
traj.vy = traj.vx + params.NT * n;
|
||||
traj.vz = traj.vy + params.NT * n;
|
||||
traj.x = (double*)xmalloc(record_steps * n * sizeof(double) * 6);
|
||||
traj.y = traj.x + record_steps * n;
|
||||
traj.z = traj.y + record_steps * n;
|
||||
traj.vx = traj.z + record_steps * n;
|
||||
traj.vy = traj.vx + record_steps * n;
|
||||
traj.vz = traj.vy + record_steps * n;
|
||||
|
||||
/* 预热 */
|
||||
for (int s = 0; s < params.warmup_steps; s++) {
|
||||
leapfrog_step(n, x, y, z, vx, vy, vz,
|
||||
atoms.masses, params.G, params.B, &bonds, atoms.fixed,
|
||||
params.box_a, params.DT);
|
||||
apply_step(params.method, n, x, y, z, vx, vy, vz,
|
||||
atoms.masses, params.G, params.B, &bonds, atoms.fixed,
|
||||
atoms.pos_0,
|
||||
params.box_a, params.DT);
|
||||
}
|
||||
|
||||
/* 记录 */
|
||||
for (int s = 0; s < params.NT; s++) {
|
||||
leapfrog_step(n, x, y, z, vx, vy, vz,
|
||||
atoms.masses, params.G, params.B, &bonds, atoms.fixed,
|
||||
params.box_a, params.DT);
|
||||
/* 保存这一帧 */
|
||||
for (int s = 0; s < record_steps; s++) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
traj.x[ s * n + i] = x[i];
|
||||
traj.y[ s * n + i] = y[i];
|
||||
@@ -482,16 +717,19 @@ int main(int argc, char **argv) {
|
||||
traj.vy[s * n + i] = vy[i];
|
||||
traj.vz[s * n + i] = vz[i];
|
||||
}
|
||||
apply_step(params.method, n, x, y, z, vx, vy, vz,
|
||||
atoms.masses, params.G, params.B, &bonds, atoms.fixed,
|
||||
atoms.pos_0,
|
||||
params.box_a, params.DT);
|
||||
}
|
||||
|
||||
/* 输出轨迹 */
|
||||
char out_path[512];
|
||||
snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir);
|
||||
write_trajectory_json(out_path, &traj, ¶ms, &atoms, &bonds);
|
||||
|
||||
clock_t t1 = clock();
|
||||
double elapsed = (double)(t1 - t0) / CLOCKS_PER_SEC;
|
||||
printf("[C-engine] 计算完成: %d 步, %.3f s\n", params.NT, elapsed);
|
||||
printf("[C-engine] 计算完成: %d 步, %.3f s\n", record_steps, elapsed);
|
||||
|
||||
/* 清理 */
|
||||
free(traj.x);
|
||||
|
||||
+681
-37
@@ -1,54 +1,633 @@
|
||||
/**
|
||||
* engines/cpp/main.cpp
|
||||
* --------------------
|
||||
* C++ 动力学模拟引擎（模板/参考实现）。
|
||||
* C++ 动力学模拟引擎。
|
||||
* 与 Python 版 (compute.py) 算法保持一致。
|
||||
*
|
||||
* 接口与 C 引擎一致：
|
||||
* 输入: param.json, <input_dir>/coord.txt, connection.txt, bond.txt
|
||||
* 输出: <output_dir>/trajectory.txt (JSON)
|
||||
* 输入: param.json, <input_dir>/coord.txt, connection.txt, bond.txt
|
||||
* 输出: <output_dir>/trajectory.txt (JSON, 与 Python 版兼容)
|
||||
*
|
||||
* 编译:
|
||||
* g++ -O3 -march=native -o build/dynamics_cpp main.cpp
|
||||
* g++ -O3 -march=native -std=c++17 -o build/dynamics_cpp main.cpp
|
||||
*
|
||||
* 用法:
|
||||
* ./build/dynamics_cpp <input_dir> <output_dir> <param_json>
|
||||
*/
|
||||
|
||||
#include <iostream>
|
||||
#include <algorithm>
|
||||
#include <chrono>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <chrono>
|
||||
#include <cstring>
|
||||
|
||||
// ========================================================================
|
||||
// 配置参数（从 param.json 读取）
|
||||
// ========================================================================
|
||||
struct SimParams {
|
||||
double box_a;
|
||||
int NT;
|
||||
double DT;
|
||||
int NSTEP;
|
||||
int warmup_steps;
|
||||
double G[3];
|
||||
double B[3];
|
||||
};
|
||||
|
||||
struct Atom {
|
||||
int id;
|
||||
double mass, radius;
|
||||
double x, y, z;
|
||||
double vx, vy, vz;
|
||||
bool fx, fy, fz; // fixed constraints
|
||||
};
|
||||
|
||||
struct Bond {
|
||||
int i, j;
|
||||
double stiffness, rest_length;
|
||||
double box_a = 10.0;
|
||||
int NT = 10000;
|
||||
double DT = 0.001;
|
||||
int NSTEP = 100;
|
||||
int warmup_steps = 0;
|
||||
std::string method = "leapfrog";
|
||||
double G[3] = {0, 0, -9.8};
|
||||
double B[3] = {0, 0, 0};
|
||||
int gravity_field = 1;
|
||||
int gravity_interaction = 0;
|
||||
int elastic_force = 1;
|
||||
int damping_force = 0;
|
||||
double gravity_strength = 1.0;
|
||||
};
|
||||
|
||||
// ========================================================================
|
||||
// 请在下方实现 Leapfrog 算法
|
||||
// 参考 engines/c/main.c 中的物理核心实现
|
||||
// 原子数据
|
||||
// ========================================================================
|
||||
struct AtomData {
|
||||
std::vector<int> ids;
|
||||
std::vector<double> masses;
|
||||
std::vector<double> radii;
|
||||
std::vector<double> pos_0; // (n_atoms * 3)
|
||||
std::vector<double> vel_0; // (n_atoms * 3)
|
||||
std::vector<int> fixed; // (n_atoms * 3), 0/1 flags
|
||||
};
|
||||
|
||||
// ========================================================================
|
||||
// 成键数据
|
||||
// ========================================================================
|
||||
struct BondData {
|
||||
std::vector<int> pairs; // (n_bonds * 2)
|
||||
std::vector<double> stiffness;
|
||||
std::vector<double> rest_lengths;
|
||||
};
|
||||
|
||||
// ========================================================================
|
||||
// 辅助函数
|
||||
// ========================================================================
|
||||
|
||||
static void die(const std::string &msg) {
|
||||
std::cerr << "[C++-engine] 错误: " << msg << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
/* 读整个文件为字符串 */
|
||||
static std::string read_file(const std::string &path) {
|
||||
std::ifstream f(path, std::ios::binary);
|
||||
if (!f) die("无法打开 " + path);
|
||||
std::ostringstream ss;
|
||||
ss << f.rdbuf();
|
||||
return ss.str();
|
||||
}
|
||||
|
||||
/* 从 JSON 中查找 key，返回冒号后的数值 */
|
||||
static double json_read_double(const std::string &json, const std::string &key) {
|
||||
auto pos = json.find("\"" + key + "\"");
|
||||
if (pos == std::string::npos) return 0.0;
|
||||
pos = json.find(':', pos);
|
||||
if (pos == std::string::npos) return 0.0;
|
||||
while (pos < json.size() && (json[pos] == ':' || json[pos] == ' ' || json[pos] == '\t' || json[pos] == '\n')) pos++;
|
||||
return std::stod(json.substr(pos));
|
||||
}
|
||||
|
||||
static int json_read_int(const std::string &json, const std::string &key) {
|
||||
return static_cast<int>(json_read_double(json, key));
|
||||
}
|
||||
|
||||
/* 从 JSON 中读取字符串值 */
|
||||
static std::string json_read_string(const std::string &json, const std::string &key) {
|
||||
auto pos = json.find("\"" + key + "\"");
|
||||
if (pos == std::string::npos) return "";
|
||||
pos = json.find(':', pos);
|
||||
if (pos == std::string::npos) return "";
|
||||
while (pos < json.size() && (json[pos] == ':' || json[pos] == ' ' || json[pos] == '\t' || json[pos] == '\n')) pos++;
|
||||
if (pos >= json.size()) return "";
|
||||
// 找到引号
|
||||
if (json[pos] != '"') return "";
|
||||
pos++;
|
||||
std::string result;
|
||||
while (pos < json.size() && json[pos] != '"') {
|
||||
result += json[pos];
|
||||
pos++;
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
/* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */
|
||||
static void json_read_double3(const std::string &json, const std::string &key, double out[3]) {
|
||||
auto pos = json.find("\"" + key + "\"");
|
||||
if (pos == std::string::npos) { out[0] = out[1] = out[2] = 0; return; }
|
||||
pos = json.find('[', pos);
|
||||
if (pos == std::string::npos) { out[0] = out[1] = out[2] = 0; return; }
|
||||
pos++;
|
||||
for (int i = 0; i < 3; i++) {
|
||||
while (pos < json.size() && (json[pos] == ' ' || json[pos] == '\t' || json[pos] == '\n' || json[pos] == ',')) pos++;
|
||||
char *end;
|
||||
out[i] = std::strtod(json.c_str() + pos, &end);
|
||||
pos = end - json.c_str();
|
||||
}
|
||||
}
|
||||
|
||||
/* 解析 param.json */
|
||||
static SimParams read_params(const std::string &path) {
|
||||
std::string buf = read_file(path);
|
||||
SimParams p;
|
||||
p.box_a = json_read_double(buf, "box_a");
|
||||
p.NT = json_read_int(buf, "NT");
|
||||
p.DT = json_read_double(buf, "DT");
|
||||
p.NSTEP = json_read_int(buf, "NSTEP");
|
||||
p.warmup_steps = json_read_int(buf, "warmup_steps");
|
||||
std::string m = json_read_string(buf, "method");
|
||||
if (!m.empty()) p.method = m;
|
||||
json_read_double3(buf, "G", p.G);
|
||||
json_read_double3(buf, "B", p.B);
|
||||
p.gravity_field = json_read_int(buf, "gravity_field");
|
||||
p.gravity_interaction = json_read_int(buf, "gravity_interaction");
|
||||
p.elastic_force = json_read_int(buf, "elastic_force");
|
||||
p.damping_force = json_read_int(buf, "damping_force");
|
||||
p.gravity_strength = json_read_double(buf, "gravity_strength");
|
||||
return p;
|
||||
}
|
||||
|
||||
/* 读取 coord.txt */
|
||||
static AtomData read_coord(const std::string &input_dir) {
|
||||
std::string path = input_dir + "/coord.txt";
|
||||
std::ifstream f(path);
|
||||
if (!f) die("无法打开 " + path);
|
||||
|
||||
std::string header;
|
||||
std::getline(f, header); // 跳过表头
|
||||
|
||||
AtomData a;
|
||||
int id, fx, fy, fz;
|
||||
double mass, rad, px, py, pz, vx, vy, vz;
|
||||
std::string line;
|
||||
|
||||
while (std::getline(f, line)) {
|
||||
if (line.empty() || line[0] == '#') continue;
|
||||
int n_parsed = std::sscanf(line.c_str(), "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %d",
|
||||
&id, &mass, &rad, &px, &py, &pz, &vx, &vy, &vz, &fx, &fy, &fz);
|
||||
if (n_parsed == 9) {
|
||||
fx = fy = fz = 0;
|
||||
} else if (n_parsed != 12) {
|
||||
continue;
|
||||
}
|
||||
a.ids.push_back(id);
|
||||
a.masses.push_back(mass);
|
||||
a.radii.push_back(rad);
|
||||
a.pos_0.push_back(px); a.pos_0.push_back(py); a.pos_0.push_back(pz);
|
||||
a.vel_0.push_back(vx); a.vel_0.push_back(vy); a.vel_0.push_back(vz);
|
||||
a.fixed.push_back(fx); a.fixed.push_back(fy); a.fixed.push_back(fz);
|
||||
}
|
||||
|
||||
if (a.ids.empty()) die("coord.txt 中没有原子数据");
|
||||
return a;
|
||||
}
|
||||
|
||||
/* 读取 connection.txt 和 bond.txt */
|
||||
static BondData read_bonds(const std::string &input_dir) {
|
||||
BondData b;
|
||||
std::string conn_path = input_dir + "/connection.txt";
|
||||
std::ifstream f(conn_path);
|
||||
if (!f) return b; // 无成键
|
||||
|
||||
std::string header;
|
||||
std::getline(f, header); // 跳过表头
|
||||
|
||||
// 先读取 bond.txt 获得键参数映射
|
||||
std::string bond_path = input_dir + "/bond.txt";
|
||||
std::ifstream fb(bond_path);
|
||||
|
||||
int a1, a2;
|
||||
std::string bond_name;
|
||||
std::vector<std::tuple<int, int, std::string>> conn_lines;
|
||||
while (f >> a1 >> a2 >> bond_name) {
|
||||
conn_lines.emplace_back(a1 - 1, a2 - 1, bond_name);
|
||||
}
|
||||
|
||||
for (auto &[i, j, name] : conn_lines) {
|
||||
double k = 1.0, r0 = 2.0;
|
||||
if (fb) {
|
||||
fb.clear();
|
||||
fb.seekg(0);
|
||||
std::string bn, header;
|
||||
double bk, br;
|
||||
std::getline(fb, header); // 跳过表头行
|
||||
while (fb >> bn >> bk >> br) {
|
||||
if (bn == name) {
|
||||
k = bk;
|
||||
r0 = br;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
b.pairs.push_back(i);
|
||||
b.pairs.push_back(j);
|
||||
b.stiffness.push_back(k);
|
||||
b.rest_lengths.push_back(r0);
|
||||
}
|
||||
|
||||
return b;
|
||||
}
|
||||
|
||||
// ========================================================================
|
||||
// 物理核心
|
||||
// ========================================================================
|
||||
|
||||
/* 加速度计算（各力独立开关控制）——与 Python compute_acceleration 一致 */
|
||||
static void compute_acceleration(
|
||||
int n,
|
||||
const double *x, const double *y, const double *z,
|
||||
const double *vx, const double *vy, const double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData &bonds,
|
||||
int gravity_field, int gravity_interaction,
|
||||
int elastic_force, int damping_force,
|
||||
double gravity_strength,
|
||||
double *ax, double *ay, double *az)
|
||||
{
|
||||
// 清零
|
||||
std::fill(ax, ax + n, 0.0);
|
||||
std::fill(ay, ay + n, 0.0);
|
||||
std::fill(az, az + n, 0.0);
|
||||
|
||||
// 均匀重力场
|
||||
if (gravity_field) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
ax[i] += G[0];
|
||||
ay[i] += G[1];
|
||||
az[i] += G[2];
|
||||
}
|
||||
}
|
||||
|
||||
// 阻尼
|
||||
if (damping_force) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
ax[i] -= B[0] * vx[i] / m[i];
|
||||
ay[i] -= B[1] * vy[i] / m[i];
|
||||
az[i] -= B[2] * vz[i] / m[i];
|
||||
}
|
||||
}
|
||||
|
||||
// 弹簧力
|
||||
if (elastic_force) {
|
||||
int nb = static_cast<int>(bonds.stiffness.size());
|
||||
for (int b = 0; b < nb; b++) {
|
||||
int i = bonds.pairs[b * 2];
|
||||
int j = bonds.pairs[b * 2 + 1];
|
||||
double dx = x[j] - x[i];
|
||||
double dy = y[j] - y[i];
|
||||
double dz = z[j] - z[i];
|
||||
double dist = std::sqrt(dx * dx + dy * dy + dz * dz);
|
||||
if (dist < 1e-12) continue;
|
||||
double stretch = dist - bonds.rest_lengths[b];
|
||||
double fmag = bonds.stiffness[b] * stretch;
|
||||
double ux = dx / dist, uy = dy / dist, uz = dz / dist;
|
||||
double fx = fmag * ux, fy = fmag * uy, fz = fmag * uz;
|
||||
ax[i] += fx / m[i]; ay[i] += fy / m[i]; az[i] += fz / m[i];
|
||||
ax[j] -= fx / m[j]; ay[j] -= fy / m[j]; az[j] -= fz / m[j];
|
||||
}
|
||||
}
|
||||
|
||||
// 万有引力（所有原子对之间）
|
||||
if (gravity_interaction) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
for (int j = i + 1; j < n; j++) {
|
||||
double dx = x[j] - x[i];
|
||||
double dy = y[j] - y[i];
|
||||
double dz = z[j] - z[i];
|
||||
double r2 = dx * dx + dy * dy + dz * dz;
|
||||
if (r2 <= 1e-12) continue;
|
||||
double r = std::sqrt(r2);
|
||||
double f_mag = gravity_strength * m[i] * m[j] / r2;
|
||||
double fx_g = f_mag * dx / r;
|
||||
double fy_g = f_mag * dy / r;
|
||||
double fz_g = f_mag * dz / r;
|
||||
ax[i] += fx_g / m[i]; ay[i] += fy_g / m[i]; az[i] += fz_g / m[i];
|
||||
ax[j] -= fx_g / m[j]; ay[j] -= fy_g / m[j]; az[j] -= fz_g / m[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */
|
||||
static void limit_in_box(double &pos, double &vel, double lo, double hi) {
|
||||
if (pos > hi) { pos = hi; vel = -vel; }
|
||||
if (pos < lo) { pos = lo; vel = -vel; }
|
||||
}
|
||||
|
||||
// ========================================================================
|
||||
// 四种积分方法（只做位置/速度更新，不含边界条件）
|
||||
// 与 Python: Explicit_Euler_Method / Implicit_Euler_Method /
|
||||
// Midpoint_Method / Leapfrog_Method 保持一致
|
||||
// ========================================================================
|
||||
|
||||
/* ── 显式欧拉法 ──────────── */
|
||||
static void explicit_euler_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData &bonds, const int *fixed, double dt,
|
||||
int gravity_field, int gravity_interaction,
|
||||
int elastic_force, int damping_force,
|
||||
double gravity_strength)
|
||||
{
|
||||
std::vector<double> ax(n), ay(n), az(n);
|
||||
compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength,
|
||||
ax.data(), ay.data(), az.data());
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
x[i] += vx[i] * dt;
|
||||
y[i] += vy[i] * dt;
|
||||
z[i] += vz[i] * dt;
|
||||
vx[i] += ax[i] * dt;
|
||||
vy[i] += ay[i] * dt;
|
||||
vz[i] += az[i] * dt;
|
||||
}
|
||||
}
|
||||
|
||||
/* ── 隐式欧拉法 ──────────── */
|
||||
static void implicit_euler_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData &bonds, const int *fixed, double dt,
|
||||
int gravity_field, int gravity_interaction,
|
||||
int elastic_force, int damping_force,
|
||||
double gravity_strength)
|
||||
{
|
||||
std::vector<double> ax(n), ay(n), az(n);
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) {
|
||||
ax[i] = ay[i] = az[i] = 0;
|
||||
continue;
|
||||
}
|
||||
double gamma_x = B[0] / m[i];
|
||||
double gamma_y = B[1] / m[i];
|
||||
double gamma_z = B[2] / m[i];
|
||||
// 隐式更新速度（重力 + 阻尼）
|
||||
double vxn = (vx[i] + G[0] * dt) / (1.0 + gamma_x * dt);
|
||||
double vyn = (vy[i] + G[1] * dt) / (1.0 + gamma_y * dt);
|
||||
double vzn = (vz[i] + G[2] * dt) / (1.0 + gamma_z * dt);
|
||||
// 用隐式速度 + 当前位置算加速度（包含各力开关）
|
||||
// 注意：Python 中 compute_acceleration(x, y, z, vx_next, ...) 用新速度+旧位置
|
||||
double tpx = x[i], tpy = y[i], tpz = z[i];
|
||||
compute_acceleration(1, &tpx, &tpy, &tpz, &vxn, &vyn, &vzn,
|
||||
&m[i], G, B, bonds,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength,
|
||||
&ax[i], &ay[i], &az[i]);
|
||||
vx[i] += ax[i] * dt;
|
||||
vy[i] += ay[i] * dt;
|
||||
vz[i] += az[i] * dt;
|
||||
x[i] += vx[i] * dt;
|
||||
y[i] += vy[i] * dt;
|
||||
z[i] += vz[i] * dt;
|
||||
}
|
||||
}
|
||||
|
||||
/* ── 中点法 ──────────── */
|
||||
static void midpoint_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData &bonds, const int *fixed, double dt,
|
||||
int gravity_field, int gravity_interaction,
|
||||
int elastic_force, int damping_force,
|
||||
double gravity_strength)
|
||||
{
|
||||
std::vector<double> ax(n), ay(n), az(n);
|
||||
compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength,
|
||||
ax.data(), ay.data(), az.data());
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
double x_mid = x[i] + 0.5 * vx[i] * dt;
|
||||
double y_mid = y[i] + 0.5 * vy[i] * dt;
|
||||
double z_mid = z[i] + 0.5 * vz[i] * dt;
|
||||
double vx_mid = vx[i] + 0.5 * ax[i] * dt;
|
||||
double vy_mid = vy[i] + 0.5 * ay[i] * dt;
|
||||
double vz_mid = vz[i] + 0.5 * az[i] * dt;
|
||||
x[i] += vx_mid * dt;
|
||||
y[i] += vy_mid * dt;
|
||||
z[i] += vz_mid * dt;
|
||||
double ax_mid, ay_mid, az_mid;
|
||||
compute_acceleration(1, &x_mid, &y_mid, &z_mid, &vx_mid, &vy_mid, &vz_mid,
|
||||
&m[i], G, B, bonds,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength,
|
||||
&ax_mid, &ay_mid, &az_mid);
|
||||
vx[i] += ax_mid * dt;
|
||||
vy[i] += ay_mid * dt;
|
||||
vz[i] += az_mid * dt;
|
||||
}
|
||||
}
|
||||
|
||||
/* ── 蛙跳法（Velocity-Verlet）——与 Python Leapfrog_Method 一致 ── */
|
||||
static void leapfrog_full_step(
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData &bonds, const int *fixed, double dt,
|
||||
int gravity_field, int gravity_interaction,
|
||||
int elastic_force, int damping_force,
|
||||
double gravity_strength)
|
||||
{
|
||||
// 第一次加速度
|
||||
std::vector<double> ax(n), ay(n), az(n);
|
||||
compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength,
|
||||
ax.data(), ay.data(), az.data());
|
||||
|
||||
// 半推速度：v_half = v + 0.5*a*dt (存入 vx, vy, vz)
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
vx[i] += ax[i] * dt * 0.5;
|
||||
vy[i] += ay[i] * dt * 0.5;
|
||||
vz[i] += az[i] * dt * 0.5;
|
||||
}
|
||||
|
||||
// 全推位置（不含边界，边界在外层统一处理）
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
x[i] += vx[i] * dt; // vx 此时是 v_half
|
||||
y[i] += vy[i] * dt;
|
||||
z[i] += vz[i] * dt;
|
||||
}
|
||||
|
||||
// 显式预测器：v_pred = v_half + 0.5*a_old*dt，用第一次加速度外推半步
|
||||
// 包含所有力的贡献（标准 Velocity-Verlet 预测步）
|
||||
std::vector<double> pred_vx(n), pred_vy(n), pred_vz(n);
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
pred_vx[i] = vx[i] + 0.5 * ax[i] * dt;
|
||||
pred_vy[i] = vy[i] + 0.5 * ay[i] * dt;
|
||||
pred_vz[i] = vz[i] + 0.5 * az[i] * dt;
|
||||
}
|
||||
|
||||
// 用新位置 + 预测速度重算加速度
|
||||
compute_acceleration(n, x, y, z, pred_vx.data(), pred_vy.data(), pred_vz.data(),
|
||||
m, G, B, bonds,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength,
|
||||
ax.data(), ay.data(), az.data());
|
||||
|
||||
// 速度后半步：v = v_half + 0.5*a_next*dt
|
||||
// vx 仍为 v_half（未被覆盖），直接加上 0.5*a_next*dt
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
vx[i] += ax[i] * dt * 0.5;
|
||||
vy[i] += ay[i] * dt * 0.5;
|
||||
vz[i] += az[i] * dt * 0.5;
|
||||
}
|
||||
}
|
||||
|
||||
/* ── 分发器：调用对应积分方法 + 边界条件（与 Python apply_motion_update 一致）── */
|
||||
static void apply_step(
|
||||
const std::string &method,
|
||||
int n, double *x, double *y, double *z,
|
||||
double *vx, double *vy, double *vz,
|
||||
const double *m, const double G[3], const double B[3],
|
||||
const BondData &bonds, const int *fixed,
|
||||
const double *pos_0,
|
||||
double box_a, double dt,
|
||||
int gravity_field, int gravity_interaction,
|
||||
int elastic_force, int damping_force,
|
||||
double gravity_strength)
|
||||
{
|
||||
// 积分
|
||||
if (method == "explicit_euler") {
|
||||
explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength);
|
||||
} else if (method == "implicit_euler") {
|
||||
implicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength);
|
||||
} else if (method == "midpoint") {
|
||||
midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength);
|
||||
} else if (method == "leapfrog") {
|
||||
leapfrog_full_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt,
|
||||
gravity_field, gravity_interaction,
|
||||
elastic_force, damping_force, gravity_strength);
|
||||
} else {
|
||||
die("未知算法: " + method);
|
||||
}
|
||||
|
||||
// 边界条件（与 Python Limit_in_box 一致）
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue;
|
||||
limit_in_box(x[i], vx[i], -box_a, box_a);
|
||||
limit_in_box(y[i], vy[i], -box_a, box_a);
|
||||
limit_in_box(z[i], vz[i], -box_a, box_a);
|
||||
}
|
||||
|
||||
// 逐自由度固定约束（与 Python apply_fixed_constraints 一致）
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (fixed[i*3+0]) { x[i] = pos_0[i*3]; vx[i] = 0.0; }
|
||||
if (fixed[i*3+1]) { y[i] = pos_0[i*3+1]; vy[i] = 0.0; }
|
||||
if (fixed[i*3+2]) { z[i] = pos_0[i*3+2]; vz[i] = 0.0; }
|
||||
}
|
||||
}
|
||||
|
||||
// ========================================================================
|
||||
// JSON 输出
|
||||
// ========================================================================
|
||||
|
||||
static void write_trajectory_json(
|
||||
const std::string &path,
|
||||
const std::vector<double> &x, const std::vector<double> &y,
|
||||
const std::vector<double> &z, const std::vector<double> &vx,
|
||||
const std::vector<double> &vy, const std::vector<double> &vz,
|
||||
int n_steps, int n_atoms,
|
||||
const SimParams ¶ms, const AtomData &atoms, const BondData &bonds)
|
||||
{
|
||||
std::ofstream f(path);
|
||||
if (!f) die("无法写入 " + path);
|
||||
f << std::setprecision(15);
|
||||
|
||||
f << "{\n";
|
||||
|
||||
// 轨迹数组
|
||||
const std::string names[] = {"traj_x","traj_y","traj_z","traj_vx","traj_vy","traj_vz"};
|
||||
const std::vector<double> *arrs[] = {&x, &y, &z, &vx, &vy, &vz};
|
||||
|
||||
for (int a = 0; a < 6; a++) {
|
||||
f << " \"" << names[a] << "\": [\n";
|
||||
const auto &data = *arrs[a];
|
||||
for (int t = 0; t < n_steps; t++) {
|
||||
f << " [";
|
||||
for (int i = 0; i < n_atoms; i++) {
|
||||
f << data[t * n_atoms + i];
|
||||
if (i < n_atoms - 1) f << ',';
|
||||
}
|
||||
f << ']';
|
||||
if (t < n_steps - 1) f << ',';
|
||||
f << '\n';
|
||||
}
|
||||
f << " ],\n";
|
||||
}
|
||||
|
||||
// 标量参数
|
||||
f << " \"NT\": " << params.NT << ",\n";
|
||||
f << " \"DT\": " << params.DT << ",\n";
|
||||
f << " \"NSTEP\": " << params.NSTEP << ",\n";
|
||||
f << " \"method\": \"" << params.method << "\",\n";
|
||||
f << " \"warmup_steps\": " << params.warmup_steps << ",\n";
|
||||
f << " \"G\": [" << params.G[0] << ", " << params.G[1] << ", " << params.G[2] << "],\n";
|
||||
f << " \"B\": [" << params.B[0] << ", " << params.B[1] << ", " << params.B[2] << "],\n";
|
||||
|
||||
// 原子信息
|
||||
f << " \"atom_ids\": [";
|
||||
for (size_t i = 0; i < atoms.ids.size(); i++) {
|
||||
if (i > 0) f << ',';
|
||||
f << atoms.ids[i];
|
||||
}
|
||||
f << "],\n";
|
||||
|
||||
f << " \"atom_masses\": [";
|
||||
for (size_t i = 0; i < atoms.masses.size(); i++) {
|
||||
if (i > 0) f << ',';
|
||||
f << atoms.masses[i];
|
||||
}
|
||||
f << "],\n";
|
||||
|
||||
// 成键
|
||||
f << " \"bond_pairs\": [";
|
||||
for (size_t b = 0; b < bonds.stiffness.size(); b++) {
|
||||
if (b > 0) f << ',';
|
||||
f << "[" << bonds.pairs[b * 2] << ", " << bonds.pairs[b * 2 + 1] << "]";
|
||||
}
|
||||
f << "],\n";
|
||||
|
||||
f << " \"bond_stiffness\": [";
|
||||
for (size_t b = 0; b < bonds.stiffness.size(); b++) {
|
||||
if (b > 0) f << ',';
|
||||
f << bonds.stiffness[b];
|
||||
}
|
||||
f << "],\n";
|
||||
|
||||
f << " \"bond_rest_lengths\": [";
|
||||
for (size_t b = 0; b < bonds.rest_lengths.size(); b++) {
|
||||
if (b > 0) f << ',';
|
||||
f << bonds.rest_lengths[b];
|
||||
}
|
||||
f << "]\n";
|
||||
|
||||
f << "}\n";
|
||||
}
|
||||
|
||||
// ========================================================================
|
||||
// 主函数
|
||||
// ========================================================================
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
@@ -63,15 +642,80 @@ int main(int argc, char **argv) {
|
||||
|
||||
auto t0 = std::chrono::high_resolution_clock::now();
|
||||
|
||||
// TODO: 读取 param.json、coord.txt、connection.txt、bond.txt
|
||||
// TODO: 执行 Leapfrog 模拟
|
||||
// TODO: 输出 trajectory.txt (JSON 格式)
|
||||
//
|
||||
// 提示：输出格式与 Python 版兼容，参考 engines/c/main.c 中的 write_trajectory_json
|
||||
// 读取参数和输入
|
||||
SimParams params = read_params(param_path);
|
||||
AtomData atoms = read_coord(input_dir);
|
||||
BondData bonds = read_bonds(input_dir);
|
||||
|
||||
std::cout << "[C++-engine] 原子数=" << atoms.ids.size()
|
||||
<< ", 键数=" << bonds.stiffness.size()
|
||||
<< ", NT=" << params.NT << ", DT=" << params.DT
|
||||
<< ", method=" << params.method << std::endl;
|
||||
|
||||
int n = static_cast<int>(atoms.ids.size());
|
||||
|
||||
// 初始化位置/速度
|
||||
std::vector<double> x(n), y(n), z(n), vx(n), vy(n), vz(n);
|
||||
for (int i = 0; i < n; i++) {
|
||||
x[i] = atoms.pos_0[i * 3];
|
||||
y[i] = atoms.pos_0[i * 3 + 1];
|
||||
z[i] = atoms.pos_0[i * 3 + 2];
|
||||
vx[i] = atoms.vel_0[i * 3];
|
||||
vy[i] = atoms.vel_0[i * 3 + 1];
|
||||
vz[i] = atoms.vel_0[i * 3 + 2];
|
||||
}
|
||||
|
||||
// 分配轨迹缓冲区
|
||||
int record_steps = params.NT - params.warmup_steps;
|
||||
std::vector<double> traj_x(record_steps * n);
|
||||
std::vector<double> traj_y(record_steps * n);
|
||||
std::vector<double> traj_z(record_steps * n);
|
||||
std::vector<double> traj_vx(record_steps * n);
|
||||
std::vector<double> traj_vy(record_steps * n);
|
||||
std::vector<double> traj_vz(record_steps * n);
|
||||
|
||||
// 预热
|
||||
for (int s = 0; s < params.warmup_steps; s++) {
|
||||
apply_step(params.method, n, x.data(), y.data(), z.data(),
|
||||
vx.data(), vy.data(), vz.data(),
|
||||
atoms.masses.data(), params.G, params.B,
|
||||
bonds, atoms.fixed.data(),
|
||||
atoms.pos_0.data(),
|
||||
params.box_a, params.DT,
|
||||
params.gravity_field, params.gravity_interaction,
|
||||
params.elastic_force, params.damping_force, params.gravity_strength);
|
||||
}
|
||||
|
||||
// 记录
|
||||
for (int s = 0; s < record_steps; s++) {
|
||||
// 保存当前帧
|
||||
for (int i = 0; i < n; i++) {
|
||||
traj_x[s * n + i] = x[i];
|
||||
traj_y[s * n + i] = y[i];
|
||||
traj_z[s * n + i] = z[i];
|
||||
traj_vx[s * n + i] = vx[i];
|
||||
traj_vy[s * n + i] = vy[i];
|
||||
traj_vz[s * n + i] = vz[i];
|
||||
}
|
||||
|
||||
apply_step(params.method, n, x.data(), y.data(), z.data(),
|
||||
vx.data(), vy.data(), vz.data(),
|
||||
atoms.masses.data(), params.G, params.B,
|
||||
bonds, atoms.fixed.data(),
|
||||
atoms.pos_0.data(),
|
||||
params.box_a, params.DT,
|
||||
params.gravity_field, params.gravity_interaction,
|
||||
params.elastic_force, params.damping_force, params.gravity_strength);
|
||||
}
|
||||
|
||||
// 输出轨迹
|
||||
std::string out_path = output_dir + "/trajectory.txt";
|
||||
write_trajectory_json(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz,
|
||||
record_steps, n, params, atoms, bonds);
|
||||
|
||||
auto t1 = std::chrono::high_resolution_clock::now();
|
||||
double elapsed = std::chrono::duration<double>(t1 - t0).count();
|
||||
std::cout << "[C++-engine] 计算完成: " << record_steps << " 步, " << elapsed << " s" << std::endl;
|
||||
|
||||
std::cout << "[C++-engine] 占位: NT=" << 0 << ", " << elapsed << " s (待实现)" << std::endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
+841
-20
@@ -1,22 +1,49 @@
|
||||
! engines/fortran/main.f90
|
||||
! ------------------------
|
||||
! Fortran 动力学模拟引擎（模板/参考实现）。
|
||||
! Fortran 动力学模拟引擎。
|
||||
! 与 Python 版 (compute.py) 算法保持一致。
|
||||
!
|
||||
! 接口与 C 引擎一致：
|
||||
! 输入: param.json, <input_dir>/coord.txt, connection.txt, bond.txt
|
||||
! 输出: <output_dir>/trajectory.txt (JSON)
|
||||
! 输入: param.json, <input_dir>/coord.txt, connection.txt, bond.txt
|
||||
! 输出: <output_dir>/trajectory.txt (JSON)
|
||||
!
|
||||
! 编译:
|
||||
! gfortran -O3 -march=native -o build/dynamics_f90 main.f90
|
||||
!
|
||||
! 用法:
|
||||
! ./build/dynamics_f90 <input_dir> <output_dir> <param_json>
|
||||
! 编译: cmake --build build --target dynamics_f90
|
||||
! 用法: ./build/dynamics_f90 <input_dir> <output_dir> <param_json>
|
||||
|
||||
program dynamics_f90
|
||||
implicit none
|
||||
character(len=256) :: input_dir, output_dir, param_path
|
||||
integer :: narg
|
||||
|
||||
character(len=256) :: input_dir, output_dir, param_path
|
||||
integer :: narg, i, n, s
|
||||
integer :: NT, NSTEP, warmup_steps, n_atoms, n_bonds
|
||||
double precision :: DT, box_a
|
||||
double precision :: G(3), B(3)
|
||||
integer :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
double precision :: gravity_strength
|
||||
character(len=32) :: method
|
||||
double precision :: t0, t1, elapsed
|
||||
|
||||
! 原子数据
|
||||
integer, allocatable :: atom_ids(:)
|
||||
double precision, allocatable :: masses(:), radii(:)
|
||||
double precision, allocatable :: pos_0(:, :), vel_0(:, :)
|
||||
integer, allocatable :: fixed(:, :)
|
||||
|
||||
! 成键数据
|
||||
integer, allocatable :: bond_pairs(:, :)
|
||||
double precision, allocatable :: bond_stiffness(:), bond_rest_lengths(:)
|
||||
|
||||
! 运行时位置/速度
|
||||
double precision, allocatable :: x(:), y(:), z(:)
|
||||
double precision, allocatable :: vx(:), vy(:), vz(:)
|
||||
|
||||
! 轨迹缓冲区
|
||||
integer :: record_steps
|
||||
double precision, allocatable :: traj_x(:, :), traj_y(:, :), traj_z(:, :)
|
||||
double precision, allocatable :: traj_vx(:, :), traj_vy(:, :), traj_vz(:, :)
|
||||
|
||||
! ========================================================================
|
||||
! 主流程
|
||||
! ========================================================================
|
||||
narg = command_argument_count()
|
||||
if (narg < 3) then
|
||||
write(*,*) "用法: dynamics_f90 <input_dir> <output_dir> <param_json>"
|
||||
@@ -27,14 +54,808 @@ program dynamics_f90
|
||||
call get_command_argument(2, output_dir)
|
||||
call get_command_argument(3, param_path)
|
||||
|
||||
! TODO: 实现与 C 引擎相同的物理模拟
|
||||
!
|
||||
! 参考实现:
|
||||
! 1. 读取 param.json → NT, DT, G, B, ...
|
||||
! 2. 读取 coord.txt → 原子数, 质量, 位置, 速度
|
||||
! 3. 读取 connection.txt + bond.txt → 成键参数
|
||||
! 4. Leapfrog 循环 NT 步
|
||||
! 5. 输出 trajectory.txt (JSON 格式)
|
||||
call cpu_time(t0)
|
||||
|
||||
! 读取 param.json
|
||||
call read_params(param_path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
|
||||
! 读取 coord.txt
|
||||
call read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed)
|
||||
|
||||
! 读取成键信息
|
||||
call read_bonds(input_dir, n_atoms, atom_ids, pos_0, n_bonds, &
|
||||
bond_pairs, bond_stiffness, bond_rest_lengths)
|
||||
|
||||
write(*, '("[Fortran-engine] 原子数=", i0, ", 键数=", i0, &
|
||||
&", NT=", i0, ", DT=", f0.6, ", method=", a)') &
|
||||
n_atoms, n_bonds, NT, DT, trim(method)
|
||||
|
||||
! 初始化位置/速度
|
||||
n = n_atoms
|
||||
allocate(x(n), y(n), z(n), vx(n), vy(n), vz(n))
|
||||
do i = 1, n
|
||||
x(i) = pos_0(i, 1); y(i) = pos_0(i, 2); z(i) = pos_0(i, 3)
|
||||
vx(i) = vel_0(i, 1); vy(i) = vel_0(i, 2); vz(i) = vel_0(i, 3)
|
||||
end do
|
||||
|
||||
! 分配轨迹缓冲区
|
||||
record_steps = NT - warmup_steps
|
||||
allocate(traj_x(record_steps, n), traj_y(record_steps, n), traj_z(record_steps, n))
|
||||
allocate(traj_vx(record_steps, n), traj_vy(record_steps, n), traj_vz(record_steps, n))
|
||||
|
||||
! 预热
|
||||
do s = 1, warmup_steps
|
||||
call apply_step(method, n, x, y, z, vx, vy, vz, masses, G, B, &
|
||||
n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths, &
|
||||
fixed, box_a, DT, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, &
|
||||
pos_0)
|
||||
end do
|
||||
|
||||
! 记录
|
||||
do s = 1, record_steps
|
||||
traj_x(s, :) = x; traj_y(s, :) = y; traj_z(s, :) = z
|
||||
traj_vx(s, :) = vx; traj_vy(s, :) = vy; traj_vz(s, :) = vz
|
||||
call apply_step(method, n, x, y, z, vx, vy, vz, masses, G, B, &
|
||||
n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths, &
|
||||
fixed, box_a, DT, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, &
|
||||
pos_0)
|
||||
end do
|
||||
|
||||
! 输出轨迹
|
||||
call write_json(output_dir, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, &
|
||||
record_steps, n_atoms, atom_ids, masses, &
|
||||
NT, DT, NSTEP, warmup_steps, method, G, B, &
|
||||
n_bonds, bond_pairs, bond_stiffness, bond_rest_lengths)
|
||||
|
||||
call cpu_time(t1)
|
||||
elapsed = t1 - t0
|
||||
write(*, '("[Fortran-engine] 计算完成: ", i0, " 步, ", f0.6, " s")') record_steps, elapsed
|
||||
|
||||
deallocate(x, y, z, vx, vy, vz)
|
||||
deallocate(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz)
|
||||
deallocate(atom_ids, masses, radii, pos_0, vel_0, fixed)
|
||||
if (allocated(bond_pairs)) deallocate(bond_pairs, bond_stiffness, bond_rest_lengths)
|
||||
|
||||
contains
|
||||
|
||||
! ========================================================================
|
||||
! 读取 param.json
|
||||
! ========================================================================
|
||||
subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
character(len=*), intent(in) :: path
|
||||
double precision, intent(out) :: box_a, DT, G(3), B(3), gravity_strength
|
||||
integer, intent(out) :: NT, NSTEP, warmup_steps
|
||||
integer, intent(out) :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
character(len=*), intent(out) :: method
|
||||
character(len=8096) :: buf
|
||||
character(len=256) :: line
|
||||
integer :: u, ios
|
||||
|
||||
box_a = 10.0d0; NT = 10000; DT = 0.001d0; NSTEP = 100; warmup_steps = 0
|
||||
method = 'leapfrog'
|
||||
G = (/ 0.0d0, 0.0d0, -9.8d0 /); B = 0.0d0
|
||||
gravity_field = 1; gravity_interaction = 0
|
||||
elastic_force = 1; damping_force = 0; gravity_strength = 1.0d0
|
||||
|
||||
open(newunit=u, file=trim(path), status='old', action='read', iostat=ios)
|
||||
if (ios /= 0) then
|
||||
write(*, '("[Fortran-engine] 警告: 无法打开 ", a, ", 使用默认值")') trim(path)
|
||||
return
|
||||
end if
|
||||
buf = ''
|
||||
do
|
||||
read(u, '(a)', iostat=ios) line
|
||||
if (ios /= 0) exit
|
||||
buf = trim(buf) // trim(line) // char(10)
|
||||
end do
|
||||
close(u)
|
||||
if (ios > 0) return
|
||||
|
||||
box_a = json_get_double(buf, 'box_a', box_a)
|
||||
NT = json_get_int(buf, 'NT', NT)
|
||||
DT = json_get_double(buf, 'DT', DT)
|
||||
NSTEP = json_get_int(buf, 'NSTEP', NSTEP)
|
||||
warmup_steps = json_get_int(buf, 'warmup_steps', warmup_steps)
|
||||
call json_get_string(buf, 'method', method)
|
||||
if (len_trim(method) == 0) method = 'leapfrog'
|
||||
call json_get_double3(buf, 'G', G)
|
||||
call json_get_double3(buf, 'B', B)
|
||||
gravity_field = json_get_int(buf, 'gravity_field', 1)
|
||||
gravity_interaction = json_get_int(buf, 'gravity_interaction', 0)
|
||||
elastic_force = json_get_int(buf, 'elastic_force', 1)
|
||||
damping_force = json_get_int(buf, 'damping_force', 0)
|
||||
gravity_strength = json_get_double(buf, 'gravity_strength', 1.0d0)
|
||||
end subroutine read_params
|
||||
|
||||
! ========================================================================
|
||||
! JSON 辅助解析
|
||||
! ========================================================================
|
||||
function json_get_double(buf, key, default) result(val)
|
||||
character(len=*), intent(in) :: buf, key
|
||||
double precision, intent(in) :: default
|
||||
double precision :: val
|
||||
integer :: p, ios
|
||||
val = default
|
||||
p = index(buf, '"' // key // '"')
|
||||
if (p == 0) return
|
||||
p = index(buf(p:), ':') + p
|
||||
do while (p <= len(buf) .and. (buf(p:p) == ':' .or. buf(p:p) == ' ' &
|
||||
.or. buf(p:p) == char(9) .or. buf(p:p) == char(10)))
|
||||
p = p + 1
|
||||
end do
|
||||
read(buf(p:), *, iostat=ios) val
|
||||
end function json_get_double
|
||||
|
||||
function json_get_int(buf, key, default) result(val)
|
||||
character(len=*), intent(in) :: buf, key
|
||||
integer, intent(in) :: default
|
||||
integer :: val
|
||||
val = nint(json_get_double(buf, key, dble(default)))
|
||||
end function json_get_int
|
||||
|
||||
subroutine json_get_double3(buf, key, vals)
|
||||
character(len=*), intent(in) :: buf, key
|
||||
double precision, intent(out) :: vals(3)
|
||||
integer :: p, i, ios
|
||||
vals = 0.0d0
|
||||
p = index(buf, '"' // key // '"')
|
||||
if (p == 0) return
|
||||
p = index(buf(p:), '[') + p
|
||||
do i = 1, 3
|
||||
do while (buf(p:p) == ' ' .or. buf(p:p) == ',' .or. &
|
||||
buf(p:p) == char(9) .or. buf(p:p) == char(10))
|
||||
p = p + 1
|
||||
end do
|
||||
read(buf(p:), *, iostat=ios) vals(i)
|
||||
if (ios /= 0) exit
|
||||
do while (p <= len(buf) .and. buf(p:p) /= ',' .and. buf(p:p) /= ']')
|
||||
p = p + 1
|
||||
end do
|
||||
end do
|
||||
end subroutine json_get_double3
|
||||
|
||||
! 从 JSON 中读取字符串值
|
||||
subroutine json_get_string(buf, key, val)
|
||||
character(len=*), intent(in) :: buf, key
|
||||
character(len=*), intent(out) :: val
|
||||
integer :: p, i, j
|
||||
val = ''
|
||||
p = index(buf, '"' // key // '"')
|
||||
if (p == 0) return
|
||||
p = index(buf(p:), ':') + p
|
||||
do while (p <= len(buf) .and. (buf(p:p) == ':' .or. buf(p:p) == ' ' &
|
||||
.or. buf(p:p) == char(9) .or. buf(p:p) == char(10)))
|
||||
p = p + 1
|
||||
end do
|
||||
if (buf(p:p) /= '"') return
|
||||
p = p + 1
|
||||
i = 1
|
||||
do while (p <= len(buf) .and. buf(p:p) /= '"' .and. i <= len(val))
|
||||
val(i:i) = buf(p:p)
|
||||
i = i + 1; p = p + 1
|
||||
end do
|
||||
end subroutine json_get_string
|
||||
|
||||
! ========================================================================
|
||||
! 读取 coord.txt
|
||||
! ========================================================================
|
||||
subroutine read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed)
|
||||
character(len=*), intent(in) :: input_dir
|
||||
integer, intent(out) :: n_atoms
|
||||
integer, allocatable, intent(out) :: atom_ids(:)
|
||||
double precision, allocatable, intent(out) :: masses(:), radii(:)
|
||||
double precision, allocatable, intent(out) :: pos_0(:, :), vel_0(:, :)
|
||||
integer, allocatable, intent(out) :: fixed(:, :)
|
||||
|
||||
character(len=512) :: path, line
|
||||
integer :: u, ios, ncols, n_cap, id, fx, fy, fz, i
|
||||
double precision :: mass, rad, px, py, pz, vx, vy, vz
|
||||
integer, parameter :: MX = 4096
|
||||
integer :: ids_tmp(MX), fix_tmp(MX, 3)
|
||||
double precision :: m_tmp(MX), r_tmp(MX), p_tmp(MX, 3), v_tmp(MX, 3)
|
||||
|
||||
n_atoms = 0
|
||||
path = trim(input_dir) // '/coord.txt'
|
||||
open(newunit=u, file=trim(path), status='old', action='read', iostat=ios)
|
||||
if (ios /= 0) then; write(*, '("[Fortran-engine] 错误: 无法打开 ", a)') trim(path); stop; end if
|
||||
|
||||
read(u, '(a)', iostat=ios) line
|
||||
ncols = 0
|
||||
do i = 1, len_trim(line)
|
||||
if (line(i:i) /= ' ' .and. (i == 1 .or. line(i-1:i-1) == ' ')) ncols = ncols + 1
|
||||
end do
|
||||
|
||||
do
|
||||
read(u, '(a)', iostat=ios) line
|
||||
if (ios /= 0) exit
|
||||
if (len_trim(line) == 0 .or. line(1:1) == '#') cycle
|
||||
n_atoms = n_atoms + 1
|
||||
if (n_atoms > MX) exit
|
||||
if (ncols >= 12) then
|
||||
read(line, *) id, mass, rad, px, py, pz, vx, vy, vz, fx, fy, fz
|
||||
else
|
||||
read(line, *) id, mass, rad, px, py, pz, vx, vy, vz
|
||||
fx = 0; fy = 0; fz = 0
|
||||
end if
|
||||
ids_tmp(n_atoms) = id
|
||||
m_tmp(n_atoms) = mass; r_tmp(n_atoms) = rad
|
||||
p_tmp(n_atoms, :) = (/ px, py, pz /)
|
||||
v_tmp(n_atoms, :) = (/ vx, vy, vz /)
|
||||
fix_tmp(n_atoms, :) = (/ fx, fy, fz /)
|
||||
end do
|
||||
close(u)
|
||||
if (n_atoms <= 0) then; write(*, '("[Fortran-engine] 错误: coord.txt 为空")') ; stop; end if
|
||||
|
||||
allocate(atom_ids(n_atoms), masses(n_atoms), radii(n_atoms))
|
||||
allocate(pos_0(n_atoms, 3), vel_0(n_atoms, 3), fixed(n_atoms, 3))
|
||||
atom_ids = ids_tmp(1:n_atoms)
|
||||
masses = m_tmp(1:n_atoms); radii = r_tmp(1:n_atoms)
|
||||
pos_0 = p_tmp(1:n_atoms, :); vel_0 = v_tmp(1:n_atoms, :)
|
||||
fixed = fix_tmp(1:n_atoms, :)
|
||||
end subroutine read_coord
|
||||
|
||||
! ========================================================================
|
||||
! 读取成键信息 (connection.txt + bond.txt)
|
||||
! ========================================================================
|
||||
subroutine read_bonds(input_dir, n_atoms, atom_ids, pos_0, n_bonds, &
|
||||
bond_pairs, bond_stiffness, bond_rest_lengths)
|
||||
character(len=*), intent(in) :: input_dir
|
||||
integer, intent(in) :: n_atoms, atom_ids(n_atoms)
|
||||
double precision, intent(in) :: pos_0(n_atoms, 3)
|
||||
integer, intent(out) :: n_bonds
|
||||
integer, allocatable, intent(out) :: bond_pairs(:, :)
|
||||
double precision, allocatable, intent(out) :: bond_stiffness(:), bond_rest_lengths(:)
|
||||
|
||||
character(len=512) :: conn_path
|
||||
integer :: u, ios, a1, a2, idx1, idx2, i
|
||||
integer, parameter :: MX = 4096
|
||||
integer :: ptmp(MX, 2), idx_map(9999)
|
||||
double precision :: ktmp(MX), r0tmp(MX), dx, dy, dz
|
||||
character(len=64) :: name
|
||||
|
||||
n_bonds = 0; idx_map = -1
|
||||
do i = 1, n_atoms
|
||||
if (atom_ids(i) > 0 .and. atom_ids(i) <= 9998) idx_map(atom_ids(i)) = i
|
||||
end do
|
||||
|
||||
conn_path = trim(input_dir) // '/connection.txt'
|
||||
open(newunit=u, file=trim(conn_path), status='old', action='read', iostat=ios)
|
||||
if (ios /= 0) return
|
||||
|
||||
read(u, *, iostat=ios)
|
||||
do
|
||||
read(u, *, iostat=ios) a1, a2, name
|
||||
if (ios /= 0) exit
|
||||
call find_bond(trim(input_dir), trim(name), ktmp(n_bonds+1), r0tmp(n_bonds+1))
|
||||
idx1 = idx_map(a1); idx2 = idx_map(a2)
|
||||
if (idx1 < 1 .or. idx2 < 1) cycle
|
||||
if (r0tmp(n_bonds+1) <= 0.0d0) then
|
||||
dx = pos_0(idx2, 1) - pos_0(idx1, 1)
|
||||
dy = pos_0(idx2, 2) - pos_0(idx1, 2)
|
||||
dz = pos_0(idx2, 3) - pos_0(idx1, 3)
|
||||
r0tmp(n_bonds+1) = sqrt(dx*dx + dy*dy + dz*dz)
|
||||
end if
|
||||
n_bonds = n_bonds + 1
|
||||
ptmp(n_bonds, :) = (/ idx1 - 1, idx2 - 1 /)
|
||||
if (n_bonds >= MX) exit
|
||||
end do
|
||||
close(u)
|
||||
|
||||
if (n_bonds <= 0) return
|
||||
allocate(bond_pairs(n_bonds, 2), bond_stiffness(n_bonds), bond_rest_lengths(n_bonds))
|
||||
bond_pairs(:, :) = ptmp(1:n_bonds, :)
|
||||
bond_stiffness = ktmp(1:n_bonds)
|
||||
bond_rest_lengths = r0tmp(1:n_bonds)
|
||||
end subroutine read_bonds
|
||||
|
||||
subroutine find_bond(bond_dir, bond_name, k, r0)
|
||||
character(len=*), intent(in) :: bond_dir, bond_name
|
||||
double precision, intent(out) :: k, r0
|
||||
character(len=512) :: path, bn
|
||||
integer :: u, ios
|
||||
double precision :: bk, br
|
||||
k = 1.0d0; r0 = -1.0d0
|
||||
path = trim(bond_dir) // '/bond.txt'
|
||||
open(newunit=u, file=trim(path), status='old', action='read', iostat=ios)
|
||||
if (ios /= 0) return
|
||||
read(u, *, iostat=ios)
|
||||
do
|
||||
read(u, *, iostat=ios) bn, bk, br
|
||||
if (ios /= 0) exit
|
||||
if (trim(bn) == trim(bond_name)) then
|
||||
k = bk
|
||||
if (br > 0.0d0) r0 = br
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
close(u)
|
||||
end subroutine find_bond
|
||||
|
||||
! ========================================================================
|
||||
! 物理核心
|
||||
! ========================================================================
|
||||
|
||||
! 加速度计算（各力独立开关控制）
|
||||
pure subroutine accel(n, x, y, z, vx, vy, vz, m, G, B, &
|
||||
nb, bp, bk, br, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, &
|
||||
ax, ay, az)
|
||||
integer, intent(in) :: n, nb, bp(nb, 2)
|
||||
integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
double precision, intent(in) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
|
||||
double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), gravity_strength
|
||||
double precision, intent(out) :: ax(n), ay(n), az(n)
|
||||
integer :: i, j, ib
|
||||
double precision :: dx, dy, dz, dist, s, fm, ux, uy, uz, r2, fg
|
||||
|
||||
! 清零
|
||||
ax = 0.0d0; ay = 0.0d0; az = 0.0d0
|
||||
|
||||
! 均匀重力场
|
||||
if (gravity_field /= 0) then
|
||||
do i = 1, n
|
||||
ax(i) = G(1); ay(i) = G(2); az(i) = G(3)
|
||||
end do
|
||||
end if
|
||||
|
||||
! 阻尼
|
||||
if (damping_force /= 0) then
|
||||
do i = 1, n
|
||||
ax(i) = ax(i) - B(1) * vx(i) / m(i)
|
||||
ay(i) = ay(i) - B(2) * vy(i) / m(i)
|
||||
az(i) = az(i) - B(3) * vz(i) / m(i)
|
||||
end do
|
||||
end if
|
||||
|
||||
! 弹簧键力
|
||||
if (elastic_force /= 0) then
|
||||
do ib = 1, nb
|
||||
i = bp(ib, 1) + 1; j = bp(ib, 2) + 1
|
||||
dx = x(j) - x(i); dy = y(j) - y(i); dz = z(j) - z(i)
|
||||
dist = sqrt(dx*dx + dy*dy + dz*dz)
|
||||
if (dist < 1.0d-12) cycle
|
||||
s = (dist - br(ib)) / dist
|
||||
fm = bk(ib) * s
|
||||
ux = fm * dx; uy = fm * dy; uz = fm * dz
|
||||
ax(i) = ax(i) + ux / m(i); ay(i) = ay(i) + uy / m(i); az(i) = az(i) + uz / m(i)
|
||||
ax(j) = ax(j) - ux / m(j); ay(j) = ay(j) - uy / m(j); az(j) = az(j) - uz / m(j)
|
||||
end do
|
||||
end if
|
||||
|
||||
! 万有引力（所有原子对之间）
|
||||
if (gravity_interaction /= 0) then
|
||||
do i = 1, n - 1
|
||||
do j = i + 1, n
|
||||
dx = x(j) - x(i); dy = y(j) - y(i); dz = z(j) - z(i)
|
||||
r2 = dx*dx + dy*dy + dz*dz
|
||||
if (r2 <= 1.0d-12) cycle
|
||||
dist = sqrt(r2)
|
||||
fg = gravity_strength * m(i) * m(j) / r2
|
||||
ux = fg * dx / dist; uy = fg * dy / dist; uz = fg * dz / dist
|
||||
ax(i) = ax(i) + ux / m(i); ay(i) = ay(i) + uy / m(i); az(i) = az(i) + uz / m(i)
|
||||
ax(j) = ax(j) - ux / m(j); ay(j) = ay(j) - uy / m(j); az(j) = az(j) - uz / m(j)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
end subroutine accel
|
||||
|
||||
! 边界条件：clamp 位置 + 速度反转
|
||||
subroutine limit_in_box(pos, vel, lo, hi)
|
||||
double precision, intent(inout) :: pos, vel
|
||||
double precision, intent(in) :: lo, hi
|
||||
if (pos > hi) then
|
||||
pos = hi; vel = -vel
|
||||
else if (pos < lo) then
|
||||
pos = lo; vel = -vel
|
||||
end if
|
||||
end subroutine limit_in_box
|
||||
|
||||
! ── 显式欧拉法 ────────────
|
||||
subroutine explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3)
|
||||
integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
|
||||
double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), dt, gravity_strength
|
||||
double precision :: ax(n), ay(n), az(n)
|
||||
integer :: i
|
||||
|
||||
call accel(n, x, y, z, vx, vy, vz, m, G, B, nb, bp, bk, br, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, ax, ay, az)
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
x(i) = x(i) + vx(i) * dt
|
||||
y(i) = y(i) + vy(i) * dt
|
||||
z(i) = z(i) + vz(i) * dt
|
||||
vx(i) = vx(i) + ax(i) * dt
|
||||
vy(i) = vy(i) + ay(i) * dt
|
||||
vz(i) = vz(i) + az(i) * dt
|
||||
end do
|
||||
end subroutine explicit_euler_step
|
||||
|
||||
! ── 隐式欧拉法 ────────────
|
||||
subroutine implicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3)
|
||||
integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
|
||||
double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), dt, gravity_strength
|
||||
double precision :: ax(n), ay(n), az(n)
|
||||
double precision :: vxn(n), vyn(n), vzn(n)
|
||||
double precision :: gx, gy, gz
|
||||
integer :: i
|
||||
|
||||
! 隐式速度（重力 + 阻尼）
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then
|
||||
vxn(i) = 0; vyn(i) = 0; vzn(i) = 0; cycle
|
||||
end if
|
||||
gx = B(1) / m(i); gy = B(2) / m(i); gz = B(3) / m(i)
|
||||
vxn(i) = (vx(i) + G(1) * dt) / (1.0d0 + gx * dt)
|
||||
vyn(i) = (vy(i) + G(2) * dt) / (1.0d0 + gy * dt)
|
||||
vzn(i) = (vz(i) + G(3) * dt) / (1.0d0 + gz * dt)
|
||||
end do
|
||||
|
||||
! 用当前位 + 隐式速度算完整加速度
|
||||
call accel(n, x, y, z, vxn, vyn, vzn, m, G, B, nb, bp, bk, br, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, ax, ay, az)
|
||||
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
vx(i) = vx(i) + ax(i) * dt
|
||||
vy(i) = vy(i) + ay(i) * dt
|
||||
vz(i) = vz(i) + az(i) * dt
|
||||
x(i) = x(i) + vx(i) * dt
|
||||
y(i) = y(i) + vy(i) * dt
|
||||
z(i) = z(i) + vz(i) * dt
|
||||
end do
|
||||
end subroutine implicit_euler_step
|
||||
|
||||
! ── 中点法 ────────────
|
||||
subroutine midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3)
|
||||
integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
|
||||
double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), dt, gravity_strength
|
||||
double precision :: ax(n), ay(n), az(n)
|
||||
double precision :: xm(n), ym(n), zm(n), vxm(n), vym(n), vzm(n)
|
||||
integer :: i
|
||||
|
||||
call accel(n, x, y, z, vx, vy, vz, m, G, B, nb, bp, bk, br, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, ax, ay, az)
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then
|
||||
xm(i)=0; ym(i)=0; zm(i)=0; vxm(i)=0; vym(i)=0; vzm(i)=0; cycle
|
||||
end if
|
||||
xm(i) = x(i) + 0.5d0 * vx(i) * dt
|
||||
ym(i) = y(i) + 0.5d0 * vy(i) * dt
|
||||
zm(i) = z(i) + 0.5d0 * vz(i) * dt
|
||||
vxm(i) = vx(i) + 0.5d0 * ax(i) * dt
|
||||
vym(i) = vy(i) + 0.5d0 * ay(i) * dt
|
||||
vzm(i) = vz(i) + 0.5d0 * az(i) * dt
|
||||
x(i) = x(i) + vxm(i) * dt
|
||||
y(i) = y(i) + vym(i) * dt
|
||||
z(i) = z(i) + vzm(i) * dt
|
||||
end do
|
||||
|
||||
call accel(n, xm, ym, zm, vxm, vym, vzm, m, G, B, nb, bp, bk, br, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, ax, ay, az)
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
vx(i) = vx(i) + ax(i) * dt
|
||||
vy(i) = vy(i) + ay(i) * dt
|
||||
vz(i) = vz(i) + az(i) * dt
|
||||
end do
|
||||
end subroutine midpoint_step
|
||||
|
||||
! ── 蛙跳法（Velocity-Verlet）──
|
||||
subroutine leapfrog_full(n, x, y, z, vx, vy, vz, m, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3)
|
||||
integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
|
||||
double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), dt, gravity_strength
|
||||
double precision :: ax(n), ay(n), az(n)
|
||||
double precision :: dmp_vx(n), dmp_vy(n), dmp_vz(n)
|
||||
double precision :: gx, gy, gz
|
||||
integer :: i
|
||||
|
||||
call accel(n, x, y, z, vx, vy, vz, m, G, B, nb, bp, bk, br, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, ax, ay, az)
|
||||
|
||||
! 速度半步推
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
vx(i) = vx(i) + ax(i) * dt * 0.5d0
|
||||
vy(i) = vy(i) + ay(i) * dt * 0.5d0
|
||||
vz(i) = vz(i) + az(i) * dt * 0.5d0
|
||||
end do
|
||||
|
||||
! 全推位置（不含边界）
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
x(i) = x(i) + vx(i) * dt
|
||||
y(i) = y(i) + vy(i) * dt
|
||||
z(i) = z(i) + vz(i) * dt
|
||||
end do
|
||||
|
||||
! 隐式阻尼处理（用临时数组 dmp_v，不覆盖 vx/vy/vz）
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) then
|
||||
dmp_vx(i) = 0; dmp_vy(i) = 0; dmp_vz(i) = 0; cycle
|
||||
end if
|
||||
gx = B(1) / m(i); gy = B(2) / m(i); gz = B(3) / m(i)
|
||||
dmp_vx(i) = (vx(i) + 0.5d0 * G(1) * dt) / (1.0d0 + 0.5d0 * gx * dt)
|
||||
dmp_vy(i) = (vy(i) + 0.5d0 * G(2) * dt) / (1.0d0 + 0.5d0 * gy * dt)
|
||||
dmp_vz(i) = (vz(i) + 0.5d0 * G(3) * dt) / (1.0d0 + 0.5d0 * gz * dt)
|
||||
end do
|
||||
|
||||
! 用新位置 + 阻尼速度重算加速度
|
||||
call accel(n, x, y, z, dmp_vx, dmp_vy, dmp_vz, m, G, B, nb, bp, bk, br, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, ax, ay, az)
|
||||
|
||||
! 速度后半步：v = v_half + 0.5*a_next*dt（vx 仍为 v_half）
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
vx(i) = vx(i) + ax(i) * dt * 0.5d0
|
||||
vy(i) = vy(i) + ay(i) * dt * 0.5d0
|
||||
vz(i) = vz(i) + az(i) * dt * 0.5d0
|
||||
end do
|
||||
end subroutine leapfrog_full
|
||||
|
||||
! ── 分发器 + 边界条件 + 自由度约束 ──
|
||||
subroutine apply_step(method, n, x, y, z, vx, vy, vz, m, G, B, &
|
||||
nb, bp, bk, br, fixed, box_a, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength, &
|
||||
pos_0)
|
||||
character(len=*), intent(in) :: method
|
||||
integer, intent(in) :: n, nb, bp(nb, 2), fixed(n, 3)
|
||||
integer, intent(in) :: gravity_field, gravity_interaction, elastic_force, damping_force
|
||||
double precision, intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
|
||||
double precision, intent(in) :: m(n), G(3), B(3), bk(nb), br(nb), box_a, dt, gravity_strength
|
||||
double precision, intent(in) :: pos_0(n, 3)
|
||||
double precision :: mm(n)
|
||||
integer :: i
|
||||
|
||||
mm = m
|
||||
|
||||
if (trim(method) == 'explicit_euler') then
|
||||
call explicit_euler_step(n, x, y, z, vx, vy, vz, mm, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
else if (trim(method) == 'implicit_euler') then
|
||||
call implicit_euler_step(n, x, y, z, vx, vy, vz, mm, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
else if (trim(method) == 'midpoint') then
|
||||
call midpoint_step(n, x, y, z, vx, vy, vz, mm, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
else if (trim(method) == 'leapfrog') then
|
||||
call leapfrog_full(n, x, y, z, vx, vy, vz, mm, G, B, &
|
||||
nb, bp, bk, br, fixed, dt, &
|
||||
gravity_field, gravity_interaction, &
|
||||
elastic_force, damping_force, gravity_strength)
|
||||
else
|
||||
write(*, '("[Fortran-engine] 未知算法: ", a)') trim(method)
|
||||
stop
|
||||
end if
|
||||
|
||||
! 边界条件
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0 .and. fixed(i,2) /= 0 .and. fixed(i,3) /= 0) cycle
|
||||
call limit_in_box(x(i), vx(i), -box_a, box_a)
|
||||
call limit_in_box(y(i), vy(i), -box_a, box_a)
|
||||
call limit_in_box(z(i), vz(i), -box_a, box_a)
|
||||
end do
|
||||
|
||||
! 逐自由度固定约束（与 Python apply_fixed_constraints 一致）
|
||||
do i = 1, n
|
||||
if (fixed(i,1) /= 0) then; x(i) = pos_0(i, 1); vx(i) = 0.0d0; end if
|
||||
if (fixed(i,2) /= 0) then; y(i) = pos_0(i, 2); vy(i) = 0.0d0; end if
|
||||
if (fixed(i,3) /= 0) then; z(i) = pos_0(i, 3); vz(i) = 0.0d0; end if
|
||||
end do
|
||||
end subroutine apply_step
|
||||
|
||||
! ========================================================================
|
||||
! JSON 输出
|
||||
! ========================================================================
|
||||
|
||||
subroutine write_json(outdir, tx, ty, tz, tvx, tvy, tvz, &
|
||||
nsteps, nat, aid, amass, &
|
||||
NT, DT, NSTEP, warmup, method, G, B, &
|
||||
nb, bp, bk, br)
|
||||
character(len=*), intent(in) :: outdir, method
|
||||
integer, intent(in) :: nsteps, nat, NT, NSTEP, warmup, nb, bp(nb, 2), aid(nat)
|
||||
double precision, intent(in) :: tx(nsteps, nat), ty(nsteps, nat), tz(nsteps, nat)
|
||||
double precision, intent(in) :: tvx(nsteps, nat), tvy(nsteps, nat), tvz(nsteps, nat)
|
||||
double precision, intent(in) :: DT, G(3), B(3), bk(nb), br(nb), amass(nat)
|
||||
|
||||
character(len=512) :: path, buf
|
||||
integer :: u, s, i, ib, ios
|
||||
|
||||
path = trim(outdir) // '/trajectory.txt'
|
||||
open(newunit=u, file=trim(path), status='replace', action='write', iostat=ios)
|
||||
if (ios /= 0) then
|
||||
write(*, '("[Fortran-engine] 错误: 无法写入 ", a)') trim(path)
|
||||
stop
|
||||
end if
|
||||
|
||||
write(u, '(a)') '{'
|
||||
|
||||
! traj_x
|
||||
write(u, '(a)') ' "traj_x": ['
|
||||
do s = 1, nsteps
|
||||
call json_arr(u, tx(s, :), nat, s < nsteps, ' ')
|
||||
end do
|
||||
write(u, '(a)') ' ],'
|
||||
|
||||
! traj_y
|
||||
write(u, '(a)') ' "traj_y": ['
|
||||
do s = 1, nsteps
|
||||
call json_arr(u, ty(s, :), nat, s < nsteps, ' ')
|
||||
end do
|
||||
write(u, '(a)') ' ],'
|
||||
|
||||
! traj_z
|
||||
write(u, '(a)') ' "traj_z": ['
|
||||
do s = 1, nsteps
|
||||
call json_arr(u, tz(s, :), nat, s < nsteps, ' ')
|
||||
end do
|
||||
write(u, '(a)') ' ],'
|
||||
|
||||
! traj_vx
|
||||
write(u, '(a)') ' "traj_vx": ['
|
||||
do s = 1, nsteps
|
||||
call json_arr(u, tvx(s, :), nat, s < nsteps, ' ')
|
||||
end do
|
||||
write(u, '(a)') ' ],'
|
||||
|
||||
! traj_vy
|
||||
write(u, '(a)') ' "traj_vy": ['
|
||||
do s = 1, nsteps
|
||||
call json_arr(u, tvy(s, :), nat, s < nsteps, ' ')
|
||||
end do
|
||||
write(u, '(a)') ' ],'
|
||||
|
||||
! traj_vz
|
||||
write(u, '(a)') ' "traj_vz": ['
|
||||
do s = 1, nsteps
|
||||
call json_arr(u, tvz(s, :), nat, s < nsteps, ' ')
|
||||
end do
|
||||
write(u, '(a)') ' ],'
|
||||
|
||||
! 标量参数
|
||||
write(buf, '(a, i0, a)') ' "NT": ', NT, ','
|
||||
write(u, '(a)') trim(buf)
|
||||
write(buf, '(a, g0, a)') ' "DT": ', DT, ','
|
||||
write(u, '(a)') trim(buf)
|
||||
write(buf, '(a, i0, a)') ' "NSTEP": ', NSTEP, ','
|
||||
write(u, '(a)') trim(buf)
|
||||
write(buf, '(a, a, a)') ' "method": "', trim(method), '",'
|
||||
write(u, '(a)') trim(buf)
|
||||
write(buf, '(a, i0, a)') ' "warmup_steps": ', warmup, ','
|
||||
write(u, '(a)') trim(buf)
|
||||
|
||||
write(buf, '(a, g0, a, g0, a, g0, a)') &
|
||||
' "G": [', G(1), ', ', G(2), ', ', G(3), '],'
|
||||
write(u, '(a)') trim(buf)
|
||||
write(buf, '(a, g0, a, g0, a, g0, a)') &
|
||||
' "B": [', B(1), ', ', B(2), ', ', B(3), '],'
|
||||
write(u, '(a)') trim(buf)
|
||||
|
||||
! 原子信息
|
||||
write(u, '(a)', advance='no') ' "atom_ids": ['
|
||||
do i = 1, nat
|
||||
if (i > 1) write(u, '(a)', advance='no') ','
|
||||
write(u, '(i0)', advance='no') aid(i)
|
||||
end do
|
||||
write(u, '(a)') '],'
|
||||
|
||||
write(u, '(a)', advance='no') ' "atom_masses": ['
|
||||
do i = 1, nat
|
||||
if (i > 1) write(u, '(a)', advance='no') ','
|
||||
write(u, '(g0)', advance='no') amass(i)
|
||||
end do
|
||||
write(u, '(a)') '],'
|
||||
|
||||
! 成键
|
||||
if (nb > 0) then
|
||||
call write_int2_arr(u, 'bond_pairs', bp, nb, .true.)
|
||||
call write_dbl_arr(u, 'bond_stiffness', bk, nb, .true.)
|
||||
call write_dbl_arr(u, 'bond_rest_lengths', br, nb, .false.)
|
||||
else
|
||||
write(u, '(a)') ' "bond_pairs": [],'
|
||||
write(u, '(a)') ' "bond_stiffness": [],'
|
||||
write(u, '(a)') ' "bond_rest_lengths": []'
|
||||
end if
|
||||
|
||||
write(u, '(a)') '}'
|
||||
close(u)
|
||||
end subroutine write_json
|
||||
|
||||
! 写出单行 JSON 数组 [v1, v2, ...]
|
||||
subroutine json_arr(u, vals, n, has_next, indent)
|
||||
integer, intent(in) :: u, n
|
||||
double precision, intent(in) :: vals(n)
|
||||
logical, intent(in) :: has_next
|
||||
character(len=*), intent(in) :: indent
|
||||
integer :: i
|
||||
write(u, '(a)', advance='no') indent // '['
|
||||
do i = 1, n
|
||||
if (i > 1) write(u, '(a)', advance='no') ','
|
||||
write(u, '(g0)', advance='no') vals(i)
|
||||
end do
|
||||
if (has_next) then
|
||||
write(u, '(a)') '],'
|
||||
else
|
||||
write(u, '(a)') ']'
|
||||
end if
|
||||
end subroutine json_arr
|
||||
|
||||
subroutine write_int2_arr(u, name, arr, n, has_next)
|
||||
integer, intent(in) :: u, n, arr(n, 2)
|
||||
character(len=*), intent(in) :: name
|
||||
logical, intent(in) :: has_next
|
||||
character(len=65536) :: buf
|
||||
integer :: i, pos
|
||||
write(u, '(a)', advance='no') ' "' // trim(name) // '": ['
|
||||
do i = 1, n
|
||||
if (i > 1) write(u, '(a)', advance='no') ','
|
||||
write(buf, '(a, i0, a, i0, a)') '[', arr(i, 1), ',', arr(i, 2), ']'
|
||||
write(u, '(a)', advance='no') trim(buf)
|
||||
end do
|
||||
if (has_next) then
|
||||
write(u, '(a)') '],'
|
||||
else
|
||||
write(u, '(a)') ']'
|
||||
end if
|
||||
end subroutine write_int2_arr
|
||||
|
||||
subroutine write_dbl_arr(u, name, arr, n, has_next)
|
||||
integer, intent(in) :: u, n
|
||||
double precision, intent(in) :: arr(n)
|
||||
character(len=*), intent(in) :: name
|
||||
logical, intent(in) :: has_next
|
||||
integer :: i
|
||||
write(u, '(a)', advance='no') ' "' // trim(name) // '": ['
|
||||
do i = 1, n
|
||||
if (i > 1) write(u, '(a)', advance='no') ','
|
||||
write(u, '(g0)', advance='no') arr(i)
|
||||
end do
|
||||
if (has_next) then
|
||||
write(u, '(a)') '],'
|
||||
else
|
||||
write(u, '(a)') ']'
|
||||
end if
|
||||
end subroutine write_dbl_arr
|
||||
|
||||
write(*,*) "[Fortran-engine] 待实现"
|
||||
end program dynamics_f90
|
||||
|
||||
@@ -1,3 +1,3 @@
|
||||
n mass radius x y z vx vy vz
|
||||
1 1 0.28 -1 0 0 0 0 0
|
||||
2 1 0.28 1 0 1 0 0 0
|
||||
n mass radius x y z vx vy vz fix_x fix_y fix_z
|
||||
1 1 0.28 -1 0 0 0 0 0 0 0 0
|
||||
2 1 0.28 1 0 1 0 0 0 0 0 0
|
||||
|
||||
@@ -7,18 +7,20 @@
|
||||
# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧
|
||||
step_simulate: 1 # 运行物理模拟 → output/trajectory.txt
|
||||
step_sample: 1 # 抽帧 → output/display.txt
|
||||
step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png
|
||||
step_plot: 1 # 绘制轨迹/能量图 → output/trajectory_plots.png
|
||||
step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
|
||||
force_calc: 0 # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
|
||||
|
||||
# ── 计算引擎 ──────────────────────────────────
|
||||
# 可选: python, c, cpp, fortran, java
|
||||
# python = Python 参考实现(compute.py)
|
||||
# c = C 引擎 (engines/c/build/dynamics_c)
|
||||
# cpp = C++ 引擎 (engines/cpp/build/dynamics_cpp)
|
||||
engine: python # 默认使用 Python 引擎
|
||||
# fortran = Fortran 引擎 (engines/fortran/build/dynamics_f90)
|
||||
engine: python # 默认使用 Python 引擎
|
||||
|
||||
# ── 盒子 ──────────────────────────────────────
|
||||
box_a: 10.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
|
||||
box_a: 20.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
|
||||
|
||||
# ── 初始构型 ──────────────────────────────────
|
||||
# 坐标文件格式：
|
||||
@@ -37,6 +39,14 @@ G: [0.0, 0.0, -9.8] # 重力场分量 (m/s²)
|
||||
# B: [0.5, 0.5, 0.5] # 阻尼分量
|
||||
B: [0.0, 0.0, 0.0] # 阻尼分量
|
||||
|
||||
# ── 力开关（0=关闭, 1=开启）──────────────────
|
||||
gravity_field: 1 # 均匀重力场 (G)
|
||||
gravity_interaction: 0 # 原子间万有引力
|
||||
elastic_force: 1 # 弹簧键力
|
||||
damping_force: 0 # 阻尼 (B)
|
||||
#
|
||||
gravity_strength: 1.0 # 万有引力强度（仅 gravity_interaction=1 时有效）
|
||||
|
||||
# ── 数值算法 ──────────────────────────────────
|
||||
# 可选：
|
||||
# explicit_euler 显式欧拉法
|
||||
@@ -51,18 +61,21 @@ method: leapfrog
|
||||
# 预热步数：模拟开始时跳过不保存的步数（用于稳定初始状态）
|
||||
warmup_steps: 0 # 默认 0（立即开始记录）
|
||||
|
||||
# 总计算步数
|
||||
NT: 10000
|
||||
# 总模拟时间（秒），程序自动计算 NT = T_total / DT
|
||||
# 如果同时指定了 NT，以 NT 为准
|
||||
T_total: 10.0
|
||||
|
||||
# 抽帧间隔（每 NSTEP 步取一帧用于动画）
|
||||
NSTEP: 100
|
||||
|
||||
# ── 时间步长 ──────────────────────────────────
|
||||
DT: 0.001 # 时间步长 (s)
|
||||
|
||||
# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧
|
||||
sample_start: null # null 表示从头开始（帧索引从 0 起）
|
||||
sample_end: null # null 表示到末尾
|
||||
|
||||
# ── 时间步长 ──────────────────────────────────
|
||||
DT: 0.001 # 时间步长 (s)
|
||||
|
||||
|
||||
# ── 显示参数 ──────────────────────────────────
|
||||
# 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序
|
||||
@@ -18,7 +18,7 @@ CASE_DIR = Path(__file__).resolve().parent
|
||||
DYNAMICS_PATH = Path("..") / ".." / "dynamics.py"
|
||||
INPUT_DIR = Path("input")
|
||||
OUTPUT_DIR = Path("output")
|
||||
CONFIG_FILE = INPUT_DIR / "parameters.yaml"
|
||||
CONFIG_FILE = INPUT_DIR / "input.txt"
|
||||
|
||||
|
||||
def load_dynamics_module(module_path: Path):
|
||||
|
||||
@@ -0,0 +1,2 @@
|
||||
bond_name k rest_length
|
||||
k1 100.0 2.0
|
||||
@@ -0,0 +1,2 @@
|
||||
n1 n2 bond_name
|
||||
|
||||
@@ -0,0 +1,3 @@
|
||||
n mass radius x y z vx vy vz fix_x fix_y fix_z
|
||||
1 1 0.28 0 0 0 0 0 0 1 1 1
|
||||
2 1 0.28 0 0 4 0 4 0 0 0 0
|
||||
@@ -0,0 +1,93 @@
|
||||
# 物理模拟参数配置
|
||||
# 格式：YAML
|
||||
# 用法：python run_dynamics.py
|
||||
|
||||
# ── 流程控制 ──────────────────────────────────
|
||||
# 每步用 0/1 单独开关，1=执行，0=跳过
|
||||
# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧
|
||||
step_simulate: 1 # 运行物理模拟 → output/trajectory.txt
|
||||
step_sample: 1 # 抽帧 → output/display.txt
|
||||
step_plot: 1 # 绘制轨迹/能量图 → output/trajectory_plots.png
|
||||
step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy）
|
||||
force_calc: 0 # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出
|
||||
|
||||
# ── 计算引擎 ──────────────────────────────────
|
||||
# 可选: python, c, cpp, fortran, java
|
||||
# python = Python 参考实现(compute.py)
|
||||
# c = C 引擎 (engines/c/build/dynamics_c)
|
||||
# cpp = C++ 引擎 (engines/cpp/build/dynamics_cpp)
|
||||
# fortran = Fortran 引擎 (engines/fortran/build/dynamics_f90)
|
||||
engine: python # 默认使用 Python 引擎
|
||||
|
||||
# ── 盒子 ──────────────────────────────────────
|
||||
box_a: 20.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
|
||||
|
||||
# ── 初始构型 ──────────────────────────────────
|
||||
# 坐标文件格式：
|
||||
# 第一行：n mass radius x y z vx vy vz
|
||||
# 后续行：原子序号 质量 半径 x y z vx vy vz
|
||||
coord_file: input/coord.txt
|
||||
connection_file: input/connection.txt
|
||||
bond_file: input/bond.txt
|
||||
|
||||
# 绘图/动画展示的原子序号（对应 coord_file 第一列 n）
|
||||
plot_atom: 1
|
||||
|
||||
# ── 物理参数 ──────────────────────────────────
|
||||
# 三个方向分量分别对应 x, y, z
|
||||
G: [0.0, 0.0, -9.8] # 重力场分量 (m/s²)
|
||||
# B: [0.5, 0.5, 0.5] # 阻尼分量
|
||||
B: [0.0, 0.0, 0.0] # 阻尼分量
|
||||
|
||||
# ── 力开关（0=关闭, 1=开启）──────────────────
|
||||
gravity_field: 0 # 均匀重力场 (G)
|
||||
gravity_interaction: 1 # 原子间万有引力
|
||||
elastic_force: 0 # 弹簧键力
|
||||
damping_force: 0 # 阻尼 (B)
|
||||
#
|
||||
gravity_strength: 100.0 # 万有引力强度（仅 gravity_interaction=1 时有效）
|
||||
|
||||
# ── 数值算法 ──────────────────────────────────
|
||||
# 可选：
|
||||
# explicit_euler 显式欧拉法
|
||||
# implicit_euler 隐式欧拉法
|
||||
# midpoint 中点法
|
||||
# leapfrog 蛙跳法
|
||||
method: leapfrog
|
||||
|
||||
# ── 步骤控制 ──────────────────────────────────
|
||||
# 以下参数控制哪些步骤被执行和保存
|
||||
|
||||
# 预热步数：模拟开始时跳过不保存的步数（用于稳定初始状态）
|
||||
warmup_steps: 0 # 默认 0（立即开始记录）
|
||||
|
||||
# 总模拟时间（秒），程序自动计算 NT = T_total / DT
|
||||
# 如果同时指定了 NT，以 NT 为准
|
||||
T_total: 10.0
|
||||
|
||||
# 抽帧间隔（每 NSTEP 步取一帧用于动画）
|
||||
NSTEP: 100
|
||||
|
||||
# ── 时间步长 ──────────────────────────────────
|
||||
DT: 0.001 # 时间步长 (s)
|
||||
|
||||
# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧
|
||||
sample_start: null # null 表示从头开始（帧索引从 0 起）
|
||||
sample_end: null # null 表示到末尾
|
||||
|
||||
|
||||
|
||||
# ── 显示参数 ──────────────────────────────────
|
||||
# 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序
|
||||
alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5]
|
||||
|
||||
# 小球颜色
|
||||
# 小球半径从 coord_file 的 radius 列读取
|
||||
ball_color_r: 0.90 # R 分量 (0~1)
|
||||
ball_color_g: 0.20 # G 分量
|
||||
ball_color_b: 0.20 # B 分量
|
||||
|
||||
# 盒子面颜色
|
||||
box_color_r: 0.80
|
||||
box_color_g: 0.80
|
||||
box_color_b: 0.85
|
||||
@@ -0,0 +1,54 @@
|
||||
"""
|
||||
Case runner for Dynamics case01.
|
||||
|
||||
This script keeps program and data separated:
|
||||
- program: ../../dynamics.py
|
||||
- input: ./input
|
||||
- output: ./output
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import importlib.util
|
||||
from pathlib import Path
|
||||
|
||||
|
||||
CASE_DIR = Path(__file__).resolve().parent
|
||||
DYNAMICS_PATH = Path("..") / ".." / "dynamics.py"
|
||||
INPUT_DIR = Path("input")
|
||||
OUTPUT_DIR = Path("output")
|
||||
CONFIG_FILE = INPUT_DIR / "input.txt"
|
||||
|
||||
|
||||
def load_dynamics_module(module_path: Path):
|
||||
spec = importlib.util.spec_from_file_location("dynamics_module", module_path)
|
||||
if spec is None or spec.loader is None:
|
||||
raise ImportError(f"无法加载 dynamics.py: {module_path}")
|
||||
module = importlib.util.module_from_spec(spec)
|
||||
spec.loader.exec_module(module)
|
||||
return module
|
||||
|
||||
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01")
|
||||
parser.add_argument("--no-plot", action="store_true", help="跳过 matplotlib 绘图")
|
||||
args = parser.parse_args()
|
||||
|
||||
dynamics_path = (CASE_DIR / DYNAMICS_PATH).resolve()
|
||||
input_dir = (CASE_DIR / INPUT_DIR).resolve()
|
||||
output_dir = (CASE_DIR / OUTPUT_DIR).resolve()
|
||||
config_path = (CASE_DIR / CONFIG_FILE).resolve()
|
||||
|
||||
module = load_dynamics_module(dynamics_path)
|
||||
module.run_case(
|
||||
config_path=config_path,
|
||||
runtime_base=CASE_DIR,
|
||||
input_dir=input_dir,
|
||||
output_dir=output_dir,
|
||||
no_plot=args.no_plot,
|
||||
)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user