GPUMDkit 约束优化与补充数据工具
实跑固定原子优化、逐帧拆分、一般晶胞转LAMMPS、原子编号、OVITO RDF、体积提取及能量力应力预测图。
阅读本篇目录
GPUMDkit除了菜单入口，还有直接调用的Python脚本。本文把常用补充流程逐一跑通，并说明输入和结果该怎样检查。下载案例资料包，解压后进入根目录：
module load gpumdkit/1.5.8
python run.py --list
每个run.py案例复制输入到独立work/名称，再运行本章列出的原生命令；原始参考数据保存在reference。计算放在申请的资源内，已有工作目录不会覆盖。
1. 固定部分原子，再优化其余原子#
本例是64原子PbTe，晶格参数6.56 Å，加入标准差0.015 Å、种子53的微小位移。固定原子位置、保持晶胞，优化自由原子至最大力0.02 eV/Å。
python run.py fixed-1
python run.py fixed-2
python run.py fixed-3
原生调用形式：
python "$GPUMDkit_path/Scripts/calculators/calc_minimize_fixed.py" \
model.xyz nep.txt 1 1 0.02 200
最后四项依次是分组方法、要固定的组编号、最大自由原子力标准、最多步数。这里的组编号按这个脚本的1起始选择顺序解释，不能直接等同于GPUMD中通常从0开始的group值。
输入model.xyz除了元素和坐标还带type、group整数列。本例显式设Te的type为0、Pb为1：
| 案例 | 参数 | 固定对象 | 固定/自由原子数 |
|---|---|---|---|
| fixed-1 | 1 1 0.02 200 |
按group选择全部Pb | 32/32 |
| fixed-2 | 2 2 0.02 200 |
按type选择全部Pb | 32/32 |
| fixed-3 | 3 2 0.02 200 |
按group与type组合选择下半区域Pb | 15/49 |
第三种输入先扰动再按实际z坐标分组，部分边界原子跨过分组面，所以固定15个，并非几何理想化的16个。固定集合要以实际输入为准。第一、第二种选择的是同一组原子，最终结构一致，不是两个不同物理条件。
输出minimized_fixed.xyz包含优化过程。实际结果：
| 条件 | 初始能量/eV | 最终能量/eV | 最终最大自由原子力/eV·Å⁻¹ |
|---|---|---|---|
| 固定全部Pb | −245.051070 | −245.084414 | 0.014824 |
| 固定下半区域Pb | −245.051070 | −245.102969 | 0.015835 |
不能把固定原子上的受力当成已优化为零，也不能把文件最后一帧自动当成收敛结构；需检查自由原子的最大力及是否达到步数上限。
2. 多帧extxyz拆成单帧文件#
python run.py split-frames
对应命令：
python "$GPUMDkit_path/Scripts/format_conversion/split_single_xyz.py" movie.xyz
输入是LLZO实际1000 K轨迹的5帧，输出model_1.xyz到model_5.xyz，各192原子。逐帧的元素、坐标和晶胞保持对应。适用于后续逐帧计算或检查；拆分不会让短轨迹自动获得更好的统计。
3. 转成LAMMPS结构文件#
python run.py cell-orthogonal
python run.py cell-general
原生转换：
python "$GPUMDkit_path/Scripts/format_conversion/exyz2lmp.py" model.xyz model.data
输出为atom_style atomic几何文件，保留元素质量、类型和坐标，不提供电荷、键、速度或力场。运行MD前还需按模型添加相应LAMMPS设置。
一般晶胞不能只把三个基矢分量抄成xhi/yhi/zhi和倾斜量；转换必须同时旋转基矢和笛卡尔坐标。当前脚本使用ASE Prism处理。本例第二个控制把实际LLZO帧仿射变成斜晶胞再绕z旋转31°，是几何转换测试，不是新的MD平衡结构。
旋转前的受限三斜基矢 / Å：
13.3 0.0 0.0
1.1 13.4 0.0
0.6 0.4 13.5
体积：2405.97 ų
转换后体积2405.97 ų、元素和质量保持，周期边界下的原子间距最大变化低于1E-8 Å。包内还实际用LAMMPS读写了一遍：
units metal
atom_style atomic
read_data model.data
write_data roundtrip.data
这里只读写几何，没有运行势能或动力学计算。
4. 用OVITO计算分元素RDF#
python run.py rdf-ovito
对应命令：
module load gpumdkit/1.5.8
module load ovito/3.16.1
python "$GPUMDkit_path/Scripts/calculators/rdf_calculator_ovito.py" \
movie.xyz 5.0 100
5.0是截断距离/Å，100是距离区间数。本例平均LLZO的5个实际轨迹帧，输出rdf.txt有100行，第一列r，后10列依次为Li–Li、Li–La、Li–Zr、Li–O、La–La、La–Zr、La–O、Zr–Zr、Zr–O、O–O；读取自己的文件时以文件头为准，不同元素或类型映射可能改变列顺序。
左图展示其中三条分量，不是总RDF；右图展示同一晶胞旋转前后的MIC原子间距。RDF峰反映短程结构，5帧仅演示操作，不能代表采样已充分。
5. LAMMPS dump重编号#
python run.py renumber
对应命令：
gpumdkit.sh -re_atoms input.dump renumbered.dump
该脚本要求原子字段第一列是id，按当前行顺序重新编号为1…N。本例Ni实际dump有2048原子，原编号本来连续，因此坐标、速度和编号均保持。对有缺口的编号可消除缺口，但不能恢复跨帧被打乱的原子身份；做扩散、位移或轨迹跟踪时，不要逐帧重编号后假定同一ID仍代表同一个原子。
6. 从温度目录提取平均体积#
python run.py volume
原生命令在含950K、1000K等子目录的位置运行：
gpumdkit.sh -get_volume
各目录放入实际GPUMD thermo.out。本例是18列格式，末9列组成逐帧晶胞矩阵，体积为行列式。输出的生产段平均体积：
| 温度/K | 平均体积/ų |
|---|---|
| 950 | 2237.832 |
| 1000 | 2246.787 |
| 1100 | 2261.301 |
| 1200 | 2277.586 |
这些数值来自LLZO短MD同一seed91序列。应先确定平衡段截取方式和热膨胀采样是否收敛，再用于物性拟合。
7. 能量、力、应力三联预测图#
python run.py prediction2-si
对应命令：
python "$GPUMDkit_path/Scripts/plt_scripts/plt_prediction2.py" save
输入为energy_train.out、force_train.out、stress_train.out，输出prediction.png。该补充脚本要求至少存在有效应力数据，不能把全部是缺失标记的应力列直接用于此图。
本例Si数据的原始RMSE：能量0.790116 eV/atom、力0.051345 eV/Å、应力0.326525 GPa。能量有明显基准偏置，不能因为力散点看起来接近对角线就宣称模型整体准确；比较能量时需核对参考零点、训练来源和标签定义。本图是训练数据回算演示，不是独立测试集精度。
更多入口见GPUMDkit功能索引、42类绘图导航、抽样与数据集对比和空位迁移NEB。