计算化学公社

标题: Gromacs+CP2K QM/MM模拟时报错 [打印本页]
作者
Author: ZetaFunction    时间: 2023-10-23 05:33
标题: Gromacs+CP2K QM/MM模拟时报错
本帖最后由 ZetaFunction 于 2023-10-23 20:44 编辑

最近尝试用Gromacs连接CP2K跑蛋白的QM/MM模拟，用的是默认设置自动生成的inp文件，遇到CP2K的报错“CPASSERT failed”。
本人对CP2K的了解约等于零，希望各位大佬能够帮忙解答，以下是我的输入文件，以及自动生成的inp文件out文件和pdb文件，QM区有74个原子，总电荷-4，目前猜测有可能是电荷太大的关系？
(, 下载次数 Times of downloads: 11) (, 下载次数 Times of downloads: 3) (, 下载次数 Times of downloads: 8) (, 下载次数 Times of downloads: 4) (, 下载次数 Times of downloads: 2) (, 下载次数 Times of downloads: 10) (, 下载次数 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下的最小化，但还是相同的报错。
(, 下载次数 Times of downloads: 2) (, 下载次数 Times of downloads: 0) (, 下载次数 Times of downloads: 0) (, 下载次数 Times of downloads: 1) (, 下载次数 Times of downloads: 1)


作者
Author: ZetaFunction    时间: 2023-10-23 14:42
重新设定了QM区域为27个原子，总电荷为0，和MM区域也没有共价键相连，但还是报错。
  1.   Translating the system in order to center the QM fragment in the QM box.

  3. *******************************************************************************
  4. *   ___                                                                       *
  5. *  /   \                                                                      *
  6. * [ABORT]                                                                     *
  7. *  \___/                             CPASSERT failed                          *
  8. *    |                                                                        *
  9. *  O/|                                                                        *
  10. * /| |                                                                        *
  11. * / \                                                  cp_ddapc_methods.F:650 *
  12. *******************************************************************************


  15. ===== Routine Calling Stack =====

  17.             8 ewald_ddapc_pot
  18.             7 cp_ddapc_create
  19.             6 cp_ddapc_init
  20.             5 qs_env_update_s_mstruct
  21.             4 qs_energies_init_hamiltonians
  22.             3 qs_energies
  23.             2 qs_forces
  24.             1 CP2K
复制代码

看输出信息是不是因为MM区域的几何结构有问题？但是我检查结构文件，发现除了部分MM区域跨过了边界其他都很正常，输入文件是Gromacs自动生成的，请问要怎么处理？
作者
Author: 啊不错的飞过海    时间: 2023-10-23 17:10
设的盒子是正交的吗？CP2K在这里应该是检查过盒子正交性，不正交的话跑不了
作者
Author: beyond    时间: 2023-10-23 18:44
&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了解不多的话，跑起来估计会有很多问题不知道如何解决的。。。
作者
Author: ZetaFunction    时间: 2023-10-23 19:22
啊不错的飞过海 发表于 2023-10-23 17:10
设的盒子是正交的吗？CP2K在这里应该是检查过盒子正交性，不正交的话跑不了

我设的是截角八面体的盒子，是必须要方形的盒子吗？
作者
Author: ZetaFunction    时间: 2023-10-23 19:27
beyond 发表于 2023-10-23 18:44
&QMMM
    &CELL
      A 29.472 0.000 0.000

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

需要把cMD最后一帧的结构最小化后重新重新加热？1fs甚至0.5fs的间隔对于QM/MM来说耗时好像有点太长了，减小间隔时间是否是为了便于SCF收敛的考量？如果SCF可以收敛是不是就没有必要了。
作者
Author: ZetaFunction    时间: 2023-10-23 20:39
beyond 发表于 2023-10-23 18:44
&QMMM
    &CELL
      A 29.472 0.000 0.000

我重新设置了周期性盒子的参数，改成了方形盒子，但还是报相同的错误，尝试调低间隔时间也不行，先做QM/MM的最小化也不行。
作者
Author: ZetaFunction    时间: 2023-10-23 20:46
beyond 发表于 2023-10-23 18:44
&QMMM
    &CELL
      A 29.472 0.000 0.000

重新设定盒子参数后产生的inp文件：
  1. &GLOBAL
  2.   PRINT_LEVEL LOW
  3.   PROJECT GROMACS
  4.   RUN_TYPE ENERGY_FORCE
  5. &END GLOBAL
  6. &FORCE_EVAL
  7.   METHOD QMMM
  8.   &DFT
  9.     CHARGE 0
  10.     MULTIPLICITY 1
  11.     BASIS_SET_FILE_NAME  BASIS_MOLOPT
  12.     POTENTIAL_FILE_NAME  POTENTIAL
  13.     &MGRID
  14.       NGRIDS 5
  15.       CUTOFF 450
  16.       REL_CUTOFF 50
  17.       COMMENSURATE
  18.     &END MGRID
  19.     &SCF
  20.       SCF_GUESS RESTART
  21.       EPS_SCF 5.0E-8
  22.       MAX_SCF 20
  23.       &OT  T
  24.         MINIMIZER  DIIS
  25.         STEPSIZE   0.15
  26.         PRECONDITIONER FULL_ALL
  27.       &END OT
  28.       &OUTER_SCF  T
  29.         MAX_SCF 20
  30.         EPS_SCF 5.0E-8
  31.       &END OUTER_SCF
  32.     &END SCF
  33.     &XC
  34.       DENSITY_CUTOFF     1.0E-12
  35.       GRADIENT_CUTOFF    1.0E-12
  36.       TAU_CUTOFF         1.0E-12
  37.       &XC_FUNCTIONAL PBE
  38.       &END XC_FUNCTIONAL
  39.     &END XC
  40.     &QS
  41.      METHOD GPW
  42.      EPS_DEFAULT 1.0E-10
  43.      EXTRAPOLATION ASPC
  44.      EXTRAPOLATION_ORDER  4
  45.     &END QS
  46.   &END DFT
  47.   &QMMM
  48.     &CELL
  49.       A 17.442 0.000 0.000
  50.       B 0.000 17.442 0.000
  51.       C 10.070 10.070 14.241
  52.       PERIODIC XYZ
  53.     &END CELL
  54.     CENTER EVERY_STEP
  55.     CENTER_GRID TRUE
  56.     &WALLS
  57.       TYPE REFLECTIVE
  58.     &END WALLS
  59.     ECOUPL GAUSS
  60.     USE_GEEP_LIB 12
  61.     &PERIODIC
  62.       GMAX     1.0E+00
  63.       &MULTIPOLE ON
  64.          RCUT     1.0E+01
  65.          EWALD_PRECISION     1.0E-06
  66.       &END
  67.     &END PERIODIC
  68.     &QM_KIND H
  69.       MM_INDEX 6849 6850 6851 6852 6854 6855 6858 6860 6864 6866 6868 6870
  70.     &END QM_KIND
  71.     &QM_KIND C
  72.       MM_INDEX 6845 6847 6853 6856 6857 6861 6862 6863 6865 6867 6869
  73.     &END QM_KIND
  74.     &QM_KIND N
  75.       MM_INDEX 6848 6859
  76.     &END QM_KIND
  77.     &QM_KIND O
  78.       MM_INDEX 6844 6846
  79.     &END QM_KIND
  80.   &END QMMM
  81.   &MM
  82.     &FORCEFIELD
  83.       DO_NONBONDED FALSE
  84.     &END FORCEFIELD
  85.     &POISSON
  86.       &EWALD
  87.         EWALD_TYPE NONE
  88.       &END EWALD
  89.     &END POISSON
  90.   &END MM
  91.   &SUBSYS
  92.     &CELL
  93.       A 86.728 0.000 0.000
  94.       B 0.000 86.728 0.000
  95.       C 43.364 43.364 61.326
  96.       PERIODIC XYZ
  97.     &END CELL
  98.     &TOPOLOGY
  99.       COORD_FILE_NAME qmmm_em/_cp2k.pdb
  100.       COORD_FILE_FORMAT PDB
  101.       CHARGE_EXTENDED TRUE
  102.       CONNECTIVITY OFF
  103.       &GENERATE
  104.          &ISOLATED_ATOMS
  105.             LIST 1..45896
  106.          &END
  107.       &END GENERATE
  108.     &END TOPOLOGY
  109.     &KIND H
  110.       ELEMENT H
  111.       BASIS_SET DZVP-MOLOPT-GTH
  112.       POTENTIAL GTH-PBE
  113.     &END KIND
  114.     &KIND C
  115.       ELEMENT C
  116.       BASIS_SET DZVP-MOLOPT-GTH
  117.       POTENTIAL GTH-PBE
  118.     &END KIND
  119.     &KIND N
  120.       ELEMENT N
  121.       BASIS_SET DZVP-MOLOPT-GTH
  122.       POTENTIAL GTH-PBE
  123.     &END KIND
  124.     &KIND O
  125.       ELEMENT O
  126.       BASIS_SET DZVP-MOLOPT-GTH
  127.       POTENTIAL GTH-PBE
  128.     &END KIND
  129.     &KIND X
  130.       ELEMENT H
  131.     &END KIND
  132.   &END SUBSYS
  133. &END FORCE_EVAL
复制代码

调低时间间隔还是同样的报错，CPASSERT failed，cp_ddapc_methods.F:650，尝试做最小化也报错
作者
Author: ZetaFunction    时间: 2023-10-26 01:53
改用了cubic盒子，并且把原先错误命名为NA的钠离子重新命名为Na+，可以正常运行，问题解决了！但是计算的速度特别慢，QM区域才27个原子，用PBE居然每帧要6-7分钟。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3