本文首发于 https://zhuanlan.zhihu.com/p/22599245408
QM/MM/MD模拟是一种研究体系反应的常用方法，于2013年获得诺奖。CP2K是最近兴起的第一性原理软件，可负责QM部分；Gromacs是常用的MD软件，且提供了CP2K接口，可负责MM和MD部分。二者的结合是一种有效的QM/MM toolchain。实际上CP2K自己也可算QM/MM/MD，但需要用到Amber的力场格式，比较复杂，还是更习惯用Gromacs。本文总结了我通过Gromacs+CP2K做QM/MM/MD的步骤和遇到的问题。
QM/MM的介绍如下。
GROMACS/CP2K interface QM/MM parameterisation
步骤
- 交叉编译Gromacs和CP2K。
- 使用Gromacs建立常规MD体系。
- 使用Gromacs将体系平衡一段时间。
- 使用Gromacs建立体系的QM部分的索引，命名为“QMatoms”。
- 在产生相的mdp文件中，加入一段QM/MM计算的关键词。
- 使用Gromacs结合CP2K生成产生相的tpr文件，并跑QM/MM/MD。
1–3与经典MD一致。
第4步中，所选QM区域原子范围，一般选择催化残基的侧链+底物。例如，PET降解酶的催化残基是SHD三联体，则将QM/MM的边界划在SHD各自与其主链的C—C键上（CA—CB键）。但这样做会使QM和MM区各自的电荷不为整数。一种解决办法是重新计算RESP电荷并限制QM区的电荷为整数。另一种做法是，注意到每个残基的电荷为整数，可将边界划在残基与前后残基相连的酯键上。有条件的话更推荐前者。
第5步可通过我的VSCode扩展Gromacs Helper来生成。输入qmmm选择模板，对电荷和多重度稍作修改即可。
第6步需要提前建立好与产生相tpr文件同名的目录（例如md/），否则无法正常跑。这是 Gromacs 2024 的 bug，2025 版已经修复了。
此外，可能遇到的错误是，因为底物使用了GAFF力场，原子类型为小写，不能被CP2K正常识别，报错信息为：Atoms 3959 does not have atomic number needed for QMMM. Check atomtypes section in your topology or forcefield.。
这里将底物拓扑用Sobtop重新生成了一遍，使用AMBER原子类型。接下来即可正常运行QM/MM/MD。
目前该方法默认只支持PBE和BLYP两种理论方法，以及DZVP一种基组。
注意：GFN1-xTB半经验方法精度较低，未必适用于QM/MM计算。谨慎使用。
如果想要使用GFN1-xTB方法的话，需要自行提供模板文件。替换掉默认生成的 _cp2k.inp（Gromacs 2025 生成的是 xxx_cp2k.inp）中的&DFT字段，使用GFN1-xTB的设置，例如我这里的：
&DFT
# WFN_RESTART_FILE_NAME HID-RESTART.wfn
CHARGE -2 #Net charge
MULTIPLICITY 1 #Spin multiplicity
&QS
EPS_DEFAULT 1.0E-10 #Set all EPS_xxx to values such that the energy will be correct up to this value
EXTRAPOLATION ASPC #Extrapolation for wavefunction during e.g. MD. ASPC is default, PS can also be used
EXTRAPOLATION_ORDER 3 #Order for PS or ASPC extrapolation. 3 is default
METHOD xTB
&xTB
DO_EWALD T
CHECK_ATOMIC_CHARGES F #xTB calculation often crashes without setting this to false
&END xTB
&END QS
&MGRID
NGRIDS 5
CUTOFF 450
REL_CUTOFF 50
COMMENSURATE
&END MGRID
&POISSON
PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
PSOLVER PERIODIC #The way to solve Poisson equation
&END POISSON
&SCF
MAX_SCF 25 #Maximum number of steps of inner SCF
EPS_SCF 1.0E-05 #Convergence threshold of density matrix of inner SCF
# SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
# IGNORE_CONVERGENCE_FAILURE #Continue calculation even if SCF not converged, works for version >= 2024.1
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS #CG is worth to consider in difficult cases
LINESEARCH 2PNT #1D line search algorithm for CG. 2PNT is default. 3PNT is more expensive but may be better. GOLD is best but very expensive
ALGORITHM STRICT #Algorithm of OT. Can be STRICT (default) or IRAC
&END OT
&OUTER_SCF
MAX_SCF 20 #Maximum number of steps of outer SCF
EPS_SCF 1.0E-05 #Convergence threshold of outer SCF
&END OUTER_SCF
&PRINT
&RESTART OFF #Do not generate wfn file to suppress meaningless I/O cost
&END RESTART
&END PRINT
&END SCF
&END DFT
并在&FORCE_EVAL-&QMMM-&PERIODIC 字段里加入如下内容：
&POISSON
&EWALD
EWALD_TYPE SPME
GMAX 25 25 25
&END EWALD
&END POISSON
此外，将 ECOUPL 选项修改为 POINT_CHARGE （xTB 没法用 GAUSS）。
将文件保存在当前目录，命名为temp.inp，并将mdp文件的默认方法改为INPUT，在grompp步骤加上选项：-qmi temp.inp即可。
注意：CP2K 会生成一个新的结构文件，名字例如 md_cp2k.pdb ，其中原子标号为从 0 开始。因此涉及到 QM 区域原子索引时（例如 COLVARS 的 atomNumbers），需要根据输出确认索引是否正确。例如，COLVARS 在输出中会给出 CV 的定义信息，可据此确认索引原子的准确性。