CP2K 静态外场与分子偶极

给分子施加静态电场，可以观察电子能量和偶极的变化。本篇固定优化后的H₂O原子位置，用同一PBE模型做五个小外场条件，先把外势、电场方向和单位说清楚。

下载完整输入、真实日志、CSV与绘图脚本，解压进入run.sh目录。采用CP2K2026.2 P100双精度构建（GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3），一个MPI、四个OpenMP线程、一张P100；分析环境为ASE3.29.0。

1 提交五个条件#

sbatch submit.slurm

脚本在同一张GPU依次执行零外势斜率、±0.0005、±0.001五个条件，再提取与画图。核时P100模板；套餐用户用yskit生成所属分区脚本，保留已分配资源中的 bash run.sh。不要在登录节点计算。inputs可修改，reference保存本文结果，work保存你的新结果；已有work时新解压。

2 用解析函数定义外势#

采用PBE、DZVP-MOLOPT-SR-GTH、O GTH-PBE-q6、H GTH-PBE-q1、12 Å非周期MT盒、600/50 Ry网格与EPS_SCF1E−8。核坐标来自分子振动教程的实际优化结果，本篇不在外场中继续放松原子。

完整input.inp：water-field-pluslarge#

&GLOBAL
 PROJECT response
 RUN_TYPE ENERGY_FORCE
 PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
 METHOD QS
 &DFT
  &EXTERNAL_POTENTIAL
   FUNCTION A*(Z-Z0)
   PARAMETERS A Z0
   VALUES [hartree*bohr^-1] 0.001 [angstrom] 6.0
  &END EXTERNAL_POTENTIAL
  &PRINT
   &MOMENTS ON
    PERIODIC .FALSE.
    REFERENCE COM
   &END MOMENTS
  &END PRINT
  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
 &END DFT
 &SUBSYS
  &CELL
   ABC 12 12 12
   PERIODIC NONE
  &END CELL
  &COORD
  O         5.9999999937        6.0000000000        5.9993806921
  H         5.9999999971        6.7653454490        6.5954871263
  H         5.9999998457        5.2346546213        6.5954872165
  &END COORD
  &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

FUNCTION A*(Z-Z0)定义沿z线性变化的静态电势，坐标变量由CP2K处理。参数A使用Hartree/Bohr，Z0使用Å并显式标注单位。Z0=6 Å只定义势零点，不是把电场限制在z=6 Å附近。

物理电场是电势的负梯度，因此本篇 E_z=−A。A=+0.001时，物理电场沿负z；把参数名叫A可以避免把电势斜率直接误当成同方向电场。1个电场原子单位约51.422 V/Å，本例幅度0.0005/0.001分别约0.02571/0.05142 V/Å。

3 取偶极与能量#

DFT/PRINT/MOMENTS采用非周期传统偶极算符，并以质心为参考点。中性分子的总偶极由电子与离子共同给出；不要只取电子部分与总能量差分比较。

grep -E 'ENERGY\| Total FORCE_EVAL' work/cases/water-field-*/output.log
grep -A 1 'Dipole moment' work/cases/water-field-*/output.log
cat work/analysis/fields.csv

原分子有约1.800482 Debye的z向永久偶极。沿其方向的物理电场使偶极增加、能量降低；正负场不必产生对称的总能量，因为永久偶极带来线性项。讨论极化率时还需要二阶能量曲率或偶极的一阶变化，不能用某一个场下的总偶极除以电场。

本例仍是孤立分子边界。这个解析线性势不能直接替代周期晶体的Berry相位有限场设置；周期体系和带电体系需要相应的边界、势参考与电荷处理。

提取和重画#

module load ase/3.29.0
python analyze.py work
python plot_results.py work
cat work/analysis/fields.csv
cat work/analysis/RESULTS.json
ls work/figures

CSV依次列出外势斜率A、物理电场、总能量Hartree和三个偶极分量Debye；JSON保存有限差分与原生线性响应结果。只查看本文图，把work换成reference。

继续看有限外场与极化率。函数和参数单位定义见CP2K2026.2 EXTERNAL_POTENTIAL手册。