CP2K的QM/MM速度与纯QM速度会相差几十倍呢？

函数与激情

我使用CP2K的QM/MM方法想跑个有机小分子（11个原子）在水盒子中的metadynamics，QM区的计算方法是PBE0/DZVP-MOLOPT-SR-GTH，力场是使用的是AMBER。结果发现加上了MM环境后（水分子总数约2700原子），速度大幅下降，SCF迭代速度几乎是气相中纯小分子的60倍。我之前做过chemshell的QM/MM，根据我的经验，MM区对整体速度的影响是十分小的，请问大家有没有在CP2K遇到过这种问题，谢谢大家啦！

1. &GLOBAL
   
2.   PROJECT MD                     ! Name of the calculation
   
3.   PRINT_LEVEL LOW                ! Verbosity of the output
   
4.   RUN_TYPE MD                    ! Calculation type: MD
   
5. &END GLOBAL
   
6. 
   
7. &FORCE_EVAL                      ! parameters needed to calculate energy and forces
   
8.   METHOD QMMM                    ! Hybrid quantum classical
   
9.   STRESS_TENSOR ANALYTICAL       ! Compute the stress tensor analytically (if available).
   
10. 
    
11.   &DFT
    
12.     BASIS_SET_FILE_NAME  BASIS_MOLOPT
    
13.     POTENTIAL_FILE_NAME  POTENTIAL
    
14.     WFN_RESTART_FILE_NAME MD-RESTART.wfn
    
15.     CHARGE    0 #Net charge
    
16.     MULTIPLICITY    1 #Spin multiplicity
    
17.     &QS
    
18.       EPS_DEFAULT 1E-10 #This is default. Set all EPS_xxx to values such that the energy will be correct up to this value
    
19.     # EPS_PGF_ORB 1E-8 #If seeing "Kohn Sham matrix not 100% occupied" when using HF, MP2 or hybrid functional, uncomment this
    
20.     &END QS
    
21.     &POISSON
    
22.       PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
    
23.       PSOLVER PERIODIC #The way to solve Poisson equation
    
24.     &END POISSON
    
25.     &XC
    
26.       &XC_FUNCTIONAL
    
27.         &PBE
    
28.           SCALE_X 0.75 #75% GGA exchange
    
29.           SCALE_C 1.0 #100% GGA correlation
    
30.         &END PBE
    
31.       &END XC_FUNCTIONAL
    
32.       &HF
    
33.         FRACTION 0.25
    
34.         &SCREENING
    
35.         FRACTION 0.25
    
36.         &SCREENING
    
37.           EPS_SCHWARZ 1E-8 #Important to improve scaling. The larger the value, the lower the cost and lower the accuracy
    
38.           SCREEN_ON_INITIAL_P T #Screening on product between maximum of density matrix elements and ERI
    
39.         &END SCREENING
    
40.         &INTERACTION_POTENTIAL
    
41.           POTENTIAL_TYPE TRUNCATED
    
42.           CUTOFF_RADIUS 6.0 #Cutoff radius for truncated 1/r potential
    
43.           T_C_G_DATA t_c_g.dat
    
44.         &END INTERACTION_POTENTIAL
    
45.         &MEMORY
    
46.           MAX_MEMORY 3000 #Memory(MB) per MPI process for calculating HF exchange
    
47.           EPS_STORAGE_SCALING 0.1
    
48.         &END MEMORY
    
49.       &END HF
    
50.     &END XC
    
51.     &MGRID
    
52.       COMMENSURATE
    
53.       CUTOFF 300
    
54.       REL_CUTOFF 40
    
55.     &END MGRID
    
56.     &SCF
    
57.       MAX_SCF 256 #Maximum number of steps of inner SCF
    
58.       EPS_SCF 1.0E-05 #Convergence threshold of density matrix of inner SCF
    
59.       SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
    
60.       &OT
    
61.         PRECONDITIONER FULL_SINGLE_INVERSE #Usually best but expensive for large system. Cheaper: FULL_SINGLE_INVERSE and FULL_KINETIC (default)
    
62.         MINIMIZER DIIS #CG is worth to consider in difficult cases
    
63.         LINESEARCH 2PNT #1D line search algorithm for CG. 2PNT is default, 3PNT is better but more costly. GOLD is best but very expensive
    
64.       &END OT
    
65.     &END SCF
    
66.   &END DFT
    
67. 
    
68. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PART2 MM
    
69. &MM                              ! Parameters to run a MM calculation
    
70.     &FORCEFIELD                  ! Set up a force_field for the classical calculations
    
71.       PARMTYPE AMBER             ! Kind of torsion potential
    
72.                                  ! Filename that contains the parameters of the FF
    
73.       PARM_FILE_NAME complex_LJ_mod.prmtop
    
74.       &SPLINE                    ! Parameters to set up the splines used in the nonboned interactions
    
75.         EMAX_SPLINE 1.0E8        ! Maximum value of the potential up to which splines will be constructed
    
76.         RCUT_NB [angstrom] 10    ! Cutoff radius for nonbonded interactions
    
77.       &END SPLINE
    
78.     &END FORCEFIELD
    
79.     &POISSON
    
80.       &EWALD
    
81.                                  ! Ewald parameters controlling electrostatic
    
82.         EWALD_TYPE SPME          ! Type of ewald
    
83.         ALPHA .40                ! Alpha parameter associated with Ewald (EWALD|PME|SPME)
    
84.         GMAX 80                  ! Number of grid points (SPME and EWALD)
    
85.       &END EWALD
    
86.     &END POISSON
    
87. &END MM
    
88. 
    
89. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PART3 QMMM
    
90.   &QMMM
    
91.     ECOUPL GAUSS
    
92.     USE_GEEP_LIB 6
    
93.     MM_POTENTIAL_FILE_NAME MM_POTENTIAL
    
94.     &CELL
    
95.       ABC 32.2 32.2 80.0
    
96.     &END CELL
    
97.    !CENTER EVERY_STEP
    
98.     CENTER NEVER
    
99.     NOCOMPATIBILITY
    
100.     &INTERPOLATOR
     
101.     &INTERPOLATOR
     
102.       EPS_R 1.0e-14
     
103.       EPS_X 1.0e-14
     
104.       MAXITER 100
     
105.     &END INTERPOLATOR
     
106.     &QM_KIND H
     
107.     MM_INDEX 2 3 4
     
108.     &END QM_KIND
     
109.     &QM_KIND O
     
110.     MM_INDEX 6 7 8 10 11
     
111.     &END QM_KIND
     
112.     &QM_KIND N
     
113.      MM_INDEX 1 5 9
     
114.     &END QM_KIND
     
115.   &END QMMM
     
116. 
     
117. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PART4 SUBSYS
     
118. &SUBSYS                          ! a subsystem: coordinates, topology, molecules and cell
     
119.     &CELL                        !Set box dimensions here
     
120.       ABC [angstrom] 32.2 32.2 80.0
     
121.       ALPHA_BETA_GAMMA 90.0 90.0 90.0
     
122.     &END CELL
     
123.     &TOPOLOGY                    ! Topology for classical runs
     
124.       CONN_FILE_FORMAT AMBER
     
125.       CONN_FILE_NAME complex_LJ_mod.prmtop
     
126.       COORD_FILE_FORMAT XYZ
     
127.       COORD_FILE_NAME 2220.xyz
     
128.     &END TOPOLOGY
     
129. &KIND O
     
130.   ELEMENT O
     
131.   BASIS_SET DZVP-MOLOPT-SR-GTH
     
132.   POTENTIAL GTH-PBE
     
133. &END KIND
     
134. &KIND H
     
135.   ELEMENT H
     
136.   BASIS_SET DZVP-MOLOPT-SR-GTH
     
137.   POTENTIAL GTH-PBE
     
138. &END KIND
     
139. &KIND N
     
140.   ELEMENT N
     
141.   BASIS_SET DZVP-MOLOPT-SR-GTH
     
142.   POTENTIAL GTH-PBE
     
143. &END KIND
     
144. 
     
145. &COLVAR
     
146.   &DISTANCE_FUNCTION
     
147.    COEFFICIENT -1
     
148.    ATOMS 1 9 1 3
     
149.   &END
     
150. &END COLVAR
     
151. 
     
152. &END SUBSYS
     
153. &END FORCE_EVAL
     
154. 
     
155. &MOTION
     
156.   &FREE_ENERGY
     
157.     &METADYN
     
158.       DO_HILLS T
     
159.      !LAGRANGE F
     
160.       NT_HILLS 50 !Specify the maximum MD step interval between spawning two hills
     
161.      !WW 1.0e-3 !Specifies the height of the gaussian to spawn. Default 0.1
     
162.       WW [kcalmol] 0.12
     
163.       &METAVAR
     
164.         SCALE [angstrom] 0.25
     
165.         COLVAR 1
     
166.         &WALL
     
167.           TYPE REFLECTIVE
     
168.           POSITION [angstrom] 1.35
     
169.           &REFLECTIVE
     
170.              DIRECTION WALL_PLUS
     
171.           &END
     
172.         &END
     
173.         &WALL
     
174.           TYPE REFLECTIVE
     
175.           POSITION [angstrom] 0.80
     
176.           &REFLECTIVE
     
177.              DIRECTION WALL_MINUS
     
178.           &END
     
179.         &END
     
180.       &END METAVAR
     
181. 
     
182.      &PRINT
     
183.        &COLVAR
     
184.           COMMON_ITERATION_LEVELS 3
     
185.          &EACH
     
186.            MD 10
     
187.          &END
     
188.        &END
     
189.        &HILLS
     
190.          COMMON_ITERATION_LEVELS 3
     
191.           &EACH
     
192.             MD 10
     
193.           &END
     
194.         &END
     
195.      &END
     
196.     &END METADYN

复制代码

Aridea

大佬，这个60差别可能是杂化泛函导致的，如果换纯泛函不知是多少倍
另外请教大佬一个问题，我看别人周期性QMMM两个区域都要定义盒子的，这个盒子尺寸怎么确定合理性呢

函数与激情

QM盒子尺寸根据你的QM区来定义，不用太大，大了会影响速度，我这个速度差异就是因为QM区盒子尺寸不合适造成。

elpa

遇到过。但是OT每一步的耗时不知道为什么增加那么多。
但是可能由于MM分子的影响，每一步结构比纯QM变化大，ASPC靠前几步波函数组合做初猜，所以要更多的OT迭代步数才收敛。

yjb

你好，有cp2k跑QM/MM MD的教程吗？