CP2K 同位素与位移步长

更换同位素通常改变核质量，从而改变质量加权振动频率；在这里的Born–Oppenheimer电子结构模型中，电子势能面不因H换成D而改变。本篇用同一优化水结构比较质量与有限位移步长。

下载完整输入、原始输出、频率表和绘图脚本，解压后进入 run.sh 所在目录。本文在鸢算科技P100上用CP2K2026.2双精度 cp2k.psmp 实算；构建采用GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3。分析环境为ASE3.29.0。

1 提交与查看结果#

sbatch submit.slurm heavy-water

模板申请四张P100、四个MPI进程，每进程一个CPU核和一个OpenMP线程。NPROC_REP 1让不同进程计算不同位移，脚本按MPI全局编号分配GPU。在已分配的计算资源里也可以执行 bash run.sh heavy-water。套餐用户用 yskit 生成对应分区的提交脚本，再保留这条运行命令；不要在登录节点直接运行。

脚本依次执行“优化结构 → 提取最终XYZ → 有限位移力、偶极与极化率 → 振动模式 → 数据提取 → 画图”。inputs为可编辑输入，reference为本文结果，work为本次新结果。每个条件请新解压一份，脚本不会覆盖已有work。

2 只改变氢质量#

重水输入把KIND H中的 MASS改为2.01410177812 u，元素仍为H，仍用H的GTH-PBE-q1赝势。普通水使用CP2K质量H=1.00794 u、O=15.9994 u。不要把D误当成一种需要新赝势的元素，也不要同时改变几何、泛函和网格后把差异全归于质量。

完整输入：d2o-ir-raman#

&GLOBAL
 PROJECT water
 RUN_TYPE VIBRATIONAL_ANALYSIS
 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 600
   REL_CUTOFF 50
  &END MGRID
  &SCF
   EPS_SCF 1.0E-8
   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
  &PRINT
   &MOMENTS ON
    PERIODIC .FALSE.
    REFERENCE COM
   &END MOMENTS
  &END PRINT
 &END DFT
 &PROPERTIES
  &LINRES
   EPS 1.0E-8
   MAX_ITER 100
   PRECONDITIONER FULL_SINGLE_INVERSE
   &POLAR ON
    DO_RAMAN .TRUE.
    PERIODIC_DIPOLE_OPERATOR .FALSE.
   &END POLAR
  &END LINRES
 &END PROPERTIES
 &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
   MASS 2.01410177812
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q1
  &END KIND
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
&END FORCE_EVAL
&VIBRATIONAL_ANALYSIS
 NPROC_REP 1
 DX 0.005
 FULLY_PERIODIC .FALSE.
 INTENSITIES .TRUE.
 THERMOCHEMISTRY .TRUE.
 TC_TEMPERATURE 298.15
 TC_PRESSURE 101325
 &PRINT
  &MOLDEN_VIB ON
  &END MOLDEN_VIB
 &END PRINT
&END VIBRATIONAL_ANALYSIS

优化仍先得到电子势能的极小点，重水质量在振动阶段参与质量加权。本例没有考虑核量子效应导致的平均键长改变。

3 比较两个DX#

另一份新目录执行 sbatch submit.slurm step，采用普通水质量，只将 DX 0.005改为 DX 0.01。DX单位是Bohr，并非Å；0.01 Bohr约0.005292 Å。更小步长会放大电子力噪声，更大步长会增加非谐有限差分误差，不能把“越小越好”作为通用规则。

grep -E 'DX|EPS_SCF|MASS' inputs/water-dx001/input.inp
grep -E 'VIB\|Frequency|VIB\|IR int|VIB\|Raman' work/cases/water-dx001/output.log

4 实际对比#

以下频率单位均为cm⁻¹。

本例两个DX的最大频率差约0.422 cm⁻¹；重水三条频率均降低，但各模式比例不同，不能把所有模式都直接除以√2。正常模式由整个分子的质量矩阵和力常数组成。

提取与重画#

module load ase/3.29.0
python analyze.py work
python plot_results.py work
ls work/analysis
ls work/figures

CSV给出模式编号、频率、IR强度和Raman活性；JSON保留质量、最终结构和热修正。原生 water-VIBRATIONS-1.mol 是Molden振动文件，可以用支持该格式的分子查看器检查模式位移。项目名前缀water是本包统一输出名，不代表CO₂输入成了水。

只重画本文结果，把上述命令中的work换成reference。脚本用实际模式的线位置与活性叠加高斯曲线，展示宽度σ=10 cm⁻¹；这是绘图展宽，不是计算得到的声子寿命。

如果出现意料外的小虚频，先确认优化结构、电子力的收敛和步长稳定性，再判断是否确有不稳定模式。两个步长的本例对照仅说明当前小分子参数下结果稳定，不替代其他体系的收敛检查。完整光谱步骤见分子红外与拉曼。