计算化学公社

标题: 求助：用gromacs+CP2K联用对蛋白-配体小分子进行QMMM时能量最小化过程无法收敛 [打印本页]

作者
Author: 李放 时间: 2025-11-24 21:27
标题: 求助：用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分析会更好吗？

作者
Author: beyond 时间: 2025-11-24 22:31
跑QM/MM之前，最好先跑classical MD优化一下

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

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

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

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

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

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

作者
Author: 李放 时间: 2025-11-25 16:29

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

作者
Author: 李放 时间: 2025-11-25 16:32

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

好的，谢谢sob老师

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

作者
Author: beyond 时间: 2025-11-25 18:12

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

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

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

作者
Author: sobereva 时间: 2025-11-26 04:52

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

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

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

作者
Author: 李放 时间: 2025-11-26 10:41

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区域，但是无法收敛

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

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

作者
Author: 李放 时间: 2025-12-4 11:25

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

好的，谢谢sob老师！

作者
Author: 李放 时间: 2025-12-4 11:39
本帖最后由李放于 2025-12-4 18:45 编辑

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

好的，谢谢

欢迎光临计算化学公社 (http://bbs.keinsci.com/)