case05 — 一维原子链驱动力学模拟 | 物理原理 & 使用文档

一、物理原理

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₀) · u_ij

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

1.3 运动方程

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

m · a_i = F_i^spring + F_i^driving

本案例中 唯一的外力 来自驱动力（仅施加于原子 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 时间步长与采样

参数	值	说明
DT	0.01 s	积分步长（远小于 1/ω ≈ 0.16 s，满足稳定性条件）
T_total	100 s	总模拟时间 → NT = 10000 步
NSTEP	50	每 NSTEP 步取一帧用于动画 → 200 帧
method	leapfrog	蛙跳法（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_z	5.0	z 方向驱动振幅
freq_z	1.0 Hz	驱动频率（周期 1 s）
phi_z	90°	驱动相位 → z(0) = 5·cos(90°) = 0
period	all	全程驱动，永不停止

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_field	0	均匀重力场（已关闭）
gravity_interaction	0	原子间万有引力（已关闭）
elastic_force	1	弹簧键力（已开启）
damping_force	0	阻尼（已关闭）
driving_force	1	驱动力开关（1=开启，需 driver.txt）
method	leapfrog	数值积分方法
DT	0.01	积分步长 (s)
T_total	100.0	总模拟时间 (s)
NSTEP	50	抽帧步数间隔
engine	python	计算引擎（python / c / cpp / fortran）
use_marker	1	渲染模式（0=Sphere 网格, 1=Marker GPU 实例化）

5.2 流程控制参数

参数	0	1
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 减少动画帧数

目录