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

CP2K跨共价键QM/MM与LINK封端

丁烷中心量子区、两端经典区的实际IMOMM输入，检查原子编号、封端氢和投影力。

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

下载完整输入、六个实际条件和三张高清图。解压后在已分配CPU资源内执行 bash run.sh，新结果在 work/cases、表格在 work/analysis、图片在 work/figures。

水与外部Ar/Na例子没有共价键跨过QM/MM边界，因此无需封端。本篇把丁烷中间片段交给QM、两端交给MM，两处C–C键跨区，需要明确处理边界。

实际丁烷原子编号、QM/MM分区与封端氢

完整模型有14个真实原子。第5–10号是中间C₂H₄片段；两端甲基属于MM。QM模型在两处边界各补一个虚拟H，变成8原子的C₂H₆，具有14个价电子。

虚拟封端H由程序生成，不要再把它们作为真实原子加进PDB。不然原子总数、电子数与原本的力场拓扑都会改变。

2 获取输入并提交#

程序为CP2K2026.2 CPU版，Intel oneAPI2024.2、Intel MPI2021.13。每个小条件4MPI×1线程，完整脚本最多同时两个条件。

module load cp2k/2026.2-8173m
ls inputs/link-base
cat inputs/link-base/butane.pdb
cat inputs/link-base/input.inp
sbatch submit.slurm

CPU模板申请56核整节点；小例配置用于对照，不作为实际大体系并行效率建议。已有计算资源时也可直接 bash run.sh。

单独运行基准输入：

mkdir my-link
cp inputs/link-base/input.inp inputs/link-base/butane.pdb inputs/link-base/butane.pot my-link/
cd my-link
mpirun -np 4 cp2k.psmp -i input.inp -o output.log
grep 'SCF run converged' output.log
grep 'ENERGY|' output.log

3 真实结构与经典参数#

以下为实际PDB，原子名/类型和B列电荷来自明确的丁烷教学模型：

ATOM      1 CTL3 BUT A   1       5.324   4.335   5.642  0.00 -0.27      BUT  C
ATOM      2 HAL3 BUT A   1       6.352   4.066   5.719  0.00  0.09      BUT  H
ATOM      3 HAL3 BUT A   1       4.880   4.006   4.772  0.00  0.09      BUT  H
ATOM      4 HAL3 BUT A   1       4.757   3.850   6.334  0.00  0.09      BUT  H
ATOM      5 CTL2 BUT A   1       5.166   5.869   5.818  0.00 -0.18      BUT  C
ATOM      6 HAL2 BUT A   1       4.272   6.135   6.421  0.00  0.09      BUT  H
ATOM      7 HAL2 BUT A   1       5.146   6.354   4.894  0.00  0.09      BUT  H
ATOM      8 CTL2 BUT A   1       6.323   6.384   6.684  0.00 -0.18      BUT  C
ATOM      9 HAL2 BUT A   1       6.983   5.652   6.996  0.00  0.09      BUT  H
ATOM     10 HAL2 BUT A   1       5.804   6.894   7.503  0.00  0.09      BUT  H
ATOM     11 CTL3 BUT A   1       7.119   7.438   5.943  0.00 -0.27      BUT  C
ATOM     12 HAL3 BUT A   1       8.149   7.278   6.094  0.00  0.09      BUT  H
ATOM     13 HAL3 BUT A   1       6.824   8.409   6.263  0.00  0.09      BUT  H
ATOM     14 HAL3 BUT A   1       6.896   7.324   4.921  0.00  0.09      BUT  H

CHARGE_BETA使CP2K从PDB的B因子列读取经典电荷。通常的实验PDB中该列是温度因子，不能未经处理就把任意蛋白PDB加上这个设置。本例的原子类型与 butane.pot对应，坐标整体移到盒内，内部结构保持。

文件 butane.pot给出键、角、Urey-Bradley、二面角与LJ参数，不能只有几个电荷就当成完整有机力场。结构与参数改编自CP2K固定v2026.2的公开丁烷例子，来源/哈希和许可在附件 UPSTREAM-SOURCES.json、LICENSE-GPL-2.txt。本例改变了计算方法、网格、边界与耦合设置，结果应与本篇数据比较。

&QMMM
 E_COUPL NONE
 CENTER NEVER
 ...
 &LINK
  ALPHA_IMOMM 1.5
  FIST_SCALE_FACTOR 1.0
  LINK_TYPE IMOMM
  MM_INDEX 1
  QM_INDEX 5
  RADIUS 0.80
 &END LINK
 &LINK
  ALPHA_IMOMM 1.5
  FIST_SCALE_FACTOR 1.0
  LINK_TYPE IMOMM
  MM_INDEX 11
  QM_INDEX 8
  RADIUS 0.80
 &END LINK
 &QM_KIND C
  MM_INDEX 5 8
 &END QM_KIND
 &QM_KIND H
  MM_INDEX 6 7 9 10
 &END QM_KIND
&END QMMM

本例的编号按完整PDB原子表填写，并从输出确认边界为QM5–MM1、QM8–MM11；不要在分区后按QM子集自行重编成1、2，再误填进原模型。

FIST_SCALE_FACTOR 1.0保留边界MM原子在经典部分的电荷。这里采用机械嵌入，QM密度不被外部电荷极化；若改成GAUSS，还要专门考虑边界电荷分布与过极化，不能只换一个关键词就认定模型完成。

完整输入如下，前面的省略号不能直接拿来计算：

展开完整可运行input.inp 
&GLOBAL
  PRINT_LEVEL MEDIUM
  PROJECT link-base
  RUN_TYPE ENERGY_FORCE
&END GLOBAL


&FORCE_EVAL
  METHOD QMMM
  &DFT
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &MGRID
      COMMENSURATE
      CUTOFF 400
      REL_CUTOFF 50
    &END MGRID
    &QS
    &END QS
    &SCF
      SCF_GUESS ATOMIC
      EPS_SCF 1E-6
      MAX_SCF 100
      &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
  &MM
    &FORCEFIELD
      PARMTYPE CHM
      PARM_FILE_NAME butane.pot
    &END FORCEFIELD
    &POISSON
      &EWALD
       EWALD_TYPE NONE
      &END EWALD
    &END POISSON
  &END MM
  &QMMM
    E_COUPL NONE
    USE_GEEP_LIB 6
    CENTER NEVER
    &CELL
      ABC 12 12 12
      PERIODIC NONE
    &END CELL
    &INTERPOLATOR
      EPS_R 1.0e-14
      EPS_X 1.0e-14
      MAXITER 100
    &END INTERPOLATOR
    &LINK
      ALPHA_IMOMM 1.5
      FIST_SCALE_FACTOR 1.0
      LINK_TYPE IMOMM
      MM_INDEX 1
      QM_INDEX 5
      RADIUS 0.80
    &END LINK
    &LINK
      ALPHA_IMOMM 1.5
      FIST_SCALE_FACTOR 1.0
      LINK_TYPE IMOMM
      MM_INDEX 11
      QM_INDEX 8
      RADIUS 0.80
    &END LINK
    &MM_KIND CTL3
      RADIUS 0.78
    &END MM_KIND
    &MM_KIND HAL3
      RADIUS 0.78
    &END MM_KIND
    &QM_KIND C
      MM_INDEX 5 8
    &END QM_KIND
    &QM_KIND H
      MM_INDEX 6 7 9 10
    &END QM_KIND
    #
    # QM_KINDS
    #
    #
    # MM_KINDS
    #
    #
    # LINKS
    #
    #
  &END QMMM
  &SUBSYS
    &CELL
      ABC 12 12 12
      PERIODIC NONE
    &END CELL
    &KIND H
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q1
    &END KIND
    &KIND C
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q4
    &END KIND
    &TOPOLOGY
      CHARGE_BETA
      COORDINATE pdb
      COORD_FILE_NAME butane.pdb
    &END TOPOLOGY
  &END SUBSYS
  &PRINT
   &FORCES ON
   &END FORCES
  &END PRINT
&END FORCE_EVAL

实际采用PBE、DZVP-MOLOPT-SR-GTH、GTH-PBE、400/50Ry网格、非周期12Å盒、最终 EPS_SCF 1E-6。这是一组固定几何单点，不是优化后的丁烷构型。

5 封端氢在哪里，alpha控制什么#

IMOMM中，封端位置满足：

rH = rQM + (rMM − rQM) / alpha

因此同一真实C–C边界上，较小alpha给出较长QM–H距离。初值应参考力场平衡QM–MM键长与QM–H键长之比；alpha不是电子精度参数，也不是为了把总能量压低而任意调节的旋钮。

两种alpha的封端距离和固定结构能量差

同一真实结构，把alpha从1.5改到1.4，当前模型的总能量改变−0.209207465eV。这说明量子封端几何已经改变，不能把能量更低直接解释成参数更正确。

grep -A8 'QM/MM LINKS' work/cases/link-base/output.log
grep 'H_LINK' work/cases/link-base/output.log
grep 'Number of electrons' work/cases/link-base/output.log
grep 'ENERGY|' work/cases/link-base/output.log
grep 'FORCES|' work/cases/link-base/output.log

应看到两条IMOMM LINK、两个H_LINK和14个QM价电子。最终输出有14个真实原子的力；虚拟H的力会按LINK关系投影回真实边界原子，不会给你变成16个可以独立移动的真实原子。

7 微移真实边界原子，检查投影力#

附件将真实MM原子1沿x移动±0.01Å、±0.005Å，分别重新做电子单点，虚拟H随边界自动移动。

grep 'ENERGY|' work/cases/link-minus/output.log
grep 'ENERGY|' work/cases/link-plus/output.log
grep 'ENERGY|' work/cases/link-halfminus/output.log
grep 'ENERGY|' work/cases/link-halfplus/output.log
module load ase/3.29.0
python -B work/analyze.py
python -B work/plot_results.py

Fx(h) = −[E(x+h)−E(x−h)]/(2h)
Fx(extrapolated) = [4Fx(h/2)−Fx(h)]/3
计算方式 MM原子1的Fx/eV·Å⁻¹
h=0.010Å中心差分 0.258474312
h=0.005Å中心差分 0.259200899
两步长Richardson外推 0.259443095
CP2K输出的投影后力 0.259448803

实际边界原子的能量差分力与投影力

位移减半后误差从约9.74×10⁻⁴降到2.48×10⁻⁴eV/Å；外推后约5.71×10⁻⁶eV/Å。有限位移本身会带来误差，因此第一组不完全相等不代表LINK力错误，应该看步长收敛。

8 换成自己的QM/MM边界#

先明确QM区要保留的化学键和电荷，再保证经典拓扑、参数和编号一致。尽量选择合理的断键位置，检查封端键长、近边界电荷与局域结构，并验证能量与力。不同LINK类型、边界电荷重分布和电嵌入需要各自的参数检查，不能把本篇机械封端例子视为完整蛋白或反应自由能方案。

参考：CP2K LINK输入、QM/MM设置。

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

查看配图

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

100%查看原图 ↗