Issue generating wave functions when atoms are in cell boundaries

[decerio]

I am working with skewed cells like the one below.

Cell vectors are: A=[-62.718, -36.210, 0.000], B=[-4.919, 144.840, 0.000], C=[0.000, 0.000, 15.000] Some of the atoms are siting in the left boundary of the cell:

Then I generate the wave-functions in the Gamma point using the sisl.physics.electron.wavefunction method as follows: Htest.nc.zip

import sisl
H = sisl.get_sile('Htest.nc').read_hamiltonian()
grid = sisl.Grid(0.2,sc=H.geometry.sc,geometry=H.geometry,dtype=np.complex128)
eigens = H.eigenstate()
eigens[1359].wavefunction(grid)

(the Htest.nc hamiltonian has a SIESTA pz orbital assigned to it)

The wavefunction is not correct in the boundaries, as it can be seen below:

This can be better seen if you take modulus square and tile the cell two times in the first lattice direction:

There is a missing feature in the left cell boundary: the features corresponding to the atoms in the middle of the ribbon backbone are missing.

This seems to be related to the atoms being in the cell boundary, and the grid not "seeing those atoms", as it can be solved by shifting the geometry to avoid having atoms sitting in the cell boundary:

Now some (red) features appear related to the atoms in the middle of the ribbon, which were missing before. This avoids the wave-function discontinuity when tiling the grid:

Do you think sisl.physics.electron.wavefunction could be changed somehow so that users don't need to worry about this issue?