dpdata读取分子量化与AMBER数据

分子量化常输出原子单位的能量和梯度，分子动力学则可能将坐标、力、能量分开保存。不能把梯度直接当作力，也不能把打印的非周期占位晶胞当成真实盒子。本篇用原始格式样本演示完整读取；样本来自固定提交的dpdata官方示例，在集群实际进行解析与转换。

1 下载和运行#

下载格式实例包，保存为dpdata-formats.zip。

python3 -m zipfile -e dpdata-formats.zip .
cd dpdata-formats
bash run.sh

最后一条列出示例名，选中一个再运行，例如bash run.sh qe-cp。脚本调用集群现有Python 3.11环境中的dpdata 1.1.0、NumPy等依赖，只读取、转换数据，不启动原计算程序。需要批作业时可用sbatch job.slurm qe-cp，只申请1个CPU。输出保存在results/。

2 Gaussian：日志、fchk与MD#

bash run.sh gaussian-log
bash run.sh gaussian-fchk
bash run.sh gaussian-md

import dpdata
ch4 = dpdata.LabeledSystem('samples/gaussian/methane.gaussianlog', fmt='gaussian/log')
water = dpdata.LabeledSystem('samples/gaussian/waterfreq.gaussianfchk', fmt='gaussian/fchk')
md = dpdata.LabeledSystem('samples/gaussian/aimd_gaussian_CH4_output', fmt='gaussian/md')
print(len(ch4), len(water), len(md))
print(ch4['energies'][0])

分别为1、1、22帧。CH₄样本读取能量约−1102.714589869 eV，5个原子；这是给定方法与结构的总能量。gaussian/log默认单点读取与gaussian/md多步读取用途不同，不能以构造对象成功推断已经读完所有优化/MD步。

图中是电子能量相对首帧的变化，不包含离子动能，不用来判断NVE总能守恒或温度是否平衡。

3 ORCA与Psi4：梯度已经转换为力#

bash run.sh orca-molecule
bash run.sh psi4-molecule

orca = dpdata.LabeledSystem('samples/orca/orca.spout', fmt='orca/spout')
psi = dpdata.LabeledSystem('samples/psi4/psi4.out', fmt='psi4/out')
print(orca['forces'].shape, psi['forces'].shape)

两个例子各1帧。Psi4原输出总能为−79.8685493165660603 Hartree，dpdata读取后转换为eV；力为负梯度，单位eV/Å。不要再对dpdata结果额外乘一次Hartree→eV，也不要再反一次符号。

4 DFTB+：输入和输出成对#

bash run.sh dftb-NH3

nh3 = dpdata.LabeledSystem(
    ('samples/dftbplus/dftb_pin.hsd', 'samples/dftbplus/detailed.out'),
    fmt='dftbplus',
)
print(nh3['energies'][0], nh3['forces'][0])

本例NH₃含4个原子，原生detailed.out的总能为−3.2963983884 Hartree，约−89.6996 eV。读取使用较多有效数字的原子单位值；日志右侧只打印四位小数的eV值存在舍入差，不人为改成完全相等。

箭头共用一个缩放系数，仅用于显示方向和相对大小；CSV保留真实eV/Å数值。原始样本为非平衡几何，较长的N—H距离按原数据保留，不用作优化平衡键长；该图不是新完成的集群DFTB计算。

5 SQM：有坐标不代表有力#

bash run.sh sqm-molecule

sqm = dpdata.LabeledSystem('samples/amber/sqm_forces.out', fmt='sqm/out')
print(sqm.get_natoms(), sqm['forces'].shape)

本例40原子，源文件确实含梯度，才可以得到力标签。只包含几何的SQM输出应按System读取；不要给不存在的力补零。

6 AMBER：四类文件缺一不可#

bash run.sh amber-md

amber = dpdata.LabeledSystem('samples/amber/02_Heat', fmt='amber/md')
print(len(amber), amber.get_natoms())

例子按共同前缀读取拓扑.parm7、坐标.nc、力.mdfrc、能量.mden，得到10帧、1912个原子。默认类型表保留AMBER力场类型（如CT、HW），它们不是化学元素；准备元素型机器学习数据时应显式转换成元素符号，见下面的命令。这里的能量是势能，不含离子动能；dpdata将相应AMBER单位转为统一eV与eV/Å。把不同任务的同名文件放在一起，即使帧数相同，也不代表数据配对正确。

# 这里的1912来自本例拓扑，其他体系先用get_natoms()读取原子数。
elements = dpdata.LabeledSystem(
    'samples/amber/02_Heat', fmt='amber/md',
    use_element_symbols=list(range(amber.get_natoms())),
)
print(elements['atom_names'])

转换类型后仍应核对原子顺序和对应坐标/力不变，不按元素重新排列单个数组。

7 非周期数据的晶胞占位#

print(ch4.data.get('nopbc', False))
ch4.to('deepmd/npy', 'ch4-data')
back = dpdata.LabeledSystem('ch4-data', fmt='deepmd/npy')
print(back.data.get('nopbc', False))

本例保持nopbc=True，坐标、能量和力在双精度导出读回后保持；读回端可能使用100 Å的占位晶胞，原始分子日志则为零晶胞。这个变化不表示分子体积膨胀，也不应计算密度或用晶胞占位值换算压力。周期性、物理盒子和标签要分别核对。

继续看对象、混合体系与键图。样本出处见官方固定提交，所用日志与输入可在包内直接检查。

图表数据与重新绘图#

bash plot.sh

图件保存在plots/，原数值另存于附件reference-data/。直接下载：AMBER-energy.csv、DFTB-NH3-geometry-forces.csv、Gaussian-CH4-trajectory.csv、Si9-reference.csv。