CP2K机械嵌入与能量拆分
实际对照QM水、经典模型和QM/MM总能量,并逐个核对原子力。
阅读本篇目录
下载完整输入、实际结果与四张高清图。解压后在已分配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边界电荷处理仍需按各自模型检查。