Example of RealspaceSI giving wrong(??) result

[AleksBL]

RealspaceSI seems to give a wrong result in the 3D case
I am currently testing a custom implementation of the 3D real-space self energy. However, I am encountered a bug(/error?) in the code where the greens function the RealspaceSI class calculates does not agree with my calculation. It also does not correspond to the one calculated by integrating with scipy.
I kind of suspect the broadening given to the electrode and device parts in the sisl-calculation might be different? Do you have any idea of why it might be giving a wrong result? I have checked the absolute sum of the greens function for convergence, and it seems the dk=200 is enough. Else it could be the weights for the Brillouin zone integration that have some obscure error in them?

Right now the evidence I have is that the smallest (tx=1, ty=1) gf integral does not agree with what e.g. scipy gives. (Bottom part of the code-snippet)

Code to reproduce problem

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import sisl
import matplotlib.pyplot as plt
import numpy as np
from AdaptiveIntegrate.surfgf import surface_self_energy, do_integrate
from AdaptiveIntegrate.twodimint_radon import grid_init
from scipy.spatial import Delaunay
# In[] Build Hamiltonian
tx_0,ty_0 = 2,2
tx, ty, tz1, tz2 = 1,1,1,3
h0_0 = sisl.Hamiltonian(sisl.geom.sc(1.0, atoms = sisl.Atom(1,R =1.1)), orthogonal=False)
h0_0.construct([[0.1, 1.11], [(0.0, 1.0), (1.0, 0.1)]])
h0_0 = h0_0.tile(tx_0, 0).tile(ty_0, 1)
it = 0
offset = np.arange(tx_0*ty_0).sum()/(tx_0*ty_0)
for i in range(tx_0):
   for j in range(ty_0):
       inp = it - offset, 1.0
       h0_0[i,j] = inp
       it+=1
h0_0 = h0_0.tile(tz1,2)
h0_s = h0_0.tile(tz2,2)
h0_s.set_nsc((3,3,1))
elec_idx = np.array([0,1,2,3])#np.arange(h0_0.na)
nk   = 200
E    = 5.0 + 0.5j
print('Building')
Test = surface_self_energy(h0_s.copy(),h0_0.copy(),
                           tx,ty,elec_idx,pivot_start=8)
# In[] Custom adaptive integrator
tol, minarea = 1e-4, 1e-7
igN = 6
init_grid = grid_init(igN)
init_grid = init_grid[init_grid[:,1]<0.5+1e-3]
tri = Delaunay(init_grid)
pnts, simpl = tri.points, tri.simplices
gi, err, nsd, feval = do_integrate(Test.gf, pnts, simpl, E, tol, minarea)
gi += gi.T
# In[] integral with sisl
se1  = sisl.physics.RecursiveSI(h0_0, '-C')
sise = sisl.physics.RealSpaceSI(se1, h0_s, (0,1),
                                dk = 200.0, unfold=(tx,ty,1))
g2   = sise.green(E)
print(np.abs(g2).sum())
diff = gi-g2
print(np.abs(diff).sum())
print(np.abs(diff).max())
print(np.abs(gi).sum())
print(np.abs(g2).sum())

# In[] Check inverse gf is built correctly
for i in range(10):
    k = np.random.random(3)
    k[2] = 0.0
    
    ig1 = sise.surface.Sk(k=k)*E - sise.surface.Hk(k=k)
    ig1[0:4,0:4] -= sise.semi.self_energy(E=E, k=k)
    ig1 = ig1.toarray()
    ig2 = Test.gf.igf(E, k) # this function is called directly in the "Test.gf.gfnp" function below
    assert np.allclose(ig1,ig2)
print(np.abs(ig1-ig2).sum())
from scipy.integrate import dblquad
print('Integrating element 0,0 of G')
def f(y,x):
    k = np.array([x, y, 0.0])
    return Test.gf.gfnp(E, k)[0,0].real
resr = dblquad(f,0, 1, 0, 1, epsabs=1e-5, epsrel=1e-5)
def f(y,x):
    k = np.array([x, y, 0.0])
    return Test.gf.gfnp(E, k)[0,0].imag
resi = dblquad(f,0, 1, 0, 1, epsabs=1e-5, epsrel=1e-5)

print((resr[0]+1j*resi[0])/ gi[0,0]) # Out : ~ (1.00431-0.004142373j)
# gi[0,0]
#Out[154]: np.complex128(0.18289848101166165-0.13902685529002906j)
#
#g2[0,0]
#Out[155]: np.complex128(0.21474171140450368-0.4308226281511309j)

Version details
Run the below code and add to bug-report:

import sisl._debug_info as sd
sd.print_debug_info()

sd.print_debug_info()
[sys]
3.10.12 (main, Nov 6 2024, 20:22:13) [GCC 11.4.0]
[sisl]
version 0.15.2.dev4+g0f589e94d
path /home/investigator/sisl/src/sisl
CC /usr/bin/x86_64-linux-gnu-gcc
CFLAGS -O3 -DNDEBUG
C version 11.4.0
FC /usr/bin/gfortran
FFLAGS -O3 -DNDEBUG -O3
FC version 11.4.0
cython build version Cython version 3.0.11
numpy build version 2.1.1
[sisl.definitions]
NPY_NO_DEPRECATED_API NPY_1_20_API_VERSION
CYTHON_NO_PYINIT_EXPORT 1
[sisl.cmake_args]
CMAKE_BUILD_TYPE Release
WITH_FORTRAN ON
F2PY_REPORT_ON_ARRAY_COPY 10
WITH_F2PY_REPORT_COPY OFF
WITH_F2PY_REPORT_EXIT OFF
WITH_COVERAGE OFF
WITH_LINE_DIRECTIVES OFF
WITH_ANNOTATE OFF
WITH_GDB OFF
NO_COMPILATION
[sisl.env]
SISL_LOG_FILE /home/investigator
SISL_LOG_LEVEL INFO
SISL_NUM_PROCS 1
SISL_PAR_CHUNKSIZE 0.1
SISL_TMP /home/investigator/.sisl_tmp
SISL_CONFIGDIR /home/investigator/.config/sisl
SISL_FILES_TESTS /home/investigator/THIS_DIRECTORY_DOES_NOT_EXIST
SISL_SHOW_PROGRESS False
SISL_IO_DEFAULT
SISL_UNIT_SIESTA codata2018
[runtime]
numpy 2.0.2
scipy 1.15.1
xarray 2024.7.0
netCDF4 1.7.1.post2
pandas 2.2.2
matplotlib 3.10.0
dill 0.3.8
pathos not found
No module named 'pathos'
skimage 0.24.0
plotly 5.24.0
ase 3.23.0
pymatgen not found
No module named 'pymatgen'
[install]
pip numpy==2.0.2 scipy==1.15.1 xarray==2024.7.0 netCDF4==1.7.1.post2 pandas==2.2.2 matplotlib==3.10.0 dill==0.3.8 scikit-image==0.24.0 plotly==5.24.0 ase==3.23.0
conda numpy==2.0.2 scipy==1.15.1 xarray==2024.7.0 netCDF4==1.7.1.post2 pandas==2.2.2 matplotlib==3.10.0 dill==0.3.8 scikit-image==0.24.0 plotly==5.24.0 ase==3.23.0