search

grimme-lab / xtb

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

Perplexing Fukui Indices #1041

Closed corinwagen closed 2 months ago

corinwagen commented 3 months ago

Describe the bug Fukui indices are bizarre and seem wrong.

For carbon monoxide, I get these indices:

     #        f(+)     f(-)     f(0)
     1C      -0.466   -0.500   -0.483
     2O      -0.534   -0.500   -0.517

This seems wrong. The SCM Fukui tutorial gives very different values:

Atom   f+     f-     f0     f(2)
C(1)   0.670  0.674  0.672  -0.004
O(2)   0.330  0.326  0.328  0.004
-------------------------------------
Total  1.000  1.000  1.000  -0.000

And negative Fukui indices are at least a bit unexpected (see discussion here).

To Reproduce Steps to reproduce the behaviour:

  1. happens with input co.xyz:
2
title
C 0.000 0.000 0.648
O 0.000 0.000 −0.486
  1. I am running xTB with xtb co.xyz --verbose --gfn2 --vfukui.
  2. here's the output:
      -----------------------------------------------------------
     |                   =====================                   |
     |                           x T B                           |
     |                   =====================                   |
     |                         S. Grimme                         |
     |          Mulliken Center for Theoretical Chemistry        |
     |                    University of Bonn                     |
      -----------------------------------------------------------

   * xtb version 6.6.1 (8d0f1dd) compiled by 'conda@dc2330fa5886' on 2023-08-01

   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:
   * C. Bannwarth, E. Caldeweyher, S. Ehlert, A. Hansen, P. Pracht,
     J. Seibert, S. Spicher, S. Grimme, WIREs Comput. Mol. Sci., 2020, 11,
     e01493. DOI: 10.1002/wcms.1493

   for GFN2-xTB:
   * C. Bannwarth, S. Ehlert and S. Grimme., J. Chem. Theory Comput., 2019,
     15, 1652-1671. DOI: 10.1021/acs.jctc.8b01176
   for GFN1-xTB:
   * S. Grimme, C. Bannwarth, P. Shushkov, J. Chem. Theory Comput., 2017,
     13, 1989-2009. DOI: 10.1021/acs.jctc.7b00118
   for GFN0-xTB:
   * 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, 59, 15665-15673.
     DOI: 10.1002/anie.202004239

   for ALPB and GBSA implicit solvation:
   * S. Ehlert, M. Stahn, S. Spicher, S. Grimme, J. Chem. Theory Comput.,
     2021, 17, 4250-4261. DOI: 10.1021/acs.jctc.1c00471

   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
   * E. Caldeweyher, J.-M. Mewes, S. Ehlert and S. Grimme, Phys. Chem. Chem. Phys.
     2020, 22, 8499-8512. DOI: 10.1039/D0CP00502A

   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

   for SPH calculations refer to:
   * S. Spicher and S. Grimme, J. Chem. Theory Comput., 2021, 17, 1701-1714
     DOI: 10.1021/acs.jctc.0c01306

   with help from (in alphabetical order)
   P. Atkinson, C. Bannwarth, F. Bohle, G. Brandenburg, E. Caldeweyher
   M. Checinski, S. Dohm, S. Ehlert, S. Ehrlich, I. Gerasimov, C. Hölzer
   A. Katbashev, J. Koopman, C. Lavigne, S. Lehtola, F. März, M. Müller,
   F. Musil, H. Neugebauer, J. Pisarek, C. Plett, P. Pracht, F. Pultar,
   J. Seibert, P. Shushkov, S. Spicher, M. Stahn, M. Steiner, T. Strunk,
   J. Stückrath, T. Rose, and J. Unsleber

 * started run on 2024/05/22 at 20:03:39.061

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

          program call               : xtb co.xyz --verbose --gfn2 --vfukui
          hostname                   : 04cf0da4d3fd
          coordinate file            : co.xyz
          xtbhome directory          : /root
          path for xtb               : /root
          xcontrol input file        : co.xyz
          omp threads                :                     4

$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 CONSTRAINTS & SCANS: DEBUG SECTION
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   ID    Z sym.   atoms
    1    6 C      1
    2    8 O      2
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

           -------------------------------------------------
          |                 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                   8          :
          :  # atomic orbitals                   8          :
          :  # shells                            4          :
          :  # electrons                        10          :
          :  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     -6.1875466 -0.618755E+01  0.266E-05    5.17       0.0  T
   2     -6.1875466 -0.611067E-12  0.237E-05    5.17    2978.6  T
   3     -6.1875466 -0.691891E-12  0.847E-06    5.17    8351.0  T

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

         #    Occupation            Energy/Eh            Energy/eV
      -------------------------------------------------------------
         1        2.0000           -0.7968127             -21.6824
         2        2.0000           -0.6561357             -17.8544
         3        2.0000           -0.6134264             -16.6922
         4        2.0000           -0.6134264             -16.6922
         5        2.0000           -0.4753299             -12.9344 (HOMO)
         6                         -0.2853273              -7.7642 (LUMO)
         7                         -0.2853273              -7.7642
         8                          0.7745546              21.0767
      -------------------------------------------------------------
                  HL-Gap            0.1900026 Eh            5.1702 eV
             Fermi-level           -0.3803286 Eh          -10.3493 eV

 SCC (total)                   0 d,  0 h,  0 min,  0.006 sec
 SCC setup                      ...        0 min,  0.000 sec (  4.015%)
 Dispersion                     ...        0 min,  0.000 sec (  6.281%)
 classical contributions        ...        0 min,  0.002 sec ( 27.765%)
 integral evaluation            ...        0 min,  0.000 sec (  6.262%)
 iterations                     ...        0 min,  0.003 sec ( 51.878%)
 molecular gradient             ...        0 min,  0.000 sec (  2.326%)
 printout                       ...        0 min,  0.000 sec (  1.213%)

         :::::::::::::::::::::::::::::::::::::::::::::::::::::
         ::                     SUMMARY                     ::
         :::::::::::::::::::::::::::::::::::::::::::::::::::::
         :: total energy              -6.121749741648 Eh    ::
         :: gradient norm              0.378677359967 Eh/a0 ::
         :: HOMO-LUMO gap              5.170234276941 eV    ::
         ::.................................................::
         :: HOMO orbital eigv.       -12.934384780631 eV    ::
         :: LUMO orbital eigv.        -7.764150503690 eV    ::
         ::.................................................::
         :: SCC energy                -6.187546588501 Eh    ::
         :: -> isotropic ES           -0.001814980214 Eh    ::
         :: -> anisotropic ES         -0.000083882471 Eh    ::
         :: -> anisotropic XC          0.002323969703 Eh    ::
         :: -> dispersion             -0.000259847058 Eh    ::
         :: repulsion energy           0.065796846852 Eh    ::
         :: add. restraining           0.000000000000 Eh    ::
         :: total charge              -0.000000000000 e     ::
         ::.................................................::
         :: atomisation energy         0.733077473384 Eh    ::
         :::::::::::::::::::::::::::::::::::::::::::::::::::::

Fukui index Calculation

 iter      E             dE          RMSdq      gap      omega  full diag
   1     -6.2209354 -0.622094E+01  0.779E+00    0.00       0.0  T
   2     -6.2371169 -0.161815E-01  0.537E+00    0.00       1.0  T
   3     -6.2397796 -0.266270E-02  0.259E+00    0.00       1.0  T
   4     -6.2400531 -0.273505E-03  0.431E-01    0.00       1.0  T
   5     -6.2398889  0.164156E-03  0.481E-01    0.00       1.0  T
   6     -6.2402367 -0.347805E-03  0.754E-02    0.00       1.0  T
   7     -6.2402415 -0.474003E-05  0.476E-02    0.00       1.5  T
   8     -6.2402446 -0.309682E-05  0.226E-03    0.00      31.3  T
   9     -6.2402446 -0.604950E-08  0.893E-04    0.00      79.1  T
  10     -6.2402446 -0.115353E-08  0.230E-05    0.00    3068.2  T
     SCC iter.                  ...        0 min,  0.000 sec
     gradient                   ...        0 min,  0.000 sec

 iter      E             dE          RMSdq      gap      omega  full diag
   1     -5.4596427 -0.545964E+01  0.686E+00    5.17       0.0  T
   2     -5.4555811  0.406161E-02  0.737E+00    5.08       1.0  T
   3     -5.5008893 -0.453082E-01  0.377E+00    4.91       1.0  T
   4     -5.4992029  0.168644E-02  0.138E+00    4.68       1.0  T
   5     -5.5014844 -0.228157E-02  0.419E-01    4.63       1.0  T
   6     -5.5018995 -0.415031E-03  0.271E-02    4.65       2.6  T
   7     -5.5019009 -0.144221E-05  0.776E-03    4.65       9.1  T
   8     -5.5019009 -0.125279E-07  0.169E-03    4.64      41.9  T
   9     -5.5019009 -0.471223E-08  0.485E-04    4.64     145.8  T
  10     -5.5019009 -0.104995E-09  0.399E-05    4.64    1772.9  T
     SCC iter.                  ...        0 min,  0.000 sec
     gradient                   ...        0 min,  0.000 sec

     #        f(+)     f(-)     f(0)
     1C      -0.466   -0.500   -0.483
     2O      -0.534   -0.500   -0.517
           -------------------------------------------------
          |                Property Printout                |
           -------------------------------------------------

    * Orbital Energies and Occupations

         #    Occupation            Energy/Eh            Energy/eV
      -------------------------------------------------------------
         1        2.0000           -0.7968127             -21.6824
         2        2.0000           -0.6561357             -17.8544
         3        2.0000           -0.6134264             -16.6922
         4        2.0000           -0.6134264             -16.6922
         5        2.0000           -0.4753299             -12.9344 (HOMO)
         6                         -0.2853273              -7.7642 (LUMO)
         7                         -0.2853273              -7.7642
         8                          0.7745546              21.0767
      -------------------------------------------------------------
                  HL-Gap            0.1900026 Eh            5.1702 eV
             Fermi-level           -0.3803286 Eh          -10.3493 eV

     #   Z          covCN         q      C6AA      α(0)
     1   6 C        0.856     0.050    33.618     9.675
     2   8 O        0.856    -0.050    15.799     5.329

 Mol. C6AA /au·bohr⁶  :         94.566487
 Mol. C8AA /au·bohr⁸  :       2381.891212
 Mol. α(0) /au        :         15.004091

Wiberg/Mayer (AO) data.
largest (>0.10) Wiberg bond orders for each atom

 ---------------------------------------------------------------------------
     #   Z sym  total        # sym  WBO       # sym  WBO       # sym  WBO
 ---------------------------------------------------------------------------
     1   6 C    2.610 --     2 O    2.610
     2   8 O    2.610 --     1 C    2.610
 ---------------------------------------------------------------------------

Topologies differ in total number of bonds
Writing topology from bond orders to xtbtopo.mol

molecular dipole:
                 x           y           z       tot (Debye)
 q only:        0.000       0.000       0.108
   full:       -0.000       0.000      -0.237       0.601
molecular quadrupole (traceless):
                xx          xy          yy          xz          yz          zz
 q only:       -0.016       0.000      -0.016       0.000       0.000       0.033
  q+dip:        0.405       0.000       0.405      -0.000       0.000      -0.809
   full:        1.144      -0.000       1.144      -0.000       0.000      -2.288

           -------------------------------------------------
          | TOTAL ENERGY               -6.121749741648 Eh   |
          | GRADIENT NORM               0.378677359967 Eh/α |
          | HOMO-LUMO GAP               5.170234276941 eV   |
           -------------------------------------------------

------------------------------------------------------------------------
 unit  open   action     filename
   10 false : read       co.xyz
   10 false : read       co.xyz
   10 false : read       co.xyz
   10 false : read       xtbrestart
   10 false : replaced   charges
   10 false : replaced   wbo
   10 false : replaced   xtbtopo.mol
   10 false : replaced   xtbrestart

------------------------------------------------------------------------
 * finished run on 2024/05/22 at 20:03:39.086
------------------------------------------------------------------------
 total:
 * wall-time:     0 d,  0 h,  0 min,  0.024 sec
 *  cpu-time:     0 d,  0 h,  0 min,  0.157 sec
 * ratio c/w:     6.428 speedup
 SCF:
 * wall-time:     0 d,  0 h,  0 min,  0.006 sec
 *  cpu-time:     0 d,  0 h,  0 min,  0.038 sec
 * ratio c/w:     6.634 speedup

normal termination of xtb

Please provide all input and output file such that we confirm your report.

Expected behaviour I expect Fukui indices more like those from SCM (vide supra).

Additional context xTB 6.6.1 installed from conda into a Docker image.

corinwagen commented 3 months ago

Note that similar issues have been seen before: #794

corinwagen commented 3 months ago

Similarly when I compute f(-) for fluorobenzene, I get ortho = -0.049, meta = -0.045, and para = -0.081.

The literature values are ortho = 0.047, meta = 0.026, and para = 0.093. Granted, these data are computed using HF/3-21G so it's not surprising that the numerical values are different, but the total flip in sign makes me think there's something pretty substantially wrong with the xTB values.

haneug commented 3 months ago

Indeed, the sign switch I introduced in #794 seems to be wrong since we are using Mulliken charges and not Mulliken populations. Therefore, the sign of the fukui functions is switched! I will submit a PR to fix this. Thank you for bringing this to our attention.

corinwagen commented 3 months ago

any time. the values are correct otherwise?

haneug commented 3 months ago

Otherwise they should be correct. They are calculated according to equations 2a, 2b, and 2c of this paper (also shown at the end of the abstract) where the term Mulliken charges can be misleading in my opinion since I believe what they are actually using is what I understand as Mulliken population. The difference with SCM I would attribute mostly to them using the Hirshfeld partitioning (probably a better idea) and secondly on the STOs they use instead of GTOs, although this difference I assume should be very small.