CP2K 有限外场与极化率

极化率描述偶极对小外场的线性响应。本篇沿用静态外场与偶极的同一固定H₂O结构，用两种步长的中心差分，与CP2K原生线性响应互相核对。

下载完整输入、真实日志、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 零场加入原生POLAR#

五个条件中，零场额外计算解析电子极化率。完整输入如下：

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

&GLOBAL
 PROJECT response
 RUN_TYPE ENERGY_FORCE
 PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
 METHOD QS
 &DFT
  &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
 &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
&END FORCE_EVAL

LINRES EPS 1E-8控制线性响应求解，POLAR ON输出极化率张量；非周期分子使用 PERIODIC_DIPOLE_OPERATOR FALSE。有限场的另外四份输入使用相同电子参数与原子位置，只加入不同A的外势。

3 中心差分的定义与单位#

E(F)=E₀−μ₀F−αF²/2+…，其中F为物理电场。本包的输入参数A=−F，因此对A有：

• 永久偶极μ₀=[E(+h)−E(−h)]/(2h)。
• α_zz=−[E(+h)+E(−h)−2E(0)]/h²。
• α_zz=−[μ_z(+h)−μ_z(−h)]/(2h)。

这些公式中的能量用Hartree、h用电场原子单位；偶极先由Debye除以2.541746473转换到e·Bohr。返回的α是原子单位，不是Debye。1个极化率原子单位约0.1481847 ų。

4 用实际CSV重算#

module load ase/3.29.0
python - <<'PY'
import numpy as np
x = np.loadtxt('work/analysis/fields.csv', delimiter=',', skiprows=1)
r = {round(row[0], 7): row for row in x}
for h in (0.0005, 0.001):
    ep, em, e0 = r[h][2], r[-h][2], r[0.0][2]
    alpha_e = -(ep + em - 2*e0) / h**2
    alpha_mu = -(r[h][5] - r[-h][5]) / (2*h*2.541746473)
    print(h, 'energy:', alpha_e, 'dipole:', alpha_mu, 'a.u.')
PY

输出应接近下表；保留原始双精度能量，不要先把总能量只留五六位小数再做二阶差分。

两个步长的能量曲率与原生LINRES均得到约 6.66876 a.u.，对应α_zz约 0.98821 ų。零场原生对角分量α_xx、α_yy、α_zz分别约4.68566、8.18683、6.66876 a.u.，与分子方向有关，不能把单个zz分量写成各向同性平均极化率。

5 什么保持不变，什么还需要收敛#

这里固定核位置，得到的是电子响应，没有加入外场引起的原子位置变化或振动贡献。两种步长与两种求导方式一致，说明当前有限模型中的数值相互吻合；紧凑基组缺少专门弥散空间，仍需做基组、盒长和网格收敛才能报告材料/分子的定量极化率。不能把内部数值一致直接解释成与实验一致。

提取和重画#

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 POLAR手册。