search

grimme-lab / xtb-python

Python API for the extended tight binding program package
https://xtb-python.readthedocs.io
GNU Lesser General Public License v3.0
102 stars 30 forks source link

Different default parameters between xtb python and xtb #22

Closed pschmidtke closed 4 years ago

pschmidtke commented 4 years ago

Describe the bug

Not sue this is really a bug, you github leaves me no choice ;) I ran xtb with default parameters (charge 0) on the input molecule here: out.sdf.zip

xtb out.sdf --charge 0

When I run xtb through python the python bindings using the sample documentation

calc = Calculator(Param.GFN2xTB, atomic_numbers, positions,charge=0)
res = calc.singlepoint()
return(res.get_energy())

I get significantly different energies.

I though that GFN2xTB is also used by default by the xtb binary so I expected similar results. Could you indicate which parameters to use in the python bindings to reproduce results from the command line version?

To Reproduce

CLI: xtb out.sdf --charge 0

Output:

xtb data/docking/out.sdf --charge 0 
      -----------------------------------------------------------      
     |                   =====================                   |     
     |                           x T B                           |     
     |                   =====================                   |     
     |                         S. Grimme                         |     
     |          Mulliken Center for Theoretical Chemistry        |     
     |                    University of Bonn                     |     
      -----------------------------------------------------------      

   * xtb version 6.3.1 (conda-forge) compiled by 'runner@Mac-1389.local' on 

   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/06/11 at 15:52:08.342     

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

          program call               : xtb data/docking/out.sdf --charge 0
          coordinate file            : data/docking/out.sdf
          omp threads                :                    12
          number of atoms            :                    23
          number of electrons        :                    99
          charge                     :                     0
          spin                       :                   0.5
          first test random number   :      0.12034868984915

   ID    Z sym.   atoms
    1    6 C      1-4, 6, 8, 9, 11, 12, 15-21, 23
    2    7 N      5, 7, 10, 13, 14
    3    8 O      22

           -------------------------------------------------
          |                 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                  92          :
          :  # atomic orbitals                  92          :
          :  # shells                           46          :
          :  # electrons                        99          :
          :  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    -54.4329000 -0.544329E+02  0.131E-04    0.02       0.0  T
   2    -54.4328991  0.933692E-06  0.437E-03    0.02       4.8  T
   3    -54.4328996 -0.532692E-06  0.302E-03    0.02       6.9  T
   4    -54.4329000 -0.401829E-06  0.257E-05    0.02     812.9  T
   5    -54.4329000  0.264180E-10  0.365E-05    0.02     570.8  T

   *** convergence criteria satisfied after 5 iterations ***

         #    Occupation            Energy/Eh            Energy/eV
      -------------------------------------------------------------
         1        2.0000           -0.8077093             -21.9789
       ...           ...                  ...                  ...
        44        2.0000           -0.4519193             -12.2974
        45        2.0000           -0.4416039             -12.0167
        46        2.0000           -0.4382805             -11.9262
        47        2.0000           -0.4364431             -11.8762
        48        1.9712           -0.4193192             -11.4103
        49        1.3939           -0.4159596             -11.3188
        50        0.9620           -0.4149744             -11.2920 (HOMO)
        51        0.6726           -0.4143188             -11.2742 (LUMO)
        52        0.0003           -0.4064085             -11.0589
        53        0.0000           -0.4013599             -10.9216
        54        0.0000           -0.3993727             -10.8675
        55                         -0.3978393             -10.8258
       ...                                ...                  ...
        92                          0.5733594              15.6019
      -------------------------------------------------------------
                  HL-Gap            0.0006557 Eh            0.0178 eV
             Fermi-level           -0.4150575 Eh          -11.2943 eV

 SCC (total)                   0 d,  0 h,  0 min,  0.075 sec
 SCC setup                      ...        0 min,  0.000 sec (  0.613%)
 Dispersion                     ...        0 min,  0.001 sec (  1.099%)
 classical contributions        ...        0 min,  0.000 sec (  0.094%)
 integral evaluation            ...        0 min,  0.014 sec ( 18.101%)
 iterations                     ...        0 min,  0.048 sec ( 64.779%)
 molecular gradient             ...        0 min,  0.009 sec ( 11.510%)
 printout                       ...        0 min,  0.003 sec (  3.768%)

         :::::::::::::::::::::::::::::::::::::::::::::::::::::
         ::                     SUMMARY                     ::
         :::::::::::::::::::::::::::::::::::::::::::::::::::::
         :: total energy             -53.817273406757 Eh    ::
         :: gradient norm              0.437211002759 Eh/a0 ::
         :: HOMO-LUMO gap              0.017841402980 eV    ::
         ::.................................................::
         :: SCC energy               -54.432899994162 Eh    ::
         :: -> isotropic ES           -0.008400535937 Eh    ::
         :: -> anisotropic ES          0.008550713598 Eh    ::
         :: -> anisotropic XC          0.042938496197 Eh    ::
         :: -> dispersion             -0.033240383834 Eh    ::
         :: repulsion energy           0.615339696451 Eh    ::
         :: add. restraining           0.000000000000 Eh    ::
         :::::::::::::::::::::::::::::::::::::::::::::::::::::

           -------------------------------------------------
          |                Property Printout                |
           -------------------------------------------------

    * Orbital Energies and Occupations

         #    Occupation            Energy/Eh            Energy/eV
      -------------------------------------------------------------
         1        2.0000           -0.8077093             -21.9789
       ...           ...                  ...                  ...
        38        2.0000           -0.4887331             -13.2991
        39        2.0000           -0.4859830             -13.2243
        40        2.0000           -0.4809836             -13.0882
        41        2.0000           -0.4761952             -12.9579
        42        2.0000           -0.4739852             -12.8978
        43        2.0000           -0.4618020             -12.5663
        44        2.0000           -0.4519193             -12.2974
        45        2.0000           -0.4416039             -12.0167
        46        2.0000           -0.4382805             -11.9262
        47        2.0000           -0.4364431             -11.8762
        48        1.9712           -0.4193192             -11.4103
        49        1.3939           -0.4159596             -11.3188
        50        0.9620           -0.4149744             -11.2920 (HOMO)
        51        0.6726           -0.4143188             -11.2742 (LUMO)
        52        0.0003           -0.4064085             -11.0589
        53        0.0000           -0.4013599             -10.9216
        54        0.0000           -0.3993727             -10.8675
        55                         -0.3978393             -10.8258
        56                         -0.3914812             -10.6527
        57                         -0.3673054              -9.9949
        58                         -0.3528530              -9.6016
        59                         -0.3265086              -8.8848
        60                         -0.3089247              -8.4063
        61                         -0.2935231              -7.9872
       ...                                ...                  ...
        92                          0.5733594              15.6019
      -------------------------------------------------------------
                  HL-Gap            0.0006557 Eh            0.0178 eV
             Fermi-level           -0.4150575 Eh          -11.2943 eV

     #   Z        covCN         q      C6AA      α(0)
     1   6 C      0.977    -0.104    39.190    10.512
     2   6 C      1.961     0.030    30.295     8.889
     3   6 C      1.962     0.038    30.109     8.861
     4   6 C      1.897    -0.072    33.089     9.290
     5   7 N      2.839     0.034    20.639     6.747
     6   6 C      3.111     0.151    24.593     8.124
     7   7 N      1.828    -0.181    25.739     7.518
     8   6 C      1.898     0.136    27.793     8.514
     9   6 C      1.962    -0.061    32.717     9.237
    10   7 N      1.833    -0.153    25.108     7.425
    11   6 C      3.177     0.075    25.926     8.314
    12   6 C      2.838     0.187    24.002     8.029
    13   7 N      0.923    -0.187    24.987     7.223
    14   7 N      1.904    -0.192    25.974     7.556
    15   6 C      3.236     0.157    23.820     7.918
    16   6 C      1.993     0.041    29.998     8.845
    17   6 C      2.984     0.080    26.200     8.396
    18   6 C      1.992    -0.087    33.432     9.338
    19   6 C      1.992    -0.006    31.203     9.021
    20   6 C      1.995    -0.007    31.226     9.024
    21   6 C      2.862     0.176    24.229     8.068
    22   8 O      1.702    -0.038    14.976     5.225
    23   6 C      0.855    -0.018    35.590     9.954

 Mol. C6AA /au·bohr⁶  :      14570.426237
 Mol. C8AA /au·bohr⁸  :     395614.902682
 Mol. α(0) /au        :        192.030328

Wiberg/Mayer (AO) data.
largest (>0.10) Wiberg bond orders for each atom
          total WBO             WBO to atom ...
     1  C   2.224        C    2 1.742    C    3 0.265    C    4 0.122
     2  C   3.774        C    3 1.826    C    1 1.742    C    4 0.155
     3  C   3.647        C    2 1.826    C    4 1.473    C    1 0.265
     4  C   2.974        C    3 1.473    N    5 1.001    C    2 0.155    C    1 0.122
     5  N   3.669        C    6 1.213    C   15 1.178    C    4 1.001
     6  C   3.862        N    7 1.301    N    5 1.213    C   11 1.106
     7  N   3.052        C    8 1.483    C    6 1.301
     8  C   3.717        N   10 1.639    N    7 1.483    N   13 0.420
     9  C   2.909        C   15 1.192    C   17 1.108    C   20 0.128    C   16 0.110
    10  N   3.109        C    8 1.639    C   12 1.124    N   13 0.208
    11  C   3.856        N   14 1.419    C    6 1.106    C   12 1.055
    12  C   3.943        N   13 1.668    N   10 1.124    C   11 1.055
    13  N   2.482        C   12 1.668    C    8 0.420    N   10 0.208
    14  N   3.067        C   11 1.419    C   15 1.313
    15  C   3.905        N   14 1.313    C    9 1.192    N    5 1.178
    16  C   3.746        C   19 1.676    C   17 1.297    C   20 0.339    C   18 0.135    C    9 0.110
    17  C   3.952        C   18 1.358    C   16 1.297    C    9 1.108
    18  C   3.607        C   21 1.363    C   17 1.358    C   19 0.420    C   20 0.136    C   16 0.135
    19  C   3.842        C   16 1.676    C   20 1.627    C   18 0.420
    20  C   3.717        C   19 1.627    C   21 1.295    C   16 0.339    C   18 0.136    C    9 0.128
    21  C   3.869        C   18 1.363    C   20 1.295    O   22 1.016
    22  O   2.621        C   23 1.535    C   21 1.016
    23  C   1.803        O   22 1.535
Topologies differ in total number of bonds
Writing corrected topology to xtbtopo.sdf

molecular dipole:
                 x           y           z       tot (Debye)
 q only:        1.231       0.461       0.259
   full:        1.222      -0.190       0.882       3.861
molecular quadrupole (traceless):
                xx          xy          yy          xz          yz          zz
 q only:      -18.992      50.577       6.724      84.932      38.401      12.269
  q+dip:      -39.411      47.020     -46.218      83.420       9.375      85.629
   full:      -42.080      46.850     -44.870      83.838       7.790      86.949

           -------------------------------------------------
          | TOTAL ENERGY              -53.817273406757 Eh   |
          | GRADIENT NORM               0.437211002759 Eh/α |
          | HOMO-LUMO GAP               0.017841402980 eV   |
           -------------------------------------------------

------------------------------------------------------------------------
 * finished run on 2020/06/11 at 15:52:08.456     
------------------------------------------------------------------------
 total:
 * wall-time:     0 d,  0 h,  0 min,  0.113 sec
 *  cpu-time:     0 d,  0 h,  0 min,  1.114 sec
 * ratio c/w:     9.884 speedup
 SCF:
 * wall-time:     0 d,  0 h,  0 min,  0.075 sec
 *  cpu-time:     0 d,  0 h,  0 min,  0.768 sec
 * ratio c/w:    10.234 speedup

normal termination of xtb
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL

real    0m0.124s
user    0m0.955s
sys     0m0.187s

Python bindings: 

from xtb.interface import Molecule
from xtb.interface import Calculator, Param
from rdkit import Chem
import numpy as np

def getSinglePointEnergy(ctab: str):
    m = Chem.MolFromMolBlock(ctab)
    n_conformers=m.GetNumConformers()
    if(n_conformers>0):
        atomic_numbers=np.array([atom.GetAtomicNum() for atom in m.GetAtoms()])
        conformer=m.GetConformers()[0]  #get the first conformer (we should only have a single one here)
        if(conformer.Is3D()):
            positions=conformer.GetPositions()
            #xtbmol = Molecule(atomic_numbers, positions)
            calc = Calculator(Param.GFN2xTB, atomic_numbers, positions,charge=0)
            res = calc.singlepoint()
            return(res.get_energy())
        else:
            raise Exception("You need to provide a 3D conformation of a molecule")
    else:
        raise Exception("You need to provide a 3D conformation of a molecule")

getSinglePointEnergy(mol)

mol="""
  Mrv1780 06112013133D          
libRbt.so/2013.1/901 2013/11/27
 23 25  0  0  0  0            999 V2000
    4.4805    8.0903   23.7408 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.2758    7.9615   24.6858 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.2357    8.5813   25.8817 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.1873   10.1075   25.9102 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.9462   10.5375   25.2461 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.7523   10.5053   25.8000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.2909   10.1293   26.9964 N   0  0  0  0  0  0  0  0  0  0  0  0
   -1.0323   10.1983   27.2440 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.0204   11.2327   23.1046 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.9366   10.7080   26.3758 N   0  0  0  0  0  0  0  0  0  0  0  0
   -0.1533   11.0667   24.8026 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5522   11.1019   25.1404 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.4072   11.5667   24.2297 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.5736   11.4000   23.7443 N   0  0  0  0  0  0  0  0  0  0  0  0
    1.8411   11.0338   24.0224 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.3520   12.0306   23.6154 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.9675   12.2618   23.6795 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.4509   13.4279   24.2087 C   0  0  0  0  0  0  0  0  0  0  0  0
    6.2183   12.9874   24.1236 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.7030   14.1479   24.6662 C   0  0  0  0  0  0  0  0  0  0  0  0
    4.3331   14.3732   24.6769 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.8034   15.5097   25.2271 O   0  0  0  0  0  0  0  0  0  0  0  0
    4.5788   16.4496   25.9793 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  2  3  1  0  0  0  0
  3  4  1  0  0  0  0
  4  5  1  0  0  0  0
  5  6  1  0  0  0  0
  6  7  2  0  0  0  0
  7  8  1  0  0  0  0
  8 10  2  0  0  0  0
  6 11  1  0  0  0  0
 10 12  1  0  0  0  0
 11 12  2  0  0  0  0
 12 13  1  0  0  0  0
 11 14  1  0  0  0  0
  5 15  1  0  0  0  0
  9 15  1  0  0  0  0
 14 15  2  0  0  0  0
  9 17  1  0  0  0  0
 16 17  2  0  0  0  0
 17 18  1  0  0  0  0
 16 19  1  0  0  0  0
 19 20  2  0  0  0  0
 18 21  2  0  0  0  0
 20 21  1  0  0  0  0
 21 22  1  0  0  0  0
 22 23  1  0  0  0  0
M  END
"""

This python script yields the following output: 

> getSinglePointEnergy(mol)
   1    -44.5513299 -0.445513E+02  0.153E+01    1.28       0.0  T
   2    -45.4179115 -0.866582E+00  0.962E+00    1.78       1.0  T
   3    -45.4024419  0.154696E-01  0.359E+00    1.11       1.0  T
   4    -45.3759017  0.265402E-01  0.196E+00    2.06       1.0  T
   5    -45.4782437 -0.102342E+00  0.971E-01    1.14       1.0  T
   6    -45.5091614 -0.309176E-01  0.289E-01    1.43       1.0  T
   7    -45.5097077 -0.546331E-03  0.154E-01    1.47       1.0  T
   8    -45.5098477 -0.139974E-03  0.469E-02    1.50       1.0  T
   9    -45.5098563 -0.859794E-05  0.143E-02    1.50       1.5  T
  10    -45.5098570 -0.753974E-06  0.503E-03    1.50       4.1  T
  11    -45.5098572 -0.167557E-06  0.166E-03    1.50      12.5  T
  12    -45.5098572 -0.552605E-08  0.530E-04    1.50      39.3  T
  13    -45.5098572 -0.124182E-08  0.175E-04    1.50     119.4  T
     SCC iter.                  ...        0 min,  0.252 sec
     gradient                   ...        0 min,  0.023 sec
-4.667743216484284

Tested on OSX and python bindings in docker container (please provide a mac version for the python bindings on conda-forge.

Expected behaviour

I'd expect the energies to be similar / identical

Additional context

awvwgk commented 4 years ago

If you provide xtb-python with the same input as xtb you will get the same result. Please keep in mind that both xtb and xtb-python are working in atomic units, rdkit is working (at least for the coordinates) in Ångström.

Running the unconverted coordinates through xtb yields the matching result:

TOTAL ENERGY               -4.667743216484 Eh
pschmidtke commented 4 years ago

Ohh that's a noob error on my side :) Thanks, that makes sense then!

awvwgk commented 4 years ago

It is indeed a quite common error to mix-up Bohr and Ångström. I guess putting a note in the documentation at a more prominent location wouldn't hurt here.

pschmidtke commented 4 years ago

maybe some example code coming from rdkit would be helpful as well in the doc

awvwgk commented 4 years ago

I'm not familiar with rdkit, but feel free to contribute an integration with rdkit or a recipe for the docs.

pschmidtke commented 4 years ago

Last question: if I specify an sdf as input xtb seems to convert automagically from Angstroms to au. When I specify as xyz from within xtb-python it expects au to be provided, is that right? If I specify a mol input it would work with Angstroms as coordinates provided as well?

awvwgk commented 4 years ago

Of course not, you get direct write access to the API objects here. The readers implemented in xtb will parse the respective input format and convert all the quantities before creating the object, it is their responsibility to translate an input following a specification to the internal representation of the object. This provenance is not necessary for the API, because the responsibility to take care of the units is shifted to the API user, which is kinda the price to pay for direct access to the API objects.