CP2K专题

CP2K SCCS腔体、盒长与介电函数

三种真实腔体模型、8/12A周期盒对照和原生介电函数CUBE可视化。

更新于 2026-10-29约 7 分钟阅读按任务与环境选择适用配置

阅读本篇目录

使用连续介质时，不仅要选择体相介电常数，还要知道溶质与介质的边界怎样定义。本篇沿用固定水分子的SCCS流程，比较盒长和三种腔体设置，再用原生介电函数网格查看它们的空间形状。

下载完整输入、真实日志、介电函数CUBE和绘图脚本，解压进入run.sh目录。采用CP2K2026.2 P100双精度构建：GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3；分析ASE3.29.0。

1 提交这一组#

sbatch submit.slurm

模板申请一张P100、一个MPI和一个CPU核，固定 OMP_NUM_THREADS=1，依次运行八个条件，再提取和画图。这里的资源设置与实际脚本一致；调整进程、线程或边界后，应重新核对结果。

这是核时P100模板，套餐用户用yskit生成所属分区脚本，保留已分配资源中的 bash run.sh。不要在登录节点计算。inputs可修改，reference为本文结果，work为你的新结果；已有work时请新解压。

2 每个盒长都用自己的参照#

12 Å条件保留分子内部几何，把所有原子一起平移2 Å，并改变三个盒长；电子参数和介电模型保持一致。分别计算12 Å中的ε=1与78.36，再求同盒差值。

盒长/Å	ε=1总能量/Hartree	ε=78.36总能量/Hartree	同盒电子能量差/eV
8	-17.219894911298	-17.226580424137	-0.181922
12	-17.219636006400	-17.226546586467	-0.188046

8与12 Å的电子能量差相差约6.12 meV。这两个点反映当前周期镜像与离散网格等数值影响，不足以外推无限稀释结果；需要定量值时继续检查盒长、网格和基组。

实际盒长对照和不同腔体模型的电子能量差

3 Fattebert–Gygi与溶剂感知模型#

Andreussi使用RHO_MIN/RHO_MAX。本例Fattebert–Gygi使用β=1.3、RHO_ZERO=0.0004，完整输入如下：

完整input.inp：water-sccs-fattebert#

&GLOBAL
 PROJECT sccs
 RUN_TYPE ENERGY
 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-8
  &END QS
  &MGRID
   CUTOFF 240
   REL_CUTOFF 40
   NGRIDS 4
   &RS_GRID
    DISTRIBUTION_TYPE DISTRIBUTED
   &END RS_GRID
  &END MGRID
  &POISSON
   PERIODIC XYZ
   POISSON_SOLVER PERIODIC
  &END POISSON
  &SCF
   EPS_SCF 1.0E-6
   MAX_SCF 300
   SCF_GUESS ATOMIC
   &DIAGONALIZATION
    ALGORITHM STANDARD
   &END DIAGONALIZATION
   &MIXING
    METHOD DIRECT_P_MIXING
    ALPHA 0.2
   &END MIXING
  &END SCF
  &XC
   &XC_FUNCTIONAL PBE
   &END XC_FUNCTIONAL
  &END XC
  &SCCS ON
   METHOD FATTEBERT-GYGI
   RELATIVE_PERMITTIVITY 78.36
   EPS_SCF 0.3
   EPS_SCCS 1.0E-6
   MAX_ITER 500
   MIXING 0.6
   DERIVATIVE_METHOD CD5
   ALPHA 0
   BETA 0
   GAMMA 0
   &FATTEBERT-GYGI
    BETA 1.3
    RHO_ZERO 0.0004
   &END FATTEBERT-GYGI
  &END SCCS
  &PRINT
   &SCCS ON
    &EACH
     QS_SCF 1
    &END EACH
    &DIELECTRIC_FUNCTION ON
     FILENAME =dielectric.cube
     APPEND FALSE
     STRIDE 1 1 1
     &EACH
      QS_SCF 1
     &END EACH
    &END DIELECTRIC_FUNCTION
    &DENSITY_GRADIENT OFF
    &END DENSITY_GRADIENT
    &POLARISATION_POTENTIAL OFF
    &END POLARISATION_POTENTIAL
   &END SCCS
  &END PRINT
 &END DFT
 &SUBSYS
  &CELL
   ABC 8.0 8.0 8.0
   PERIODIC XYZ
  &END CELL
  &COORD
O         4.9999870401        4.5606769781        4.5606769833
      H         5.0000000491        4.4536015938        5.5327585878
      H         5.0000000214        5.5327585972        4.4536017371
  &END COORD
  &KIND H
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q1
  &END KIND
  &KIND O
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q6
  &END KIND
 &END SUBSYS
&END FORCE_EVAL

溶剂感知算法选 METHOD SAA_ANDREUSSI，并设置：

&SAA_ANDREUSSI
 RHO_MIN 0.0001
 RHO_MAX 0.001
 R_SOLV 2.6
 ALPHA_ZETA 2.0
 DELTA_ZETA 0.5
 DELTA_ETA 0.02
 F0 0.65
&END SAA_ANDREUSSI

Fattebert–Gygi子块的BETA=1.3控制介面平滑宽度，与SCCS层BETA=0的色散体积参数是不同的量，不能混用。

R_SOLV与DELTA_ZETA的长度单位为Bohr；其他参数控制溶剂感知平滑。不要把R_SOLV当成输入原子的半径。完整文件为inputs/water-sccs-solvent-aware/input.inp，仍是同一几何、8 Å周期盒和ε=78.36。

腔体模型	8 Å总能量/Hartree	相对同盒ε=1/eV
Andreussi	-17.226580424137	-0.181922
Fattebert–Gygi	-17.236289848430	-0.446129
Solvent-aware	-17.226580427637	-0.181922

这组小分子中Andreussi与溶剂感知结果很接近；Fattebert–Gygi的腔体参数不同，能量变化也不同。不能按照哪一个总能量更低，就把该模型认定为更准确；模型与参数需要根据目标问题选择。

4 输出和读取真正的介电函数#

三种模型的8 Å、ε=78.36条件都会输出dielectric.cube：

&DFT
 &PRINT
  &SCCS ON
   &DIELECTRIC_FUNCTION ON
    FILENAME =dielectric.cube
    APPEND FALSE
    STRIDE 1 1 1
    &EACH
     QS_SCF 1
    &END EACH
   &END DIELECTRIC_FUNCTION
  &END SCCS
 &END PRINT
&END DFT

该段放在已有DFT/PRINT中，不要在同一输入末尾重复新建DFT。固定文件名和APPEND FALSE使循环中覆盖同名文件，完成后保留最后一次写出的场；子块EACH显式指定输出频率。STRIDE 1保留完整网格，便于明确坐标和周期盒长度。

这个CUBE的数值是无量纲介电常数ε(r)，不是电子密度，也不是电势。ASE读取时将CUBE坐标转换成Å，标量ε不乘密度单位换算因子。

module load ase/3.29.0
python - <<'PY'
from pathlib import Path
from ase.io.cube import read_cube
p = Path('work/cases/water-periodic-sccs-78p36/dielectric.cube')
with p.open() as f:
    c = read_cube(f)
print('grid:', c['data'].shape)
print('epsilon range:', c['data'].min(), c['data'].max())
print('cell in angstrom:', c['atoms'].cell.lengths())
PY

本例得到75×75×75网格，ε范围1至78.36，盒长接近8 Å；CUBE文本头的有限小数位会带来微小长度舍入。

三种模型实际dielectric.cube的完整周期盒截面与分子位置

切片取最接近原子平面的实际网格x≈5.0133 Å。原子附近是低介电区域，远处趋近体相值，白线给出相同介电常数等值线。图来自这三个真实CUBE文件，不是按照原子半径画出来的示意球。

提取与重画#

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

CSV用条件名、盒长、腔体模型、介电常数和各项能量标识每个结果，能量单位Hartree。JSON保留分项与CUBE网格信息；作图时转换为eV。只看本文参考，把work换成reference。

平滑函数和单位见Andreussi、Fattebert–Gygi与溶剂感知模型。

需要根据你的体系选择资源？联系讨论计算需求。实际可用软件与权限，以当前账户和运行环境为准。