Calorine读取GPUMD数据与qNEP响应

下载模型、真实数据、完整脚本与参考图片。相关教程：偶极与极化率 · 红外与拉曼谱 · NEP89微调 · 模型编辑 · GPUMD数据读写 · 完整接口索引

1 准备与运行#

先解压并重画已经完成的参考结果：

unzip calorine-response-workflows-wave266-*.zip
cd calorine-response-workflows-wave266
bash run.sh

六张参考图在figures，CSV在analysis。若要自己重新计算，在新解压目录以普通账号执行：

bash prepare.sh
bash calculate.sh
bash status.sh

calculate通过Slurm分配一张P100，自动完成本包所有计算；新结果在work，参考图与新算结果分开。已提交时继续查看status，不要删除记录重复提交。结束后：

cat work/COMPLETE.json
ls work/analysis

本例环境为Calorine 4.0、Python 3.11、GPUMD 5.9.1；GPUMD采用GCC 13.3/CUDA 12.9系列环境。可以用module show calorine/4.0和module show gpumd/5.9.1-p100查看当前模块。下面的Python片段在计算完成后的work目录执行；MD和训练由calculate提交。

2 15个读写函数#

资料包native-reader-data保存集群已经完成的真实输出。自己新算的分子、MCMD和dpdt输出会替换到work/Reader-API/native，历史热输运样本保留来源。

3 先看列名，再导出#

cd work
python - <<'CALORINE_READERS'
from calorine.gpumd import read_thermo, read_shc, read_mcmd
root = 'Reader-API/native/'
df = read_thermo(root+'thermo.out', normalize=False)
print(df.columns)
print(df.head())
df.to_csv('my-thermo.csv', index=False)
corr, spec = read_shc(root+'shc.out')
corr.to_csv('my-SHC-correlation.csv', index=False)
spec.to_csv('my-SHC-spectrum.csv', index=False)
print(read_mcmd(root+'mcmd.out').head())
CALORINE_READERS

thermo能量可明确选择总量或每原子量；温度和压力不再按原子数除一次。本例分子极化率模型的标签是整个分子的ų，dump_polarizability保存整体预测，读取使用normalize=False。只有确实需要每原子量时才按3除；不要看到num_atoms就自动归一化。nep预测文件的每原子输出约定与MD的dump文件也要分开看。

历史HAC样本是512原子Si、300 K，晶胞边长21.724 Å；get_rtc_from_hac示例用同一体积和温度作矩形积分。它与原生程序的积分离散方式可能不同，短样本不代表收敛热导。

使用同一HAC样本的实际晶胞与300 K输入，核心积分命令是：

import numpy as np
from ase.io import read
from calorine.gpumd import read_hac
from calorine.tools import get_rtc_from_hac
h = read_hac('Reader-API/native/hac.out',
             exclude_currents=False, exclude_in_out=False)
V = read('Reader-API/native/HAC-model.xyz').get_volume()
dt = float(h.time.iloc[1] - h.time.iloc[0])
kx = get_rtc_from_hac(h.jjx_tot.to_numpy(), V=V, T=300, dt=dt)
np.savetxt('my-running-kappa.csv', np.c_[h.time, kx],
           delimiter=',', header='time_ps,kappa_W_mK', comments='')

4 速度与分组的读写#

read_xyz把GPUMD的Å/fs转换成ASE速度。当前接口会同时保留旧vel数组；再次写入前去掉重复表示，ASE的momenta仍在：

from calorine.gpumd import read_xyz, write_xyz
a = read_xyz('Reader-API/native/model.xyz')
a.arrays.pop('vel', None)
write_xyz('my-model.xyz', a, groupings=[[[0],[1,2]]])

这只修改自己的Atoms对象，不会删除输入文件。参考脚本同时比较了位置、速度和分组，避免把Å/fs与ASE内部单位混用。

5 qNEP的边界与长程设置#

qNEP的这里两类模型需要周期性边界；孤立分子的非周期设置不能直接套用。CPU/GPU精细数值比较要使用相同长程求和方式：GPU默认PPPM 1 Å网格与CPU的Ewald不能当成同一设置。本例显式添加：

kspace ewald
time_step 0.1
ensemble nve
compute_dpdt 5
dump_thermo 5
dump_xyz 5 movie.xyz precision double velocity charge bec force
run 200

这是本包Slurm作业中自动运行的20 fs示例，40个输出时刻。盒子是24 Å周期立方体，不是液态水采样。带速度快照的thermo应力包含动能项，与CPUNEP纯构型应力比较时需分清二者；原始速度不会因为再写一个很小的velocity参数就自动消失。

6 dpdt.out里的P是什么意思#

read_dpdt返回time，单位ps；dPx/dPy/dPz单位eÅ/fs。P列由BEC×速度从零开始按采样步长积分，不能当成绝对sum(q*r)。本包直接用每帧真实输出的BEC和速度重新收缩、积分，与原生列对照。

head -n 9 work/qNEP-dpdt/dpdt.out
cat work/qNEP-dpdt/RESULTS.json
ls work/Reader-API/*.csv

参考Calorine GPUMD读写接口。