CP2K核量子效应与同位素对照

在相同外部势中改变核质量，可以直观看到经典与量子分布的区别。本篇使用已实际运行的14条CP2K PIMD轨迹，比较珠数、质量、两个种子与时间分块，说明怎样从结果文件提取可解释的数据。

1 准备已完成的操作包#

下载原生PIMD、珠数与完整续算操作包，保存为cp2k-pimd.zip。解压后可用sbatch job.slurm提交普通CPU作业；以下命令用于已经通过Slurm分配的CPU计算环境。

CPU分区按56核整节点申请与计费。包内job.slurm已设置单节点、1个Slurm任务、每任务56核；案例脚本在该分配中按各自所需的少量核心运行，每个CP2K进程仍为1线程。

python3 -m zipfile -e cp2k-pimd.zip .
cd cp2k-pimd
bash setup.sh

bash run.sh
bash verify.sh
bash plot.sh

操作包先运行14个珠数、质量与种子条件，再做连续2000步与800＋1200步的完整续算对照。每个独立条件为1个MPI进程、1线程，外层按已分配的物理核心并行，最多12个条件同时运行。结果在results/，图片在plots/。同一输入已完成时会复用结果；修改物理设置后请建立新的工作目录。

环境为CP2K2026.2 CPU模块cp2k/2026.2-8173m，Intel ifx2024.2.1、IntelMPI2021.13、MKL2024.2，配套分析使用已有Python/NumPy环境。

2 对照条件保持哪些量相同#

温度均为300 K，谐振力常数均为0.05 Ha/Bohr²，时间步长0.25 fs，每条10 ps、去掉前2 ps。质量1 u采用P=1、4、8、16、32，质量2 u采用P=1、32；每个条件种子101、202各一条。改变质量时保持势能函数相同。

这里采用自定义质量的谐振粒子，质量对照用于解释核量子效应。真实H/D或其他同位素研究还需采用相应体系的真实质量与势能面。

3 从终端提取本次数据#

source env.sh
"$YS_PY" - <<'PYCODE'
import json
from pathlib import Path
v = json.loads(Path('results/PIMD-VALIDATION.json').read_text())
for r in v['results']:
    print(r['name'], 'P=', r['P'], 'mass/u=', r['mass_u'],
          'mean-square/Bohr^2=', r['variance_Bohr2'],
          'block-SE=', r['block_mean_stderr_Bohr2'])
print('量子/经典展宽比:', v['pooled_quantum_classical_variance_ratio'])
print('轻/重质量展宽比:', v['pooled_H_D_variance_ratio'])
PYCODE

variance_Bohr2是相对势阱中心的逐坐标均方位置，单位Bohr²；换为Ų需乘Bohr到Å的换算系数的平方。gyration_radius_A则是同一时刻珠子相对质心的回转半径，二者含义不同。

4 实际珠数与质量结果#

本次两个种子合并后，质量1 u的P=32位置均方展宽约为P=1的2.84倍；在P=32时，1 u约为2 u的1.35倍。经典极限的同一谐振势位置分布由温度和K控制，质量变化主要体现在动力学；量子分布则与质量有关。

图中的误差条采用每条生产段8个1 ps时间块、两个种子合并后的分块标准误差。帧数多不代表独立样本同样多，正式研究应继续检查块长、自相关与独立种子数量。

5 解析参考为什么有帮助#

一维谐振方向的量子均方位置为ħ/(2mω) coth(ħω/2kBT)，其中ω=sqrt(K/m)；经典极限为kBT/K。三维各向同性模型可以分别对三个方向统计。

有限P还会带来虚时间离散误差。因此图中同时给出有限P配分函数的精确参考和P趋于无穷的量子参考，再叠加真实PIMD结果与采样误差。P=16、32的变化应结合误差条看，不能把有限时间下某个点更接近解析线当作任意真实体系都已收敛的依据。

6 迁移到自己的体系#

先确定温度、势能方法、核质量和目标观测量，再用小规模珠数序列考察收敛。结构分布、氢键、平均能量和量子动能可能需要不同采样时间；对照同位素时保持势能方法、温度、边界条件与分析定义一致。完整输入与文件含义见PIMD珠数操作，长任务接续见完整检查点续算。

参考：CP2K PINT、原生virial估计量实现。