GPUMDkit 按能量、力和virial筛选构型

本篇接着标签范围检查演示怎样保留符合条件的构型。阈值只是九组Si标签的操作示例，不是通用的数据质量标准；按自己的体系、能量参考和目标训练范围选择阈值，保留原文件及被排除构型的来源。

1. 先得到可回查的结果#

下载输入、命令和参考结果，解压进入附件根目录：

bash run.sh

程序会将各条件放在不同的work/子目录，避免每次生成的filtered.xyz互相覆盖。输入共九帧、每帧两个Si原子，来源见PROVENANCE.json；reference/native/保存原始终端输出与真实筛选文件。这里只有标签处理，不生成新的电子结构计算。

2. 原生命令究竟比较什么#

等于阈值时保留。任意所需标签缺失、非有限或重复定义互相矛盾都会停止。不要因为命令成功生成文件就直接开始训练，先检查保留数量、结构、元素和标签。

3. 以力为例，自己敲一遍#

先在附件根目录新建独立工作副本，再运行短命令：

module load gpumdkit/1.5.8
mkdir inspect-filter
cp inputs/canonical.xyz inspect-filter/train.xyz
cd inspect-filter
gpumdkit.sh -filter_value train.xyz force 1.991418

当前目录生成filtered.xyz，原始train.xyz保留。实际保留5帧，原始ID依次为0, 1, 2, 5, 6，仍按输入顺序排列；这些帧的最大绝对力分量不超过1.991418 eV/Å。下面直接读出结果，能量和力通过ASE的正式接口获取：

python - <<'ASE'
from ase.io import read
frames = read('filtered.xyz', index=':')
print('frames:', len(frames))
for a in frames:
    print(a.info['source_frame_id'], a.get_potential_energy(),
          abs(a.get_forces()).max())
ASE

4. 能量与virial的实际结果#

附件分别在独立目录运行以下条件；需要手动尝试时也请为每个条件复制一份输入，避免固定名称的输出覆盖上一例：

能量与virial的原生命令如下：

gpumdkit.sh -filter_value train.xyz energy -10.68471052
gpumdkit.sh -filter_value train.xyz virial 0.7583178060226703

上面两行在同一目录连续执行时，最终filtered.xyz是第二条命令的结果；附件自动运行时各放独立目录。它们不是串联筛选，如需串联，先给前一步输出另取文件名，再作为下一步输入。

第7帧的virial最大有符号分量约0.758318 eV，而最大绝对分量约2.357193 eV，仍能通过本例有符号上限。这说明原生命令不能代替“所有分量绝对值都小于某个阈值”的筛选。压力还需要体积、符号与单位换算，不能把这里的eV直接写作GPa。

5. 为什么力模长与分量上限不同#

附件inputs/rotated.xyz将同一批结构的坐标、晶胞、力和virial共同绕z轴旋转45°，没有再做DFT。能量和每个原子的力模长保持不变，而Cartesian分量会改变，因此原生分量筛选具有坐标基底依赖。

module load gpumdkit/1.5.8
python -B rotation_check.py

该脚本核对完整坐标/晶胞、力旋转、virial的R W Rᵀ以及模长，结果写入work/ROTATION.json。例如ID2的力模长约1.991418 eV/Å；旋转后最大绝对分量约1.408145 eV/Å。附件的旋转筛选使用1.40814517 eV/Å，实际保留ID 0, 1, 2, 5, 6。这不是物理力变小，而是坐标分量改变。

如果自己的标准是“所有原子力模长都小于阈值”，请用ASE明确计算模长并筛选，不要把-filter_value ... force误认为这一标准。下面把自己的模长阈值填入第一行；示例仍以1.991418 eV/Å演示，输出单独保存：

python - <<'ASE'
import numpy as np
from ase.io import read, write
threshold = 1.991418
frames = read('train.xyz', index=':')
kept = [a for a in frames
        if np.all(np.linalg.norm(a.get_forces(), axis=1) <= threshold)]
if not kept:
    raise SystemExit('No frames selected; review the threshold.')
write('filtered_by_norm.xyz', kept, format='extxyz')
print('kept:', len(kept))
ASE

该模长示例实际保留4帧，ID为0, 1, 5, 6。ID2除了边界上的x分量，还有很小的横向力，所以其模长严格大于1.991418；用四舍五入的终端数字判断边界会忽略这一点。实际筛选时明确使用的量和精度，不要随意把边界附近的差异归为坏数据。

这个ASE补充操作针对当前规范的forces标签文件；旧字段force需先规范，不能指望ASE自动当作calculator力结果。

6. 输出检查与下一步#

filtered.xyz保留入选帧的原始顺序、坐标、晶胞和标签；旧字段输入仍保留旧字段，兼容读取不会替你完成NEP训练格式规范化。空选择会生成0字节文件并提示保留0帧，这不是有效训练集。原始文件若恰好就是当前目录的filtered.xyz，请先更换工作副本名称，程序会拒绝输入输出同一路径。

附件的work/labels.csv汇总总能量、最大力模长、最大绝对力分量和两种virial范围；work/ANALYSIS.json记录各条件实际选出的ID，便于回查。先确认为什么某帧被选中，再进行抽样与训练数据诊断。