CP2K相关方法力的差分核对

检查解析力时，可以把原子向正、负方向各移动一点，计算总能量斜率。本篇实际使用0.002和0.001 Å两个步长，核对同一水分子H1的z分量。

下载解析力、两组位移输入、实际日志与工具脚本，解压进入run.sh目录。CP2K2026.2 P100双精度：GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3；分析ASE3.29.0。

1 从输入生成到提交#

python3 make_displacements.py
sbatch submit.slurm

make_displacements.py从两个基准输入生成8个正负位移输入，只移动指定原子，不重新居中。run.sh也会自动调用它；模板默认完成10次运行：RI-MP2与Minimax RI-RPA各一次解析力，加上每种方法两种步长的正负位移单点。

模板为核时P100分区，一GPU一CPU，每个条件一个MPI与一个OpenMP线程。套餐用户用yskit生成所属分区脚本，在已分配资源中保留bash run.sh。多卡仅使用已经分配的可见GPU/CPU亲和范围，每卡一个任务，空闲卡接下一项；不自动申请资源。不要在登录节点计算，已有work重新解压。

inputs是新计算输入；reference是本文实际原始输出；work是新结果。只运行一个条件可在sbatch后附CASES.json中的名称，例如：

sbatch submit.slurm water-mp2-force

单条件只输出自己的数据；完整力图和差分对照需要整组。

2 位移前先固定整体位置#

本例基准XYZ已经写入固定的居中坐标，TOPOLOGY中不再使用CENTER_COORDINATES。CP2K这个自动居中功能依据包围盒中点；改变一个原子后，包围盒中点可能改变，其他原子也随之平移。

要核对“只移动H1-z”的导数，就应只改变这一坐标，保持其他原子和晶胞不变。make_displacements.py从固定基准复制结构、修改指定原子分量，并保留同一方法的其余数值设置。

ASE/脚本索引从0开始，H1是atom_index=1，z是axis=2。CASES.json记录每个目录的正负位移；修改原子/分量时先检查XYZ顺序，再重新生成整套输入。

3 一个正位移单点的完整输入#

完整input.inp：water-mp2-plus#

&GLOBAL
  PRINT_LEVEL MEDIUM
  PROJECT RI_MP2_H2O
  RUN_TYPE ENERGY
  &TIMINGS
    THRESHOLD 0.01
  &END TIMINGS
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep
  &DFT
    BASIS_SET_FILE_NAME ../DATA/HFX_BASIS
    POTENTIAL_FILE_NAME ../DATA/POTENTIAL
    &MGRID
      CUTOFF 300
      REL_CUTOFF 50
    &END MGRID
    &POISSON
      PERIODIC NONE
      POISSON_SOLVER WAVELET
    &END POISSON
    &QS
      EPS_DEFAULT 1.0E-12
      EPS_PGF_ORB 1.0E-18
      METHOD GPW
    &END QS
    &SCF
      EPS_SCF 1.0E-7
      MAX_SCF 100
      SCF_GUESS ATOMIC
    &END SCF
    &XC
      &HF
        FRACTION 1.0000000
        &SCREENING
          EPS_SCHWARZ 1.0E-8
              EPS_SCHWARZ_FORCES 1.0E-8
          SCREEN_ON_INITIAL_P FALSE
        &END SCREENING
      &END HF
      &WF_CORRELATION
        &CANONICAL_GRADIENTS
          EPS_CANONICAL 1.0E-7
          FREE_HFX_BUFFER FALSE
          &CPHF
            EPS_CONV 1.0E-8
            MAX_ITER 50
          &END CPHF
        &END CANONICAL_GRADIENTS
        MEMORY 200.
        NUMBER_PROC 1
        &INTEGRALS
          &WFC_GPW
            CUTOFF 300
            EPS_FILTER 1.0E-8
            EPS_GRID 1.0E-8
            REL_CUTOFF 50
          &END WFC_GPW
        &END INTEGRALS
        &RI_MP2
        &END RI_MP2
      &END WF_CORRELATION
      &XC_FUNCTIONAL NONE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC [angstrom] 14.0 14.0 14.0
      PERIODIC NONE
    &END CELL
    &KIND H
      BASIS_SET DZVP-GTH
      BASIS_SET RI_AUX RI_DZVP-GTH
      POTENTIAL GTH-HF-q1
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH
      BASIS_SET RI_AUX RI_DZVP-GTH
      POTENTIAL GTH-HF-q6
    &END KIND
    &TOPOLOGY
      COORD_FILE_FORMAT xyz
      COORD_FILE_NAME H2O_gas.xyz
    &END TOPOLOGY
  &END SUBSYS
  &PRINT
    &FORCES ON
    &END FORCES
  &END PRINT
&END FORCE_EVAL

基准force目录求解析力；位移目录用RUN_TYPE ENERGY求总能，用于构造斜率。这里保留相同的方法参数，梯度求解设置只在请求力时使用。

4 中心差分与单位#

Fz≈−[E(z+δ)−E(z−δ)]/(2δ)。本例能量原生为Hartree，位移为Å，所以计算能量差分时乘Hartree到eV即可得到eV/Å；对比的解析力原生为Hartree/Bohr，则需要同时换算长度。

grep 'ENERGY|' work/cases/water-mp2-minus/output.log
grep 'ENERGY|' work/cases/water-mp2-plus/output.log
grep 'ENERGY|' work/cases/water-rpa-minus-001/output.log
grep 'ENERGY|' work/cases/water-rpa-plus-001/output.log

分析脚本读取完整总能，包含该方法的相关项；不只用SCF参考能来对比完整相关力。

5 实际结果#

本例两个步长下均吻合到约3E-5 eV/Å以内；0.001 Å反而比0.002 Å的数值差略大。位移变小时，有限电子求解和能量打印等噪声也更容易被除法放大，不能只认为步长越小必然越好。

这是当前固定输入中一个力分量的一致性检查，不替代轨道/RI基组、积分网格、盒长、电子收敛和其他力分量的检查，也不表示新的方法已经得到优化结构。

6 换成自己的体系#

先更新两个方法的基准XYZ与KIND、赝势和电荷设置，确认整体位置固定；再检查CASES.json中的原子编号、分量与步长，运行make_displacements.py，最后提交。每对正负位移都用相同方法与数值参数，不在位移之间重新优化其余原子。

提取与重画#

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

CSV包含方法、位移步长、解析力、差分力和两者之差，力单位eV/Å；JSON还保存三个原子的完整力向量。看本文数据时将work改成reference。

基准输入与解析力读取见MP2与RPA解析力。