鸢算科技 yuansuan Technology
用户门户
CP2K专题

CP2K电嵌入与QM极化实操

Na离子GAUSS嵌入的输入、同网格参照、QM偶极及能量差分力对照。

更新于 2026-10-28约 8 分钟阅读按任务与环境选择适用配置
阅读本篇目录

下载完整输入、实际结果与四张高清图。解压后在已分配CPU资源中执行 bash run.sh；输入在 inputs，新的结果在 work/cases，分析在 work/analysis，图片在 work/figures。

1 让外部电荷影响QM电子密度#

机械嵌入中，Na⁺与水的相互作用来自经典参数，QM水不被其电场极化。这里改用 E_COUPL GAUSS，把MM电荷产生的势加入QM计算。

结构仍为QM水＋距O4Å的MM Na⁺。Na没有加入量子区，不能把这个结果解释成“水和Na的全量子计算”。

2 输入要同时设置分区、网格和电荷#

&DFT
 CHARGE 0
 &MGRID
  COMMENSURATE .TRUE.
  CUTOFF 400
  REL_CUTOFF 50
 &END MGRID
 ...
&END DFT
&QMMM
 E_COUPL GAUSS
 CENTER NEVER
 USE_GEEP_LIB 6
 &CELL
  ABC 12 12 12
  PERIODIC NONE
 &END CELL
 &QM_KIND O
  MM_INDEX 1
 &END QM_KIND
 &QM_KIND H
  MM_INDEX 2 3
 &END QM_KIND
 &MM_KIND Na
  RADIUS [angstrom] 1.0
 &END MM_KIND
&END QMMM

这里省略号只用于说明关键设置，完整可运行文件如下：

展开Na电嵌入完整input.inp 
&GLOBAL
 PROJECT na-qmmm
 RUN_TYPE ENERGY_FORCE
 PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
 METHOD QMMM
 &DFT
  BASIS_SET_FILE_NAME BASIS_MOLOPT
  POTENTIAL_FILE_NAME GTH_POTENTIALS
  CHARGE 0
  &MGRID
   COMMENSURATE .TRUE.
   CUTOFF 400
   REL_CUTOFF 50
  &END MGRID
  &SCF
   EPS_SCF 1E-6
   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
    REFERENCE COM
   &END MOMENTS
  &END PRINT
 &END DFT
 &MM
  &FORCEFIELD
   &SPLINE
    EPS_SPLINE 1E-10
   &END SPLINE
   SHIFT_CUTOFF .TRUE.
   &BOND
    ATOMS O H
    K [kcalmol*angstrom^-2] 450.0
    KIND CHARMM
    R0 [angstrom] 0.9572
   &END BOND
   &BEND
    ATOMS H O H
    K [kcalmol*rad^-2] 55.0
    KIND CHARMM
    THETA0 [deg] 104.52
   &END BEND
   &CHARGE
    ATOM O
    CHARGE -0.834
   &END CHARGE
   &CHARGE
    ATOM H
    CHARGE 0.417
   &END CHARGE
   &CHARGE
    ATOM Na
    CHARGE 1.0
   &END CHARGE
   &NONBONDED
    &LENNARD-JONES
     ATOMS O O
     EPSILON [kcalmol] 0.152100000000
     SIGMA [angstrom] 3.150610000000
     RCUT [angstrom] 5.5
    &END LENNARD-JONES
    &LENNARD-JONES
     ATOMS O H
     EPSILON [kcalmol] 0.000000000000
     SIGMA [angstrom] 2.075305000000
     RCUT [angstrom] 5.5
    &END LENNARD-JONES
    &LENNARD-JONES
     ATOMS O Na
     EPSILON [kcalmol] 0.084459990528
     SIGMA [angstrom] 2.790268126865
     RCUT [angstrom] 5.5
    &END LENNARD-JONES
    &LENNARD-JONES
     ATOMS H H
     EPSILON [kcalmol] 0.000000000000
     SIGMA [angstrom] 1.000000000000
     RCUT [angstrom] 5.5
    &END LENNARD-JONES
    &LENNARD-JONES
     ATOMS H Na
     EPSILON [kcalmol] 0.000000000000
     SIGMA [angstrom] 1.714963126865
     RCUT [angstrom] 5.5
    &END LENNARD-JONES
    &LENNARD-JONES
     ATOMS Na Na
     EPSILON [kcalmol] 0.046900000000
     SIGMA [angstrom] 2.429926253730
     RCUT [angstrom] 5.5
    &END LENNARD-JONES
   &END NONBONDED
  &END FORCEFIELD
  &POISSON
   &EWALD
    EWALD_TYPE NONE
   &END EWALD
  &END POISSON
 &END MM
 &QMMM
  E_COUPL GAUSS
  CENTER NEVER
  USE_GEEP_LIB 6
  &CELL
   ABC 12 12 12
   PERIODIC NONE
  &END CELL
  &QM_KIND O
   MM_INDEX 1
  &END QM_KIND
  &QM_KIND H
   MM_INDEX 2 3
  &END QM_KIND
  &MM_KIND Na
   RADIUS [angstrom] 1.0
  &END MM_KIND
 &END QMMM
 &SUBSYS
  &CELL
   ABC 12 12 12
   PERIODIC NONE
  &END CELL
  &COORD
   O 6.000000000000 6.000000000000 6.000000000000
   H 6.000000000000 6.756950327264 6.585882276618
   H 6.000000000000 5.243049672736 6.585882276618
   Na 10.000000000000 6.000000000000 6.000000000000
  &END COORD
  &TOPOLOGY
   CONN_FILE_FORMAT PSF
   CONN_FILE_NAME topology.psf
  &END TOPOLOGY
  &KIND O
   ELEMENT O
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q6
  &END KIND
  &KIND H
   ELEMENT H
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q1
  &END KIND
  &KIND Na
   ELEMENT Na
  &END KIND
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
&END FORCE_EVAL

GAUSS要求DFT多重网格设为 COMMENSURATE。这个设置会影响子网格，因此比较能量时，孤立水参照也必须使用同样网格；附件对应目录是 water-qm-grid，不能拿另一个网格的 water-qm随意替代。

RADIUS 1Å是这个MM粒子的高斯电荷模型参数，不是水合离子半径或QM截断半径；真实模型要检查参数和结果对它的敏感性。QM水用 CHARGE 0，整个模型因MM Na⁺带+1电荷。

3 按顺序运行与取数#

cat inputs/na-qmmm/input.inp
cat inputs/na-qmmm/topology.psf
sbatch submit.slurm

作业结束后：

grep 'SCF run converged' work/cases/na-qmmm/output.log
grep 'ENERGY|' work/cases/na-qmmm/output.log
grep -A1 'Dipole moment' work/cases/na-qmmm/output.log
cat work/analysis/energies.csv
条件 实际总能量/Hartree
同网格孤立QM水 −17.219478567653560
Na机械嵌入 −17.222554011388251
Na GAUSS电嵌入 −17.224358064081070

相同几何下，机械嵌入相互作用能−0.083687080eV，GAUSS为−0.132777850eV。这是当前固定几何、基组、网格和电荷模型的差异，不是结合自由能，不能仅凭能量更低就认定所有模型参数都合适。

GAUSS输出有很大的电子与核静电分项，两者会相消；不能只取负的电子项作为结合能，也不能把经典跨区固定电荷静电再加一遍。使用完整 ENERGY| ... ( QMMM )总量，并配合同模型参照。

4 看QM水是否真的发生极化#

QM水子系统实际偶极分量与变化

条件 QM水μx/D QM水μz/D 偶极大小/D
孤立QM，同网格 −0.00000147 1.832321 1.832321
Na机械嵌入 −0.00000147 1.832321 1.832321
Na GAUSS −0.19848250 1.856061 1.866644

Na位于水的+x方向，QM电子密度发生响应，偶极x分量明显改变。这些偶极来自中性QM水子系统。不要把它们写成总带电组合的原点无关偶极矩。

5 用微小位移检查力#

附件还做了Na沿x方向±0.01Å的两个独立单点：na-xminus与 na-xplus。

grep 'ENERGY|' work/cases/na-xminus/output.log
grep 'ENERGY|' work/cases/na-xplus/output.log
grep 'FORCES|' work/cases/na-qmmm/output.log
module load ase/3.29.0
python -B work/analyze.py
python -B work/plot_results.py

Fx ≈ −[E(x+0.01Å)−E(x−0.01Å)] / 0.02Å

实际微移能量与CP2K输出力的对照

能量差分得到−0.099837941eV/Å，程序输出−0.099836445eV/Å，相差约1.50×10⁻⁶eV/Å。这比只看到程序正常退出更能帮助判断输入、单位和力的方向是否一致。

6 扩展到自己的QM/MM体系#

先保证MM拓扑和力场在纯MM模型中可用，再检查QM区的电荷、自旋、边界和基组。明确是否跨共价键、边界如何处理、是否周期及长程静电如何求解，再做QM/MM单点与力检查，最后开展优化或动力学。

当前四原子模型没有跨区共价键，也没有液体采样；不能据此直接给出溶剂化自由能、反应势垒或大体系平衡性质。参考：CP2K QM/MM输入说明。

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

查看配图

放大后可滚动查看细节，按 Esc 关闭。

100%查看原图 ↗