CP2K液体扩散实操：MSD拟合与统计
从四条完整氩轨迹解缠坐标、多起点平均MSD，比较两个拟合窗口，报告真实扩散值与样本标准差。
阅读本篇目录
完成预平衡与生产轨迹后，本篇说明怎样从自己的CP2K输出得到扩散参考值，以及哪些误差仍需继续检查。
1. 先复算参考曲线#
python3 -m zipfile -e cp2k-Ar108-liquid-diffusion-wave307-*.zip .
cd cp2k-Ar108-liquid-diffusion-wave307
bash run.sh
cat results/DATA.json
head results/seed1-MSD.csv
全部分析以实际1201帧108原子坐标和完整能量表为输入。bash run.sh不重新运行MD，也不从图片上手工取点。
2. 读取自己的CP2K轨迹#
若已用上一篇脚本在work/完成四组计算：
module purge
module load pymatgen/2026.9.24
python read_cp2k_liquid.py work --output my-raw
python analyze.py --data my-raw --output my-results
python plot.py --data my-results --output my-plots
cat my-results/DATA.json
读取器读取seedN-prod-pos-1.xyz、seedN-prod-1.ener以及预平衡能量表。当前分析器对应108个Ar、18 Å盒、2 fs步长、每20 fs输出和24 ps生产；换体系时同步修改这些范围、原子选择与拟合窗口，不将其他长度轨迹硬改成相同形状。
3. 周期坐标先解缠#
粒子穿过周期边界时，如果直接用包裹坐标相减，MSD会出现虚假的大跳跃。本例利用相邻20 fs帧的最小镜像位移累加得到连续坐标，再逐帧扣除整个体系的质心平移。
delta = positions[1:] - positions[:-1]
delta -= box_length * np.rint(delta / box_length)
unwrapped = np.concatenate([
positions[:1],
positions[:1] + np.cumsum(delta, axis=0)
])
relative = unwrapped - unwrapped.mean(axis=1)[:, None, :]
本例是固定正交盒，box_length=18.0 Å；倾斜盒、NPT变胞或刚体整体移动要采用相应处理，不能直接逐笛卡尔分量套用。保存坐标过稀会使相邻帧移动超过半盒长，从而无法唯一恢复跨越次数。
4. 用多个时间起点求MSD#
对每个时间间隔τ，对所有能组成该间隔的时间起点和全部108个原子平均：
MSD(τ) = < |r_i(t+τ) − r_i(t)|² >_(i,t)
D = (MSD对时间的线性斜率) / 6
分母6对应三维各向同性扩散，二维通常为4，不能直接照抄。短时间弹道区也不适合用同一线性扩散公式。本例主拟合2–6 ps，另用1–4 ps检查窗口依赖，拟合保留自由截距，不强迫直线过原点。
5. 本次实际数值#
| 初速度样本 | 生产平均温度 / K | D，2–6 ps / Ų·ps⁻¹ | D，1–4 ps / Ų·ps⁻¹ |
|---|---|---|---|
| 1 | 124.2129 | 0.522918 | 0.512977 |
| 2 | 122.1886 | 0.448230 | 0.455608 |
| 3 | 123.3966 | 0.591913 | 0.573710 |
| 4 | 122.9525 | 0.437628 | 0.438566 |
四个主窗口结果的平均值为0.500172 Ų/ps，即5.00172×10⁻⁹ m²/s；样本标准差为0.071980 Ų/ps，四样本平均值的标准误为0.035990 Ų/ps。标准差与标准误描述不同量，不要都标成同一种“误差”。
1 Ų/ps = 10⁻⁸ m²/s
本组四个样本的相对标准差约14.4%。换拟合窗口所得变化与独立样本差异都应保留，不挑其中一个值当成所有液体的精确答案。MSD时间起点互相重叠，因此不能把所有起点视为独立重复样本，再据此声称误差极小。
6. RDF怎样帮助判断结构#
cat results/seed1-RDF.csv | head
ls plots/Ar108-liquid-structure-RDF.png
RDF用生产后半段12 ps的实际原子对最小镜像距离统计。附件采用有限N的N(N−1)归一，不包含原子与自己的配对；若与采用N²归一的工具比较，应先统一约定。RDF的宽峰和MSD的扩散区间应结合起来看，不能只因某一条曲线有峰就宣布液体已充分平衡。
7. 怎样把这个流程用于自己的体系#
- 选择适合材料与相态的势、密度、温度和边界条件，先得到合理预平衡状态。
- 保存足够密的轨迹，并在生产阶段检查能量及温度趋势。
- 先画MSD，再选择扩散线性区；比较多个窗口和不同起始速度。
- 增加生产时长、独立样本和体系尺寸，观察D及误差是否稳定。
- 含固体基底、固定原子或多组分时，分别选择真正移动的物种，并按所需物理方向计算。
本次没有做体系尺寸外推、黏度相关有限尺寸修正或充分长时间极限收敛，也没有把LJ模型的结果标成第一性原理或实验氩扩散系数。需要晶体振动、短轨迹与续算的入门对照，可回看续算、RDF与MSD。