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?

[zerothi]

Thanks for this issue! Really should be solved, yes!

[zerothi]

Could you add a complete minimal example, I can't read the basis information from the supplied file?

[decerio]

Sure. See the jupyter notebook attached. It should be enough with that. test_example.zip

[zerothi]

how many times did you tile this NPG? 20? I am trying things out, but a smaller geometry would be ideal, i.e. say 5 times?

Also, after this reduction, which eigenstate would it be?

I tried something like this:

UT = 5
g = g.untile(UT, 1)
... rest of code
eigens[len(g)//2 - 1]

but didn't get something that looked similar?

Which index of the wavefunction would it then be?

[zerothi]

I could reproduce with a minimal graphene lattice, no need for anything, thanks.

[decerio]

Sorry! I guess the problem should be reproducible with any geometry, as far as some of the atoms sit in the cell boundary. I sent you that strange geometry because it is the one that gave the issue, but now I realize it could be any other too. The eigenstate could also be any that has a nonzero value in the orbital sitting in the cell boundary. In my particular example I choose the VB.

[zerothi]

Ok, the problem is because the code translates to the primary unit-cell. I can see that this can be problematic when coordinates are close to the boundary, as they may get shifted to the other side.

Now the coordinates won't be shifted, but the code will inform you if they were to be shifted. Then it should be up to the user to decide what to do. As you say, a rigid shift of everything into the unit-cell should be the easiest and proper solution.

[zerothi]

Could you have a look with the latest commit?

[decerio]

I had at look at it and I have some comments.

1. SISL only warns you if the atoms are outside the cell, but having them in the boundary is also problematic, and SISL doesn't rise the warning in that case.

2. In fact I believe having the atoms strictly outside the cell wasn't problematic before the commit and only having them in the boundary was. My example geometry has the atoms exactly in the left boundary:

That's problematic before and after the commit:

If we translate the geometry 0.1 Å to the left (not even shifting it inside the cell!) so that no atoms lie in the boundary:

In that case, before the commit the wave-function is correctly generated:

But after the commit it is not because the atoms aren't translated to the primary unit-cell and then you just loose part of the wave-function:

That's even more clearly visible if the shift is larger, 5 Å, in this case:

So, are you sure there is any problem in translating the atoms to the primary unit-cell? I think it is even better to have them translated. In fact I believe the problem can be that the atoms in the boundary are not translated to the origin in the primary unit cell (are they?). E.g., if the cell vectors are: v1 = [a1,0,0] v2 = [0,a2,0], and the origin is [0,0,0], then an atom sitting in the boundary at [a1,0,0] is in the next cell, not in the primary unit-cell, and should be translated to [0,0,0].

In any case, I don't know why should this be a problem. I would need to have a close look to the code to learn how is the wavefunction method implemented. I'll try to do it.

[decerio]

A small update on my comment. This "boundary" atoms aren't really in the boundary but very close to it, and still inside the primary unit-cell:

Look at the scale in the x-axis... I wonder how close to the boundary does this start to be problematic. In the new commit you do:

dxyz = geometry.lattice.cell2length(1e-6).sum(0)
dxyz = geometry.move(dxyz).translate2uc(axes=(0, 1, 2)).move(-dxyz).xyz - geometry.xyz

Such small shift is not enough to translate those atoms outside the primary unit-cell. I guess that's why SISL doesn't rise the warning. But still that atomic arrangement is definitely problematic.

[zerothi]

Thanks for your detailed tests! There has always been some obscurity in selecting supercells connections, and this is what is causing this. I'll have to check again and see what can be done, in the end the neighbouring algorithm should be changed to fix this in a simpler manner. AFAIK.