Gromacs+CP2K QM/MM模拟时报错

ZetaFunction · 发表于 Post on 2023-10-23 05:33:34

本帖最后由 ZetaFunction 于 2023-10-23 20:44 编辑

最近尝试用Gromacs连接CP2K跑蛋白的QM/MM模拟，用的是默认设置自动生成的inp文件，遇到CP2K的报错“CPASSERT failed”。
本人对CP2K的了解约等于零，希望各位大佬能够帮忙解答，以下是我的输入文件，以及自动生成的inp文件out文件和pdb文件，QM区有74个原子，总电荷-4，目前猜测有可能是电荷太大的关系？

_cp2k.inp (3.11 KB, 下载次数 Times of downloads: 11)

_cp2k.out (22.43 KB, 下载次数 Times of downloads: 3)

mdout.mdp (10.77 KB, 下载次数 Times of downloads: 8)

QMatoms.ndx (796.74 KB, 下载次数 Times of downloads: 4)

qmmm.log (17.47 KB, 下载次数 Times of downloads: 2)

qmmm.mdp (2.7 KB, 下载次数 Times of downloads: 10)

_cp2k.pdb.rar (437.54 KB, 下载次数 Times of downloads: 9)

gmx_mpi_d grompp -f qmmm.mdp -c md_20.gro -t md_20.cpt -p FgaPT2_IRP_LSN174.top -n QMatoms.ndx -o qmmm.tpr -maxwarn 5
mpirun -np 1 gmx_mpi_d mdrun -deffnm qmmm -dlb yes

万分感谢各位大佬！

更新：原来由amber转换而来的盒子信息好像有点问题，我重新生成了正方形的盒子，并且在楼下的建议下尝试先做QM/MM下的最小化，但还是相同的报错。

_cp2k.inp (2.6 KB, 下载次数 Times of downloads: 2)

_cp2k.out (11.21 KB, 下载次数 Times of downloads: 0)

_cp2k.pdb.rar (298.12 KB, 下载次数 Times of downloads: 0)

qmmm_min.mdp (1.22 KB, 下载次数 Times of downloads: 1)

qmmm.mdp (2.7 KB, 下载次数 Times of downloads: 1)

ZetaFunction · 发表于 Post on 2023-10-23 14:42:56

重新设定了QM区域为27个原子，总电荷为0，和MM区域也没有共价键相连，但还是报错。

Translating the system in order to center the QM fragment in the QM box.
*******************************************************************************
* ___ *
* / \ *
* [ABORT] *
* \___/ CPASSERT failed *
* | *
* O/| *
* /| | *
* / \ cp_ddapc_methods.F:650 *
*******************************************************************************
===== Routine Calling Stack =====
8 ewald_ddapc_pot
7 cp_ddapc_create
6 cp_ddapc_init
5 qs_env_update_s_mstruct
4 qs_energies_init_hamiltonians
3 qs_energies
2 qs_forces
1 CP2K

复制代码

看输出信息是不是因为MM区域的几何结构有问题？但是我检查结构文件，发现除了部分MM区域跨过了边界其他都很正常，输入文件是Gromacs自动生成的，请问要怎么处理？

啊不错的飞过海 · 发表于 Post on 2023-10-23 17:10:35

设的盒子是正交的吗？CP2K在这里应该是检查过盒子正交性，不正交的话跑不了

beyond · 发表于 Post on 2023-10-23 18:44:07

&QMMM
&CELL
   A 29.472 0.000 0.000
   B -9.824 27.786 0.000
   C -9.824 -13.893 24.063
   PERIODIC XYZ
&END CELL

有点看不懂这个QM的box，不是长方体或者立方体的，建议检查一下。
另外跑QM/MM建议用NVT的系综，步长2fs也过大了，1fs或者0.5fs
另外从cMD到QMMM，最好先做QMMM的minimizaiton和equilibration。。
如果对于CP2K了解不多的话，跑起来估计会有很多问题不知道如何解决的。。。

ZetaFunction · 发表于 Post on 2023-10-23 19:22:07

啊不错的飞过海发表于 2023-10-23 17:10
设的盒子是正交的吗？CP2K在这里应该是检查过盒子正交性，不正交的话跑不了

我设的是截角八面体的盒子，是必须要方形的盒子吗？

ZetaFunction · 发表于 Post on 2023-10-23 19:27:53

beyond 发表于 2023-10-23 18:44
&QMMM
&CELL
A 29.472 0.000 0.000

感谢解答，刚刚检查了一下盒子确实有问题，因为我的拓扑文件是从amber格式转换的，本来设置的是截角八面体，好像出了点问题没有读取这块信息，模拟的时候又变成方形了，我重新设置一下试试看。

需要把cMD最后一帧的结构最小化后重新重新加热？1fs甚至0.5fs的间隔对于QM/MM来说耗时好像有点太长了，减小间隔时间是否是为了便于SCF收敛的考量？如果SCF可以收敛是不是就没有必要了。

ZetaFunction · 发表于 Post on 2023-10-23 20:39:53

beyond 发表于 2023-10-23 18:44
&QMMM
&CELL
A 29.472 0.000 0.000

我重新设置了周期性盒子的参数，改成了方形盒子，但还是报相同的错误，尝试调低间隔时间也不行，先做QM/MM的最小化也不行。

ZetaFunction · 发表于 Post on 2023-10-23 20:46:09

beyond 发表于 2023-10-23 18:44
&QMMM
&CELL
A 29.472 0.000 0.000

重新设定盒子参数后产生的inp文件：

&GLOBAL
PRINT_LEVEL LOW
PROJECT GROMACS
RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
METHOD QMMM
&DFT
CHARGE 0
MULTIPLICITY 1
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
NGRIDS 5
CUTOFF 450
REL_CUTOFF 50
COMMENSURATE
&END MGRID
&SCF
SCF_GUESS RESTART
EPS_SCF 5.0E-8
MAX_SCF 20
&OT T
MINIMIZER DIIS
STEPSIZE 0.15
PRECONDITIONER FULL_ALL
&END OT
&OUTER_SCF T
MAX_SCF 20
EPS_SCF 5.0E-8
&END OUTER_SCF
&END SCF
&XC
DENSITY_CUTOFF 1.0E-12
GRADIENT_CUTOFF 1.0E-12
TAU_CUTOFF 1.0E-12
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&QS
METHOD GPW
EPS_DEFAULT 1.0E-10
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 4
&END QS
&END DFT
&QMMM
&CELL
A 17.442 0.000 0.000
B 0.000 17.442 0.000
C 10.070 10.070 14.241
PERIODIC XYZ
&END CELL
CENTER EVERY_STEP
CENTER_GRID TRUE
&WALLS
TYPE REFLECTIVE
&END WALLS
ECOUPL GAUSS
USE_GEEP_LIB 12
&PERIODIC
GMAX 1.0E+00
&MULTIPOLE ON
RCUT 1.0E+01
EWALD_PRECISION 1.0E-06
&END
&END PERIODIC
&QM_KIND H
MM_INDEX 6849 6850 6851 6852 6854 6855 6858 6860 6864 6866 6868 6870
&END QM_KIND
&QM_KIND C
MM_INDEX 6845 6847 6853 6856 6857 6861 6862 6863 6865 6867 6869
&END QM_KIND
&QM_KIND N
MM_INDEX 6848 6859
&END QM_KIND
&QM_KIND O
MM_INDEX 6844 6846
&END QM_KIND
&END QMMM
&MM
&FORCEFIELD
DO_NONBONDED FALSE
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE NONE
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
A 86.728 0.000 0.000
B 0.000 86.728 0.000
C 43.364 43.364 61.326
PERIODIC XYZ
&END CELL
&TOPOLOGY
COORD_FILE_NAME qmmm_em/_cp2k.pdb
COORD_FILE_FORMAT PDB
CHARGE_EXTENDED TRUE
CONNECTIVITY OFF
&GENERATE
&ISOLATED_ATOMS
LIST 1..45896
&END
&END GENERATE
&END TOPOLOGY
&KIND H
ELEMENT H
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND C
ELEMENT C
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND N
ELEMENT N
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND O
ELEMENT O
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND X
ELEMENT H
&END KIND
&END SUBSYS
&END FORCE_EVAL

复制代码

调低时间间隔还是同样的报错，CPASSERT failed，cp_ddapc_methods.F:650，尝试做最小化也报错

ZetaFunction · 发表于 Post on 2023-10-26 01:53:38

改用了cubic盒子，并且把原先错误命名为NA的钠离子重新命名为Na+，可以正常运行，问题解决了！但是计算的速度特别慢，QM区域才27个原子，用PBE居然每帧要6-7分钟。

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[GROMACS] Gromacs+CP2K QM/MM模拟时报错

浏览过的版块