CP2K原生PIMD与珠数选择
从真实路径积分轨迹学习PINT、虚时间珠子、PILE恒温、输出文件与珠数对照，给出完整输入和终端步骤。
阅读本篇目录
研究轻核、低温振动或同位素效应时，核的量子分布可能影响结果。CP2K原生PINT可以用环聚合物珠子表示量子粒子。本篇先用具有解析参考的三维谐振子模型，把输入、运行、输出与珠数选择串起来。
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 珠子代表什么#
P=1给出该模型的经典极限；P=8或32表示同一个粒子在虚时间方向由8或32个珠子表示。它们共同用于采样量子平衡分布，粒子的代表位置可取珠子质心；珠子围绕质心的展宽反映量子离域的一部分。
本例是一个质量1 u或2 u的粒子，处在U=½K(x²+y²+z²)中，K=0.05 Ha/Bohr²。模型在空间中由外部谐振势约束，用于理解路径积分。自己的材料或分子计算需要替换为相应的DFT或可靠势能方法。可选取一份input.inp复制到新目录，在作业脚本加载同一模块后，用cp2k.psmp -i input.inp -o output.log运行。配套run/verify脚本用于本文固定对照，新的体系应按原子数和目标量调整分析。
3 一份实际执行的完整输入#
下面为质量1 u、P=8、种子101的inputs/mass1-P8-s101/input.inp：
&GLOBAL
PROJECT mass1-P8-s101
RUN_TYPE PINT
PRINT_LEVEL MEDIUM
SEED 20261109
&END GLOBAL
&MOTION
&PINT
DT 0.25
HARM_INT EXACT
NRESPA 1
NUM_STEPS 40000
P 8
PROC_PER_REPLICA 1
TEMP 300
TRANSFORMATION NORMAL
&PILE
TAU [fs] 10
&END PILE
&INIT
LEVY_POS_SAMPLE TRUE
LEVY_SEED 101
THERMOSTAT_SEED 101
&END INIT
&PRINT
&ENERGY
FILENAME =pint-energy.dat
&EACH
PINT 1
&END EACH
&END ENERGY
&CENTROID_POS
FILENAME =centroid.xyz
UNIT angstrom
&EACH
PINT 10
&END EACH
&END CENTROID_POS
&CENTROID_GYR
FILENAME =gyration.xyz
UNIT angstrom
&EACH
PINT 10
&END EACH
&END CENTROID_GYR
&END PRINT
&END PINT
&PRINT
&RESTART_HISTORY OFF
&END RESTART_HISTORY
&RESTART
BACKUP_COPIES 0
&EACH
PINT 10000
&END EACH
&END RESTART
&TRAJECTORY
FORMAT XYZ
FILENAME bead-pos
&EACH
PINT 10
&END EACH
&END TRAJECTORY
&FORCES
FORMAT XYZ
FILENAME bead-force
UNIT hartree/bohr
&EACH
PINT 10
&END EACH
&END FORCES
&END PRINT
&END MOTION
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
DO_ELECTROSTATICS FALSE
&CHARGE
ATOM H
CHARGE 0
&END CHARGE
&NONBONDED
&LENNARD-JONES
ATOMS H H
EPSILON 0.0
SIGMA [angstrom] 1
RCUT [angstrom] 4
&END LENNARD-JONES
&END NONBONDED
&END FORCEFIELD
&POISSON
PERIODIC NONE
&EWALD
EWALD_TYPE NONE
&END EWALD
&END POISSON
&END MM
&EXTERNAL_POTENTIAL
ATOMS_LIST 1
FUNCTION 0.5*K*(X^2+Y^2+Z^2)
PARAMETERS K
VALUES 0.05
&END EXTERNAL_POTENTIAL
&SUBSYS
&CELL
ABC 20 20 20
PERIODIC NONE
&END CELL
&TOPOLOGY
CONNECTIVITY OFF
&END TOPOLOGY
&KIND H
ELEMENT H
MASS 1.0
&END KIND
&COORD
H 0.1 0.1 0.1
&END COORD
&END SUBSYS
&END FORCE_EVAL
RUN_TYPE PINT进入路径积分；P决定珠数；TEMP=300设置物理温度；DT=0.25以fs为单位，40000步为10 ps。前2 ps用于平衡，后8 ps分为8个1 ps时间块做统计。
TRANSFORMATION NORMAL采用正规模坐标，HARM_INT EXACT用于珠子间谐和弹簧的积分；本例配套PILE恒温，TAU=10 fs设置质心恒温时间常数。LEVY_POS_SAMPLE生成初始珠子分布，两个随机种子101、202分别采样。RESTART_HISTORY OFF控制历史文件数量，周期性完整检查点仍按RESTART设置输出。
PROC_PER_REPLICA=1在这里配合每个条件一个MPI进程，珠子的力计算由该进程依次完成。P决定虚时间离散精度，申请CPU核心数则决定并行方式；配套脚本在不同参数条件之间并行。
4 输出文件怎么看#
source env.sh
ls results/cases/mass1-P8-s101/*bead-pos*.xyz
head -n 6 results/cases/mass1-P8-s101/pint-energy.dat
"$YS_PY" - <<'PYCODE'
import numpy as np
x = np.load('results/cases/mass1-P8-s101/fields.npz')
print('记录帧数:', len(x['steps']))
print('位置数组:', x['positions_A'].shape)
print('最后时间 / fs:', x['time_fs'][-1])
PYCODE
bead-pos-pos-1-1.xyz、...-2-1.xyz等分别对应不同珠子；centroid.xyz给出质心，gyration.xyz给出珠子围绕质心的回转半径。本例位置数组形状为(4001,P,3)，第二轴为珠子序号，末轴为xyz分量。位置和力每10步记录一次，0–40000步共4001帧；能量每一步输出，共40001行。
能量文件列为步数、时间fs、热力学动能估计量、virial动能估计量、扩展体系温度、物理势能、守恒量和每步耗时。量子动能应选择相应估计量，并结合误差分析；不要只把扩展珠子速度的经典动能当作核的量子动能。
5 怎样确定珠数#
本例对质量1 u计算P=1、4、8、16、32，每个P使用两个独立种子。应同时看珠数离散误差与有限采样误差，而不是只看某一条轨迹最后一帧。珠数和质量的实际结果见核量子效应与同位素对照。
PIMD在扩展相空间采样静态量子分布，积分时间用于采样过程。讨论真实时间的扩散、振动关联等动力学量时，应选择并验证相应的动力学方法。完整检查点接续见PIMD续算，不同温度之间的交换见温度副本交换。