计算化学公社

标题: 求助：ORCA联用XTB跑QMMM时SCC不收敛 [打印本页]

作者
Author: 霜晨月 时间: 2021-10-26 11:00
标题: 求助：ORCA联用XTB跑QMMM时SCC不收敛
本帖最后由霜晨月于 2021-10-27 16:51 编辑

（标题原来误写成SCF不收敛，已修改，抱歉）

我用Amber 18 + Orca 5.01 + xtb 6.4 跑一个酶的QM/MM动力学，QM区域是个有机小分子配体和几个氨基酸残基。
经常出现SCF不收敛的情况，而且不是一上来就不收敛，而是QM/MM MD正常跑了几百步、几千步之后出现不收敛。
（我仔细检查了MD轨迹和临时结构文件inpfile.xyz，未见结构异常。）

感觉XTB计算这么小的有机分子体系，SCF不收敛有些奇怪。请教各位老师问题可能出在何处？

这是Amber参数文件的相关设置：
&qmmm
qmmask=' XXXXXX ' ! These atoms should be treated using QM
qmshake=0, ! DO NOT SHAKE
qmcharge=1, ! Charge on QM region is +1
qm_theory='EXTERN',
qmcut=12.0 ! Use 12 angstrom cut off for QM region
qm_ewald=0 ! Required for qm_theory='EXTERN'
writepdb = 1
/
&orc
use_template=1,
/

这是ORCA template文件：
! XTB2
%pal nprocs 1 end
%maxcore 2000

错误信息：（orc_job.dat）
================================================================================
                                    WARNINGS
                     Please study these warnings very carefully!
================================================================================

WARNING: Found dipole moment calculation with XTB calculation
  ===> : Switching off dipole moment calculation

WARNING: TRAH-SCF for XTB is not implemented!
  ===> : Turning TRAH off!

================================================================================
                                    INPUT FILE
================================================================================
NAME = orc_job.inp
|  1> ! XTB2
|  2> %pal nprocs 1 end
|  3> %maxcore 2000
|  4> ! ENGRAD
|  5> ! Angs NoUseSym
|  6> %pointcharges "ptchrg.xyz"
|  7> *xyzfile 1 1 inpfile.xyz
|  8>
|  9>                         ****END OF INPUT****
================================================================================

         -------------------------------------------------
      |             Calculation Setup             |
         -------------------------------------------------

      program call             : /opt/software/orca5.01/orca_5_0_1/otool_xtb orc_job_XTB.xyz --grad -c 1 -u 0 -P 1 --namespace orc_job --input orc_job_XTB.input.tmp --acc 0.100000
      hostname                : localhost.localdomain
      calculation namespace    : orc_job
      coordinate file          : orc_job_XTB.xyz
      omp threads             :                   1
      number of atoms          :                   63
      number of electrons       :                168
      charge                   :                   1
      spin                      :                0.0
      first test random number :    0.46629982145293

ID Z sym. atoms
1 6 C    1, 4, 6, 10, 12, 14, 20-29, 49-52
2 1 H    2, 3, 7, 9, 11, 13, 15-17, 19, 36-47, 53-63
3 7 N    5, 8, 30, 31
4 8 O    18, 32-35, 48

molecular fragmentation (1/2 indicates fragments):
111111111112222222222222222222222222222222222222222222222212222
# atoms in fragment 1/2: 12 51
fragment masses (1/2) :    82.10    343.40
CMA distance (Bohr) : 9.559
constraining FC (au) :  0.0500

         -------------------------------------------------
      |                G F N 2 - x T B                |
         -------------------------------------------------

      Reference                   10.1021/acs.jctc.8b01176
   * Hamiltonian:
      H0-scaling (s, p, d)          1.850000 2.230000 2.230000
      zeta-weighting                0.500000
   * Dispersion:
      s8                            2.700000
      a1                            0.520000
      a2                            5.000000
      s9                            5.000000
   * Repulsion:
      kExp                         1.500000 1.000000
      rExp                         1.000000
   * Coulomb:
      alpha                         2.000000
      third order                   shell-resolved
      anisotropic                   true
      a3                            3.000000
      a5                            4.000000
      cn-shift                      1.200000
      cn-exp                      4.000000
      max-rad                      5.000000

q/qsh data taken from xtbrestart
CAMM data taken from xtbrestart

      ...................................................
      :                   SETUP                   :
      :.................................................:
      :  # basis functions                153       :
      :  # atomic orbitals                153       :
      :  # shells                         93       :
      :  # electrons                      168       :
      :  max. iterations                250       :
      :  Hamiltonian                GFN2-xTB       :
      :  restarted?                      true       :
      :  GBSA solvation                false       :
      :  PC potential                   true       :
      :  -> # point charges             1978       :
      :  -> sum of PC             -0.3358700    e :
      :  electronic temp.       300.0000000    K :
      :  accuracy                   0.1000000       :
      :  -> integral cutoff       0.3500000E+02    :
      :  -> integral neglect       0.1000000E-08    :
      :  -> SCF convergence       0.1000000E-06 Eh :
      :  -> wf. convergence       0.1000000E-04 e :
      :  Broyden damping          0.4000000       :
      ...................................................

iter    E          dE       RMSdq    gap    omega  full diag
1 -95.7691403 -0.957691E+02  0.223E+00 0.12    0.0  T
2 -96.1462537 -0.377113E+00  0.218E+00 0.26    1.0  T
3 -96.4393084 -0.293055E+00  0.927E-01 2.95    1.0  T
4 -96.4752596 -0.359512E-01  0.629E-01 2.57    1.0  T
5 -96.4865898 -0.113303E-01  0.478E-01 2.26    1.0  T
6 -96.4912654 -0.467562E-02  0.386E-01 1.71    1.0  T
7 -96.4969152 -0.564977E-02  0.224E-01 0.70    1.0  T
8 -95.6909752  0.805940E+00  0.233E+00 0.17    1.0  T
9 -96.4949687 -0.803993E+00  0.285E-01 1.11    1.0  T
  10 -96.4952918 -0.323127E-03  0.275E-01 1.05    1.0  T
…………………………………………
221 -96.4979424 -0.158422E-02  0.176E-01 0.44    1.0  T
222 -96.4988412 -0.898797E-03  0.135E-01 0.11    1.0  T
223 -96.4984675  0.373624E-03  0.155E-01 0.33    1.0  T
224 -96.4775389  0.209286E-01  0.403E-01 0.07    1.0  T
225 -96.4957224 -0.181834E-01  0.261E-01 0.96    1.0  T
226 -96.4986269 -0.290448E-02  0.149E-01 0.31    1.0  T
227 -96.4754178  0.232091E-01  0.420E-01 0.07    1.0  T
228 -96.4955178 -0.201000E-01  0.268E-01 1.00    1.0  T
229 -96.4985244 -0.300666E-02  0.153E-01 0.32    1.0  T
230 -96.4818816  0.166429E-01  0.365E-01 0.07    1.0  T
231 -96.4962476 -0.143660E-01  0.244E-01 0.85    1.0  T
232 -96.4991079 -0.286029E-02  0.134E-01 0.27    1.0  T
233 -96.4447522  0.543557E-01  0.619E-01 0.05    1.0  T
234 -96.4921655 -0.474133E-01  0.360E-01 1.58    1.0  T
235 -96.4975483 -0.538278E-02  0.193E-01 0.54    1.0  T
236 -96.5011361 -0.358787E-02  0.502E-02 0.17    1.0  T
237 -96.4567729  0.443632E-01  0.550E-01 0.05    1.0  T
238 -96.4934069 -0.366340E-01  0.329E-01 1.38    1.0  T
239 -96.4978144 -0.440748E-02  0.182E-01 0.47    1.0  T
240 -96.5014928 -0.367838E-02  0.368E-03 0.15    3.4  T
241 -96.4989076  0.258520E-02  0.133E-01 0.11    1.0  T
242 -96.4984142  0.493423E-03  0.157E-01 0.34    1.0  T
243 -96.4903642  0.804994E-02  0.275E-01 0.09    1.0  T
244 -96.4971622 -0.679792E-02  0.209E-01 0.64    1.0  T
245 -96.5005857 -0.342356E-02  0.804E-02 0.19    1.0  T
246 -96.4363976  0.641882E-01  0.663E-01 0.04    1.0  T
247 -96.4911151 -0.547175E-01  0.384E-01 1.73    1.0  T
248 -96.4972676 -0.615250E-02  0.205E-01 0.61    1.0  T
249 -96.5001448 -0.287719E-02  0.985E-02 0.21    1.0  T
250 -96.4243660  0.757788E-01  0.721E-01 0.04    1.0  T

*** convergence criteria cannot be satisfied within 250 iterations ***

      # Occupation          Energy/Eh          Energy/eV
   -------------------------------------------------------------
      1       2.0000          -0.8070918          -21.9621
   ...          ...                ...                ...
      78       2.0000          -0.4511214          -12.2756
      79       2.0000          -0.4500486          -12.2464
      80       2.0000          -0.4448855          -12.1059
      81       2.0000          -0.4372995          -11.8995
      82       2.0000          -0.4308839          -11.7249
      83       2.0000          -0.3971413          -10.8068
      84       1.3381          -0.3761309          -10.2350 (HOMO)
      85       0.6619          -0.3747934          -10.1986 (LUMO)
      86                      -0.3591114             -9.7719
      87                      -0.3027914             -8.2394
      88                      -0.2896251             -7.8811
      89                      -0.2210509             -6.0151
   ...                               ...                ...
   153                         0.6289752             17.1153
   -------------------------------------------------------------
               HL-Gap          0.0013375 Eh          0.0364 eV
         Fermi-level          -0.3754621 Eh       -10.2168 eV

SCC (total)                0 d,  0 h,  0 min,  1.483 sec
SCC setup                   ...       0 min,  0.003 sec (  0.227%)
Dispersion                   ...       0 min,  0.002 sec (  0.108%)
classical contributions       ...       0 min,  0.000 sec (  0.019%)
integral evaluation          ...       0 min,  0.016 sec (  1.054%)
iterations                   ...       0 min,  1.411 sec ( 95.115%)
molecular gradient          ...       0 min,  0.051 sec (  3.441%)
printout                      ...       0 min,  0.000 sec (  0.034%)

########################################################################
[ERROR] Program stopped due to fatal error
-3- Single point calculation terminated
-2- xtb_calculator_singlepoint: Electronic structure method terminated
-1- scf: Self consistent charge iterator did not converge
########################################################################
[file orca_main/run.cpp, line 12281]: ORCA finished with error return - aborting the run

作者
Author: wzkchem5 时间: 2021-10-26 15:56
如果确认没有结构异常的话，有没有奇怪或者特殊（但是正确）的结构？比如自由基、卡宾、非休克尔共轭体系之类的电子结构比较特殊的结构？
此外可以把最后一帧拿出来算个QM/MM单点，看看是否收敛；如果不收敛，只算QM区，看看是否收敛。逐次缩小问题范围，我们才好回答

作者
Author: 霜晨月 时间: 2021-10-26 17:22
本帖最后由霜晨月于 2021-10-26 17:25 编辑

wzkchem5 发表于 2021-10-26 15:56
如果确认没有结构异常的话，有没有奇怪或者特殊（但是正确）的结构？比如自由基、卡宾、非休克尔共轭体系之 ...

谢谢您。这些特殊结构都没有，两个有机小分子+两个氨基酸残基，请见附件。（未配对电子数为0，电荷为+1，吡啶部分带一个正电荷）

QM区域单独用xtb算，很容易就收敛了，比如以下命令很容易收敛：
xtb orc_job_XTB.xyz -c 1 -u 0 > xtb-singlepoint-test.log
xtb orc_job_XTB.xyz --grad -c 1 -u 0 > xtb-grad-test.log
xtb orc_job_XTB.xyz --grad -c 1 -u 0 --input orc_job_XTB.input.tmp > xtb-grad-input-test.log

但是如果用ORCA调用xtb的命令（见下），就不收敛：
/opt/soft/orca5.01/orca501/otool_xtb orc_job_XTB.xyz --grad -c 1 -u 0 -P 1 --namespace orc_job --input orc_job_XTB.input.tmp --acc 0.100000 > otool_xtb-test.log

附件是一个不收敛的构象，以及相关的输入输出文件，烦请老师指点。（这个构象前面已经跑了1700步，都没有出错，而且构象也没有显著变化。）
(, 下载次数 Times of downloads: 14)

作者
Author: wzkchem5 时间: 2021-10-26 18:23

霜晨月发表于 2021-10-26 10:22
谢谢您。这些特殊结构都没有，两个有机小分子+两个氨基酸残基，请见附件。（未配对电子数为0，电荷为+1， ...

有没有试过用otool_xtb调用但不加点电荷？怀疑是某个点电荷和某个QM区原子跑得太近了，导致polarization catastrophe

作者
Author: 霜晨月 时间: 2021-10-26 20:28

wzkchem5 发表于 2021-10-26 18:23
有没有试过用otool_xtb调用但不加点电荷？怀疑是某个点电荷和某个QM区原子跑得太近了，导致polarization ...

尝试了otool_xtb orc_job_XTB.xyz --grad -c 1 -u 0 --input orc_job_XTB.input.tmp
正常结束，3步收敛。

我测试的结果是：
QM/MM调用命令otool_xtb orc_job_XTB.xyz --grad -c 1 -u 0 -P 1 --namespace orc_job --input orc_job_XTB.input.tmp --acc 0.100000 不收敛
但只要把 --namespace orc_job 去掉，即otool_xtb orc_job_XTB.xyz --grad -c 1 -u 0 -P 1 --input orc_job_XTB.input.tmp --acc 0.100000 就正常收敛

比较一下输出结果，没有--namespace orc_job 的时候：
      ...................................................
      :                   SETUP                   :
      :.................................................:
      :  # basis functions                153       :
      :  # atomic orbitals                153       :
      :  # shells                         93       :
      :  # electrons                      168       :
      :  max. iterations                250       :
      :  Hamiltonian                GFN2-xTB       :
      :  restarted?                      true       :
      :  GBSA solvation                false       :
      :  PC potential                   false       :
      :  electronic temp.       300.0000000    K :
      :  accuracy                   1.0000000       :
      :  -> integral cutoff       0.2500000E+02    :
      :  -> integral neglect       0.1000000E-07    :

      :  -> SCF convergence       0.1000000E-05 Eh :
      :  -> wf. convergence       0.1000000E-03 e :
      :  Broyden damping          0.4000000       :
      ...................................................

iter    E          dE       RMSdq    gap    omega  full diag
1 -96.3712391 -0.963712E+02  0.109E-05 0.53    0.0  T
2 -96.3712391 -0.852651E-12  0.773E-06 0.53 1629.8  T
3 -96.3712391 -0.139266E-11  0.467E-06 0.53 2695.0  T

*** convergence criteria satisfied after 3 iterations ***

有--namespace orc_job 的时候：
      ...................................................
      :                   SETUP                   :
      :.................................................:
      :  # basis functions                153       :
      :  # atomic orbitals                153       :
      :  # shells                         93       :
      :  # electrons                      168       :
      :  max. iterations                250       :
      :  Hamiltonian                GFN2-xTB       :
      :  restarted?                      true       :
      :  GBSA solvation                false       :
        :  PC potential                   true       :
      :  -> # point charges             1983       :
      :  -> sum of PC             -2.7782068    e :
      :  electronic temp.       300.0000000    K :
      :  accuracy                   0.1000000       :

      :  -> integral cutoff       0.3500000E+02    :
      :  -> integral neglect       0.1000000E-08    :
      :  -> SCF convergence       0.1000000E-06 Eh :
      :  -> wf. convergence       0.1000000E-04 e :
      :  Broyden damping          0.4000000       :
      ...................................................

iter    E          dE       RMSdq    gap    omega  full diag
1 -95.7711306 -0.957711E+02  0.247E+00 0.34    0.0  T
2 -96.1628485 -0.391718E+00  0.237E+00 0.50    1.0  T
3 -96.5172788 -0.354430E+00  0.992E-01 3.44    1.0  T
4 -96.5562828 -0.390039E-01  0.681E-01 2.96    1.0  T
5 -96.5696220 -0.133392E-01  0.504E-01 2.08    1.0  T
6 -96.5761702 -0.654818E-02  0.363E-01 1.27    1.0  T

……………………
248 -96.5846674 -0.317981E-02  0.655E-02 0.18    1.0  T
249 -96.5657256  0.189418E-01  0.380E-01 0.08    1.0  T
250 -96.5800780 -0.143524E-01  0.233E-01 0.61    1.0  T

*** convergence criteria cannot be satisfied within 250 iterations ***

这是什么道理？如果我没理解错的话，--namespace orc_job 仅仅是为xtb指定一个命名前缀而已，为什么会影响PC potential的设定？
   --namespace STRING give this xtb(1) run a namespace.
                  All files, even temporary ones, will be named accordingly

@wuzhiyi

作者
Author: ylcong 时间: 2024-1-24 22:58
请问这个问题你解决吗？我也遇到了相同的问题

欢迎光临计算化学公社 (http://bbs.keinsci.com/)