diff --git a/compute.py b/compute.py index 68d7fa8..81f9cc1 100644 --- a/compute.py +++ b/compute.py @@ -69,8 +69,12 @@ GRAVITY_FIELD = 1 # 均匀重力场 GRAVITY_INTERACTION = 0 # 原子间万有引力 ELASTIC_FORCE = 1 # 弹簧键力 DAMPING_FORCE = 0 # 阻尼 +DRIVING_FORCE = 0 # 驱动力 GRAVITY_STRENGTH = 1.0 +# 驱动力数据 +DRIVER_DATA = None # 加载自 driver.txt + # 派生边界（根据 box_a 计算） X_MIN = None X_MAX = None @@ -328,6 +332,135 @@ def load_bond_connections(connection_path, atom_ids, atom_positions, bond_map): ) +# =========================================================================== +# 驱动力加载 & 应用 +# =========================================================================== + +def load_driver_file(driver_path, atom_ids): + """从 driver.txt 加载驱动力定义。 + + 格式： + n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period + 数值 0 0 0 0 0 10 0 0 90 all + + 其中： + position = A * cos(2π f t + φ), φ 为角度制（自动转弧度） + period = all | 数值（周期数，结束后原子静止） + """ + if not os.path.exists(driver_path): + print(f"[compute] 警告: 驱动力文件不存在: {driver_path}") + return None + + atom_index = {int(aid): idx for idx, aid in enumerate(atom_ids)} + drivers = [] + ncols = 11 + + with open(driver_path, "r", encoding="utf-8") as f: + header = f.readline().strip().split() + expected = ["n", "amp_x", "amp_y", "amp_z", + "freq_x", "freq_y", "freq_z", + "phi_x", "phi_y", "phi_z", "period"] + if header != expected: + raise ValueError( + f"driver.txt 表头应为: {' '.join(expected)}，实际为: {' '.join(header)}") + + for line_no, line in enumerate(f, start=2): + stripped = line.strip() + if not stripped or stripped.startswith("#"): + continue + parts = stripped.split() + if len(parts) != ncols: + raise ValueError( + f"{driver_path}:{line_no} 应有 {ncols} 列，实际为 {len(parts)} 列") + + n = int(parts[0]) + if n not in atom_index: + raise ValueError(f"{driver_path}:{line_no} 原子序号 {n} 不在 coord.txt 中") + + amp = np.array([float(parts[1]), float(parts[2]), float(parts[3])], + dtype=np.float64) + freq = np.array([float(parts[4]), float(parts[5]), float(parts[6])], + dtype=np.float64) + phi_deg = np.array([float(parts[7]), float(parts[8]), float(parts[9])], + dtype=np.float64) + phi_rad = phi_deg * np.pi / 180.0 + + period_str = parts[10] + + drivers.append({ + "atom_index": atom_index[n], + "atom_id": n, + "amp": amp, + "freq": freq, + "phi": phi_rad, + "period_str": period_str, + "period_cycles": None if period_str == "all" else float(period_str), + # 在模拟中动态记录：冻结步数索引、冻结位置 + "freeze_step": None, + "freeze_pos": None, + }) + + if not drivers: + print(f"[compute] 警告: driver.txt 中没有有效数据") + return None + + print(f"[compute] 已加载驱动力: {len(drivers)} 条定义") + for d in drivers: + print(f" 原子 {d['atom_id']}: " + f"A=({d['amp'][0]},{d['amp'][1]},{d['amp'][2]}), " + f"f=({d['freq'][0]},{d['freq'][1]},{d['freq'][2]}), " + f"φ=({phi_deg[d['amp'].tolist().index(max(d['amp']))]}° 等), " + f"period={d['period_str']}") + return drivers + + +def apply_driving_force(x, y, z, vx, vy, vz, t, step, drivers, dt): + """对受驱原子按驱动力函数覆盖位置/速度。 + + 驱动力公式：pos = A * cos(2π f t + φ) + vel = -A * 2π f * sin(2π f t + φ) + """ + if drivers is None: + return x, y, z, vx, vy, vz + + for d in drivers: + idx = d["atom_index"] + + # 确定驱动力持续到哪一步 + if d["period_cycles"] is not None: + max_freq = np.max(np.abs(d["freq"])) + if max_freq > 1e-12: + period_duration = d["period_cycles"] / max_freq + period_steps = int(period_duration / dt) + else: + period_steps = 0 + + if step > period_steps: + # 驱动力已结束：原子静止（保持最后位置，速度归零） + if d["freeze_pos"] is not None: + x[idx], y[idx], z[idx] = d["freeze_pos"] + vx[idx] = vy[idx] = vz[idx] = 0.0 + continue + else: + period_steps = None # 全程驱动 + + # 当前驱动力下的位置 / 速度 + t_vec = np.array([t, t, t], dtype=np.float64) + pos_drive = d["amp"] * np.cos(2.0 * np.pi * d["freq"] * t_vec + d["phi"]) + vel_drive = -d["amp"] * 2.0 * np.pi * d["freq"] * np.sin(2.0 * np.pi * d["freq"] * t_vec + d["phi"]) + + x[idx] = pos_drive[0] + y[idx] = pos_drive[1] + z[idx] = pos_drive[2] + vx[idx] = vel_drive[0] + vy[idx] = vel_drive[1] + vz[idx] = vel_drive[2] + + # 若周期有限，记录冻结位置（驱动力最后一帧的位置） + if period_steps is not None and step == period_steps: + d["freeze_pos"] = (float(pos_drive[0]), float(pos_drive[1]), float(pos_drive[2])) + + return x, y, z, vx, vy, vz def parse_gravity_vector(value): """Parse gravity into a 3D acceleration vector. @@ -379,6 +512,7 @@ def run_from_config(config, out_dir=None): 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 + global DRIVING_FORCE, DRIVER_DATA box_a = float(config.get("box_a", 10.0)) raw_alpha = config.get("alpha", 0.2) @@ -439,11 +573,22 @@ def run_from_config(config, out_dir=None): # 力开关 global GRAVITY_FIELD, GRAVITY_INTERACTION, ELASTIC_FORCE, DAMPING_FORCE, GRAVITY_STRENGTH + global DRIVING_FORCE, DRIVER_DATA 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)) + DRIVING_FORCE = int(config.get("driving_force", 0)) + + # 加载驱动力定义 + DRIVER_DATA = None + if DRIVING_FORCE: + driver_rel = str(config.get("driver_file", os.path.join("input", "driver.txt"))) + driver_path = driver_rel + if out_dir is not None and not os.path.isabs(driver_rel): + driver_path = os.path.join(out_dir, driver_rel) + DRIVER_DATA = load_driver_file(driver_path, ATOM_IDS) print(f"[compute] 使用算法: {METHOD}") print(f"[compute] 已加载成键信息: {len(BOND_PAIRS)} 条键") @@ -1062,11 +1207,14 @@ def run_simulation(): vy = ATOM_VELOCITIES[:, 1].copy() vz = ATOM_VELOCITIES[:, 2].copy() x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz) + # 初始时刻驱动力（t=0 时原子 1 的位置由驱动力决定而非 coord.txt） + x, y, z, vx, vy, vz = apply_driving_force(x, y, z, vx, vy, vz, 0.0, 0, DRIVER_DATA, DT) if warmup_steps is not None and warmup_steps > 0: print(f"[compute] 预热阶段: 前 {warmup_steps} 步不记录") for step in trange(warmup_steps, desc="[compute] 预热"): t = (step + 1) * DT + x, y, z, vx, vy, vz = apply_driving_force(x, y, z, vx, vy, vz, t, step, DRIVER_DATA, DT) x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B) x, y, z = wrap_position(x, y, z) x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz) @@ -1086,6 +1234,10 @@ def run_simulation(): traj_vz = np.zeros((record_steps, n_atoms), dtype=np.float64) for step in trange(record_steps, desc="[compute] 计算中"): + t = (step + (warmup_steps or 0)) * DT + # 先施加驱动力（受驱原子的位置覆盖初始/固定约束，为弹簧力提供正确参考） + x, y, z, vx, vy, vz = apply_driving_force(x, y, z, vx, vy, vz, t, step, DRIVER_DATA, DT) + traj_x[step] = x traj_y[step] = y traj_z[step] = z diff --git a/draw.py b/draw.py index 530d221..904f7d4 100644 --- a/draw.py +++ b/draw.py @@ -72,6 +72,9 @@ PLOT_ATOM_ROW = int(disp_data.get("plot_atom_row", 0)) PLOT_ATOM_ID = int(disp_data.get("plot_atom_id", ATOM_IDS[0])) BOND_PAIRS = disp_data.get("bond_pairs", []) +# 渲染方式：0=Sphere(网格球体), 1=Marker(GPU点精灵) +USE_MARKER = int(disp_data.get("use_marker", 0)) + if N_FRAMES <= 0: raise ValueError( "output/display.txt 中没有可播放的帧，请检查 sample_start/sample_end/NSTEP 配置。") @@ -171,7 +174,10 @@ axes_group.append(scene.visuals.Text(text="y", color=(0.2, 1.0, 0.2, 1.0), font_ axes_group.append(scene.visuals.Text(text="z", color=(0.3, 0.6, 1.0, 1.0), font_size=18, pos=(0, 0, axis_length + 0.2), anchor_x="left", anchor_y="bottom", parent=view.scene)) -# 所有小球（每个原子一个球，不同颜色） +# ── 原子渲染 ────────────────────────────────── +# 两种渲染模式： +# Sphere = 独立网格球体（精细，适合少量原子） +# Marker = GPU 实例化点精灵（高效，适合大量原子） TAB10_RGB = np.array([ [0.1216, 0.4667, 0.7059], # 蓝 [0.8902, 0.4667, 0.1137], # 橙 @@ -184,23 +190,41 @@ TAB10_RGB = np.array([ [0.7373, 0.7412, 0.1333], # 黄绿 [0.0902, 0.7451, 0.8118], # 青 ]) -balls = [] +# 每个原子的颜色（循环使用 tab10 色板） +atom_colors = np.zeros((N_ATOMS, 4), dtype=np.float32) for i in range(N_ATOMS): r, g, b = TAB10_RGB[i % len(TAB10_RGB)] - s = scene.visuals.Sphere( - radius=float(ATOM_RADII[i]), method="latitude", - color=(r, g, b, 1.0), edge_color=None, parent=view.scene) - s.mesh.shading = "smooth" - balls.append(s) + atom_colors[i] = [r, g, b, 1.0] + +if USE_MARKER: + # ── Marker 模式：GPU 实例化，一次 draw call ── + marker_pos = np.zeros((N_ATOMS, 3), dtype=np.float32) + marker_sizes = np.array([float(ATOM_RADII[i]) * 40 for i in range(N_ATOMS)], dtype=np.float32) + balls = scene.visuals.Markers(parent=view.scene) + balls.set_data( + pos=marker_pos, face_color=atom_colors, + size=marker_sizes, symbol='disc', edge_width=0, + ) +else: + # ── Sphere 模式：每个原子一个独立网格球体 ── + balls = [] + for i in range(N_ATOMS): + r, g, b, _ = atom_colors[i] + s = scene.visuals.Sphere( + radius=float(ATOM_RADII[i]), method="latitude", + color=(r, g, b, 1.0), edge_color=None, parent=view.scene) + s.mesh.shading = "smooth" + balls.append(s) # 成键线（原子之间的连接） if len(BOND_PAIRS) > 0: n_bonds = len(BOND_PAIRS) - bond_pos = np.zeros((n_bonds * 2, 3), dtype=np.float32) + bond_pos_buffer = np.zeros((n_bonds * 2, 3), dtype=np.float32) # 预分配，后续原地修改 bond_lines = scene.visuals.Line( - pos=bond_pos, color=(0.7, 0.7, 0.7, 0.8), width=3, + pos=bond_pos_buffer, color=(0.7, 0.7, 0.7, 0.8), width=3, connect="segments", parent=view.scene) else: + bond_pos_buffer = None bond_lines = None # 六个面形成立方体边界（每个面独立透明度，alpha<=0 时隐藏该面） @@ -361,19 +385,10 @@ def handle_key_press(event): def reset_camera_view(): global frame_idx frame_idx = 0 - # 立即复位所有小球到第 0 帧 - for i in range(N_ATOMS): - balls[i].transform = STTransform(translate=( - float(DISP_ALL_X[frame_idx, i]), - float(DISP_ALL_Y[frame_idx, i]), - float(DISP_ALL_Z[frame_idx, i]), - )) + # 立即复位所有原子到第 0 帧 + _update_atom_positions(frame_idx) if bond_lines is not None and len(BOND_PAIRS) > 0: - pos = np.zeros((len(BOND_PAIRS) * 2, 3), dtype=np.float32) - for b_idx, (i, j) in enumerate(BOND_PAIRS): - pos[b_idx * 2] = [DISP_ALL_X[frame_idx, i], DISP_ALL_Y[frame_idx, i], DISP_ALL_Z[frame_idx, i]] - pos[b_idx * 2 + 1] = [DISP_ALL_X[frame_idx, j], DISP_ALL_Y[frame_idx, j], DISP_ALL_Z[frame_idx, j]] - bond_lines.set_data(pos=pos) + _update_bond_positions(frame_idx) # 复位相机 view.camera.distance = initial_camera["distance"] view.camera.elevation = initial_camera["elevation"] @@ -413,31 +428,51 @@ def handle_mouse_press(event): axis.visible = axes_visible +# =========================================================================== +# 位置/键线更新辅助函数（两种渲染模式共用） +# =========================================================================== + +def _update_atom_positions(f_idx): + """更新所有原子到第 f_idx 帧的位置。""" + if USE_MARKER: + for i in range(N_ATOMS): + marker_pos[i] = [DISP_ALL_X[f_idx, i], DISP_ALL_Y[f_idx, i], DISP_ALL_Z[f_idx, i]] + balls.set_data(pos=marker_pos) + else: + for i in range(N_ATOMS): + balls[i].transform = STTransform(translate=( + float(DISP_ALL_X[f_idx, i]), + float(DISP_ALL_Y[f_idx, i]), + float(DISP_ALL_Z[f_idx, i]), + )) + + +def _update_bond_positions(f_idx): + """原地更新键线顶点位置，避免重新分配内存。""" + for b_idx, (i, j) in enumerate(BOND_PAIRS): + b2 = b_idx * 2 + bond_pos_buffer[b2] = [DISP_ALL_X[f_idx, i], DISP_ALL_Y[f_idx, i], DISP_ALL_Z[f_idx, i]] + bond_pos_buffer[b2 + 1] = [DISP_ALL_X[f_idx, j], DISP_ALL_Y[f_idx, j], DISP_ALL_Z[f_idx, j]] + bond_lines.set_data(pos=bond_pos_buffer) + + # =========================================================================== # 动画初始化 # =========================================================================== frame_idx = 0 -# 初始帧：摆放所有小球并刷新 UI -for i in range(N_ATOMS): - balls[i].transform = STTransform(translate=( - float(DISP_ALL_X[frame_idx, i]), - float(DISP_ALL_Y[frame_idx, i]), - float(DISP_ALL_Z[frame_idx, i]), - )) -# 初始帧：更新成键线 +# 初始帧：摆放所有原子并刷新 UI +_update_atom_positions(frame_idx) if bond_lines is not None and len(BOND_PAIRS) > 0: - pos = np.zeros((len(BOND_PAIRS) * 2, 3), dtype=np.float32) - for b_idx, (i, j) in enumerate(BOND_PAIRS): - pos[b_idx * 2] = [DISP_ALL_X[frame_idx, i], DISP_ALL_Y[frame_idx, i], DISP_ALL_Z[frame_idx, i]] - pos[b_idx * 2 + 1] = [DISP_ALL_X[frame_idx, j], DISP_ALL_Y[frame_idx, j], DISP_ALL_Z[frame_idx, j]] - bond_lines.set_data(pos=pos) + _update_bond_positions(frame_idx) reposition_camera_info() update_ball_info(frame_idx, float(DISP_X[frame_idx]), float(DISP_Y[frame_idx]), float(DISP_Z[frame_idx]), float(DISP_VX[frame_idx]), float(DISP_VY[frame_idx]), float(DISP_VZ[frame_idx])) +mode_str = "Marker (GPU 实例化)" if USE_MARKER else "Sphere (独立网格)" print(f"[draw] 加载 output/display.txt: {N_FRAMES} 帧, {N_ATOMS} 个原子, NT={NT}, DT={DT}, NSTEP={NSTEP}") +print(f"[draw] 渲染方式: {mode_str}") print(f"[draw] 绘图参数: ball_radius={ball_radius}, box_color=({box_color_r:.2f},{box_color_g:.2f},{box_color_b:.2f}), alpha={alpha_list}") @@ -448,20 +483,12 @@ def update(event): global frame_idx frame_idx = (frame_idx + 1) % N_FRAMES # 循环播放 - # 更新所有小球位置 - for i in range(N_ATOMS): - x = float(DISP_ALL_X[frame_idx, i]) - y = float(DISP_ALL_Y[frame_idx, i]) - z = float(DISP_ALL_Z[frame_idx, i]) - balls[i].transform = STTransform(translate=(x, y, z)) + # 更新所有原子位置（两种模式共用逻辑） + _update_atom_positions(frame_idx) - # 更新成键线 + # 更新成键线（原地修改，无内存分配） if bond_lines is not None and len(BOND_PAIRS) > 0: - pos = np.zeros((len(BOND_PAIRS) * 2, 3), dtype=np.float32) - for b_idx, (i, j) in enumerate(BOND_PAIRS): - pos[b_idx * 2] = [DISP_ALL_X[frame_idx, i], DISP_ALL_Y[frame_idx, i], DISP_ALL_Z[frame_idx, i]] - pos[b_idx * 2 + 1] = [DISP_ALL_X[frame_idx, j], DISP_ALL_Y[frame_idx, j], DISP_ALL_Z[frame_idx, j]] - bond_lines.set_data(pos=pos) + _update_bond_positions(frame_idx) # 信息面板显示 plot_atom 的数据 x = float(DISP_X[frame_idx]) diff --git a/dynamics.py b/dynamics.py index 0cd3504..3e04d01 100644 --- a/dynamics.py +++ b/dynamics.py @@ -12,6 +12,7 @@ dynamics.py import os import sys import subprocess +import time import argparse from contextlib import contextmanager from pathlib import Path @@ -316,6 +317,8 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", "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)), + "driving_force": int(config.get("driving_force", 0)), + "use_marker": int(config.get("use_marker", 0)), } save_display_txt(disp_data, str(runtime_base)) print(f"[run] 抽帧完成: {sample_end - sample_start} 步 -> {n_frames} 帧") @@ -331,7 +334,7 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'WenQuanYi Micro Hei', 'DejaVu Sans'] plt.rcParams['axes.unicode_minus'] = False - time = np.arange(NT) * DT + time_arr = np.arange(NT) * DT n_atoms = all_x.shape[1] atom_ids_list = data.get("atom_ids", np.arange(n_atoms) + 1) @@ -353,7 +356,7 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", for ax, data_arr, title in plot_configs: for i in range(n_atoms): atom_id = int(atom_ids_list[i]) - ax.plot(time, data_arr[:, i], color=colors[i], linewidth=1.5, label=f"原子 {atom_id}") + ax.plot(time_arr, data_arr[:, i], color=colors[i], linewidth=1.5, label=f"原子 {atom_id}") ax.set_title(title) ax.set_xlabel("时间 (s)") ax.grid(True, alpha=0.3) @@ -423,7 +426,7 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", ax_e = axes[2, 0] for i in range(n_atoms): aid = int(atom_ids_list[i]) - ax_e.plot(time, e_total[:, i], color=colors[i], linewidth=1.5, label=f"原子 {aid}") + ax_e.plot(time_arr, e_total[:, i], color=colors[i], linewidth=1.5, label=f"原子 {aid}") ax_e.set_title("各原子总能量") ax_e.set_xlabel("时间 (s)") ax_e.set_ylabel("能量") @@ -432,13 +435,13 @@ 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="均匀重力势能") + ax_sys.plot(time_arr, ek_sys, 'b-', linewidth=1.5, label="系统动能") + ax_sys.plot(time_arr, 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="系统弹性势能") + ax_sys.plot(time_arr, 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.plot(time_arr, ug_grav_sys, color='purple', linewidth=1.5, label="万有引力势能") + ax_sys.plot(time_arr, e_sys, 'r--', linewidth=1.5, label="系统总能量") ax_sys.set_title("系统总能量") ax_sys.set_xlabel("时间 (s)") ax_sys.set_ylabel("能量") @@ -461,19 +464,39 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", # 5. 自动播放动画（可选） if config.get("step_animation", 0): draw_script = os.path.join(os.path.dirname(os.path.abspath(__file__)), "draw.py") - if os.path.exists(draw_script): - try: - print("[run] 正在启动 VisPy 3D 动画窗口…") - subprocess.Popen( - [sys.executable, draw_script], - cwd=runtime_base, - stdout=subprocess.DEVNULL, - stderr=subprocess.DEVNULL, - ) - except Exception as e: - print(f"[run] 启动动画失败: {e}") - else: + if not os.path.exists(draw_script): print(f"[run] 未找到动画脚本: {draw_script}") + else: + # 检查 display.txt 是否存在（step_sample=0 时可能没有） + disp_path = os.path.join(output_dir_abs, "display.txt") + if not os.path.exists(disp_path): + print(f"[run] 错误: 找不到 {disp_path}") + print(f"[run] 启动动画需要先运行抽帧（step_sample: 1），或手动保留 output/display.txt") + else: + try: + print("[run] 正在启动 VisPy 3D 动画窗口…") + ansi_log = os.path.join(output_dir_abs, "animation.log") + if sys.platform == "win32": + # Windows 上给子进程独立控制台，避免父进程退出时连带关闭 + creation_flags = subprocess.CREATE_NEW_CONSOLE + else: + creation_flags = 0 + proc = subprocess.Popen( + [sys.executable, draw_script, output_dir_abs], + cwd=runtime_base, + stdout=subprocess.DEVNULL, + stderr=open(ansi_log, "w", encoding="utf-8"), + creationflags=creation_flags, + ) + # 等待半秒检查子进程是否成功启动（未立即崩溃） + time.sleep(0.5) + if proc.poll() is not None: + print(f"[run] ⚠ 动画进程已退出，返回码={proc.returncode}") + print(f"[run] 请查看错误日志: {ansi_log}") + else: + print(f"[run] VisPy 动画窗口已启动（PID={proc.pid}）") + except Exception as e: + print(f"[run] 启动动画失败: {e}") else: print("[run] 运行 python draw.py 查看动画。") diff --git a/examples/case01/Readme.md b/examples/case01/Readme.md new file mode 100644 index 0000000..2dc8a66 --- /dev/null +++ b/examples/case01/Readme.md @@ -0,0 +1 @@ +本案例描述两个粒子同一个键连接，受到重力作用下的动力学过程。 \ No newline at end of file diff --git a/examples/case01/input/input.txt b/examples/case01/input/input.txt index 8be2c8d..549da91 100644 --- a/examples/case01/input/input.txt +++ b/examples/case01/input/input.txt @@ -29,6 +29,7 @@ box_a: 20.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ coord_file: input/coord.txt connection_file: input/connection.txt bond_file: input/bond.txt +driver_file: input/driver.txt # 驱动力定义文件（driving_force=1 时生效） # 绘图/动画展示的原子序号（对应 coord_file 第一列 n） plot_atom: 1 @@ -44,6 +45,7 @@ gravity_field: 1 # 均匀重力场 (G) gravity_interaction: 0 # 原子间万有引力 elastic_force: 1 # 弹簧键力 damping_force: 0 # 阻尼 (B) +driving_force: 0 # 驱动力（需 driver_file 定义） # gravity_strength: 1.0 # 万有引力强度（仅 gravity_interaction=1 时有效） @@ -77,6 +79,12 @@ sample_end: null # null 表示到末尾 +# ── 渲染方式 ────────────────────────────────── +# 3D 动画中原子渲染方式： +# 0 = Sphere （网格球体，效果精细，原子数少时推荐） +# 1 = Marker （GPU 实例化点，原子数多时性能更佳） +use_marker: 0 + # ── 显示参数 ────────────────────────────────── # 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序 alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5] diff --git a/examples/case02/Readme.md b/examples/case02/Readme.md new file mode 100644 index 0000000..c0cf993 --- /dev/null +++ b/examples/case02/Readme.md @@ -0,0 +1 @@ +本案例描述地球围绕太阳做椭圆运动 \ No newline at end of file diff --git a/examples/case02/input/coord.txt b/examples/case02/input/coord.txt index 417b237..573905e 100644 --- a/examples/case02/input/coord.txt +++ b/examples/case02/input/coord.txt @@ -1,3 +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 +1 1 0.28 0 0 0 0 0 0 1 1 1 +2 1 0.28 4 0 0 0 0 4 0 0 0 diff --git a/examples/case02/input/input.txt b/examples/case02/input/input.txt index 222af60..a5bf5de 100644 --- a/examples/case02/input/input.txt +++ b/examples/case02/input/input.txt @@ -29,6 +29,7 @@ box_a: 20.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ coord_file: input/coord.txt connection_file: input/connection.txt bond_file: input/bond.txt +driver_file: input/driver.txt # 驱动力定义文件（driving_force=1 时生效） # 绘图/动画展示的原子序号（对应 coord_file 第一列 n） plot_atom: 1 @@ -44,6 +45,7 @@ gravity_field: 0 # 均匀重力场 (G) gravity_interaction: 1 # 原子间万有引力 elastic_force: 0 # 弹簧键力 damping_force: 0 # 阻尼 (B) +driving_force: 0 # 驱动力（需 driver_file 定义） # gravity_strength: 100.0 # 万有引力强度（仅 gravity_interaction=1 时有效） @@ -77,6 +79,12 @@ sample_end: null # null 表示到末尾 +# ── 渲染方式 ────────────────────────────────── +# 3D 动画中原子渲染方式： +# 0 = Sphere （网格球体，效果精细，原子数少时推荐） +# 1 = Marker （GPU 实例化点，原子数多时性能更佳） +use_marker: 0 + # ── 显示参数 ────────────────────────────────── # 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序 alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5] diff --git a/examples/case03/Readme.md b/examples/case03/Readme.md new file mode 100644 index 0000000..3c1fbf6 --- /dev/null +++ b/examples/case03/Readme.md @@ -0,0 +1 @@ +地球围绕太阳转，月球围绕地球转，失败案例 \ No newline at end of file diff --git a/examples/case03/input/bond.txt b/examples/case03/input/bond.txt new file mode 100644 index 0000000..2b183d5 --- /dev/null +++ b/examples/case03/input/bond.txt @@ -0,0 +1,2 @@ +bond_name k rest_length +k1 100.0 2.0 diff --git a/examples/case03/input/connection.txt b/examples/case03/input/connection.txt new file mode 100644 index 0000000..c845a64 --- /dev/null +++ b/examples/case03/input/connection.txt @@ -0,0 +1,2 @@ +n1 n2 bond_name + diff --git a/examples/case03/input/coord.txt b/examples/case03/input/coord.txt new file mode 100644 index 0000000..7b5d01a --- /dev/null +++ b/examples/case03/input/coord.txt @@ -0,0 +1,4 @@ +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 4 0 0 0 0 4 0 0 0 +3 0.1 0.18 5 0 0 0 0 6 0 0 0 \ No newline at end of file diff --git a/examples/case03/input/input.txt b/examples/case03/input/input.txt new file mode 100644 index 0000000..fbe8bcb --- /dev/null +++ b/examples/case03/input/input.txt @@ -0,0 +1,101 @@ +# 物理模拟参数配置 +# 格式：YAML +# 用法：python run_dynamics.py + +# ── 流程控制 ────────────────────────────────── +# 每步用 0/1 单独开关，1=执行，0=跳过 +# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧 +step_simulate: 1 # 运行物理模拟 → output/trajectory.txt +step_sample: 1 # 抽帧 → output/display.txt +step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy） +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 +driver_file: input/driver.txt # 驱动力定义文件（driving_force=1 时生效） + +# 绘图/动画展示的原子序号（对应 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) +driving_force: 0 # 驱动力（需 driver_file 定义） +# +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: 2 + +# ── 时间步长 ────────────────────────────────── +DT: 0.001 # 时间步长 (s) + +# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧 +sample_start: null # null 表示从头开始（帧索引从 0 起） +sample_end: null # null 表示到末尾 + + + +# ── 渲染方式 ────────────────────────────────── +# 3D 动画中原子渲染方式： +# 0 = Sphere （网格球体，效果精细，原子数少时推荐） +# 1 = Marker （GPU 实例化点，原子数多时性能更佳） +use_marker: 0 + +# ── 显示参数 ────────────────────────────────── +# 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序 +alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5] + +# 小球颜色 +# 小球半径从 coord_file 的 radius 列读取 +ball_color_r: 0.90 # R 分量 (0~1) +ball_color_g: 0.20 # G 分量 +ball_color_b: 0.20 # B 分量 + +# 盒子面颜色 +box_color_r: 0.80 +box_color_g: 0.80 +box_color_b: 0.85 diff --git a/examples/case03/run_dynamics.py b/examples/case03/run_dynamics.py new file mode 100644 index 0000000..d092d49 --- /dev/null +++ b/examples/case03/run_dynamics.py @@ -0,0 +1,54 @@ +""" +Case runner for Dynamics case01. + +This script keeps program and data separated: + - program: ../../dynamics.py + - input: ./input + - output: ./output +""" + +from __future__ import annotations + +import argparse +import importlib.util +from pathlib import Path + + +CASE_DIR = Path(__file__).resolve().parent +DYNAMICS_PATH = Path("..") / ".." / "dynamics.py" +INPUT_DIR = Path("input") +OUTPUT_DIR = Path("output") +CONFIG_FILE = INPUT_DIR / "input.txt" + + +def load_dynamics_module(module_path: Path): + spec = importlib.util.spec_from_file_location("dynamics_module", module_path) + if spec is None or spec.loader is None: + raise ImportError(f"无法加载 dynamics.py: {module_path}") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def main(): + parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01") + parser.add_argument("--no-plot", action="store_true", help="跳过 matplotlib 绘图") + args = parser.parse_args() + + dynamics_path = (CASE_DIR / DYNAMICS_PATH).resolve() + input_dir = (CASE_DIR / INPUT_DIR).resolve() + output_dir = (CASE_DIR / OUTPUT_DIR).resolve() + config_path = (CASE_DIR / CONFIG_FILE).resolve() + + module = load_dynamics_module(dynamics_path) + module.run_case( + config_path=config_path, + runtime_base=CASE_DIR, + input_dir=input_dir, + output_dir=output_dir, + no_plot=args.no_plot, + ) + + +if __name__ == "__main__": + main() diff --git a/examples/case04/Readme.md b/examples/case04/Readme.md new file mode 100644 index 0000000..dc86f6c --- /dev/null +++ b/examples/case04/Readme.md @@ -0,0 +1,15 @@ +地球围绕太阳转，月球围绕地球转，成功案例 +| 项目 | 数值 | 比例 | +| --------- | --------------------: | ----------------: | +| 太阳—地球平均距离 | (1.496\times 10^8) km | 约为地月距离的 **389 倍** | +| 地球—月球平均距离 | (3.844\times 10^5) km | 1 | + +R日地? : R地?月 ≈ 389 + +| 天体 | 质量/kg | 与地球质量之比 | +| -- | --------------------: | -------: | +| 太阳 | (1.989\times 10^{30}) | (333000) | +| 地球 | (5.972\times 10^{24}) | (1) | +| 月球 | (7.35\times 10^{22}) | (0.0123) | + +M太阳?:M地球?:M月球?≈27100000:81.3:1 \ No newline at end of file diff --git a/examples/case04/input/bond.txt b/examples/case04/input/bond.txt new file mode 100644 index 0000000..2b183d5 --- /dev/null +++ b/examples/case04/input/bond.txt @@ -0,0 +1,2 @@ +bond_name k rest_length +k1 100.0 2.0 diff --git a/examples/case04/input/connection.txt b/examples/case04/input/connection.txt new file mode 100644 index 0000000..c845a64 --- /dev/null +++ b/examples/case04/input/connection.txt @@ -0,0 +1,2 @@ +n1 n2 bond_name + diff --git a/examples/case04/input/coord.txt b/examples/case04/input/coord.txt new file mode 100644 index 0000000..7b5d01a --- /dev/null +++ b/examples/case04/input/coord.txt @@ -0,0 +1,4 @@ +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 4 0 0 0 0 4 0 0 0 +3 0.1 0.18 5 0 0 0 0 6 0 0 0 \ No newline at end of file diff --git a/examples/case04/input/input.txt b/examples/case04/input/input.txt new file mode 100644 index 0000000..fbe8bcb --- /dev/null +++ b/examples/case04/input/input.txt @@ -0,0 +1,101 @@ +# 物理模拟参数配置 +# 格式：YAML +# 用法：python run_dynamics.py + +# ── 流程控制 ────────────────────────────────── +# 每步用 0/1 单独开关，1=执行，0=跳过 +# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧 +step_simulate: 1 # 运行物理模拟 → output/trajectory.txt +step_sample: 1 # 抽帧 → output/display.txt +step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy） +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 +driver_file: input/driver.txt # 驱动力定义文件（driving_force=1 时生效） + +# 绘图/动画展示的原子序号（对应 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) +driving_force: 0 # 驱动力（需 driver_file 定义） +# +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: 2 + +# ── 时间步长 ────────────────────────────────── +DT: 0.001 # 时间步长 (s) + +# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧 +sample_start: null # null 表示从头开始（帧索引从 0 起） +sample_end: null # null 表示到末尾 + + + +# ── 渲染方式 ────────────────────────────────── +# 3D 动画中原子渲染方式： +# 0 = Sphere （网格球体，效果精细，原子数少时推荐） +# 1 = Marker （GPU 实例化点，原子数多时性能更佳） +use_marker: 0 + +# ── 显示参数 ────────────────────────────────── +# 盒子透明度：单个数值（统一）或 6 个数的数组，按 [-x,+x,-y,+y,-z,+z] 顺序 +alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5] + +# 小球颜色 +# 小球半径从 coord_file 的 radius 列读取 +ball_color_r: 0.90 # R 分量 (0~1) +ball_color_g: 0.20 # G 分量 +ball_color_b: 0.20 # B 分量 + +# 盒子面颜色 +box_color_r: 0.80 +box_color_g: 0.80 +box_color_b: 0.85 diff --git a/examples/case04/run_dynamics.py b/examples/case04/run_dynamics.py new file mode 100644 index 0000000..d092d49 --- /dev/null +++ b/examples/case04/run_dynamics.py @@ -0,0 +1,54 @@ +""" +Case runner for Dynamics case01. + +This script keeps program and data separated: + - program: ../../dynamics.py + - input: ./input + - output: ./output +""" + +from __future__ import annotations + +import argparse +import importlib.util +from pathlib import Path + + +CASE_DIR = Path(__file__).resolve().parent +DYNAMICS_PATH = Path("..") / ".." / "dynamics.py" +INPUT_DIR = Path("input") +OUTPUT_DIR = Path("output") +CONFIG_FILE = INPUT_DIR / "input.txt" + + +def load_dynamics_module(module_path: Path): + spec = importlib.util.spec_from_file_location("dynamics_module", module_path) + if spec is None or spec.loader is None: + raise ImportError(f"无法加载 dynamics.py: {module_path}") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def main(): + parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01") + parser.add_argument("--no-plot", action="store_true", help="跳过 matplotlib 绘图") + args = parser.parse_args() + + dynamics_path = (CASE_DIR / DYNAMICS_PATH).resolve() + input_dir = (CASE_DIR / INPUT_DIR).resolve() + output_dir = (CASE_DIR / OUTPUT_DIR).resolve() + config_path = (CASE_DIR / CONFIG_FILE).resolve() + + module = load_dynamics_module(dynamics_path) + module.run_case( + config_path=config_path, + runtime_base=CASE_DIR, + input_dir=input_dir, + output_dir=output_dir, + no_plot=args.no_plot, + ) + + +if __name__ == "__main__": + main() diff --git a/examples/case05/Readme.md b/examples/case05/Readme.md new file mode 100644 index 0000000..af29871 --- /dev/null +++ b/examples/case05/Readme.md @@ -0,0 +1,28 @@ +# case05: 一维原子链（60原子） + +60个原子沿 x 轴排列，相邻原子用弹簧连接。 + +## 物理设定 + +| 参数 | 值 | +|---|---| +| 原子数 | 60 | +| 排列 | 沿 x 轴等间距排列，间距为 1 | +| 约束 | 原子**只沿 z 方向**振动（fix_x=1, fix_y=1, fix_z=0） | +| 弹簧 | 劲度系数 k=1.0，原长 L₀=1.0 | +| 重力 | 无 | +| 万有引力 | 无 | +| 阻尼 | 无 | +| 算法 | leapfrog（蛙跳法，能量守恒） | + +## 初始条件 + +- 所有原子初始速度为零 +- **第1个原子**在 z 方向有初始位移，位移量 = 1 +- 其余原子初始 z=0 + +## 动力学行为 + +初始时刻第1个原子的 z 位移拉伸了它和第2个原子之间的弹簧，产生一个沿 z 方向的扰动。该扰动将以波的形式沿一维原子链传播，在链的两端反射，形成驻波叠加。 + +由于无阻尼，系统总能量守恒。 diff --git a/examples/case05/doc/index.html b/examples/case05/doc/index.html new file mode 100644 index 0000000..07267ac --- /dev/null +++ b/examples/case05/doc/index.html @@ -0,0 +1,471 @@ + + + + + +case05 — 一维原子链驱动力学模拟 | 物理原理 & 使用文档 + + + + + + + +
+

一维原子链驱动力学模拟

+

60 个原子沿 x 轴排列 · 弹簧连接 · z 方向受迫振动

+ case05 · examples/case05 +
+ +
+ + + + +
+

目录

+
    +
  1. 物理原理
    2. +
  2. 数值算法
    3. +
  3. 驱动力模型
    4. +
  4. 使用方法
    5. +
  5. 参数参考
    6. +
  6. 文件结构
    7. +
  7. 常见问题
    8. +
+
+ + + + +
+

一、物理原理

+ +
+

1.1 一维原子链

+

60 个原子沿 x 轴 等间距排列，原子间距为 1。相邻原子之间用 理想弹簧 连接，弹簧的劲度系数 k = 1.0，原长 L₀ = 1.0（与原子间距一致，初始状态弹簧无拉伸）。

+

每个原子被限制在 z 方向 自由振动，x 和 y 方向锁定（fix_x=1, fix_y=1, fix_z=0）。

+
+ +
+

1.2 弹簧力（胡克定律）

+

当原子 i 和 j 之间有弹簧连接时，原子 i 受到的弹簧力为：

+
+ F = −k · (d − L₀) · uij +
+

其中 d = |rj − ri| 为两原子间距离，uij 为从 i 指向 j 的单位向量。由于原子只在 z 方向振动，弹簧在 z 方向的分量是 几何非线性 的——对于小振幅近似，z 方向等效于一个三次方恢复力（FPU 型非线性）。

+
+ +
+

1.3 运动方程

+

对于第 i 个自由原子（非受驱），牛顿第二定律给出：

+
+ m · ai = Fispring + Fidriving +
+

本案例中 唯一的外力 来自驱动力（仅施加于原子 1）。无重力、无万有引力、无阻尼，系统总能量守恒。

+
+ +
+

1.4 波传播

+

原子 1 的受迫振动通过弹簧逐次传递给相邻原子，形成沿链传播的 横波。由于横向振动的几何非线性（弹簧大部分张力在 x 方向，z 方向的有效刚度远小于 1），波的传播速度较慢，且高阶频率成分会在链中产生复杂的非线性动力学行为（类似 FPU 回波现象）。

+
+
+ + + + +
+

二、数值算法

+ +
+

2.1 蛙跳法（Leapfrog / Velocity-Verlet）

+

采用能量守恒特性优异的 蛙跳法（二阶辛积分器），更新公式为：

+
+ v(t + ½Δt) = v(t) + ½ a(t) · Δt
+ r(t + Δt) = r(t) + v(t + ½Δt) · Δt
+ a(t + Δt) = F(r(t + Δt), v(t + ½Δt)) / m
+ v(t + Δt) = v(t + ½Δt) + ½ a(t + Δt) · Δt +
+

蛙跳法在长时间模拟中能量漂移极小（本案例验证 < 0.004%），适合无阻尼的保守系统。

+
+ +
+

2.2 时间步长与采样

+ + + + + + +
参数值说明
DT0.01 s积分步长（远小于 1/ω ≈ 0.16 s，满足稳定性条件）
T_total100 s总模拟时间 → NT = 10000 步
NSTEP50每 NSTEP 步取一帧用于动画 → 200 帧
methodleapfrog蛙跳法（Velocity-Verlet）
+
+ +
+

2.3 计算流程

+
+ 读入 coord.txt
connection.txt
bond.txt + → + 施加驱动力
（驱动原子 1） + → + 记录轨迹 + → + 蛙跳法
更新位置/速度 + → + 固定约束
（x, y 锁定） + → + 循环
NT 次 +
+

注意：驱动力在 每次积分前 施加，确保受驱原子的位置正确传递给弹簧力计算。

+
+
+ + + + +
+

三、驱动力模型

+ +
+

3.1 定义文件

+

驱动力由 input/driver.txt 定义，格式如下：

+ 
n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period
+1     0     0     5      0      0      1     0     0    90    all
+
+ +
+

3.2 数学公式

+

受驱原子的位置由下式决定（完全替换 coord.txt 中的初始坐标和固定约束）：

+
+ r(t) = A · cos(2πf · t + φ) +
+

速度由解析导数给出：

+
+ v(t) = −A · 2πf · sin(2πf · t + φ) +
+

其中 A = (amp_x, amp_y, amp_z)，f = (freq_x, freq_y, freq_z) 为不同方向的驱动频率，φ = (phi_x, phi_y, phi_z) 为相位（角度制，代码自动转换为弧度）。

+
+ +
+

3.3 本案例驱动参数

+ + + + + + +
参数值含义
amp_z5.0z 方向驱动振幅
freq_z1.0 Hz驱动频率（周期 1 s）
phi_z90°驱动相位 → z(0) = 5·cos(90°) = 0
periodall全程驱动，永不停止
+
+ z(t) = 5.0 · cos(2π · 1.0 · t + 90°) +
+
+ +
+

3.4 有限周期驱动

+

period 参数支持三种模式：

+
    +
  • all — 全程驱动
    • +
  • 数值 — 驱动指定周期数后 静止（冻结在最终位置，速度归零）。例如 period: 1 表示驱动 1 个完整周期后停止。
    • +
+
+ +
+

3.5 驱动与固定约束的关系

+

对于受驱原子（driver.txt 中 n 指定的原子），其在 coord.txt 中的初始坐标和 fix_x/fix_y/fix_z 约束被 完全忽略。原子的位置和速度完全由驱动力公式决定。

+
+
+ + + + +
+

四、使用方法

+ +
+

4.1 完整运行（模拟 + 动画）

+ 
cd examples/case05
+python run_dynamics.py
+

这步会依次执行：物理模拟 → 抽帧 → 打开 VisPy 3D 动画窗口。

+
+ +
+

4.2 仅查看已有结果

+

如果已经跑完模拟且生成了 output/display.txt，可以通过修改 input.txt 跳过计算，只开动画：

+ 
step_simulate:  0   # 跳过模拟
+step_sample:    0   # 跳过抽帧
+step_animation: 1   # 播放动画
+

然后运行：python run_dynamics.py

+
+ +
+

4.3 手动 3D 动画

+

也可以单独启动 VisPy 窗口：

+ 
python ../../draw.py output/
+
+ +
+

4.4 强制重新计算

+

修改参数后需要重新运行模拟时，设置：

+ 
force_calc: 1        # 忽略缓存，强制重新计算
+
+ +
+

4.5 动画交互

+ + + + + + + + +
操作效果
鼠标拖动旋转视角
滚轮缩放
左上角 reset 按钮复位视角到初始位置
左上角 info 按钮切换信息面板显示/隐藏
左上角 axes 按钮切换坐标轴显示/隐藏
Q / E 键画面绕视向旋转 -90° / +90°
+
+
+ + + + +
+

五、参数参考

+ +
+

5.1 input.txt 关键参数

+ + + + + + + + + + + + + +
参数默认值说明
gravity_field0均匀重力场（已关闭）
gravity_interaction0原子间万有引力（已关闭）
elastic_force1弹簧键力（已开启）
damping_force0阻尼（已关闭）
driving_force1驱动力开关（1=开启，需 driver.txt）
methodleapfrog数值积分方法
DT0.01积分步长 (s)
T_total100.0总模拟时间 (s)
NSTEP50抽帧步数间隔
enginepython计算引擎（python / c / cpp / fortran）
use_marker1渲染模式（0=Sphere 网格, 1=Marker GPU 实例化）
+
+ +
+

5.2 流程控制参数

+ + + + + + + +
参数01
step_simulate跳过模拟（加载已有轨迹）运行物理模拟
step_sample跳过抽帧从轨迹抽取显示帧
step_plot不生成图表生成轨迹/能量图
step_animation不启动动画自动打开 VisPy 3D 窗口
force_calc自动检测缓存强制重新计算
+
+
+ + + + +
+

六、文件结构

+ + 
case05/
+├── input/
+│   ├── input.txt          # 主配置文件（YAML 格式）
+│   ├── coord.txt          # 原子坐标（60 个原子）
+│   ├── connection.txt     # 弹簧连接关系（59 条键）
+│   ├── bond.txt           # 弹簧参数（k=1.0, L₀=1.0）
+│   └── driver.txt       # 驱动力定义（本案例新增）
+├── output/
+│   ├── trajectory.txt     # 全量轨迹数据（10000 步 × 60 原子）
+│   ├── display.txt        # 抽帧后的动画数据（200 帧 × 60 原子）
+│   ├── dynamics.log       # 计算日志
+│   └── animation.log      # 动画启动日志（闪退时排查用）
+├── doc/
+│   └── index.html         # 本文档
+├── Readme.md              # 案例简介
+└── run_dynamics.py        # 案例运行入口
+
+ + + + +
+

七、常见问题

+ +
+

7.1 动画窗口闪退

+

如果 VisPy 窗口一闪就消失，请检查：

+
    +
  • output/animation.log 中是否有错误信息
    • +
  • output/display.txt 是否存在（需先跑 step_sample: 1）
    • +
+
+ +
+

7.2 原子不振动

+

可能原因：

+
    +
  • NSTEP 过大：抽帧间隔大于驱动周期的一半时，动画会丢失振动细节。建议 NSTEP ≤ 1/(freq · DT · 10)
    • +
  • 相位 φ 使采样点落在零值：试试 phi_z: 0 让原子在 t=0 处于振幅峰值
    • +
  • 确认 driving_force: 1 且 driver.txt 中 amp_z 不为 0
    • +
+
+ +
+

7.3 渲染性能慢

+

原子数多时动画卡顿：

+
    +
  • 设置 use_marker: 1（使用 GPU 实例化渲染替代独立网格球体）
    • +
  • 增大 NSTEP 减少动画帧数
    • +
+
+
+ +
+ + + +
+ + diff --git a/examples/case05/input/bond.txt b/examples/case05/input/bond.txt new file mode 100644 index 0000000..f33e17f --- /dev/null +++ b/examples/case05/input/bond.txt @@ -0,0 +1,2 @@ +bond_name k rest_length +k1 50.0 1.0 diff --git a/examples/case05/input/connection.txt b/examples/case05/input/connection.txt new file mode 100644 index 0000000..c0d1094 --- /dev/null +++ b/examples/case05/input/connection.txt @@ -0,0 +1,60 @@ +n1 n2 bond_name +1 2 k1 +2 3 k1 +3 4 k1 +4 5 k1 +5 6 k1 +6 7 k1 +7 8 k1 +8 9 k1 +9 10 k1 +10 11 k1 +11 12 k1 +12 13 k1 +13 14 k1 +14 15 k1 +15 16 k1 +16 17 k1 +17 18 k1 +18 19 k1 +19 20 k1 +20 21 k1 +21 22 k1 +22 23 k1 +23 24 k1 +24 25 k1 +25 26 k1 +26 27 k1 +27 28 k1 +28 29 k1 +29 30 k1 +30 31 k1 +31 32 k1 +32 33 k1 +33 34 k1 +34 35 k1 +35 36 k1 +36 37 k1 +37 38 k1 +38 39 k1 +39 40 k1 +40 41 k1 +41 42 k1 +42 43 k1 +43 44 k1 +44 45 k1 +45 46 k1 +46 47 k1 +47 48 k1 +48 49 k1 +49 50 k1 +50 51 k1 +51 52 k1 +52 53 k1 +53 54 k1 +54 55 k1 +55 56 k1 +56 57 k1 +57 58 k1 +58 59 k1 +59 60 k1 diff --git a/examples/case05/input/coord.txt b/examples/case05/input/coord.txt new file mode 100644 index 0000000..f34b022 --- /dev/null +++ b/examples/case05/input/coord.txt @@ -0,0 +1,61 @@ +n mass radius x y z vx vy vz fix_x fix_y fix_z +1 1 0.1 0 0 1 0 0 0 1 1 0 +2 1 0.1 1 0 0 0 0 0 1 1 0 +3 1 0.1 2 0 0 0 0 0 1 1 0 +4 1 0.1 3 0 0 0 0 0 1 1 0 +5 1 0.1 4 0 0 0 0 0 1 1 0 +6 1 0.1 5 0 0 0 0 0 1 1 0 +7 1 0.1 6 0 0 0 0 0 1 1 0 +8 1 0.1 7 0 0 0 0 0 1 1 0 +9 1 0.1 8 0 0 0 0 0 1 1 0 +10 1 0.1 9 0 0 0 0 0 1 1 0 +11 1 0.1 10 0 0 0 0 0 1 1 0 +12 1 0.1 11 0 0 0 0 0 1 1 0 +13 1 0.1 12 0 0 0 0 0 1 1 0 +14 1 0.1 13 0 0 0 0 0 1 1 0 +15 1 0.1 14 0 0 0 0 0 1 1 0 +16 1 0.1 15 0 0 0 0 0 1 1 0 +17 1 0.1 16 0 0 0 0 0 1 1 0 +18 1 0.1 17 0 0 0 0 0 1 1 0 +19 1 0.1 18 0 0 0 0 0 1 1 0 +20 1 0.1 19 0 0 0 0 0 1 1 0 +21 1 0.1 20 0 0 0 0 0 1 1 0 +22 1 0.1 21 0 0 0 0 0 1 1 0 +23 1 0.1 22 0 0 0 0 0 1 1 0 +24 1 0.1 23 0 0 0 0 0 1 1 0 +25 1 0.1 24 0 0 0 0 0 1 1 0 +26 1 0.1 25 0 0 0 0 0 1 1 0 +27 1 0.1 26 0 0 0 0 0 1 1 0 +28 1 0.1 27 0 0 0 0 0 1 1 0 +29 1 0.1 28 0 0 0 0 0 1 1 0 +30 1 0.1 29 0 0 0 0 0 1 1 0 +31 1 0.1 30 0 0 0 0 0 1 1 0 +32 1 0.1 31 0 0 0 0 0 1 1 0 +33 1 0.1 32 0 0 0 0 0 1 1 0 +34 1 0.1 33 0 0 0 0 0 1 1 0 +35 1 0.1 34 0 0 0 0 0 1 1 0 +36 1 0.1 35 0 0 0 0 0 1 1 0 +37 1 0.1 36 0 0 0 0 0 1 1 0 +38 1 0.1 37 0 0 0 0 0 1 1 0 +39 1 0.1 38 0 0 0 0 0 1 1 0 +40 1 0.1 39 0 0 0 0 0 1 1 0 +41 1 0.1 40 0 0 0 0 0 1 1 0 +42 1 0.1 41 0 0 0 0 0 1 1 0 +43 1 0.1 42 0 0 0 0 0 1 1 0 +44 1 0.1 43 0 0 0 0 0 1 1 0 +45 1 0.1 44 0 0 0 0 0 1 1 0 +46 1 0.1 45 0 0 0 0 0 1 1 0 +47 1 0.1 46 0 0 0 0 0 1 1 0 +48 1 0.1 47 0 0 0 0 0 1 1 0 +49 1 0.1 48 0 0 0 0 0 1 1 0 +50 1 0.1 49 0 0 0 0 0 1 1 0 +51 1 0.1 50 0 0 0 0 0 1 1 0 +52 1 0.1 51 0 0 0 0 0 1 1 0 +53 1 0.1 52 0 0 0 0 0 1 1 0 +54 1 0.1 53 0 0 0 0 0 1 1 0 +55 1 0.1 54 0 0 0 0 0 1 1 0 +56 1 0.1 55 0 0 0 0 0 1 1 0 +57 1 0.1 56 0 0 0 0 0 1 1 0 +58 1 0.1 57 0 0 0 0 0 1 1 0 +59 1 0.1 58 0 0 0 0 0 1 1 0 +60 1 0.1 59 0 0 0 0 0 1 1 0 diff --git a/examples/case05/input/driver.txt b/examples/case05/input/driver.txt new file mode 100644 index 0000000..79857fc --- /dev/null +++ b/examples/case05/input/driver.txt @@ -0,0 +1,2 @@ +n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period +1 0 0 1 0 0 0.1 0 0 90 all \ No newline at end of file diff --git a/examples/case05/input/input.txt b/examples/case05/input/input.txt new file mode 100644 index 0000000..039356c --- /dev/null +++ b/examples/case05/input/input.txt @@ -0,0 +1,96 @@ +# 物理模拟参数配置 +# 格式：YAML +# 用法：python run_dynamics.py + +# ── 流程控制 ────────────────────────────────── +# 每步用 0/1 单独开关，1=执行，0=跳过 +# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧 +step_simulate: 1 # 运行物理模拟 → output/trajectory.txt +step_sample: 1 # 抽帧 → output/display.txt +step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_animation: 1 # 自动播放 VisPy 3D 动画窗口（需安装 vispy） +force_calc: 0 # 强制重新计算：1=跳过缓存强算，0=自动使用已有输出 + +# ── 计算引擎 ────────────────────────────────── +# 可选: python, c, cpp, fortran, java +engine: python # 默认使用 Python 引擎 + +# ── 盒子 ────────────────────────────────────── +box_a: 80.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内 + +# ── 初始构型 ────────────────────────────────── +# 坐标文件格式： +# 第一行：n mass radius x y z vx vy vz fix_x fix_y fix_z +# 后续行：原子序号 质量 半径 x y z vx vy vz fix_x fix_y fix_z +coord_file: input/coord.txt +connection_file: input/connection.txt +bond_file: input/bond.txt +driver_file: input/driver.txt # 驱动力定义文件（driving_force=1 时生效） + +# 绘图/动画展示的原子序号（对应 coord_file 第一列 n） +plot_atom: 1 + +# ── 物理参数 ────────────────────────────────── +# 三个方向分量分别对应 x, y, z +G: [0.0, 0.0, 0.0] # 重力场分量 (m/s²) +B: [0.0, 0.0, 0.0] # 阻尼分量 + +# ── 力开关（0=关闭, 1=开启）────────────────── +gravity_field: 0 # 均匀重力场 (G) +gravity_interaction: 0 # 原子间万有引力 +elastic_force: 1 # 弹簧键力 +damping_force: 0 # 阻尼 (B) +driving_force: 1 # 驱动力（需 driver_file 定义） +# +gravity_strength: 1.0 # 万有引力强度（仅 gravity_interaction=1 时有效） + +# ── 数值算法 ────────────────────────────────── +# 可选： +# explicit_euler 显式欧拉法 +# implicit_euler 隐式欧拉法 +# midpoint 中点法 +# leapfrog 蛙跳法 +method: leapfrog + +# ── 步骤控制 ────────────────────────────────── +# 以下参数控制哪些步骤被执行和保存 + +# 预热步数：模拟开始时跳过不保存的步数（用于稳定初始状态） +warmup_steps: 0 # 默认 0（立即开始记录） + +# 总模拟时间（秒），程序自动计算 NT = T_total / DT +# 如果同时指定了 NT，以 NT 为准 +T_total: 100.0 + +# 抽帧间隔（每 NSTEP 步取一帧用于动画） +NSTEP: 50 + +# ── 时间步长 ────────────────────────────────── +DT: 0.001 # 时间步长 (s) + +# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧 +sample_start: null # null 表示从头开始（帧索引从 0 起） +sample_end: null # null 表示到末尾 + + + +# ── 渲染方式 ────────────────────────────────── +# 3D 动画中原子渲染方式： +# 0 = Sphere （网格球体，效果精细，原子数少时推荐） +# 1 = Marker （GPU 实例化点，原子数多时性能更佳） +use_marker: 1 + +# ── 显示参数 ────────────────────────────────── +# 盒子透明度：单个数值（统一）或 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.20 # R 分量 (0~1) +ball_color_g: 0.60 # G 分量 +ball_color_b: 0.90 # B 分量 + +# 盒子面颜色 +box_color_r: 0.80 +box_color_g: 0.80 +box_color_b: 0.85 diff --git a/examples/case05/run_dynamics.py b/examples/case05/run_dynamics.py new file mode 100644 index 0000000..37c77d2 --- /dev/null +++ b/examples/case05/run_dynamics.py @@ -0,0 +1,54 @@ +""" +Case runner for Dynamics case05 — 1D atomic chain. + +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()