求助：用gromacs+CP2K联用对蛋白-配体小分子进行QMMM时能量最小化过程无法收敛

李放

研究内容：酶（含金属CU）降解真菌毒素（蛋白-配体小分子）的机理研究，想通过gromacs+CP2K进行QMMM分析降解机理，有以下几个问题想请教大家
问题1、QM区域：T1 Cu其配位残基（HIS-427、HIS-489、CYS-484）、底物ZEN及作用残基ASP239，进行能量最小化时报错，是否和QM区域中各部分之间没有成键有关？因为当QM区域仅含底物ZEN时，可以正常运行，但会出现问题2

问题2、gromacs+CP2K联用进行QMMM分析时能量最小化一直震荡，无法收敛，CPU中运行2000步大概花了两三天（命令：文件生成gmx grompp -f em-qmmm.mdp -p topol.top -c input.gro -o em.tpr -n index.ndx；跑能量最小化OMP_NUM_THREADS=48 mpirun -np 4 gmx mdrun -v -s em.tpr -deffnm em&，服务器有两个物理 CPU，每个 CPU 96 个核心，只有该任务运行时也无法满核调动），增加步数是否可行或者有意义


以下是em-qmmm.mdp文件
; em-qmmm.mdp - used as input into grompp to generate stilbene-sol-opt.tpr
integrator    = steep ; Algorithm (steep = steepest descent minimization)
emtol        = 100.0  ; Stop minimization when the maximum force < 100.0 kJ/mol/nm
emstep          = 0.02  ; Energy step size
nsteps        = 2000   ; Maximum number of (minimization) steps to perform

; Set output frequency to each step
nstxout                  = 1 ; Coordinates to trr
nstlog                   = 1 ; Energies to em.log
nstcalcenergy            = 1 ; Energies to ener.edr
nstenergy                = 1 ; Energies to ener.edr

; Set cut-offs
rlist                    = 1.2
coulombtype              = PME
coulomb-modifier         = Potential-shift-Verlet
rcoulomb-switch          = 1.0
rcoulomb                 = 1.2
vdwtype                  = Cut-off
vdw-modifier             = Force-switch
rvdw-switch              = 1.0
rvdw                     = 1.2

; CP2K QMMM parameters
qmmm-cp2k-active              = true   ; Activate QMMM MdModule
qmmm-cp2k-qmgroup             = QMatoms; Index group of QM atoms
qmmm-cp2k-qmmethod            = PBE    ; Method to use
qmmm-cp2k-qmcharge            = 0      ; Charge of QM system
qmmm-cp2k-qmmultiplicity      = 1      ; Multiplicity of QM system

问题3：针对该体系和反应机理研究，用其他软件如orca分析会更好吗？

beyond

跑QM/MM之前，最好先跑classical MD优化一下

Amber还是CHARMM力场呢？Cu+的参数呢？

Gromacs-CP2K interface默认的QM box 偏大，QM计算的cutoff也偏大，
如果想更快一点，可以考虑使用-qmi,  自己修改一下cp2k的输入文件，
特别是你的静电荷为零。。。

sobereva

如果做能量极小化的目的是跑MD不崩，没必要收敛得太充分
绝对不要盲目设大步数上限，都震荡了设大还有什么意义

gromacs目前版本不支持和ORCA联用

如果你的目的仅仅是做静态的研究，根本没必要用CP2K结合gromacs，Gaussian做ONIOM(QM:MM)在这方面极为成熟、应用十分广泛，北京科音中级量子化学培训班（http://www.keinsci.com/KBQC）里专门有一节详细讲ONIOM，包括算生物分子体系。

另外，很多情况也可以简化为簇模型计算，省得折腾，看
要善用簇模型，不要盲目用ONIOM算蛋白质-小分子相互作用问题
http://sobereva.com/597

李放

beyond 发表于 2025-11-24 22:31
跑QM/MM之前，最好先跑classical MD优化一下

Amber还是CHARMM力场呢？Cu+的参数呢？

好的谢谢，那我先跑一下classical MD，然后试一下修改cp2k.inp文件

我用的是CHARMM力场，以下是Cu2+的参数

;; GROMACS topology file for CU2P
;;


[ moleculetype ]
; name        nrexcl
CU2P             1

[ atoms ]
; nr        type        resnr        residu        atom        cgnr        charge        mass
     1       CU2P      4     CU2P   CU2P      1     2.000000    63.5460   ; qtot  2.000

李放

sobereva 发表于 2025-11-25 02:04
如果做能量极小化的目的是跑MD不崩，没必要收敛得太充分
绝对不要盲目设大步数上限，都震荡了设大还有什么 ...

好的，谢谢sob老师

如果要做动态研究的话，单纯用ORCA可以做吗？您觉得用哪种方法对于新手是最方便的呢？

beyond

李放 发表于 2025-11-25 16:29
好的谢谢，那我先跑一下classical MD，然后试一下修改cp2k.inp文件

我用的是CHARMM力场，以下是Cu2+的 ...

如果是Cu2+，multiplicity应该不是1吧

sobereva

李放 发表于 2025-11-25 16:32
好的，谢谢sob老师

如果要做动态研究的话，单纯用ORCA可以做吗？您觉得用哪种方法对于新手是最方便的 ...

能做QM/MM MD，但不适合新手，远没GROMACS联用CP2K简单

李放

beyond 发表于 2025-11-25 18:12
如果是Cu2+，multiplicity应该不是1吧

请问multiplicity是对应QM区域的吗？因为我的QM区域只纳入了底物小分子，所以设置的1。
按照我的研究目的的话，QM区域应该为：T1 Cu（也就是1个Cu2+）其配位残基（HIS-427、HIS-489、CYS-484）、底物ZEN及作用残基ASP239，但是这样的话直接就报错了，所以就想着先只将底物纳入QM区域，但是无法收敛

beyond

如果CU2+参与了反应，那肯定要把Cu2+放在QM区域的，暂时不放，优化也不一定会解决问题，后期你还是要把Cu2+放进去的。
一般先跑 分子动力学，优化整个体系，在选取一帧跑QMMM。
况且你也需要跑基于力场的分子动力学，以了解你的体系。。。

上来就跑QM/MM，报错很正常。。。可以读读一些QM/MM的文章，看看他们的protocol

李放

sobereva 发表于 2025-11-26 04:52
能做QM/MM MD，但不适合新手，远没GROMACS联用CP2K简单

好的，谢谢sob老师！

李放

beyond 发表于 2025-11-26 19:49
如果CU2+参与了反应，那肯定要把Cu2+放在QM区域的，暂时不放，优化也不一定会解决问题，后期你还是要把Cu2+ ...

好的，谢谢