第7章 第一性原理分子动力学（Ab Initio Molecular Dynamics with OpenMX）

本章讲解如何用 OpenMX 4.0 做第一性原理分子动力学（first-principles / ab initio MD, AIMD）：NVE 微正则系综、速度标度（NVT_VS、NVT_VS2、NVT_VS4）与 Nosé–Hoover（NVT_NH）恒温，以及 4.0 新增的恒温恒压 NPT_VS_PR / NPT_VS_WV。我们会逐个解释 MD.TimeStep、<MD.TempControl、NH.Mass.HeatBath、NPT.Mass.Barostat、<MD.CellPressureControl、<MD.Init.Velocity 等关键词，读懂 .ene、.md、.md2、.dat# 输出，学会续算，并用 Python 检查能量守恒、温度，计算径向分布函数（RDF）与均方位移（MSD），最后用 dpdata 把 OpenMX 轨迹转换成 DeePMD 等机器学习势的训练集。

---

1. 本章目标与适用场景

学完本章，你应当能够：

• 根据体系选择系综（NVE / NVT / NPT）与 MD.Type，并按最快振动周期确定 MD.TimeStep；
• 写出两种格式的 <MD.TempControl（速度标度与 Nosé–Hoover 格式不同），设计升温—保温方案；
• 理解 MD.maxIter 是"绝对步数"，能用自动生成的 System.Name.dat# 正确续算；
• 从 System.Name.ene 判断能量是否守恒、温度是否达到目标；从 System.Name.npt.ene 读取体积与压力；
• 用 Python 解析 System.Name.md（4.0 格式含坐标、力、速度、净电荷、磁矩），计算 g(r)、配位数、MSD；
• 用 dpdata 的 openmx/md 格式把 AIMD 数据导出为 DeePMD 训练数据，或交给 NEP 等势函数的训练流程。

需求	推荐 MD.Type	本章算例
检验能量守恒、计算振动谱、小分子短时间动力学	NVE	§4 H₂O 分子
快速把体系加热 / 冷却到目标温度（平衡阶段）	NVT_VS	§5 Si64
需要正则系综统计（生产阶段）	NVT_NH	§6 Si64
不同区域设定不同温度（如热流、局域加热）	NVT_VS4（多热浴）	§5.4
有限温度 + 外压下的密度、相变、热膨胀	NPT_VS_PR / NPT_VS_WV（4.0 新增）	§7 Si8

Note

本章默认你已掌握 第5章 的 SCF 与 k 点设置和 第6章 的结构优化。MD 本质上是"成千上万次单点计算"：每一步的截断能、k 点、SCF 收敛精度都会直接影响轨迹质量，而计算量 ≈ 步数 × 每步 SCF 时间，所以输入参数要在精度与成本之间认真权衡。

---

2. 原理速览

2.1 Born–Oppenheimer MD 与速度 Verlet 积分

第一性原理 MD 在每个时间步都求解 Kohn–Sham 方程，由 Hellmann–Feynman 定理（加上 OpenMX 中基函数随原子移动带来的 Pulay 修正）得到原子受力 \(\mathbf{F}_I=-\partial E/\partial \mathbf{R}_I\)，再按牛顿方程推进原子：

\[ \mathbf{R}_I(t+\Delta t)=\mathbf{R}_I(t)+\mathbf{v}_I(t)\Delta t+\frac{\mathbf{F}_I(t)}{2M_I}\Delta t^2,\qquad \mathbf{v}_I(t+\Delta t)=\mathbf{v}_I(t)+\frac{\mathbf{F}_I(t)+\mathbf{F}_I(t+\Delta t)}{2M_I}\Delta t \]

这就是速度 Verlet 算法。与经验势 MD 相比，AIMD 的力来自电子结构计算，因此能自然描述化学键的断裂与生成（官方手册第 20 章开头也强调了这一点）。

2.2 温度与动能

OpenMX 用全部原子的动能定义瞬时温度：

\[ T=\frac{2U_{\mathrm{kc}}}{3Nk_\mathrm{B}},\qquad U_{\mathrm{kc}}=\sum_I \tfrac12 M_I v_I^2 \]

\(N\) 个原子体系的瞬时温度存在固有涨落，相对涨落约为 \(\sqrt{2/(3N)}\)：3 原子分子约 ±47%，64 原子约 ±10%。所以小体系的温度曲线"抖得厉害"是正常的，要看时间平均值。

2.3 恒温方法

速度标度（velocity scaling, NVT_VS）：每隔若干步把所有原子速度乘以同一个因子

\[ s=\sqrt{\frac{T_\mathrm{given}+(T_\mathrm{calc}-T_\mathrm{given})\,\alpha}{T_\mathrm{calc}}} \]

其中 \(\alpha\) 是 <MD.TempControl 第 4 列的"缩放参数"：\(\alpha=0\) 表示一次性把温度拉到目标值，\(\alpha\to 1\) 表示几乎不缩放。速度标度简单、稳健，适合升温、降温和平衡阶段，但它不产生严格的正则分布，涨落相关的物理量（比热等）不要用它来算。

Nosé–Hoover（NVT_NH）：引入一个热浴自由度 \(\zeta\) 及其"质量" \(Q\)（NH.Mass.HeatBath），运动方程变为 \(M_I\ddot{\mathbf{R}}_I=\mathbf{F}_I-\zeta M_I\dot{\mathbf{R}}_I\)。对遍历性良好的体系，它给出正则系综；扩展体系的 Nosé–Hoover 哈密顿量守恒，可用于检验积分质量。

2.4 恒压方法（4.0 新增）

NPT_VS_PR 和 NPT_VS_WV 都用速度标度控温，晶胞自由度分别按 Parrinello–Rahman 和 Wentzcovitch 的运动方程演化：晶格矢量受"内压（由应力张量计算） − 外压"驱动，惯性由活塞质量 \(W\)（NPT.Mass.Barostat）决定。\(W\) 越小晶胞响应越快、体积振荡越剧烈；\(W\) 越大响应越慢。两种方法的 \(W\) 定义和量级不同，不能互相套用。

2.5 时间步长的选择

时间步长应为体系最快振动周期的 1/10–1/20。OpenMX 默认 MD.TimeStep 0.5（fs），经验值如下：

体系	推荐 MD.TimeStep（fs）	说明
含 H（O–H、C–H 伸缩周期约 9–11 fs）	0.25–0.5	官方甘氨酸示例用 0.2 fs
把 H 质量改为氘（2.0 amu，见 §3.3）	0.5–1.0	动力学变成氘代体系的动力学，但结构/热力学平衡性质基本不变
C、N、O 等轻元素，不含 H	0.5–1.0	—
Si、Al 等第三周期元素	1.0–2.0	官方 Si_NPT_VS_PR.dat 用 1.0 fs
过渡金属、重元素	2.0–3.0	高温下适当减小

检验方法：先用 NVE 跑几百步，看守恒量 Utot + Ukc 是否漂移（§4.5）。

---

3. MD 相关关键词总览

3.1 系综与步数

关键词	取值 / 默认	说明
MD.Type	NOMD（默认）、NVE、NVT_VS、NVT_VS2、NVT_VS4、NVT_NH、NPT_VS_PR、NPT_VS_WV	后两者为 OpenMX 4.0 新增并有文档；大小写不敏感
MD.maxIter	默认 1	MD 的总步数（绝对值）。续算时 .dat# 中会写入 MD.Current.Iter，程序在"当前步 > MD.maxIter"时停止，所以续算必须把 MD.maxIter 改大（§8.3）
MD.TimeStep	默认 0.5	时间步长，单位 fs
scf.ExtCharge.History	手册推荐 2	MD/优化中用前几步的电荷密度外推作为下一步的初猜；手册说明更大的值反而可能不稳定

Note

在 4.0 的源码（Input_std.c）里还能看到 NVT_Langevin、NPT_NH_PR、NPT_NH_WV 等选项名，但 4.0 用户手册没有为它们写文档。本教程只使用手册中有说明的类型；如需使用未文档化的功能，请先阅读源码并做充分测试。

3.2 温度、压力与热浴

关键词	适用	默认	说明
<MD.TempControl	NVT_VS/VS2/VS4、NPT_VS_*	—	4 列：步数 间隔 温度(K) α，见 §5.2
<MD.TempControl	NVT_NH	—	2 列：步数 温度(K)，见 §6.2
NH.Mass.HeatBath	NVT_NH	20.0	热浴质量 \(Q\)，量纲为 长度²×质量，单位 Bohr²·amu（C 主同位素质量取 12.0）
MD.num.AtomGroup / <MD.AtomGroup	NVT_VS4	1	原子分组数与"原子序号 组号"映射，组号从 1 开始连续编号
NPT.Mass.Barostat	NPT_VS_*	自动	活塞质量 \(W\)（amu）；不写时由程序自动计算并在标准输出打印 A default NPT.Mass.Barostat is set to ...
<MD.CellPressureControl	NPT_VS_*	—	2 列：步数 外压(GPa)，步间线性插值；正值表示压缩
NPT.LatticeRestriction	NPT_VS_*	none	none / cubic / ortho：限制晶胞形状
MD.TempTolerance	NPT_VS_*	100	单位 K；$
MD.applied.pressure.flag	NPT_VS_*	1 1 1	各晶轴是否施加压力；正交晶胞下可做单轴加压，六方可设 0 0 1

3.3 初速度、质量与约束

初速度（单位 m/s，每个原子一行，可与 <MD.Fixed.XYZ 同时使用）：

<MD.Init.Velocity
   1    3000.000    0.0   0.0
   2   -3000.000    0.0   0.0
MD.Init.Velocity>

不给初速度时，所有原子从静止开始，体系只能靠"偏离平衡位置的受力"获得动能。

修改原子质量：在 <Definition.of.Atomic.Species 的第 4 列写质量（amu），不写则用程序内置值（定义在 SetPara_DFT.c 的 Set_Atom_Weight()）：

<Definition.of.Atomic.Species
  H   H6.0-s2p1     H_PBE19    2.0     # treat H as deuterium
  O   O6.0-s2p2d1   O_PBE19
Definition.of.Atomic.Species>

固定原子：<MD.Fixed.XYZ 对优化和 MD 都有效（格式见第6章），称为"约束 MD"。注意：约束 MD 与 NPT 不兼容。

---

4. 算例一：H₂O 分子 NVE（检验能量守恒）

完整文件见 ../examples/07-md-h2o-nve/README.md。

4.1 为什么从 NVE 开始

NVE 中没有任何外加控制，守恒量 \(E_\mathrm{cons}=U_\mathrm{tot}+U_\mathrm{kc}\) 是否守恒直接反映了"时间步长 + SCF 精度 + 实空间网格"组合的质量。任何正式的 NVT/NPT 研究之前，都值得先用 NVE 跑几百步做这个体检。我们让一个 O–H 键略微拉长（1.02 Å，另一个 0.97 Å），从静止出发，分子会自发振动。

4.2 完整输入 H2O_NVE.dat

# H2O molecule, NVE molecular dynamics (OpenMX 4.0)
System.CurrrentDirectory    ./
System.Name                 H2O_NVE
DATA.PATH                   /opt/openmx4.0/DFT_DATA19   # EDIT: dir containing PAO/ and VPS/
level.of.stdout             1
level.of.fileout            0          # do not write cube files

Species.Number              2
<Definition.of.Atomic.Species
  H   H6.0-s2p1      H_PBE19
  O   O6.0-s2p2d1    O_PBE19
Definition.of.Atomic.Species>

Atoms.Number                3
Atoms.SpeciesAndCoordinates.Unit   Ang
<Atoms.SpeciesAndCoordinates
  1  O   5.00000   5.00000   5.30000   3.0  3.0
  2  H   5.80650   5.00000   4.67554   0.5  0.5
  3  H   4.23303   5.00000   4.70615   0.5  0.5
Atoms.SpeciesAndCoordinates>

Atoms.UnitVectors.Unit      Ang
<Atoms.UnitVectors
  10.0   0.0   0.0
   0.0  10.0   0.0
   0.0   0.0  10.0
Atoms.UnitVectors>

scf.XcType                  GGA-PBE
scf.SpinPolarization        off
scf.ElectronicTemperature   300.0
scf.energycutoff            250.0
scf.maxIter                 100
scf.EigenvalueSolver        Cluster
scf.Kgrid                   1 1 1
scf.Mixing.Type             RMM-DIIS
scf.Init.Mixing.Weight      0.30
scf.Min.Mixing.Weight       0.001
scf.Max.Mixing.Weight       0.40
scf.Mixing.History          7
scf.Mixing.StartPulay       5
scf.criterion               1.0e-9

MD.Type                     NVE
MD.maxIter                  1000
MD.TimeStep                 0.5
scf.ExtCharge.History       2

4.3 关键词解释

关键词	本例取值	为什么这样设
System.Name	H2O_NVE	与文件名主干一致，后面 dpdata 读取时要求 <前缀>.dat 与 <前缀>.md 成对存在
level.of.fileout	0	MD 每步都写 cube 文件会迅速占满磁盘；0 只输出 .out/.md/.ene/.xyz 等
物种名	H、O	物种名直接用元素符号：.md 文件第 1 列写的是元素符号，dpdata 按物种名匹配
H6.0-s2p1 / O6.0-s2p2d1	DFT_DATA19 的 Standard 基组	见第4章；MD 也可以先用 Quick 基组试跑
上/下自旋电子数	O：3.0 3.0；H：0.5 0.5	两列之和等于 VPS 价电子数（O_PBE19 为 6，H_PBE19 为 1）
scf.EigenvalueSolver	Cluster	孤立分子，用 10 Å 盒子 + Γ 点
scf.energycutoff	250 Ry	实空间网格越粗，"egg-box 效应"（能量随原子相对网格位置的人为起伏）越明显，会破坏能量守恒
scf.criterion	1.0e-9 Hartree	SCF 不够收敛时力有噪声，NVE 中表现为能量漂移
MD.TimeStep	0.5 fs	O–H 伸缩周期约 9 fs，每个周期约 18 步
MD.maxIter	1000	500 fs，FFT 频率分辨率约 1/(500 fs) ≈ 67 cm⁻¹

4.4 运行

cd OpenMX/examples/07-md-h2o-nve
rm -f H2O_NVE.ene H2O_NVE.md      # .ene/.md are opened in append mode, remove old ones
mpirun -np 4 openmx H2O_NVE.dat -nt 1 | tee H2O_NVE.std

Warning

OpenMX 以追加模式写 System.Name.ene 和 System.Name.md（源码 iterout_md.c、iterout.c 中用 fopen(..., "a")）。用同一个 System.Name 从头重跑时，新数据会接在旧数据后面，分析脚本会得到"时间倒流"的轨迹。从头重跑前请先删除或改名旧文件；续算时这一特性反而很方便（§8.3）。

运行时监控：

grep -c "time=" H2O_NVE.md                       # number of finished MD steps
tail -n 3 H2O_NVE.ene | awk '{print $1, $2, $15, $16, $19}'   # step, time, Utot, Utot+Ukc, T
grep "scf_conv=0" H2O_NVE.md | head              # any MD step whose SCF did not converge?

4.5 System.Name.ene 的列与能量守恒判据

.ene 每个 MD 步一行、无表头。官方手册（第 20.1–20.3 节）列出的常用列如下（完整定义见源码 iterout_md.c）：

列	含义	单位
1	MD 步	—
2	MD 时间	fs
14	原子动能 \(U_\mathrm{kc}\)	Hartree
15	DFT 总能 \(U_\mathrm{tot}\)	Hartree
16	\(U_\mathrm{tot}+U_\mathrm{kc}\)（NVE 守恒量）	Hartree
17	费米能（化学势）	Hartree
18	目标温度	K
19	计算温度	K
22	Nosé–Hoover 哈密顿量（NVT_NH 时有意义）	Hartree
30	最大原子力（据源码 iterout_md.c）	Hartree/Bohr

用一行 awk 粗看漂移和涨落：

awk 'NR==1{e0=$16} {d=$16-e0; if(d>mx)mx=d; if(d<mn)mn=d} END{printf "E_cons span = %.3e Ha (%.3f meV)\n", mx-mn, (mx-mn)*27211.386}' H2O_NVE.ene

更完整的检查用示例目录中的 nve_check.py：线性拟合守恒量得到漂移速率，并对 O–H 键长做傅里叶变换得到振动频率：

python3 nve_check.py H2O_NVE --skip 0

（示意输出，数值仅供参考）

steps = 1000, dt = 0.500 fs
<T> = 1xx.x K
E_cons drift = x.xxx meV/ps, peak-to-peak = x.xxx meV
Ukc peak-to-peak = xxx.x meV
bond 1-2       mean =   0.98xx  spectrum peaks near 3xxx cm^-1
bond 1-3       mean =   0.98xx  spectrum peaks near 3xxx cm^-1
angle 2-1-3    mean = 10x.xxxx  spectrum peaks near 1xxx cm^-1

判据（经验性）：

• 守恒量的峰峰值应远小于动能的峰峰值（例如不到其 1–2%），且没有单调漂移；
• 若漂移明显：先减小 MD.TimeStep，再收紧 scf.criterion、提高 scf.energycutoff，并检查 .md 中是否有 scf_conv=0；
• O–H 伸缩峰应落在约 3500–3900 cm⁻¹ 范围（实验气相 H₂O 的对称/反对称伸缩基频约 3657/3756 cm⁻¹；PBE 与实验的差异、较大振幅带来的非谐性以及 500 fs 轨迹约 67 cm⁻¹ 的分辨率都会造成偏差，只需量级吻合）；H–O–H 弯曲峰在约 1500–1650 cm⁻¹（实验约 1595 cm⁻¹）。

4.6 时间步长测试与氘代技巧

for dt in 0.25 0.5 1.0; do
  sed -e "s/^System.Name .*/System.Name   H2O_dt${dt}/" \
      -e "s/^MD.TimeStep .*/MD.TimeStep   ${dt}/" \
      -e "s/^MD.maxIter .*/MD.maxIter   $(python3 -c "print(int(500/${dt}))")/" \
      H2O_NVE.dat > H2O_dt${dt}.dat
  mpirun -np 4 openmx H2O_dt${dt}.dat -nt 1 > H2O_dt${dt}.std
done

比较三条守恒量曲线：通常 1.0 fs 对含 H 体系已明显漂移。若想用 1.0 fs，可以在物种定义第 4 列把 H 质量设为 2.0（§3.3），这在研究结构、密度、RDF 等平衡性质时很常用，但扩散系数、振动频率等动力学量就对应氘代体系了。

4.7 给定初始温度：<MD.Init.Velocity

OpenMX 不会按麦克斯韦分布自动生成初速度（不给 <MD.Init.Velocity 时原子从静止出发）。示例目录提供 gen_init_velocity.py：读取 .dat 中的物种与坐标，按目标温度抽取麦克斯韦–玻尔兹曼速度、去掉质心速度，并精确缩放到 \(T=2U_\mathrm{kc}/(3Nk_\mathrm{B})\)，输出可直接粘贴进输入文件的块：

python3 gen_init_velocity.py H2O_NVE.dat --temp 300 --seed 1 >> H2O_NVE.dat

（示意输出，数值仅供参考）

# generated by gen_init_velocity.py: T = 300.0 K, seed = 1, unit = m/s
<MD.Init.Velocity
    1     -301.2345     120.9876     -88.1234
    ...
MD.Init.Velocity>

Tip

对含 3 个原子的分子，即使初速度对应 300 K，平均温度也会因为动能与势能之间的均分而下降到约一半（初始结构在势能极小附近时）。这是 NVE 的正常行为；需要真正"恒温"的统计请用 NVT。

---

5. 算例二：64 原子 Si 的 NVT 速度标度（NVT_VS）

完整文件见 ../examples/07-md-si64-nvt-npt/README.md。

5.1 输入要点

64 原子 Si 是 2×2×2 个金刚石惯用胞（边长 10.86 Å）。坐标较长，下面只列出全部关键词和前 3 行坐标，完整的 64 行坐标在示例文件 Si64_NVT_VS.dat 中：

System.CurrrentDirectory    ./
System.Name                 Si64_NVT_VS
DATA.PATH                   /opt/openmx4.0/DFT_DATA19   # EDIT
level.of.stdout             1
level.of.fileout            0

Species.Number              1
<Definition.of.Atomic.Species
  Si   Si7.0-s2p1d1   Si_PBE19
Definition.of.Atomic.Species>

Atoms.Number                64
Atoms.SpeciesAndCoordinates.Unit   Ang
<Atoms.SpeciesAndCoordinates
    1  Si   -0.00069   0.01381   0.01109   2.0 2.0
    2  Si   -0.00623   2.76255   2.58901   2.0 2.0
    3  Si    2.71293  -0.00140   2.69739   2.0 2.0
  ... (64 lines in total, see the example file)
Atoms.SpeciesAndCoordinates>

Atoms.UnitVectors.Unit      Ang
<Atoms.UnitVectors
  10.86   0.00   0.00
   0.00  10.86   0.00
   0.00   0.00  10.86
Atoms.UnitVectors>

scf.XcType                  GGA-PBE
scf.SpinPolarization        off
scf.ElectronicTemperature   300.0
scf.energycutoff            200.0
scf.maxIter                 100
scf.EigenvalueSolver        Band
scf.Kgrid                   2 2 2
scf.Mixing.Type             RMM-DIISK
scf.Init.Mixing.Weight      0.10
scf.Min.Mixing.Weight       0.001
scf.Max.Mixing.Weight       0.30
scf.Mixing.History          20
scf.Mixing.StartPulay       8
scf.criterion               1.0e-7

MD.Type                     NVT_VS
MD.maxIter                  1500
MD.TimeStep                 2.0
<MD.TempControl
  3
    10    1    300.0   0.0
   400    5   1000.0   0.0
  1500   10   1000.0   0.3
MD.TempControl>

关键词	取值	说明
坐标	理想金刚石格点 + 约 0.05 Å 的随机位移	重要：理想晶格受力为零、又没有初速度时，动能恒为 0，速度标度因子 \(s\) 的分母为零，体系"永远不动"甚至得到 NaN。随机位移（或 <MD.Init.Velocity）给了体系初始动能来源
Si7.0-s2p1d1	Quick 基组	每原子 10 个基函数，MD 中性价比高；正式研究前应与 Standard 基组比较 RDF/压力
scf.Kgrid	2 2 2	10.86 Å 的立方胞；4.0 默认 scf.Generation.Kpoint GCenter，偶数网格也包含 Γ 与 TRIM 点。只想快速试跑可用 1 1 1
scf.criterion	1.0e-7	NVT 下的守恒要求比 NVE 宽松，1e-7 是常用折中；做 NVE 或需要精确力时收紧到 1e-8–1e-9
MD.TimeStep	2.0 fs	Si 最高光学声子频率约 15.5 THz（周期约 65 fs），2 fs 足够

5.2 <MD.TempControl（速度标度格式）逐列解释

<MD.TempControl
  3                          # number of following lines
    10    1    300.0   0.0   # until step 10: scale every 1 step,  T = 300 K,  alpha = 0.0
   400    5   1000.0   0.0   # steps 11-400: scale every 5 steps, T ramps 300 -> 1000 K
  1500   10   1000.0   0.3   # steps 401-1500: every 10 steps, hold 1000 K, softer scaling
MD.TempControl>

• 第 1 列：该段结束的 MD 步；第 2 列：该段内每隔多少步做一次速度标度；第 3 列：该段结束时的目标温度；第 4 列：缩放参数 \(\alpha\)。
• 段内的目标温度按"上一段末温度 → 本段末温度"线性插值，因此第 2 行实现了 300 → 1000 K 的线性升温。第一段之前的温度取第一行的值。
• 速度标度在 (当前步 − 上一段末步) 能被"间隔"整除的那一步执行，其余步纯 NVE 演化。
• 到达最后一行的步数之后，NVT 自动切换为 NVE（手册 20.2 节）。所以最后一行的步数应 ≥ MD.maxIter，否则后面的轨迹不再控温。

Tip

平衡阶段常用 \(\alpha=0\)、间隔 1–5 步，迅速把温度拉到目标；保温阶段可以用 \(\alpha\approx0.3\)–0.7、间隔 10–40 步，让体系有更自然的涨落。速度标度不产生正则分布，生产阶段建议改用 NVT_NH（§6）。

5.3 运行与检查温度

cd OpenMX/examples/07-md-si64-nvt-npt
rm -f Si64_NVT_VS.ene Si64_NVT_VS.md
mpirun -np 16 openmx Si64_NVT_VS.dat -nt 1 | tee Si64_NVT_VS.std
# step, time(fs), target T, calculated T
awk '{printf "%5d %8.1f %8.1f %8.1f\n", $1, $2, $18, $19}' Si64_NVT_VS.ene | tail
# mean temperature of the last 500 steps
tail -n 500 Si64_NVT_VS.ene | awk '{s+=$19; n++} END{print "<T> =", s/n, "K"}'

判据：保温段 \(\langle T\rangle\) 应在目标值附近（64 原子的瞬时涨落约 ±10%）；.std 中每一步都应出现 SCF 收敛（或 .md 中没有 scf_conv=0）。

5.4 NVT_VS2 与多热浴 NVT_VS4

• NVT_VS2：用每个原子自身速度计算"局域温度"，逐原子做速度标度；<MD.TempControl 格式与 NVT_VS 相同。适合局部过热区域需要快速退火的情况。
• NVT_VS4（4.0 文档化的多热浴 MD）：把原子分组，每组由各自的速度标度热浴控温。

MD.Type              NVT_VS4
MD.num.AtomGroup     2
<MD.AtomGroup
   1  1
   2  1
   3  1
   4  2
   5  2
MD.AtomGroup>

每行是"原子序号 组号"，原子数有多少行就写多少行，组号必须从 1 开始连续编号。手册指出 4.0 中所有组的温度曲线都由同一个 <MD.TempControl（速度标度格式）控制。.md 的注释行会分组打印 Temperature[1]=、Temperature[2]=。

---

6. 算例三：Nosé–Hoover 恒温（NVT_NH）

6.1 与速度标度的差别

在 Si64_NVT_VS.dat 的基础上只改动 MD 部分（完整文件 Si64_NVT_NH.dat）：

MD.Type                     NVT_NH
MD.maxIter                  2000
MD.TimeStep                 2.0
<MD.TempControl
  2
     1    1000.0
  2000    1000.0
MD.TempControl>
NH.Mass.HeatBath            30.0

6.2 <MD.TempControl（Nosé–Hoover 格式）

注意格式不同：只有两列 步数 温度(K)，温度在相邻行之间线性插值。下面是手册中的升温—降温写法：

<MD.TempControl
  4
    1   1000.0
  100   1000.0
  400    700.0
  700    600.0
MD.TempControl>

若把速度标度的 4 列格式原样用于 NVT_NH，程序按 2 列读取，行数对不上时会报 Format error for MD.TempControl，或者更糟——读成错误的温度表。改 MD.Type 时务必同时改这个块。

6.3 热浴质量 NH.Mass.HeatBath 的选择

\(Q\) 决定热浴与体系交换能量的快慢。按 Nosé 的标准关系，热浴固有周期约为

\[ \tau_\mathrm{NH}\approx 2\pi\sqrt{\frac{Q}{g\,k_\mathrm{B}T}},\qquad g\approx 3N \]

（OpenMX 内部实现的常数可能略有不同，此式只作量级估计。）以 64 个 Si、1000 K、\(Q=30\) Bohr²·amu 代入，\(\tau_\mathrm{NH}\) 约为几十 fs，与 Si 声子周期（约 65–200 fs）同一量级，是合理的起点。经验：

• \(Q\) 太小：温度高频振荡、与声子共振，扰动动力学；
• \(Q\) 太大：温度很久才达到目标，表现为长周期大幅振荡；
• 实操：用 \(Q=10, 30, 100\) 各跑 500 步，看温度曲线的振荡周期和衰减，选能在 1 ps 内平稳进入目标温度的值。体系原子数或温度变化很大时要重新估计。

检查 Nosé–Hoover 哈密顿量（.ene 第 22 列）是否守恒，是判断积分与热浴设置是否合理的好办法：

awk '{print $2, $19, $22}' Si64_NVT_NH.ene > nh_check.txt   # time, T, H_NH

6.4 推荐流程：先 VS 平衡，再 NH 生产

1. NVT_VS 升温并平衡 1–2 ps；
2. 把 Si64_NVT_VS.dat#（最后一步的坐标和速度，见 §8.3）复制成新输入，MD.Type 改为 NVT_NH、把 <MD.TempControl 改为 2 列格式、加上 NH.Mass.HeatBath；若同时把 System.Name 改为新名字（推荐，轨迹文件分开），还要删除 MD.Current.Iter 行（步数从 1 重新计）和 scf.restart on 行（_rst 目录属于旧名字）；
3. 生产运行 5–10 ps 以上，统计时丢掉开头 0.5–1 ps。

---

7. 算例四：恒温恒压 NPT_VS_PR / NPT_VS_WV（OpenMX ≥ 4.0）

7.1 输入

8 原子 Si 惯用胞，600 K，外压在前 300 步保持 0 GPa，然后线性升到 5 GPa（完整文件 Si8_NPT_VS_PR.dat）：

System.CurrrentDirectory    ./
System.Name                 Si8_NPT_VS_PR
DATA.PATH                   /opt/openmx4.0/DFT_DATA19   # EDIT
level.of.stdout             1
level.of.fileout            0

Species.Number              1
<Definition.of.Atomic.Species
  Si   Si7.0-s2p1d1   Si_PBE19
Definition.of.Atomic.Species>

Atoms.Number                8
Atoms.SpeciesAndCoordinates.Unit   Ang
<Atoms.SpeciesAndCoordinates
    1  Si    0.00005   0.01195  -0.01097   2.0 2.0
    ... (8 lines, see the example file)
Atoms.SpeciesAndCoordinates>

Atoms.UnitVectors.Unit      Ang
<Atoms.UnitVectors
  5.43   0.00   0.00
  0.00   5.43   0.00
  0.00   0.00   5.43
Atoms.UnitVectors>

scf.XcType                  GGA-PBE
scf.SpinPolarization        off
scf.ElectronicTemperature   300.0
scf.energycutoff            250.0
scf.maxIter                 100
scf.EigenvalueSolver        Band
scf.Kgrid                   4 4 4
scf.Mixing.Type             RMM-DIISK
scf.Init.Mixing.Weight      0.10
scf.Min.Mixing.Weight       0.001
scf.Max.Mixing.Weight       0.30
scf.Mixing.History          20
scf.Mixing.StartPulay       8
scf.criterion               1.0e-8

MD.Type                     NPT_VS_PR
MD.maxIter                  1000
MD.TimeStep                 1.0
NPT.LatticeRestriction      none
MD.TempTolerance            100
<MD.TempControl
  1
  1000   10   600.0   0.5
MD.TempControl>
<MD.CellPressureControl
  2
   300    0.0
  1000    5.0
MD.CellPressureControl>

7.2 关键点

关键词 / 事项	说明
应力	NPT 时程序自动打开应力张量计算（源码中对 NPT_VS_* 强制 scf_stress_flag = 1），无需另设 scf.stress.tensor
scf.energycutoff	应力对实空间网格比力更敏感，NPT 建议比同体系 NVT 再提高 20–50 Ry，并做一次"固定晶胞的应力 vs 截断能"测试
NPT.Mass.Barostat	本例不写，用默认值。PR 方法默认 \(W=\frac{3}{4\pi^2}\sum_a m_a\)（amu）；对 8 个 Si，手册示例打印的默认值约 17.0。WV 方法的默认 \(W\) 还与晶胞体积有关、量级小得多（手册 Si8 示例约 1.5×10⁻³），不要把 PR 的值用于 WV
<MD.TempControl	与 NVT_VS 相同的 4 列格式
<MD.CellPressureControl	步数 压力(GPa)；第一行之前取第一行的压力，段间线性插值；正值压缩晶胞
MD.TempTolerance	100 K：温差小于 100 K 时不做速度标度，避免 NPT 中频繁扰动
NPT.LatticeRestriction	cubic 保持立方，ortho 保持正交；研究可能发生剪切或相变时用 none
不兼容	非共线 scf.SpinPolarization nc（报 Non-collinear calculation is not compatible with stress calculation.）、基函数优化、DFT+U 的 full 占据投影（'full' local propjector in DFT+U is not compatible with stress calculation.），以及 <MD.Fixed.XYZ 约束 MD
DFT-D	手册没有说明 NPT 与 scf.dftD 的兼容性（晶胞优化已明确不支持 DFT-D），使用前请先做小测试

单轴加压（正交晶胞，晶格矢量只有对角元）：

MD.applied.pressure.flag    1 0 0     # pressure only along a

7.3 System.Name.npt.ene

NPT 额外输出 System.Name.npt.ene（手册 20.6 节）：

列	含义	单位
1	MD 步	—
2	MD 时间	fs
3	体积	Bohr³
4	计算温度	K
5	内压 \(=-\tfrac13\mathrm{Tr}(\sigma)/V\)	GPa
6	\(U_\mathrm{tot}\)	Hartree
7	\(U_\mathrm{kc}\)	Hartree
8	\(U_\mathrm{baro}\)（晶胞动能）	Hartree
9	\(W_\mathrm{baro}\)（外压 × 体积）	Hartree
10	\(U_\mathrm{tot}+U_\mathrm{kc}+U_\mathrm{baro}+W_\mathrm{baro}\)	Hartree

由体积求平均晶格常数（立方胞，1 Bohr = 0.529177 Å）：

# average over steps 200-300 (0 GPa segment, after 200 fs equilibration)
awk '$1>200 && $1<=300 {v+=$3; n++} END{a=(v/n)^(1/3)*0.529177; printf "<V>=%.2f Bohr^3, a=%.4f A\n", v/n, a}' Si8_NPT_VS_PR.npt.ene

（示意输出，数值仅供参考）<V>=1xxx.xx Bohr^3, a=5.4xxx A。PBE 的 Si 晶格常数约 5.46–5.47 Å（比实验 5.431 Å 大约 0.5–0.7%），600 K 下还应再加上热膨胀。8 原子体系的体积涨落很大，要得到可靠的密度，需要更大的超胞和更长的轨迹。也可以与第6章的 0 K 状态方程结果对照。

Note

在本文核对的 4.0 系源码（iterout.c）中，NPT 的 .md 轨迹速度列写的是 0，需要速度做分析（如速度自相关）时请在 NVT/NVE 下进行。晶胞随时间的变化记录在 .md 每帧注释行的 Cell_Vectors= 中。

本章引用的源码细节（Input_std.c、iterout.c、iterout_md.c、Make_InputFile_with_FinalCoord.c、MD_pac.c）核对自 OpenMX 4.0 源码及基于 4.0 的 GPU 移植版 dc1394/openmx4.0_gpu；补丁版本之间可能有细微差异，以你实际使用的源码为准。

---

8. 输出文件、可视化与续算

8.1 文件一览

文件	内容	用途
System.Name.ene	每步一行的能量、温度等（§4.5）	温度、守恒量、费米能随时间变化
System.Name.md	每步一帧的扩展 xyz 轨迹（追加写入）	动画、RDF/MSD、dpdata
System.Name.md2	最后一步的坐标，带用户定义的物种名	直接粘贴为新输入的坐标块
System.Name.dat#	用最后一步坐标（和速度等）重写的完整输入文件	续算
System.Name.npt.ene	NPT 的体积/压力/能量（§7.3）	密度、状态方程
System.Name.out	最后一步的 SCF、能量分解、Mulliken 电荷等	检查电子结构
System.Name.xyz	最终结构	快速查看
System.Name_rst/	重启用的哈密顿量与电荷密度	scf.restart on

8.2 4.0 的 .md 格式

OpenMX 4.0 手册（第 18.6 节）首次详细记录了 .md 格式。每一帧：

Natom
time= 0.000 (fs) Energy= -xx.xxxxx (Hartree) Temperature= ... Cell_Vectors= ax ay az bx by bz cx cy cz  scf_conv=1
Elem  x  y  z  Fx  Fy  Fz  vx  vy  vz  Q  m  theta  phi
...

• 坐标 Å，力 Hartree/Bohr，速度 m/s，Q 为 Mulliken 净电荷（电子电荷记为负），m 为局域磁矩（μB），theta、phi 为自旋方向欧拉角（度）；
• 注释行中 Cell_Vectors= 后 9 个数是晶格矢量（Å），NPT 中逐帧变化；scf_conv=1/0 表示该步 SCF 是否收敛；
• 第 1 列是元素符号（不是你在输入中定义的物种名）。用 3.9 生成的 .md 请先 head 核对列数再分析。

8.3 用 .dat# 续算

每一步 OpenMX 都会写出 System.Name.dat#。对 MD，它包含（据 4.0 源码 Make_InputFile_with_FinalCoord.c）：

• 最后一步的坐标（NPT 时还有新的晶格矢量）；
• <MD.Init.Velocity（当前速度）和 <MD.Init.Velocity.Prev；
• Nosé–Hoover 热浴状态 NH.R、NH.nzeta、NH.czeta；
• MD.Current.Iter（已完成的步数）、scf.fixed.grid（网格原点）以及 scf.restart on。

续算步骤：

cp Si64_NVT_NH.dat# Si64_NVT_NH_cont.dat
# MD.maxIter is the ABSOLUTE last step: e.g. 2000 done -> run to 4000
sed -i 's/^MD.maxIter .*/MD.maxIter   4000/' Si64_NVT_NH_cont.dat
# the <MD.TempControl table must also cover steps up to 4000
mpirun -np 16 openmx Si64_NVT_NH_cont.dat -nt 1 | tee Si64_NVT_NH_cont.std

Important

• 程序在"MD.Current.Iter + 本次步数 > MD.maxIter"时停止。忘记增大 MD.maxIter，续算会一步不走就结束。
• .dat# 中保留了原来的 System.Name，新的帧会追加到同一个 .ene/.md 末尾，时间列继续累加，正好得到连续轨迹；scf.restart on 读取 System.Name_rst/ 加速第一步 SCF，所以不要删掉该目录。
• 速度标度的 <MD.TempControl 最后一行步数也要 ≥ 新的 MD.maxIter，否则超出部分按 NVE 演化。

8.4 MD 中的 SCF 外推与重启

MD 的每一步会自动利用前几步的电荷密度外推得到初猜（scf.ExtCharge.History，手册推荐 2），通常第 2 步以后每步只需 5–15 次 SCF 迭代。若同一 System.Name 已有收敛的 SCF（例如先做过单点），可在 MD 输入中加 scf.restart on，第 1 步也能从该密度出发。若某一步 SCF 未收敛，.md 中该帧的 scf_conv=0：偶发一两次通常无害；频繁出现说明 scf.maxIter 太小或混合参数不合适（见第5章），这种帧不应进入机器学习训练集。

8.5 可视化

• OpenMX Viewer（浏览器，https://www.openmx-square.org/viewer/ ）：直接拖入 System.Name.md，可播放动画、按净电荷/磁矩/力/速度着色，并能正确显示周期性晶胞（读取 Cell_Vectors=）。
• XCrySDen：能读 xyz 格式的 .md，但不能画周期扩展；需要时先用 md2axsf 转换：

cd openmx4.0/source && gcc md2axsf.c -lm -o md2axsf && cp md2axsf ~/bin/
md2axsf Si64_NVT_NH.md Si64_NVT_NH.axsf
xcrysden --axsf Si64_NVT_NH.axsf

• ASE / OVITO：仓库脚本 ../scripts/omx_traj.py 可把 .md 转为 extxyz / ASE Trajectory，并绘制能量、温度随步数的变化，用 python ../scripts/omx_traj.py -h 查看用法。

---

9. 轨迹分析：温度、能量、RDF 与 MSD

示例目录中的 md_analysis.py（仅依赖 numpy，matplotlib 可选）实现了以下分析，读者可以直接使用或改写：

python3 md_analysis.py Si64_NVT_NH --skip 250 --rmax 5.0 --plot

9.1 读取 .ene

import numpy as np
ene = np.loadtxt("Si64_NVT_NH.ene")
time_fs, T, Econs = ene[:, 1], ene[:, 18], ene[:, 15]   # columns 2, 19, 16 (0-based 1, 18, 15)
print("<T> =", T[250:].mean(), "+-", T[250:].std(), "K")

9.2 解析 .md

def read_md(fname):
    frames = []
    with open(fname) as f:
        while True:
            line = f.readline()
            if not line.strip():
                break
            nat = int(line.split()[0])
            head = f.readline()
            cell = None
            if "Cell_Vectors=" in head:
                v = head.split("Cell_Vectors=")[1].split()[:9]
                cell = np.array(v, float).reshape(3, 3)
            sym, pos = [], []
            for _ in range(nat):
                p = f.readline().split()
                sym.append(p[0]); pos.append([float(x) for x in p[1:4]])
            frames.append((sym, np.array(pos), cell, head))
    return frames

9.3 径向分布函数 g(r) 与配位数

对一般三斜晶胞，用分数坐标取最小像（要求 \(r_\mathrm{max}\) 小于晶胞最短"面间距"的一半）：

\[ g(r)=\frac{V}{N^2}\,\frac{\langle n(r,r+\Delta r)\rangle}{4\pi r^2\Delta r},\qquad n_\mathrm{c}(r_\mathrm{min})=4\pi\rho\int_0^{r_\mathrm{min}} g(r)\,r^2\,dr \]

晶体 Si 的第一峰在约 2.35 Å，配位数 4；1000 K 下峰明显展宽但第二峰（约 3.84 Å）仍清晰。若看到第一峰之后几乎平坦、配位数明显大于 4，说明体系已熔化（PBE 下 Si 熔点的计算值依赖于方法，这里不给具体数值）。

9.4 均方位移与扩散

\[ \mathrm{MSD}(t)=\frac1N\sum_I\left\langle|\mathbf{R}_I(t_0+t)-\mathbf{R}_I(t_0)|^2\right\rangle_{t_0},\qquad D=\lim_{t\to\infty}\frac{\mathrm{MSD}(t)}{6t} \]

md_analysis.py 先用相邻帧的最小像位移"展开"坐标（不依赖 OpenMX 是否把原子折回晶胞），再做多时间原点平均。固体中 MSD 很快饱和（等于热振动幅度的平方），液体中随时间线性增长。

（示意输出，数值仅供参考）

ene: 1750 steps used (skip 250)
<T> = 99x.x +- 9x.x K (target 1000.0 K at last step)
drift of Utot+Ukc     =   xxx.xxx meV/ps (meaningful for NVE / NVT_NH respectively)
drift of H_NoseHoover =     x.xxx meV/ps (meaningful for NVE / NVT_NH respectively)
md: 1750 frames used, dt = 2.000 fs, 64 atoms
g(r) Si-Si: first peak 2.3xx A, first minimum 2.9xx A, coordination 4.0x
MSD at 1750 fs = 0.0xxx A^2; D from linear fit (second half) = x.xxxe-0x cm^2/s

脚本同时写出 前缀_rdf.dat（r, g(r), n(r)）与 前缀_msd.dat，加 --plot 生成三联图（温度、g(r)、MSD）。对 NPT 运行，它还会读取 .npt.ene 给出平均体积、温度、内压和等效立方晶格常数。

---

10. 从 AIMD 到机器学习势：dpdata 与 DeePMD / NEP

OpenMX 4.0 手册（第 18.6 节）说明 .md 格式与 dpdata 的 openmx/md 格式一致。dpdata 读取时需要同一前缀的 <前缀>.dat 与 <前缀>.md，并从 .md 每帧读取能量（Hartree → eV）、坐标与力（Hartree/Bohr → eV/Å）、晶胞：

import dpdata

# Si64_NVT_NH.dat (input) and Si64_NVT_NH.md (trajectory) must be in the cwd
ls = dpdata.LabeledSystem("Si64_NVT_NH", fmt="openmx/md")
print(ls)                         # number of frames, atoms, element types
ls = ls[250::5]                   # drop equilibration, take every 5th frame
ls.to("deepmd/npy", "deepmd_si64", set_size=200)
ls.to("extxyz", "si64_train.xyz")  # e.g. for NEP / GPUMD or MACE workflows

注意事项：

• openmx/md 读取器和 extxyz 输出别名都需要较新的 dpdata（本文写作时最新为 1.1.0）；旧版本没有 extxyz 别名时可改用 quip/gap/xyz；
• 输入文件的物种名要与元素符号一致（dpdata 用物种名在 .md 行中匹配元素符号）；若用了 Si1、Fe_up 这类物种名，转换会失败或漏原子；
• dpdata 遇到 scf_conv=0 的帧会给出 SCF not converged! 警告，训练前应剔除这些帧；
• 相邻 MD 帧高度相关，训练集应稀疏采样（每隔 5–20 帧）并覆盖多个温度/压力；
• 机器学习势的训练流程（NEP/GPUMD 数据集脚本）见仓库 ../../Scripts/README.md，更多工作流见 第6章 工作流自动化 与 第7章 学习路线。

---

参数选择与收敛建议

参数	推荐起点	收敛 / 检查方法
MD.TimeStep	含 H：0.25–0.5 fs；Si：1–2 fs	NVE 下守恒量无漂移
scf.criterion	NVE：1e-8–1e-9；NVT/NPT：1e-7	收紧一个数量级后 RDF、温度、压力不变
scf.energycutoff	200–250 Ry；NPT 再加 20–50 Ry	固定构型下力/应力随截断能收敛；NVE 中 egg-box 引起的能量起伏可忽略
基组	先 Quick 试跑，正式用 Standard	比较若干帧的力（与 Standard/Precise 的力差 < 0.05 eV/Å 量级）
k 点	超胞 > 10 Å 时 Γ 或 2×2×2	抽几帧做更密 k 点单点，比较能量与力
平衡时间	0.5–2 ps	温度/能量/体积的滑动平均不再有趋势
生产时间	≥ 5–10 ps（RDF）、更长（扩散）	分块平均的误差棒足够小
NH.Mass.HeatBath	20–50（几十原子）	温度振荡周期合理、第 22 列守恒
NPT.Mass.Barostat	程序默认值	体积振荡周期约 0.1–1 ps；太快则增大，太慢则减小
超胞大小	≥ 50–100 原子（液体/扩散）	有限尺寸效应：比较两个尺寸的 g(r)、D

---

常见问题与报错

现象 / 报错	原因	解决
温度一直为 0、原子不动，或 .ene 中出现 nan	从理想晶格出发且没有初速度：受力为零、动能为零，速度标度除以零	给坐标加 0.02–0.05 Å 的随机位移，或用 <MD.Init.Velocity（gen_init_velocity.py）
Format error for MD.TempControl	速度标度（4 列）与 Nosé–Hoover（2 列）格式混用，或首行行数与实际行数不符	按 MD.Type 改写块，核对行数
跑到某一步后温度不再受控、开始漂移	速度标度在 <MD.TempControl 最后一行步数之后自动切换为 NVE	最后一行步数 ≥ MD.maxIter
续算立即结束，没有新步	.dat# 里有 MD.Current.Iter，而 MD.maxIter 未增大	把 MD.maxIter 改成新的绝对终止步
.ene / .md 中出现重复步号或时间倒流	同一 System.Name 从头重跑，文件以追加模式写入	重跑前删除旧的 .ene/.md，或换 System.Name
NVE 守恒量明显漂移	时间步太大、SCF 不够收敛、网格太粗（egg-box）、偶发 SCF 不收敛	依次减小步长、收紧 scf.criterion、提高 scf.energycutoff；检查 scf_conv=0
NVE 中温度"自己下降一半"	从势能极小附近 + 初速度出发，动能向势能均分	正常现象；需要固定温度请用 NVT
Please check your input for atom group!!	NVT_VS4 的组号超出 MD.num.AtomGroup 或不连续	组号从 1 连续编号，行数等于原子数
Non-collinear calculation is not compatible with stress calculation.	NPT 需要应力，而非共线计算不支持	NPT 只用共线（off/on）计算
'full' local propjector in DFT+U is not compatible with stress calculation.（源码原文拼写）	NPT + DFT+U 用了 scf.Hubbard.Occupation full	改用默认的 dual 占据
NPT 中晶胞剧烈振荡或塌缩	活塞质量太小、截断能不足导致应力误差大、外压设置过大	提高 scf.energycutoff，增大 NPT.Mass.Barostat，先在 0 GPa 下平衡
大体系 MD 中途段错误（segmentation fault）	栈空间不足或内存不够	作业脚本中 ulimit -s unlimited、export OMP_STACKSIZE=512M，减少每节点进程数；见第18章
dpdata 读取后原子数为 0 或报维度错误	物种名不是元素符号；.dat 与 .md 前缀不一致	物种名用元素符号；保证 前缀.dat 与 前缀.md 同名

---

练习

1. 用 H2O_NVE.dat 分别以 0.25、0.5、1.0 fs 跑 500 fs，画出三条守恒量曲线，找出能量开始明显漂移的步长；再把 H 质量改为 2.0 重复 1.0 fs 的计算。
2. 把 Si64_NVT_VS.dat 中 <MD.TempControl 的 \(\alpha\) 从 0.0 改为 0.7，比较保温段的温度涨落。
3. 用 NH.Mass.HeatBath = 5、30、300 各跑 1 ps 的 Si64_NVT_NH，比较温度曲线与第 22 列的守恒情况。
4. 把 Si8_NPT_VS_PR.dat 改为 NPT_VS_WV（删除 NPT.Mass.Barostat 让程序给默认值），比较两种恒压方法的体积曲线。
5. 把 Si64 加热到 2000 K 以上并保温，观察 g(r) 与 MSD 的变化，判断是否熔化；用 dpdata 把 1000 K 与高温两段轨迹合并导出为一个训练集。

---

参考资料

官方文档

• OpenMX 4.0 手册 第 20 章 Molecular dynamics：20.1 NVE、20.2 NVT_VS、20.3 NVT_NH、20.4 NVT_VS4、20.6 NPT_VS_PR、20.7 NPT_VS_WV、20.8 初速度、20.9 原子质量、20.10 md2axsf
• 手册 17.2 MD 与优化中的电荷外推、17.3 续算输入 .dat#
• 手册 18.6 md 文件格式（与 dpdata 的关系）
• 手册 8.2 关键词（MD.Type、MD.TempControl、NH.Mass.HeatBath 等）
• OpenMX What's new（4.0 新增 NPT MD、md 格式与 DeePMD-kit 联动）
• OpenMX Viewer
• dpdata 文档

相关 GitHub 仓库（仅链接，未复制其代码）

• deepmodeling/dpdata — 读取 OpenMX .dat + .md（openmx/md 格式），转换为 DeePMD 等训练数据
• deepmodeling/deepmd-kit — 用 AIMD 数据训练深度势能
• issp-center-dev/abICS — 构型采样工具，可生成 OpenMX 输入作为训练数据
• Tomoki-YAMASHITA/CrySPY — 晶体结构预测，可调用 OpenMX 做弛豫
• cmsi/MateriAppsLive — 预装 OpenMX 的 Linux 环境，适合先体验 MD
• 本仓库相关：../scripts/omx_traj.py（轨迹转换与作图）、../scripts/job_omx_slurm.sh（作业模板）、HPC 作业调度、后处理与可视化、相关仓库总表

关键文献

• L. Verlet, Phys. Rev. 159, 98 (1967). DOI: 10.1103/PhysRev.159.98；W. C. Swope, H. C. Andersen, P. H. Berens, K. R. Wilson, J. Chem. Phys. 76, 637 (1982)（速度 Verlet）. DOI: 10.1063/1.442716
• L. V. Woodcock, Chem. Phys. Lett. 10, 257 (1971)（速度标度）. DOI: 10.1016/0009-2614(71)80281-6
• S. Nosé, J. Chem. Phys. 81, 511 (1984). DOI: 10.1063/1.447334；W. G. Hoover, Phys. Rev. A 31, 1695 (1985). DOI: 10.1103/PhysRevA.31.1695
• M. Parrinello, A. Rahman, J. Appl. Phys. 52, 7182 (1981). DOI: 10.1063/1.328693
• R. M. Wentzcovitch, Phys. Rev. B 44, 2358 (1991)（变胞 MD）. DOI: 10.1103/PhysRevB.44.2358
• T. Ozaki, Phys. Rev. B 67, 155108 (2003)（OpenMX 的变分优化 PAO 基组）. DOI: 10.1103/PhysRevB.67.155108
• L. Zhang, J. Han, H. Wang, R. Car, W. E, Phys. Rev. Lett. 120, 143001 (2018)（DeePMD）. DOI: 10.1103/PhysRevLett.120.143001
• D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press（RDF、MSD、分块平均等分析方法）

---

← 上一章 | 返回目录 | 下一章 →