CP2K核磁响应：从水分子计算屏蔽张量

CP2K可以在收敛基态上求解外场微扰响应。这里用原生LINEAR_RESPONSE / LINRES / NMR计算一个水分子的磁屏蔽张量，重点是输入放在哪里、怎样提交和如何读取结果。

1. 下载与实际结果#

下载完整输入、原生结果和绘图包，保存为example.zip：

python3 -m zipfile -e example.zip .
cd cp2k-linear-response-nmr
bash results.sh
bash read.sh
bash plot.sh

这些命令直接读取随包实算结果；reference-results/存原生输出，data/存CSV，figures/存图片。read.sh和plot.sh默认使用集群ASE环境的Python，无需重新计算。自己的Python须安装NumPy、ASE和Matplotlib，可用YS_ANALYSIS_PYTHON指定解释器。

2. 体系与方法#

本例3原子水置于5 Å周期立方盒，BLYP、DZVP-GTH基组、GTH-BLYP赝势，GPW截断150 Ry。基态和线性响应收敛阈值均1E-6，最大100次。它是低成本操作演示，5 Å周期盒与150 Ry未做收敛扫描，不能视作孤立水高精度NMR数据。

GPW赝势结果不包含完整全电子核区响应，尤其O的屏蔽值不能直接与全电子实验数据比较。

3. 完整输入#

cat inputs/water-GPW/input.inp

&GLOBAL
  PRINT_LEVEL LOW
  PROJECT H2O-NMR-1
  RUN_TYPE LINEAR_RESPONSE
&END GLOBAL

&FORCE_EVAL
  &DFT
    BASIS_SET_FILE_NAME GTH_BASIS_SETS
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &MGRID
      CUTOFF 150
    &END MGRID
    &QS
    &END QS
    &SCF
      EPS_SCF 1E-6
      MAX_SCF 100
      SCF_GUESS ATOMIC
    &END SCF
    &XC
      &XC_FUNCTIONAL BLYP
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &PROPERTIES
    &LINRES
      EPS 1E-6
      MAX_ITER 100
      PRECONDITIONER FULL_ALL
      RESTART_EVERY 5
      &CURRENT
        GAUGE R_AND_STEP_FUNCTION
      &END CURRENT
      ##       RESTART
      &LOCALIZE
        MAX_ITER 100
        &PRINT
          &LOC_RESTART
          &END LOC_RESTART
        &END PRINT
        ##           RESTART
      &END LOCALIZE
      &NMR
      &END NMR
      &PRINT
        &RESTART
        &END RESTART
      &END PRINT
    &END LINRES
  &END PROPERTIES
  &SUBSYS
    &CELL
      ABC 5.0 5.0 5.0
    &END CELL
    &COORD
      H       .000689482500      -.765854689700      1.551514984100
      H       .000592940100       .767987328300      1.553204324000
      O       .046239583400       .001896697900      0.9564530746
    &END COORD
    &KIND H
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-BLYP-q1
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-BLYP-q6
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

关键是RUN_TYPE LINEAR_RESPONSE，以及FORCE_EVAL/PROPERTIES/LINRES下的NMR。DFT/SCF/EPS_SCF控制基态，LINRES/EPS控制响应残差；两者各自需要收敛。LOCALIZE准备局域轨道，PRINT/RESTART保存后续可读回的响应文件。

4. 提交#

本例使用cp2k/2026.2-8173m CPU版，Intel oneAPI2024.2/Intel MPI构建。新计算通过Slurm提交：

sbatch run.slurm

脚本申请56核整节点，实际2个MPI进程、每进程1线程，顺序运行5个小条件；不要直接在登录节点进行这些电子结构计算。结束后：

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

5. 提取屏蔽常数#

grep 'SCF run converged' reference-results/water-GPW/native.log
grep 'linear solver converged' reference-results/water-GPW/native.log
cat reference-results/water-GPW/H2O-NMR-1-SHIFT-1.data
cat data/CP2K-NMR-isotropy-anisotropy.csv

SIGMA TOTAL是3×3屏蔽张量；各向同性值为迹除以3。原生TOTAL还包含宏观磁化率贡献，不能只把SOFT J与LOCAL J相加作为完整总值。张量不必逐项对称；本包保留全部9个原生分量。

该版本输出的各向异性取对称部分本征值排序后λ₃−(λ₂+λ₁)/2，不要不分情况直接用笛卡尔σzz−(σxx+σyy)/2替代。数据文件头中的核位置单位为Bohr，张量单位为ppm。

屏蔽常数σ与实验化学位移δ不是同一个量；常见近似δ≈σ_ref−σ_sample还需要相同方法下的参照体系及适当校准。本例没有参照计算，图中只标屏蔽常数。

6. 更换体系前#

先优化结构，再选择GPW/GAPW与合适的核区处理；检查基组、网格、盒子/周期边界、泛函和响应收敛。全电子基组对照见H₂ GAPW例子，文件复用见响应续算。

输入来源为CP2K官方NMR小例，这里收紧基态/响应要求并保留实际输出。