CP2K Ehrenfest动力学：电子与原子核一起运动

RT_PROPAGATION可以在固定原子核条件下传播电子；EHRENFEST_DYN则将电子和经典原子核耦合推进。本篇用两个H原子的短任务展示后者，重点是实际操作与输出读取。

1. 下载与看真实轨迹#

下载实际输入、原生结果与数据绘图包，保存为example.zip：

python3 -m zipfile -e example.zip .
cd cp2k-ehrenfest-H2
bash results.sh
bash read.sh
bash plot.sh

上述命令先读随包真实结果；CSV在data/，原生输出在reference-results/，图在figures/。默认使用集群ASE环境中的Python，可用YS_ANALYSIS_PYTHON指定含NumPy、ASE和Matplotlib的解释器。

本例初始H–H距离0.9 Å，故意取拉伸结构，初速度为0。没有外场时，核也会因初始受力而运动；加短脉冲后电子响应与核轨迹改变。三项计算仅传播0.483777 fs，用于演示操作，不表示得到真实光解、寿命或收敛吸收光谱。

2. 方法与完整输入#

PBE、DZVP-MOLOPT-GTH/GTH-PBE-q1，8 Å非周期盒、MT Poisson，200 Ry网格与40 Ry相对截断。基态EPS_SCF 1E-6，实时迭代EPS_ITER 1E-8。以有场0.5原子时间单位步长为例：

cat inputs/H2-pulse/input.inp

&GLOBAL
  PROJECT H2-EMD
  RUN_TYPE EHRENFEST_DYN
  PRINT_LEVEL LOW
  SEED 42
&END GLOBAL
&MOTION
  &MD
    ENSEMBLE NVE
    STEPS 40
    TEMPERATURE 0
    TIMESTEP [au_t] 0.5
  &END MD
  &PRINT
    &TRAJECTORY
      &EACH
        MD 1
      &END EACH
    &END TRAJECTORY
    &VELOCITIES
      &EACH
        MD 1
      &END EACH
    &END VELOCITIES
  &END PRINT
&END MOTION
&FORCE_EVAL
  METHOD QUICKSTEP
  &DFT
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &EFIELD
      ENVELOP GAUSSIAN
      INTENSITY 1E12
      POLARISATION 0 0 1
      WAVELENGTH 210
      &GAUSSIAN_ENV
        SIGMA [au_t] 2.0
        T0 [au_t] 5.0
      &END GAUSSIAN_ENV
    &END EFIELD
    &MGRID
      CUTOFF 200
      REL_CUTOFF 40
    &END MGRID
    &QS
      EPS_DEFAULT 1E-10
    &END QS
    &POISSON
      PERIODIC NONE
      POISSON_SOLVER MT
    &END POISSON
    &SCF
      EPS_SCF 1E-6
      MAX_SCF 100
      SCF_GUESS ATOMIC
    &END SCF
    &REAL_TIME_PROPAGATION
      EPS_ITER 1E-8
      INITIAL_WFN SCF_WFN
      MAT_EXP TAYLOR
      MAX_ITER 100
      PROPAGATOR ETRS
    &END REAL_TIME_PROPAGATION
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC 8 8 8
      PERIODIC NONE
    &END CELL
    &COORD
      H 4 4 3.55
      H 4 4 4.45
    &END COORD
    &VELOCITY
      0 0 0
      0 0 0
    &END VELOCITY
    &KIND H
      BASIS_SET DZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q1
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

GLOBAL/RUN_TYPE EHRENFEST_DYN进入耦合动力学；MOTION/MD定义核步数、时间步、NVE；DFT/REAL_TIME_PROPAGATION定义电子传播器。两部分缺一不可。INITIAL_WFN SCF_WFN先求基态再作为实时传播起点，每一步并不是重新优化电子到瞬时基态的BOMD。

SUBSYS/VELOCITY显式给两原子零初速度，TEMPERATURE 0不引入随机热运动；后续出现非零核动能来自受力运动。两原子沿z放置，脉冲也沿z；此例INTENSITY 1E12是原生光强参数，采用W/cm²，WAVELENGTH 210采用nm。高斯包络σ=2 au、中心5 au是很短的数学脉冲，未按某种实验激光或真实分子反应标定。

3. 提交#

本例使用cp2k/2026.2-8173m，Intel Fortran2024.2.1与Intel MPI构建。新解压目录提交：

sbatch run.slurm

脚本申请56核整节点，以2个MPI进程、每进程1线程顺序运行三个小条件。真正电子结构与动力学计算放到Slurm任务，不在登录节点直接运行。结束后提取自己的输出：

bash read.sh work analysis-new
bash plot.sh analysis-new figures-new work

原生底层调用形式为：

module load cp2k/2026.2-8173m
export OMP_NUM_THREADS=1
mpirun -np 2 cp2k.psmp -i input.inp > native.log

本包自动准备H2-zero、H2-pulse和H2-pulse-halfstep三个工作目录。

4. 先确认基态，再读实时与轨迹#

grep 'SCF run converged' reference-results/H2-pulse/native.log
grep 'MD STEP:' reference-results/H2-pulse/native.log
head -n 12 reference-results/H2-pulse/H2-EMD-pos-1.xyz
head reference-results/H2-pulse/H2-EMD-1.ener
cat data/CP2K-Ehrenfest-H2-conditions.csv

pos-1.xyz坐标为Å，ener时间列为fs，势能和动能为Hartree。本例40个核步加初始帧，共41帧；半时间步80步加初始帧，共81帧。两个核坐标都有实际变化，不是固定核电子传播。

5. 怎样理解能量#

无场任务的E_pot+E_kin总跨度约2.4×10⁻⁸ Eh。加场时外场可以做功，不要求这两项之和保持常数；原子核质心也不能单独按无场守恒条件判断。不要把有场能量变化直接叫作数值漂移。

实际脉冲与无场的最大键长差约2.88×10⁻⁴ Å，只说明这组输入产生可辨认响应，不能据此推断反应已发生。

下一步见时间步与结果读取。输入参考CP2K官方Ehrenfest小例，本篇重新设置几何、边界、初速度、收敛和脉冲，并实际完成上述对照。