CP2K专题

CP2K机械嵌入与能量拆分

实际对照QM水、经典模型和QM/MM总能量,并逐个核对原子力。

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

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

1 用Ar例子理解机械嵌入#

先完成原子分区与拓扑。这里的 E_COUPL NONE关闭外部MM电荷对QM哈密顿量的静电嵌入,保留经典相互作用;它不是“QM和MM完全互不作用”。

Ar是中性MM粒子,与水之间仍有LJ作用。QM水的电子结构则与相同坐标、相同网格的孤立QM水相同。

2 完整可运行输入#

cat inputs/ar-qmmm/input.inp
cat inputs/ar-qmmm/topology.psf
bash run.sh

下面完整输入与前篇的PSF放在同一目录即可运行:

展开Ar机械嵌入完整input.inp
&GLOBAL
 PROJECT ar-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
   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 Ar
    CHARGE 0.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 Ar
     EPSILON [kcalmol] 0.190262450315
     SIGMA [angstrom] 3.277805000000
     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 Ar
     EPSILON [kcalmol] 0.000000000000
     SIGMA [angstrom] 2.202500000000
     RCUT [angstrom] 5.5
    &END LENNARD-JONES
    &LENNARD-JONES
     ATOMS Ar Ar
     EPSILON [kcalmol] 0.238000000000
     SIGMA [angstrom] 3.405000000000
     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 NONE
  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 Ar
   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
   Ar 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 Ar
   ELEMENT Ar
  &END KIND
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
&END FORCE_EVAL

本例非周期,MM静电采用 EWALD_TYPE NONE,QM使用MT孤立边界;Ar电荷为0。LJ截断5.5Å、截断能量平移、样条阈值1E−10Hartree。它们是明确的教学模型设置,不是通用水体系长程静电方案。

3 总能量读取哪一行#

grep 'ENERGY|' work/cases/water-qm/output.log
grep 'ENERGY|' work/cases/water-mm/output.log
grep 'ENERGY|' work/cases/ar-mm/output.log
grep 'ENERGY|' work/cases/ar-qmmm/output.log
条件 实际总能量/Hartree
水,纯QM −17.219538220205369
水,纯MM 0.000000000000000
水+Ar,纯MM −0.000204128432457
水+Ar,QM/MM −17.219742348637826

这里相互作用能为:

ΔE = E(QM/MM水+Ar) − E(QM水)
   = E(MM水+Ar) − E(MM水)
   = −0.005554617 eV

机械嵌入与电嵌入的实际相互作用能

水的MM能量为零,因为这个构型恰好取经典水模型的平衡键长和键角,不是所有水结构的经典内能都为零。QM水和MM水能量零点不同,不能把两者绝对总能量直接相减当作方法误差。

输出还会出现 QMMM Total Energy (QM + QMMM electronic + QMMM nuclear),这行没有包含完整MM能量。读最终结果应使用 ENERGY| Total FORCE_EVAL ( QMMM ) energy。

4 力也能拆分,不能只核对能量#

grep 'FORCES|' work/cases/ar-qmmm/output.log
module load ase/3.29.0
python -B work/analyze.py
cat work/analysis/ar-qmmm-forces.csv

Ar的Fx为约−0.00591203eV/Å,方向指向水的O,表示这个距离处的吸引。附件同时核对:

F(QM/MM) = F(QM水) + F(MM水+Ar) − F(MM水)

前三个原子按相同编号相加,第四个Ar的纯QM参考力为零。输出原始力单位是Hartree/Bohr,CSV已换成eV/Å。QM实空间网格会带来小的平移误差,因此需要比较同一网格、同一坐标的相互作用增量。

5 换成Na⁺,机械嵌入还保留什么#

附件中的 na-mechanical保持同样分区,把Ar换成有经典电荷的Na⁺。在这组实际输入中,机械嵌入仍保留经典水电荷与Na的相互作用,总相互作用能−0.083687080eV,与纯MM对照一致;QM水的偶极保持不变。

所以 E_COUPL NONE不能解释成“所有跨区静电都为零”。要让QM密度响应Na⁺,再看GAUSS电嵌入。

有共价键跨区时的下一步#

本篇水与外部粒子之间没有跨区共价键。丁烷LINK封端实操进一步展示两处C–C边界、IMOMM虚拟H、完整PDB与力场输入,并用两种微移步长核对投影力。机械边界与GAUSS边界电荷处理仍需按各自模型检查。

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

查看配图

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

100%查看原图 ↗