Calorine轨迹分析与GPUMD数据读取

下载本专题结构、势文件、真实参考数据和全部Python脚本。以下命令在解压目录执行；CPU计算使用已分配的CPU资源，GPU通过Slurm分配。原始结果放在 reference/，复跑结果写入当前目录。

1 Si短NVE轨迹#

在NEP计算器上使用ASE积分器，64个Si、300 K初始速度、随机种子28、1 fs步长、1000步；去掉整体平动。温度未加恒温器持续控制。

import numpy as np
from ase.build import bulk
from ase.units import fs
from ase.md.verlet import VelocityVerlet
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary
from calorine.calculators import CPUNEP
a = bulk('Si','diamond',a=5.431,cubic=True).repeat((2,2,2))
a.calc = CPUNEP('si-nep.txt')
MaxwellBoltzmannDistribution(a, temperature_K=300,
                            rng=np.random.default_rng(28))
Stationary(a)
md = VelocityVerlet(a, timestep=fs, logfile='nve-new.log')
md.run(1000)

python cpu_workflows.py
python extra_workflows.py
python plot.py
head -n 5 nve.log

完整脚本保存101个结构帧，每10步一帧；热力学日志每一步记录。最大总能量漂移3.95E-5 eV/atom。本例用于理解接口与短时能量检查，不提供收敛扩散或热导率。

2 相关性与误差估计#

from calorine.tools import (analyze_data, get_autocorrelation_function,
                            get_correlation_length, get_error_estimate)
# E为单段日志中逐步总能量的一维数组
print(analyze_data(E, max_lag=50))
acf = get_autocorrelation_function(E, max_lag=50)
print(get_correlation_length(E), get_error_estimate(E))

函数按样本间隔给出相关长度；必须乘实际记录间隔才能换成时间。短NVE的振荡和守恒约束不同于平衡涨落采样，相关估计不应直接宣称为严格独立热力学误差。若相关函数未降到软件判据以内，返回None或NaN有明确含义，不能用0代替。

3 读取已有GPUMD真实输出#

包内real-gpumd/选自已发布的Si MD、液态Si输运和微型NEP训练原始输出，不额外编造曲线。脚本实际运行8类读取器：

from calorine.gpumd import read_thermo, read_msd, read_shc
thermo = read_thermo('real-gpumd/thermo.out')
msd = read_msd('real-gpumd/msd.out')
correlation, spectrum = read_shc('real-gpumd/shc.out')
print(thermo.columns)
print(msd.head())
correlation.to_csv('shc-time.csv', index=False)
spectrum.to_csv('shc-frequency.csv', index=False)

不同读取器的返回类型不完全一致；read_shc返回两个DataFrame，不能当成单表。GPUMD多段run会追加部分相关输出，读取器可能要求单段块，不能把不同段相关函数直接拼接解释。

4 读写run.in与重绘#

from calorine.gpumd import read_runfile, write_runfile
settings = read_runfile('real-gpumd/run.in')
write_runfile('new-run.in', settings)
print(settings)

python extra_workflows.py
ls parsed-*.csv parity-*.csv nve-acf.csv
python plot.py

本教程覆盖本次明确实跑的API；TNEP偶极/极化率、IR/Raman光谱及NEP89专用微调是另一些输入与模型条件，不能由能量型Si势自动获得。