CP2K第一性原理MD与时间步

从静态结构计算切换到分子动力学，关键是每一步都得到可用的电子基态与原子力，并检查积分步长。本篇用一个水分子的短轨迹演示Born–Oppenheimer MD：电子按PBE自洽，原子核按经典动力学运动。

下载全部输入、实际轨迹、能量表与分析脚本，解压后进入run.sh目录。使用CP2K2026.2 P100双精度，GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3；分析用ASE3.29.0。

1 从提交开始#

sbatch submit.slurm

模板申请一张P100和一个CPU核，每个计算为一个MPI进程、一个OpenMP线程。整组包含三个NVE时间步、一个NVT，以及半程/续算两段，共五个教学条件、六次程序运行。默认依次计算，最后提取画图；在已经分配的多卡资源里，脚本按可见GPU和CPU亲和范围每卡安排一个任务，续算等待半程完成。

核时模板为gpu-8p100-hs；套餐用户用yskit生成所属分区脚本，在已分配资源中保留bash run.sh。不要在登录节点计算。已有work时重新解压；reference是本文结果，work是你的新结果。

只先跑0.5 fs条件，可以提交：

sbatch submit.slurm water-nve-dt05

2 使用同一初始结构和速度#

起点来自水分子振动教程的实际优化结构。12 Å非周期MT盒，PBE/GTH-PBE、DZVP-MOLOPT-SR-GTH、400/50 Ry网格、EPS_SCF=1E-7。GLOBAL固定SEED=20261029，初始速度温度300 K；三种NVE步长使用相同初始结构、动能和势能。

这里是一个气相水分子，不是液态水盒。20 fs足以看到快的O—H振动与积分差异，不用于估计液体扩散、热导或平衡分布，也不是路径积分核量子动力学。

完整input.inp：water-nve-dt05#

&GLOBAL
 PROJECT water
 RUN_TYPE MD
 SEED 20261029
 PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
 METHOD QS
 &DFT
  BASIS_SET_FILE_NAME ../DATA/BASIS_MOLOPT
  POTENTIAL_FILE_NAME ../DATA/GTH_POTENTIALS
  &QS
   EPS_DEFAULT 1.0E-12
  &END QS
  &MGRID
   CUTOFF 400
   REL_CUTOFF 50
  &END MGRID
  &SCF
   EPS_SCF 1.0E-7
   MAX_SCF 100
   SCF_GUESS ATOMIC
   &OT
    MINIMIZER DIIS
    PRECONDITIONER FULL_SINGLE_INVERSE
   &END OT
  &END SCF
  &XC
   &XC_FUNCTIONAL PBE
   &END XC_FUNCTIONAL
  &END XC
  &POISSON
   PERIODIC NONE
   POISSON_SOLVER MT
  &END POISSON
 &END DFT
 &SUBSYS
  &CELL
   ABC 12 12 12
   PERIODIC NONE
  &END CELL
  &TOPOLOGY
   COORD_FILE_NAME optimized.xyz
   COORD_FILE_FORMAT XYZ
  &END TOPOLOGY
  &KIND O
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q6
  &END KIND
  &KIND H
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q1
  &END KIND
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
&END FORCE_EVAL

&MOTION
 &MD
  ENSEMBLE NVE
  STEPS 40
  TIMESTEP [fs] 0.5
  TEMPERATURE 300
 &END MD
 &PRINT
  &TRAJECTORY
   FORMAT XYZ
   &EACH
    MD 1
   &END EACH
  &END TRAJECTORY
  &VELOCITIES
   &EACH
    MD 1
   &END EACH
  &END VELOCITIES
  &RESTART
   BACKUP_COPIES 0
   &EACH
    MD 10
   &END EACH
  &END RESTART
 &END PRINT
&END MOTION

RUN_TYPE MD启用动力学；ENSEMBLE NVE不连接温控器；TIMESTEP以fs为单位。改变步长时同时改变步数，本组0.25 fs×80、0.5 fs×40、1 fs×20都对应20 fs，避免把不同总时长直接混在一起比较。

3 先检查每一步电子自洽#

grep -E 'SCF run converged|SCF run NOT converged|PROGRAM ENDED' work/cases/water-nve-dt05/output.log
head -3 work/cases/water-nve-dt05/water-1.ener
tail -3 work/cases/water-nve-dt05/water-1.ener
head -7 work/cases/water-nve-dt05/water-pos-1.xyz

water-1.ener依次列步号、时间fs、动能Hartree、瞬时温度K、势能Hartree、守恒量Hartree、单步耗时s。water-pos-1.xyz是逐帧坐标；不能把ENERGY行末尾的势能误当成NVE应守恒的总能。

NVE里检查势能+动能。分子振动时两者相互转换，瞬时温度随动能变化；没有温控器的NVE不会强制每一步保持300 K。

4 实际时间步对照#

0.5 fs能量波动约为0.25 fs的四倍，1 fs又约为0.5 fs的四倍，符合这组短轨迹中二阶积分误差的趋势。这里报告最大值−最小值，不把这种有界振荡称为已经测准的长期能量漂移率。

键长曲线也来自实际XYZ帧。对于自己的体系，需要根据最快振动、约束、温度、电子精度和目标时长检查步长；本水分子的结果不自动等于金属、重元素或含约束体系的推荐步长。

5 用ASE提取某一根键#

module load ase/3.29.0
python - <<'PY'
from ase.io import read
frames = read('work/cases/water-nve-dt05/water-pos-1.xyz', index=':')
print('frames:', len(frames))
print('first O-H / angstrom:', frames[0].get_distance(0, 1))
print('last O-H / angstrom:', frames[-1].get_distance(0, 1))
PY

ASE原子编号从0开始；本XYZ顺序为O、H、H，因此0和1对应第一根O—H。换结构后先检查原子顺序，不要机械沿用编号。

提取与重画#

module load ase/3.29.0
python analyze.py work
python plot_results.py work
python check_results.py work
cat work/analysis/aimd.csv
ls work/figures

CSV列出步长、步数、系综、结束时间、守恒量波动与瞬时温度范围；JSON保留续算差值。只查看本文数据时将work改为reference。

继续看温控与真正续算。MD参数见CP2K2026.2 MD手册。