cp2k QMMM PBE+def2- 震荡不收敛

yinxing

[size=0.92857em]@sobereva 卢老师您好，我参考https://www.cp2k.org/howto:biochem_qmmm 的例子，通过QMMM研究镧离子对蛋白酶催化活性的影响。复合物系统采用Ambertools18中的MCPB.py进行建模。一直顺利进行到cp2k教程中如下图片标记的步骤，初步采取PBE纯泛函与小基组也长时间震荡无法收敛。输入文件protein_solv_monitor_pbe.inp 中采取PBE相关参数采用Multiwfn生成，QMMM参数参考cp2k8.2相关资料设置。鉴于我理论知识浅薄描述不够详细，初步猜测自行设置的相关参数有问题。想请您百忙之中给予指导，非常感谢！系统的所有文件见附件。

#Generated by Multiwfn
&GLOBAL
  PROJECT Protein_solv_Monitor_PBE
  PRINT_LEVEL LOW
  RUN_TYPE MD
&END GLOBAL

&FORCE_EVAL
  METHOD QMMM
  STRESS_TENSOR ANALYTICAL
  &DFT
    BASIS_SET_FILE_NAME  EMSL_BASIS_SETS
    POTENTIAL_FILE_NAME  POTENTIAL
    WFN_RESTART_FILE_NAME 5NMY_solv_monitor_pbe-RESTART.wfn
    CHARGE 0 #Net charge
    MULTIPLICITY    1 #Spin multiplicity
    &QS
      EPS_DEFAULT 1E-10 #This is default. Set all EPS_xxx to values such that the energy will be correct up to this value
      EXTRAPOLATION PS #Extrapolation for wavefunction during e.g. MD. ASPC is default, PS also be used
      EXTRAPOLATION_ORDER 2 #Order for PS or ASPC extrapolation. 3 is default
      METHOD GAPW
    &END QS
    &POISSON
      PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
      PSOLVER PERIODIC #The way to solve Poisson equation
    &END POISSON
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
      &VDW_POTENTIAL
        POTENTIAL_TYPE PAIR_POTENTIAL
        &PAIR_POTENTIAL
          PARAMETER_FILE_NAME dftd3.dat
          TYPE DFTD3
          REFERENCE_FUNCTIONAL PBE
        &END PAIR_POTENTIAL
      &END VDW_POTENTIAL
    &END XC
    &MGRID
      CUTOFF 400
      REL_CUTOFF 55
      COMMENSURATE
    &END MGRID
    &SCF
      MAX_SCF 128
      EPS_SCF 3.0E-05 #Convergence threshold of density matrix during SCF
      SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
      &DIAGONALIZATION
        ALGORITHM STANDARD #Algorithm for diagonalization. DAVIDSON is faster for large systems
      &END DIAGONALIZATION
      &MIXING #How to mix old and new density matrices
        METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
        ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one
        NBROYDEN 32 #Default is 4. Number of previous steps stored for the actual mixing scheme
      &END MIXING
      &OUTER_SCF
        MAX_SCF 64 #Maximum number of steps of outer SCF
        EPS_SCF 1.0E-05 #Convergence threshold of outer SCF
      &END OUTER_SCF
      &SMEAR
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE 300 #Electronic temperature of Fermi-Dirac smearing in K
      &END SMEAR
      ADDED_MOS 100 #Number of virtual MOs to be solved
      &PRINT
        &RESTART #Use "&RESTART OFF" can prevent generating wfn file
          BACKUP_COPIES 0 #Maximum number of backup copies of wfn file
        &END RESTART
      &END PRINT
    &END SCF
    &PRINT
      &MO_MOLDEN
        NDIGITS 8
        GTO_KIND SPHERICAL
      &END MO_MOLDEN
      &MULLIKEN
         FILENAME=mulliken
         COMMON_ITERATION_LEVELS 10
         &EACH
           MD 1
         &END EACH
      &END MULLIKEN
      &PDOS
          NLUMO -1
      &END PDOS
    &END PRINT
  &END DFT

  &MM
    &FORCEFIELD
      PARMTYPE AMBER
      PARM_FILE_NAME 5NMY_solv_LJ_mod.prmtop
      &SPLINE
        EMAX_SPLINE 1.0E3
        RCUT_NB [angstrom] 10
      &END SPLINE
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE SPME
        ALPHA .40
        GMAX 80
      &END EWALD
    &END POISSON
  &END MM
  &SUBSYS
    &CELL
      A    68.76848667     0.00000000     0.00000000
      B     0.00000000   66.604543937     0.00000000
      C     0.00000000     0.00000000   86.978811596
      PERIODIC XYZ #Direction of applied PBC (geometry aspect)
    &END CELL
    &TOPOLOGY
      CONN_FILE_FORMAT AMBER
      CONN_FILE_NAME 5NMY_solv_LJ_mod.prmtop
    &END TOPOLOGY
#   &VELOCITY #You can set initial atomic velocities in this section
#   &END VELOCITY
    &KIND N   
      ELEMENT N
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND H   
      ELEMENT H
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND HX   
      ELEMENT H
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND C   
      ELEMENT C
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND O   
      ELEMENT O
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND He   
      ELEMENT He
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND S   
      ELEMENT S
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND La   
      ELEMENT La
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND M1   
      ELEMENT La
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &KIND Cl-   
      ELEMENT Cl
      BASIS_SET Ahlrichs-def2-SVP
      POTENTIAL ALL
    &END KIND
    &COLVAR
      &DISTANCE
        ATOMS 3739 3740
      &END DISTANCE
    &END COLVAR
    &COLVAR
      &DISTANCE
            ATOMS 3733 3734
      &END DISTANCE
    &END COLVAR
    &COLVAR
      &DISTANCE
        ATOMS 3725 3727
      &END DISTANCE
    &END COLVAR
  &END SUBSYS
  &QMMM
    ECOUPL GAUSS
    CENTER EVERY_STEP
    NOCOMPATIBILITY
    USE_GEEP_LIB 8
    PARALLEL_SCHEME GRID
    &CELL
      ABC 40 40 40
      ALPHA_BETA_GAMMA 90 90 90
    &END CELL
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 444
      QM_INDEX 446
    &END LINK
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 463
      QM_INDEX 461
    &END LINK

    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 1216
      QM_INDEX 1214
    &END LINK
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 1197
      QM_INDEX 1199
    &END LINK

    #ASP
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 505
      QM_INDEX 507
    &END LINK
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 519
      QM_INDEX 509
    &END LINK

    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 1692
      QM_INDEX 1694
    &END LINK
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 1711
      QM_INDEX 1696
    &END LINK

    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 1675
      QM_INDEX 1673
    &END LINK
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 1656
      QM_INDEX 1658
    &END LINK

    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 3728
      QM_INDEX 3727
    &END LINK
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 3736
      QM_INDEX 3734
    &END LINK

    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 3741
      QM_INDEX 3740
    &END LINK
    &LINK
      ALPHA 1.50
      LINK_TYPE IMOMM
      MM_INDEX 3738
      QM_INDEX 3739
    &END LINK

    # ASP283(35)  H362(114) H360(112) H329(81) H279(31) PEP
    &QM_KIND O
      MM_INDEX 1215 510 513 514 1697 3726 3730 3735 3744
    &END QM_KIND
    &QM_KIND N
      MM_INDEX 446 454 458 461 1199 1207 1210 507 1658 1666 1669 1694 1701 1703 3727 3733 3739
    &END QM_KIND
    &QM_KIND C
      MM_INDEX 448 450 453 456 459 1201 1203 1206 1208 1212 1214 508 509 511 512 1660 1662 1665 1667 1671 1673 1695 1696 1698 1699 1700 1702 3724 3725 3734 3740
    &END QM_KIND
    &QM_KIND H
      MM_INDEX  449 451 452 455 457 460 1200 1202 1204 1205 1209 1211 1213 515 516 517 518 1659 1661 1663 1664 1668 1670 1672 1704 1705 1706 1707 1708 1709 1710 3753 3754 3760 3766
    &END QM_KIND
    &QM_KIND LA
      MM_INDEX  3723
    &END QM_KIND
    &WALLS
      WALL_SKIN [angstrom] 1.5 1.5 1.5
      TYPE QUADRATIC
      K .1
    &END WALLS
  &END QMMM
&END FORCE_EVAL

&MOTION
  &MD
    ENSEMBLE NVT
    STEPS 10000
    TIMESTEP [fs] 0.5 #fs. Decrease it properly for high temperature simulation
    TEMPERATURE 310.15 #Initial and maintained temperature (K)
    &THERMOSTAT
      TYPE CSVR
      &CSVR
        TIMECON 200 #Time constant in fs. Smaller/larger results in stronger/weaker temperature coupling
      &END CSVR
    &END THERMOSTAT
  &END MD
  &FREE_ENERGY
    METHOD METADYN
    &METADYN
      DO_HILLS  .FALSE.
      NT_HILLS 100
      WW [kcalmol] 1.5
      &METAVAR
        WIDTH 0.5 !Also known as scale
        COLVAR 1
      &END METAVAR
      &METAVAR
        WIDTH 0.5 !Also known as scale
        COLVAR 2
      &END METAVAR
      &METAVAR
        WIDTH 0.5 !Also known as scale
        COLVAR 3
      &END METAVAR
      &PRINT
       &COLVAR
          COMMON_ITERATION_LEVELS 3
          &EACH
            MD 10
          &END EACH
       &END COLVAR
       &HILLS
        &EACH
          METADYNAMICS 1
        &END EACH
        COMMON_ITERATION_LEVELS 10
       &END HILLS
      &END
    &END METADYN
  &END FREE_ENERGY

  &PRINT
    &RESTART_HISTORY                            ! This section controls dumping of unique restart files during the run keeping all of them.Most useful if recovery is needed at a later point.
      &EACH                                     ! A new restart file will be printed every 10000 md steps
        MD 500
      &END
    &END
    &TRAJECTORY
      &EACH
        MD  100 #Output frequency of geometry
      &END EACH
      FORMAT xyz
    &END TRAJECTORY
    &VELOCITIES
      &EACH
        MD   100 #Output frequency of velocity
      &END EACH
    &END VELOCITIES
    &RESTART
      BACKUP_COPIES 0 #Maximum number of backing up restart file
      &EACH
        MD 1000 #Frequency of updating last restart file
      &END EACH
    &END RESTART
  &END PRINT
&END MOTION

&EXT_RESTART
  RESTART_FILE_NAME 5NMY_solv_NPT-1.restart
  RESTART_DEFAULT .FALSE.
  RESTART_POS true
&END

beyond

QMMM的equilibration你是怎么做的啊， minimization 还是annealing？
def2- 是全电子基组，可以检查下设置是否正确。。。
权限太高，下载不了附件。。。

话说，pdb id 都给出来了，就没有必要加这么高的权限限制吧。。。

yinxing

beyond 发表于 2022-4-25 18:30
QMMM的equilibration你是怎么做的啊， minimization 还是annealing？
def2- 是全电子基组，可以检查下设置 ...

非常感谢提醒，鉴于项目要求。equilibration用的是minimization.基组为def2-SVP。能量从开始运行始终在正负4.0左右震荡。如果乐意帮忙的话可以加我微信18810617756，将不胜感激。