search

grimme-lab / xtb

Semiempirical Extended Tight-Binding Program Package
https://xtb-docs.readthedocs.io/
GNU Lesser General Public License v3.0
595 stars 148 forks source link

Error termination for charged molecules using GFN-FF (xTB version 6.3.2) #322

Closed NicolaiRee closed 4 years ago

NicolaiRee commented 4 years ago

xTB error termination when performing calculations on charged molecules using GFN-FF (xTB version 6.3.2).

This works with GFN-X (X=0,1,2)

xTB version:

xtb version 6.3.2 (954f15c) compiled by 'conda@5a45a0871d67' on 2020-07-02

Command-line input:

xtb --gfnff file.xyz --opt --chrg 1 --uhf 0 or xtb --gfnff file.sdf --opt --chrg 1 --uhf 0

(file.xyz and file.sdf are given below)

xTB terminations:

[ERROR] Program stopped due to fatal error
-3- Single point calculation terminated
-2- gfnff_calculator_singlepoint: Force-field method terminated
-1- gfnff_eg: EEQ charge constrain error

However, if I e.g. start by running a GFN-1 calculation using file.sdf as input and then run a GFN-FF calculation also with file.sdf as input, then everything works because the program uses some of the files created by the GFN-1 calculation.

Additionally

file.xyz

13

O          1.46380       -2.28250       -0.03700
C          0.83720       -1.22990       -0.01210
C          1.48320        0.05370        0.23580
C          0.78320        1.22590        0.25700
C         -0.58340        1.21710        0.04530
C         -1.29780       -0.01700       -0.19940
N         -0.49180       -1.19440       -0.21430
H          2.55640        0.03220        0.40160
H          1.29720        2.17180        0.44200
H         -1.15950        2.13480        0.05850
H         -1.83570        0.07420       -1.15500
H         -2.08310       -0.11920        0.56530
H         -0.96940       -2.06670       -0.38760

file.sdf


  xtb     08072014173D
xtb: 6.3.2 (954f15c)
 13 13  0     0  0            999 V2000
    1.4896   -2.2438   -0.0275 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.8475   -1.2058   -0.0075 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.4981    0.0466    0.2360 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7744    1.2240    0.2564 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.5681    1.2354    0.0512 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.3469   -0.0099   -0.2052 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.5125   -1.2225   -0.2193 N   0  0  0  0  0  0  0  0  0  0  0  0
    2.5680    0.0344    0.3998 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.2958    2.1567    0.4406 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.0960    2.1819    0.0742 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8599    0.0606   -1.1755 H   0  0  0  0  0  0  0  0  0  0  0  0
   -2.1168   -0.1374    0.5699 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.9728   -2.1202   -0.3930 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0  0  0  0
  2  3  1  0  0  0  0
  3  4  2  0  0  0  0
  3  8  1  0  0  0  0
  4  5  1  0  0  0  0
  4  9  1  0  0  0  0
  5  6  1  0  0  0  0
  5 10  1  0  0  0  0
  6  7  1  0  0  0  0
  6 11  1  0  0  0  0
  6 12  1  0  0  0  0
  7  2  1  0  0  0  0
  7 13  1  0  0  0  0
M  CHG  1   5   1
M  END
>  <total energy / Eh>
-18.421705869411

>  <gradient norm / Eh/a0>
0.000695317397

>  <total energy / Eh>
-18.421706027730

>  <gradient norm / Eh/a0>
0.000176190486

$$$$

The complete output of the GFN-FF calculation:

      -----------------------------------------------------------
     |                   =====================                   |
     |                           x T B                           |
     |                   =====================                   |
     |                         S. Grimme                         |
     |          Mulliken Center for Theoretical Chemistry        |
     |                    University of Bonn                     |
      -----------------------------------------------------------

   * xtb version 6.3.2 (954f15c) compiled by 'conda@5a45a0871d67' on 2020-07-02

   xtb is free software: you can redistribute it and/or modify it under
   the terms of the GNU Lesser General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   xtb is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   Cite this work as:
   * S. Grimme, C. Bannwarth, P. Shushkov, J. Chem. Theory Comput., 2017,
     13, 1989-2009. DOI: 10.1021/acs.jctc.7b00118
   * C. Bannwarth, S. Ehlert and S. Grimme., J. Chem. Theory Comput., 2019,
     15, 1652-1671. DOI: 10.1021/acs.jctc.8b01176
   * P. Pracht, E. Caldeweyher, S. Ehlert, S. Grimme, ChemRxiv, 2019, preprint.
     DOI: 10.26434/chemrxiv.8326202.v1

   for GFN-FF:
   * S. Spicher and S. Grimme, Angew. Chem. Int. Ed., 2020,
     DOI: 10.1002/anie.202004239

   for DFT-D4:
   * E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017,
     147, 034112. DOI: 10.1063/1.4993215
   * E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher,
     C. Bannwarth and S. Grimme, J. Chem. Phys., 2019, 150, 154122.
     DOI: 10.1063/1.5090222

   for sTDA-xTB:
   * S. Grimme and C. Bannwarth, J. Chem. Phys., 2016, 145, 054103.
     DOI: 10.1063/1.4959605

   in the mass-spec context:
   * V. Asgeirsson, C. Bauer and S. Grimme, Chem. Sci., 2017, 8, 4879.
     DOI: 10.1039/c7sc00601b
   * J. Koopman and S. Grimme, ACS Omega 2019, 4, 12, 15120-15133.
     DOI: 10.1021/acsomega.9b02011

   for metadynamics refer to:
   * S. Grimme, J. Chem. Theory Comput., 2019, 155, 2847-2862
     DOI: 10.1021/acs.jctc.9b00143

   with help from (in alphabetical order)
   C. Bannwarth, F. Bohle, G. Brandenburg, E. Caldeweyher, M. Checinski,
   S. Dohm, S. Ehlert, S. Ehrlich, F. März, H. Neugebauer, J. Pisarek,
   P. Pracht, P. Shushkov, and S. Spicher.

 * started run on 2020/08/07 at 15:04:13.783

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

          program call               : /home/Ree/anaconda/envs/rdkit-env/bin/xtb --gfnff file.sdf --opt --chrg 1 --uhf 0
          coordinate file            : file.sdf
          omp threads                :                     1
          number of atoms            :                    13
          number of electrons        :                    36
          charge                     :                     1
          spin                       :                   0.0
          first test random number   :      0.02158246236620

   ID    Z sym.   atoms
    1    8 O      1
    2    6 C      2-6
    3    7 N      7
    4    1 H      8-13
########################################################################
[WARNING] Please study the warnings concerning your input carefully
-1- prog_main: Charge in sdf/mol input was overwritten
########################################################################

          ==================== Thresholds ====================
          CN  :   150.00000
          rep :   500.00000
          disp:  2500.00000
          HB1 :   250.00000
          HB2 :   450.00000

          Pauling EN used:
          Z : 1  EN :  2.20
          Z : 6  EN :  2.55
          Z : 7  EN :  3.04
          Z : 8  EN :  3.44
          electric field strengths (au): 0.000

           -------------------------------------------------
          |           Force Field Initialization            |
           -------------------------------------------------

          distances ...
          ----------------------------------------
          generating topology and atomic info file ...
          pair mat ...
          computing topology distances matrix with Floyd-Warshall algo ...
          making topology EEQ charges ...
          #fragments for EEQ constrain: 1
          ----------------------------------------
          generating topology and atomic info file ...
          pair mat ...
          computing topology distances matrix with Floyd-Warshall algo ...
          making topology EEQ charges ...
          #fragments for EEQ constrain: 1
          rings ...
          # BATM   93
          # H in HB   6
          maxhb123   240 60 0
          doing iterative Hueckel for 1 subsystem(s) ...

  atom   neighbors  erfCN metchar sp-hybrid imet pi  qest     coordinates
    1  O       1    0.99   0.00         2    0    1  -0.382    2.814936   -4.240167   -0.051967
    2  C       3    2.80   0.00         2    0    1   0.165    1.601543   -2.278632   -0.014173
    3  C       3    2.81   0.00         2    0    1  -0.014    2.830998    0.088061    0.445975
    4  C       3    2.80   0.00         2    0    1  -0.020    1.463404    2.313025    0.484526
    5  C       3    2.80   0.00         2    0    1  -0.024   -1.073553    2.334567    0.096754
    6  C       4    3.49   0.00         3    0    0   0.009   -2.545272   -0.018708   -0.387772
    7  N       3    2.79   0.00        -3    0    1  -0.188   -0.968485   -2.310190   -0.414417
    8  H       1    0.97   0.00         0    0    0   0.072    4.852816    0.065007    0.755512
    9  H       1    0.97   0.00         0    0    0   0.065    2.448707    4.075572    0.832613
   10  H       1    0.97   0.00         0    0    0   0.063   -2.071140    4.123193    0.140218
   11  H       1    0.97   0.01         0    0    0   0.061   -3.514701    0.114517   -2.221373
   12  H       1    0.97   0.01         0    0    0   0.061   -4.000172   -0.259648    1.076955
   13  H       1    0.98   0.00         0    0    0   0.131   -1.838325   -4.006597   -0.742662

          #atoms :   13
          #bonds :   13
          #angl  :   21
          #tors  :   33
          #nmol  :   1
          #optfrag :   1
           -------------------------------------------------
          |                   G F N - F F                   |
          |          A general generic force-field          |
          |                  Version 1.0.1                  |
           -------------------------------------------------

 E+G (total)                   0 d,  0 h,  0 min,  0.000 sec
 distance/D3 list               ...        0 min,  0.000 sec (  0.748%)
 non bonded repulsion           ...        0 min,  0.000 sec (  8.265%)
 dCN                            ...        0 min,  0.000 sec (  8.810%)
 EEQ energy and q               ...        0 min,  0.000 sec ( 13.161%)
 D3                             ...        0 min,  0.000 sec ( 16.304%)
 EEQ gradient                   ...        0 min,  0.000 sec (  3.832%)
 bonds                          ...        0 min,  0.000 sec (  4.383%)
 bend and torsion               ...        0 min,  0.000 sec ( 18.418%)
 bonded ATM                     ...        0 min,  0.000 sec ( 14.311%)
 HB/XB (incl list setup)        ...        0 min,  0.000 sec (  9.166%)

 -0.45621720105099950       0.30717927970483128       -4.2949395042293281E-002  -2.2084787564739111E-002  -2.9970971709208372E-002   7.0511699971714629E-003 -0.29055300120902655       0.10094647190368422        6.4743877434214947E-002   6.0449446533440088E-002   6.6716200341658149E-002   6.6797978786567680E-002  0.16789093187469839
  -6.6613381477509392E-016           1
########################################################################
[ERROR] Program stopped due to fatal error
-3- Single point calculation terminated
-2- gfnff_calculator_singlepoint: Force-field method terminated
-1- gfnff_eg: EEQ charge constrain error
########################################################################
abnormal termination of xtb
ERROR STOP

Error termination. Backtrace:
#0  0x7f56c799da47 in ???
#1  0x7f56c7944ea5 in ???
#2  0x7f56c7940035 in ???
#3  0x7f56c5354b44 in ???
#4  0x7f56c79400ad in ???
    at ../sysdeps/x86_64/elf/start.S:103
awvwgk commented 4 years ago

Can confirm, only happens for sdf and mol input with GFN-FF.

The GFN-FF single point calculation uses its own storage in the force field topology field qfrag to have access to fragment charges:

https://github.com/grimme-lab/xtb/blob/c3f203135a419a5f4c4ee40e442ab120c7ba2e5d/src/gfnff/gfnff_eg.f90#L1288-L1290

The sanity check for the actual total charge constraint will fail when the sum of the resulting charges is not equal to the total system charge ichrg:

https://github.com/grimme-lab/xtb/blob/c3f203135a419a5f4c4ee40e442ab120c7ba2e5d/src/gfnff/gfnff_eg.f90#L555-L558

The charge is correctly stored in the molecular structure in the chrg field as well as in the ichrg global which is used in the GFN-FF setup routines. Apparently the GFN-FF setup routine fails to populate the qfrag field from ichrg in for the sdf and mol file case. The correct initalisation is only present in the default case (xyz, coord, ...):

https://github.com/grimme-lab/xtb/blob/c3f203135a419a5f4c4ee40e442ab120c7ba2e5d/src/gfnff/gfnff_setup.f90#L211-L212

So this is certainly a bug in the GFN-FF setup.