Bug in gradient calculation

[hanimustafa]

Describe the bug
There is a bug in the gradient calculation in XTB. Enabling the debug numerical differentiation code confirms it. Bizarrely, the error seems to be in the Z component of only 2 of the carbons (!!). This bug has been sitting with me for a while, but after procrastinating the bug hunting for too long I decided to file it instead.

To Reproduce
Run XTB against the structure provided. The gradient norm should be around 0.0059248, instead XTB calculates it as 0.0065176. Enable the debug numerical differentiation in order to confirm it. I attach the output with the debug code enabled.

output.log
NEW2-vstrict-valinomycin-xtb.xyz.log

[foxtran]

There is a much smaller structure with the same bug

13
symmetry c1
C       -0.495237121      0.340532816      3.198195224
C        0.114762879     -2.074467184      1.869195224
O        0.709762879      0.098532816      2.696195224
C        1.151762879     -1.279467184      2.696195224
O       -1.178237121     -0.490467184      3.742195224
O       -0.422237121     -1.512467184      0.937195224
C        1.527762879     -1.673467184      4.120195224
H        2.058762879     -1.268467184      2.077195224
H        2.259762879     -2.478467184      4.107195224
H        1.991762879     -0.812467184      4.602195224
H        0.658762879     -1.963467184      4.707195224
H       -0.887516562      1.343243926      3.114009292
H       -0.124377242     -3.098899245      2.113635305

[foxtran]

And then for...

2
symmetry c1
O        0.709762879      0.098532816      2.696195224
C        1.151762879     -1.279467184      2.696195224

Printout:

analytical gradient
    1    -0.087614239593840     0.273150276380796     0.000000000000000
    2     0.087614239593840    -0.273150276380796    -0.000000000000000

numerical gradient
    1    -0.087614242083589     0.273150279284096     0.000000000000000
    2     0.087614238530875    -0.273150274843204    -0.000000004884981

difference gradient
    1     0.000000002489749    -0.000000002903299     0.000000000000000
    2     0.000000001062965    -0.000000001537593     0.000000004884981

[foxtran]

Even water molecule is affected :|

So, Water number 1:

3
Water number 1
O          0.00000        0.00000        0.11779
H          0.00000        0.75545       -0.47116
H          0.00000       -0.75545       -0.47116

And gradients difference:

difference gradient
    1     0.000000000000000     0.000000000444089     0.000000001103382
    2     0.000000000000000    -0.000000000921435    -0.000000004770539
    3     0.000000000000000     0.000000001365524     0.000000000114443

And let's rotate this water to Water number 2:

3
Water number 2
O        0.000000000      0.000000000      0.957897074
H        0.000000000      0.000000000      0.000000000
H        0.928956335      0.000000000      1.191578339

And gradients difference:

difference gradient
    1     0.005063288638830    -0.000000000000000     0.000000000962910
    2    -0.005063287881972     0.000000000000000    -0.000000002126616
    3    -0.000000000312769    -0.000000000000000    -0.000000000612652

The problem happens somewhere in xtb hamiltonian. But it seems to me that Pew variable is not properly constructed by somewhat.

The problematic code where I stopped:

xtb/src/xtb/hamiltonian.f90

Lines 702 to 705 in 0da2327

	g_xyz(:) = g_xyz + 2*HPij*S(jao,iao)*dshpoly/shpoly &
	& + sdqg(:,1,jj,ii)*(2*HPij - 2*Pew(jao, iao) &
	& - Pij*(ves(ish,iat)+ves(jsh,jat)) &
	& + Pij*(vs(iat)+vs(jat)))

P.S.: numbers may slightly be different for you since I've already disabled some terms both for grads and energies.

[foxtran]

GFN0, GFNFF works fine.

[foxtran]

I caught it :)
Will fix at evening :)

[thfroitzheim]

#1208 seems to have fixed most of the problems. Only the rotated water example still seems to be wrong:

analytical gradient (xtb): 
 1   1.8820131034495E-03   7.0544083123213E-17   2.4801291048848E-03
 2  -1.4111403512881E-03  -6.4571316037434E-17   1.4042508759805E-03
 3  -4.7087275216139E-04  -5.9727670857791E-18  -3.8843799808653E-03

numerical gradient (xtb): 
  1   -3.18114866705858E-03    0.00000000000000E+00    2.48025365172532E-03
  2    3.65211642373162E-03    0.00000000000000E+00    1.40417153981430E-03
  3   -4.70703144545509E-04    0.00000000000000E+00   -3.88442519153963E-03

Here it is strangely only the x-component. tblite gets this and the other cases right (the remaining differences are numerical settings and can be converged):

numerical gradient (xtb): 
  1   -3.1812802113483E-03   0.0000000000000E+00   2.4801302338442E-03
  2    3.6521521826287E-03   0.0000000000000E+00   1.4042502386185E-03
  3   -4.7087197128034E-04   0.0000000000000E+00  -3.8843804724627E-03

So I think this is a further small bug. From looking a bit, I don't quite see the problem. The good thing is that it appears to be a very special case and tblite fixes it (so by moving to tblite for the Hamiltonians in 7.0.0 should fix this).