CP2K中QM/MM使用GFN1-xTB报错：EWALD模块缺少GMAX

ZetaFunction

最近在尝试使用Gromacs+CP2K做酶体系的QM/MM模拟，QM区域187个原子，打算使用GFN1-xTB描述。根据网上的教程，首先用Gromacs内置的PBE方法生成CP2K输入文件_cp2k.inp，然后再编辑修改这个文件以使用自定义的方法。GFN1-xTB部分的输入参数我使用了从Multiwfn生成的参数，复制替换了原输入文件的&DFT部分，同时删去了原文件用来设定基组和XC的部分，SCF迭代精度也参照Multiwfn生成的文件做了修改。

运行模拟后CP2K报错：Value requested, but no value set getting value from keyword GMAX of section EWALD，关闭&xTB中DO_EWALD F后依然得到相同的报错。在网上查阅了一些信息后，对照CP2K手册发现ECOUPL GAUSS只有GPW/GAPW可以试用，DFTB需要用ECOUPL COULOMB，和DFTB原理类似的GFN1-xTB应该也差不多，但是我在修改了耦合方式为ECOUPL COULOMB后还是得到一样的报错。我又发现网上改用ECOUPL COULOMB的输入文件对比ECOUPL GAUSS删去了&QMMM中&PERIODIC的部分，我也删去这一部分后报错变成了：Cholesky decompose failed: the matrix is not positive definite or  ill-conditioned.

我对CP2K并不熟悉，请教各位大佬，若要在QM区域使用GFN1-xTB还需要做哪些修改？感谢各位大佬！

以下是使用PBE方法时Gromacs产生的输入文件，无论是直接模拟还是复制后用-qmi指定为自定义文件都可以正常运行没有问题：

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 34.417 0.000 0.000
    
50.     B 0.000 34.417 0.000
    
51.     C 0.000 0.000 34.417
    
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 1416 1417 1419 1420 1582 1583 1585 1588 1589 1591 1592 2784 2785 2787 2788 2790 2791 2998 2999 3001 3002 3004 3005 3006 3028 3029 3032 3034 3037 3039 3041 3934 3935 3937 3940 3941 3943 3944 4171 4172 4174 4175 4177 4178 4179 4206 4207 4210 4212 4215 4217 4219 5548 5549 5552 5554 5557 5559 5561 6397 6398 6401 6403 6406 6408 6410 6655 6656 6659 6661 6664 6666 6668 6831 6832 6834 6837 6838 6839 6841 6842 6843 6849 6850 6851 6852 6854 6855 6858 6860 6864 6866 6868 6870
    
70.   &END QM_KIND
    
71.   &QM_KIND C
    
72.     MM_INDEX 1415 1418 1421 1581 1586 2783 2786 2997 3000 3027 3030 3031 3033 3035 3038 3040 3933 3938 4170 4173 4205 4208 4209 4211 4213 4216 4218 5547 5550 5551 5553 5555 5558 5560 6396 6399 6400 6402 6404 6407 6409 6654 6657 6658 6660 6662 6665 6667 6830 6833 6835 6836 6840 6845 6847 6853 6856 6857 6861 6862 6863 6865 6867 6869
    
73.   &END QM_KIND
    
74.   &QM_KIND N
    
75.     MM_INDEX 1584 1587 1590 2789 3003 3936 3939 3942 4176 6848 6859
    
76.   &END QM_KIND
    
77.   &QM_KIND O
    
78.     MM_INDEX 1422 1423 3036 4214 5556 6405 6663 6821 6823 6824 6825 6827 6828 6829 6844 6846
    
79.   &END QM_KIND
    
80.   &QM_KIND P
    
81.     MM_INDEX 6822 6826
    
82.   &END QM_KIND
    
83.   &LINK
    
84.     QM_INDEX 1415
    
85.     MM_INDEX 1413
    
86.   &END LINK
    
87.   &LINK
    
88.     QM_INDEX 1581
    
89.     MM_INDEX 1578
    
90.   &END LINK
    
91.   &LINK
    
92.     QM_INDEX 2783
    
93.     MM_INDEX 2780
    
94.   &END LINK
    
95.   &LINK
    
96.     QM_INDEX 2997
    
97.     MM_INDEX 2994
    
98.   &END LINK
    
99.   &LINK
    
100.     QM_INDEX 3027
     
101.     MM_INDEX 3025
     
102.   &END LINK
     
103.   &LINK
     
104.     QM_INDEX 3933
     
105.     MM_INDEX 3930
     
106.   &END LINK
     
107.   &LINK
     
108.     QM_INDEX 4170
     
109.     MM_INDEX 4167
     
110.   &END LINK
     
111.   &LINK
     
112.     QM_INDEX 4205
     
113.     MM_INDEX 4203
     
114.   &END LINK
     
115.   &LINK
     
116.     QM_INDEX 5547
     
117.     MM_INDEX 5545
     
118.   &END LINK
     
119.   &LINK
     
120.     QM_INDEX 6396
     
121.     MM_INDEX 6394
     
122.   &END LINK
     
123.   &LINK
     
124.     QM_INDEX 6654
     
125.     MM_INDEX 6652
     
126.   &END LINK
     
127.   &END QMMM
     
128.   &MM
     
129.     &FORCEFIELD
     
130.       DO_NONBONDED FALSE
     
131.     &END FORCEFIELD
     
132.     &POISSON
     
133.       &EWALD
     
134.         EWALD_TYPE NONE
     
135.       &END EWALD
     
136.     &END POISSON
     
137.   &END MM
     
138.   &SUBSYS
     
139.     &CELL
     
140.       A 82.706 0.000 0.000
     
141.       B 0.000 82.706 0.000
     
142.       C 0.000 0.000 82.706
     
143.       PERIODIC XYZ
     
144.     &END CELL
     
145.     &TOPOLOGY
     
146.       COORD_FILE_NAME qmmm_xTB/_cp2k.pdb
     
147.       COORD_FILE_FORMAT PDB
     
148.       CHARGE_EXTENDED TRUE
     
149.       CONNECTIVITY OFF
     
150.       &GENERATE
     
151.         &ISOLATED_ATOMS
     
152.           LIST 1..56453
     
153.         &END
     
154.       &END GENERATE
     
155.     &END TOPOLOGY
     
156.     &KIND H
     
157.       ELEMENT H
     
158.       BASIS_SET DZVP-MOLOPT-GTH
     
159.       POTENTIAL GTH-PBE
     
160.     &END KIND
     
161.     &KIND C
     
162.        ELEMENT C
     
163.        BASIS_SET DZVP-MOLOPT-GTH
     
164.        POTENTIAL GTH-PBE
     
165.     &END KIND
     
166.     &KIND N
     
167.       ELEMENT N
     
168.      BASIS_SET DZVP-MOLOPT-GTH
     
169.      POTENTIAL GTH-PBE
     
170.    &END KIND
     
171.    &KIND O
     
172.      ELEMENT O
     
173.      BASIS_SET DZVP-MOLOPT-GTH
     
174.      POTENTIAL GTH-PBE
     
175.    &END KIND
     
176.    &KIND P
     
177.      ELEMENT P
     
178.      BASIS_SET DZVP-MOLOPT-GTH
     
179.      POTENTIAL GTH-PBE
     
180.    &END KIND
     
181.    &KIND X
     
182.      ELEMENT H
     
183.    &END KIND
     
184. &END SUBSYS
     
185. &END FORCE_EVAL

复制代码

以下是修改后的xTB_cp2k.inp，修改部分有注释：

1. <blockquote>&GLOBAL

复制代码

ZetaFunction

好像是帖子太长了。。。
以下是修改后的xTB_cp2k.inp，修改部分有注释：
&GLOBAL
  PRINT_LEVEL LOW
  PROJECT GROMACS
  RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
  METHOD QMMM
  &DFT    #删去了定义基组和XC的部分
    CHARGE 0
    MULTIPLICITY 1
    &MGRID
      NGRIDS 5
      CUTOFF 450
      REL_CUTOFF 50
      COMMENSURATE
    &END MGRID
    &SCF    #根据Multiwfn的输入文件修改了迭代精度
      SCF_GUESS RESTART
      EPS_SCF 1.0E-5
      MAX_SCF 25
      &OT T    #根据Multiwfn的输入文件修改了方法
        MINIMIZER  DIIS
        LINESEARCH 2PNT
        ALGORITHM STRICT
        PRECONDITIONER FULL_SINGLE_INVERSE
      &END OT
      &OUTER_SCF T
        MAX_SCF 25
        EPS_SCF 1.0E-05
      &END OUTER_SCF
      &PRINT    #根据Multiwfn的输入文件关闭无意义的输出节约I/O时间
        &RESTART OFF
        &END RESTART
      &END PRINT
    &END SCF
    &QS    #主要修改部分
      EPS_DEFAULT 1.0E-10
      EXTRAPOLATION ASPC
      EXTRAPOLATION_ORDER 3
      METHOD XTB
      &xTB    #使用GFN1-xTB
        DO_EWALD T    #尝试过设为F，但还是报错
        CHECK_ATOMIC_CHARGES F
        &PARAMETER
          DISPERSION_PARAMETER_FILE dftd3.dat
          PARAM_FILE_NAME xTB_parameters
        &END PARAMETER
      &END xTB
      #&SE    #参照网上找到的一些输入文件添加的部分，添加了还是报错
      #  &COULOMB
      #    CUTOFF [angstrom] 10.0
      #  &END
      #  &EXCHANGE
      #    CUTOFF [angstrom] 10.0
      #  &END
      #&END
    &END QS
  &END DFT
  &QMMM
    &CELL
      A 34.417 0.000 0.000
      B 0.000 34.417 0.000
      C 0.000 0.000 34.417
      PERIODIC XYZ
    &END CELL
    CENTER EVERY_STEP
    CENTER_GRID TRUE
    &WALLS
      TYPE REFLECTIVE
    &END WALLS
    ECOUPL COULOMB
    #USE_GEEP_LIB 12    #注释掉还是报错
    &PERIODIC    #尝试过删去，然后报错：Cholesky decompose failed
      GMAX     1.0E+00
      &MULTIPOLE ON
        RCUT     1.0E+01
        EWALD_PRECISION     1.0E-06
      &END
    &END PERIODIC
    &QM_KIND H
      MM_INDEX 1416 1417 1419 1420 1582 1583 1585 1588 1589 1591 1592 2784 2785 2787 2788 2790 2791 2998 2999 3001 3002 3004 3005 3006 3028 3029 3032 3034 3037 3039 3041 3934 3935 3937 3940 3941 3943 3944 4171 4172 4174 4175 4177 4178 4179 4206 4207 4210 4212 4215 4217 4219 5548 5549 5552 5554 5557 5559 5561 6397 6398 6401 6403 6406 6408 6410 6655 6656 6659 6661 6664 6666 6668 6831 6832 6834 6837 6838 6839 6841 6842 6843 6849 6850 6851 6852 6854 6855 6858 6860 6864 6866 6868 6870
    &END QM_KIND
    &QM_KIND C
      MM_INDEX 1415 1418 1421 1581 1586 2783 2786 2997 3000 3027 3030 3031 3033 3035 3038 3040 3933 3938 4170 4173 4205 4208 4209 4211 4213 4216 4218 5547 5550 5551 5553 5555 5558 5560 6396 6399 6400 6402 6404 6407 6409 6654 6657 6658 6660 6662 6665 6667 6830 6833 6835 6836 6840 6845 6847 6853 6856 6857 6861 6862 6863 6865 6867 6869
    &END QM_KIND
    &QM_KIND N
      MM_INDEX 1584 1587 1590 2789 3003 3936 3939 3942 4176 6848 6859
    &END QM_KIND
    &QM_KIND O
      MM_INDEX 1422 1423 3036 4214 5556 6405 6663 6821 6823 6824 6825 6827 6828 6829 6844 6846
    &END QM_KIND
    &QM_KIND P
      MM_INDEX 6822 6826
    &END QM_KIND
    &LINK
      QM_INDEX 1415
      MM_INDEX 1413
    &END LINK
    &LINK
      QM_INDEX 1581
      MM_INDEX 1578
    &END LINK
    &LINK
      QM_INDEX 2783
      MM_INDEX 2780
    &END LINK
    &LINK
      QM_INDEX 2997
      MM_INDEX 2994
    &END LINK
    &LINK
      QM_INDEX 3027
      MM_INDEX 3025
    &END LINK
    &LINK
      QM_INDEX 3933
      MM_INDEX 3930
    &END LINK
    &LINK
      QM_INDEX 4170
      MM_INDEX 4167
    &END LINK
    &LINK
      QM_INDEX 4205
      MM_INDEX 4203
    &END LINK
    &LINK
      QM_INDEX 5547
      MM_INDEX 5545
    &END LINK
    &LINK
      QM_INDEX 6396
      MM_INDEX 6394
    &END LINK
    &LINK
      QM_INDEX 6654
      MM_INDEX 6652
    &END LINK
  &END QMMM
  &MM
    &FORCEFIELD
      DO_NONBONDED FALSE
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE NONE
      &END EWALD
    &END POISSON
  &END MM
  &SUBSYS
    &CELL
      A 82.706 0.000 0.000
      B 0.000 82.706 0.000
      C 0.000 0.000 82.706
      PERIODIC XYZ
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME qmmm_xTB/_cp2k.pdb
      COORD_FILE_FORMAT PDB
      CHARGE_EXTENDED TRUE
      CONNECTIVITY OFF
      &GENERATE
        &ISOLATED_ATOMS
          LIST 1..56453
        &END
      &END GENERATE
    &END TOPOLOGY
    &KIND H
      ELEMENT H
    &END KIND
    &KIND C
      ELEMENT C
    &END KIND
    &KIND N
      ELEMENT N
    &END KIND
    &KIND O
      ELEMENT O
    &END KIND
    &KIND P
      ELEMENT P
    &END KIND
    &KIND X
      ELEMENT H
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

啊不错的飞过海

报错讲的不是&EWALD段落缺参数吗？关它干什么
GMAX这个参数描述时EWALD/SPME时的最大格点数，应当出现在在你给出的输入文件134行，加上就是。没想法的话可以先参照输入参考里的25写，后面再调大/小些看这个格点规模够不够描述这个体系。
参考：
https://manual.cp2k.org/trunk/CP ... /POISSON/EWALD.html

ZetaFunction

啊不错的飞过海 发表于 2023-10-29 20:39
报错讲的不是&EWALD段落缺参数吗？关它干什么
GMAX这个参数描述时EWALD/SPME时的最大格点数，应当出现在在 ...

感谢回复。
请问你指的是在&MM部分加上GMAX吗？我尝试了以下，还是同样的报错。&MM部分我复制的是Gromacs自动生成的输入文件，这个输入文件可以正常运行，而且根据手册，设定了EWALD_TYPE NONE后无需GMAX参数。

啊不错的飞过海

ZetaFunction 发表于 2023-10-29 23:08
感谢回复。
请问你指的是在&MM部分加上GMAX吗？我尝试了以下，还是同样的报错。&MM部分我复制的是Gromac ...

是&EWALD部分。

ZetaFunction

啊不错的飞过海 发表于 2023-10-30 20:56
是&EWALD部分。

我就是在&MM部分添加了&EWALD，还是相同的报错。

啊不错的飞过海

ZetaFunction 发表于 2023-10-30 21:34
我就是在&MM部分添加了&EWALD，还是相同的报错。

不知道你加的时候设没设gmax，总之按下面注释的说法cp2k的xtb模块确实没设gmax的默认值，需要手动设。xtb调用ewald的部分不待要找了，楼主或许可以找找看。[url=]https://github.com/cp2k/cp2k/blob/c03a13580630548167f5a509f6a2d79eca29359f/src/xtb_types.F#L13[/url]

we000000

要在&FORCE_EVAL-&QMMM-&PERIODIC-&POISSON-&EWALD这个部分设定GMAX参数，就是在你之前删掉的&PERIODIC里面。不是&xTB下面的DO_EWALD和%MM下面的GMAX。具体细节可以参考CP2K手册：https://manual.cp2k.org/trunk/CP ... /POISSON/EWALD.html