CP2K专题

CP2K第一性原理MD温控与续算

Nosé-Hoover短轨迹、扩展守恒量与物理能量，原生restart加新生成波函数的连续轨迹对照。

更新于 2026-10-29约 8 分钟阅读按任务与环境选择适用配置

阅读本篇目录

分子动力学的温控和续算都不只是修改一个关键词。本篇继续使用同一水分子，实际运行Nosé–Hoover温控，并把0.5 fs NVE分成两段，核对续算与连续40步的差别。

下载全部输入、实际轨迹、能量表与分析脚本，解压后进入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是你的新结果。

2 NVT的完整输入#

完整input.inp：water-nvt-nose50#

&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 NVT
  STEPS 50
  TIMESTEP [fs] 0.5
  TEMPERATURE 300
  &THERMOSTAT
   TYPE NOSE
   &NOSE
    TIMECON 20
   &END NOSE
  &END THERMOSTAT
 &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

TYPE NOSE采用Nosé–Hoover温控，TIMECON=20 fs是温控时间尺度，目标温度300 K。这里仅50步×0.5 fs=25 fs，且只有一个水分子，是如何设置和查看温控量的操作例子，不是充分平衡的正则系综采样。

3 区分物理能量和扩展守恒量#

head -3 work/cases/water-nvt-nose50/water-1.ener
tail -3 work/cases/water-nvt-nose50/water-1.ener
cat work/analysis/aimd.csv

本次瞬时温度范围约44.3至695.5 K。一个小分子的瞬时动能涨落很大，目标温度不等于每步打印的温度，25 fs也不能据此评价长期温控分布。

NVT中温控器与物理体系交换能量，因此势能+动能可以变化。water-1.ener的Cons Qty包含温控扩展体系的相关项；本例扩展守恒量波动约0.820668 meV。

真实短时温度轨迹及扩展守恒量与物理能量变化

不能拿本例的瞬时温度最低/最高值当作液态水的温控表现。实际热力学统计应准备合理周期体系、充分平衡并延长采样，检查自相关和统计误差。

4 怎样真正续算，而不是重新给速度#

本组water-half先运行20步，即10 fs；程序自己生成water-1.restart文本状态与water-RESTART.wfn电子检查点。water-resume再读这两个文件，继续20步，得到最终40步/20 fs。

要只跑这一流程：

sbatch submit.slurm water-resume

脚本自动先运行water-half，再运行water-resume，不会直接对空目录提交后半段。文本restart保存MD状态；电子文件在你的计算中生成，下载包不附带二进制波函数。

完整input.inp：water-resume#

&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
  WFN_RESTART_FILE_NAME ../water-half/water-RESTART.wfn
 &SCF
   EPS_SCF 1.0E-7
   MAX_SCF 100
   SCF_GUESS RESTART
   &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

&EXT_RESTART
 RESTART_FILE_NAME ../water-half/water-1.restart
 RESTART_DEFAULT TRUE
&END EXT_RESTART

&MOTION
 &MD
  ENSEMBLE NVE
  STEPS 20
  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

EXT_RESTART/RESTART_DEFAULT读取原子位置、速度、步号及适用的动力学状态；SCF_GUESS RESTART与WFN_RESTART_FILE_NAME读取前半段生成的电子波函数。只拿最后一帧XYZ并重新指定TEMPERATURE，会重新初始化速度，不能当作同一轨迹的续算。

这里STEPS=20表示后半段再运行20步，输出步号继续到40；不要把它误解为读取第20步后不再往前算。

5 检查起止步号和拼接边界#

head -3 work/cases/water-half/water-1.ener
tail -2 work/cases/water-half/water-1.ener
head -3 work/cases/water-resume/water-1.ener
tail -2 work/cases/water-resume/water-1.ener
ls work/cases/water-half/*restart work/cases/water-half/*wfn

与连续40步条件比较，本次分段续算的能量最大差约1.00E-09 Hartree，最终坐标最大差约4.30E-09 Å，处于本次输出舍入和数值精度范围。

实际连续40步与20步加续算20步的能量及差值

分析脚本按步号对齐，并移除前后两段重复的边界行；没有简单把两个文件直接cat起来就当成更长时间序列。更改时间步、系综或方法后，不再假设它与原来的连续轨迹逐点相同。

6 换成实际体系时#

更换结构、元素和电荷后先验证电子自洽与力，再选择步长和温控尺度。NPT需要合理周期体相模型与应力、压强控制，不能直接对这个带真空的单水分子开晶胞涨落。经典周期体系的NVE/NVT/NPT、RDF和MSD流程可结合专题中的氩晶体教程使用；不要把气相短轨迹当成液体扩散数据。

提取与重画#

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。

重启字段见CP2K2026.2 EXT_RESTART手册。

需要根据你的体系选择资源？联系讨论计算需求。实际可用软件与权限，以当前账户和运行环境为准。