zerothi / sisl

Electronic structure Python package for post analysis and large scale tight-binding DFT/NEGF calculations
https://zerothi.github.io/sisl
Mozilla Public License 2.0
173 stars 57 forks source link

sisl.viz.FatBandsPlot error in fatbands_mode = 'line/scatter' #799

Open rreho opened 1 week ago

rreho commented 1 week ago

Describe the bug From issue #790 I created a test with PbTe trying to test if sisl.viz.FatBandsPlot works with fatbands_mode different than area_line I share here a zip folder with input files and databases. You will also find a jupyter notebook : ``pp_pbte_test.ipynb' with comments and questions inside. test.zip

Briefly, I reproduced the bug related to fatbands_mode. Moreover, I tried to visualize the projected (fatbands) band structure starting from the fatbands.nodes['bands_data'] and it seems to be work

ing. It would be great if you can have a look at the notebook to check that I do it properly and that it reproduce the expected results.

I am not sure I interpreted correctly the values stored in data.norm2. Are those normalized to 1? I find the maximum value to be around 0.92

Version details 3.12.2 | packaged by conda-forge | (main, Feb 16 2024, 20:50:58) [GCC 12.3.0] 0.1.dev4752+g1fa01ad

pfebrer commented 1 week ago

Thanks! I will try to look at it during this week.

I think that all norm2 from an orbital should sum to 2 (if the calculation is spin unpolarized)

zerothi commented 1 week ago

Thanks! I will try to look at it during this week.

I think that all norm2 from an orbital should sum to 2 (if the calculation is spin unpolarized)

No, norm2 should also give 1 for eigenstates, it just doesn't take the square root of the inner product. It corresponds to $\langle \psi | \mathbf S | \psi \rangle$.

pfebrer commented 1 week ago

Ok, so both:

norm2.sum("band")
# and
norm2.sum("orb")

should give an array of ones for a spin unpolarized calculation, no?

I recall that in one case summing over bands didn't result into 1. With non-colinear spin maybe?

rreho commented 1 week ago

In my case:

 norm2.sum("band") = all 2
 norm2.sum("orb") =all 1

I also tried: np.sum(np.sum(data.norm2[:,:,:],axis=2),axis=1) = all 72 where 72 is two times the number of orbitals

zerothi commented 1 week ago

norm2 is equivalent to $\langle \psi | \mathbf S | \psi \rangle$, so for $\psi$ being eigenstates it should always be 1 (otherwise there is an error! ;)). Even for non-colinear.

In @rreho case it seems like you are doing a non-colinear calculation, since you have 2 times number of orbitals.

pfebrer commented 1 week ago

Ok, so the problem here is that norm2 is sometimes negative, and therefore can't be used as the size directly. Does a negative norm2 make sense? In that case we just would need to take the absolute value to determine the size I guess. But I don't understand why norm2 should be negative.

By the way, in the example that you passed, setting the number of processors to 4 makes the calculation slower than if you just do it with 1.

pfebrer commented 1 week ago

I have also seen that eigenstate.norm2() of an eigenstate coming from spin-orbit returns an array of complex numbers. Is the imaginary part meaningful or just an artifact?

By the way, I retract what I said about using plots being as fast as doing it yourself. For now, with spin polarized hamiltonians, the spin moments and the ipr are also computed even if you don't need them, which slows the process down.

rreho commented 1 week ago

Regarding norm2, if as explained by @zerothi it represent the projection of the overlap matrix S then it can give negative values, as the inner product between orbitals with different phase can be negative (in a non-orthogonal basis set). I also find the naming norm2 quite confusing, as I would expect a positive definite quantity.

Now, my question would be: should we always plot the absolute value of norm2 to get information on the projection?

zerothi commented 1 week ago

Ok, so the problem here is that norm2 is sometimes negative, and therefore can't be used as the size directly. Does a negative norm2 make sense? In that case we just would need to take the absolute value to determine the size I guess. But I don't understand why norm2 should be negative.

By the way, in the example that you passed, setting the number of processors to 4 makes the calculation slower than if you just do it with 1.

Could you share an example of when this happens? I guess this only happens if one does some kind of state mangling. I don't think it should ever happen because the states originate from an eigenvalue calculation where this is the output values... I am a little inclined to say this might be a bug, but I would really like this to

Something that could change states is the degenerate argument in derivative. In fact, in #797 I will by default disable degenerate arguments for things that shouldn't touch it if not requested. A soon PR will completely remove that from the derivative argument, I don't think it should be used there!

So, does your test files contain stuff like this for direct comparison? I would like to figure that part out!

zerothi commented 1 week ago

I have also seen that eigenstate.norm2() of an eigenstate coming from spin-orbit returns an array of complex numbers. Is the imaginary part meaningful or just an artifact?

The imaginary part is an artifact because it can be used on State objects (which doesn't necessarily needs to be an eigenstate.

rreho commented 1 week ago

@zerothi what would you like to see in the test files precisely? It's not clear to me from the above message. Is it about which state give negative norm2 ?

Moreover, this issue was opened for a bug in fatbands_mode. @pfebrer , were you able to reproduce the bug? Is it related to norm2?

pfebrer commented 1 week ago

Yes, the bug is because we are trying to set the size of scatter points to norm2, which has some negative values, and a negative size doesn't make sense.

@zerothi in the zip file that was shared at the issue description, if you read the Hamiltonian, compute an eigenstate and then its norm2 projected on orbitals you will see that there are negative values.

zerothi commented 1 week ago

Ah, yes, I see. I think these relations are useful:

norm = norm2(projection="orbital")
assert np.allclose(norm.sum(-1), 1)
# I *think* the below will also hold:
assert np.all(np.absolute(norm) <= 1)
pfebrer commented 1 week ago

Ok, so in the fatbands does it make sense to show the absolute value?

What about when summing orbitals? Do we take abs first and then sum the absolute values or sum first and then take the abs value of the sum?

zerothi commented 1 week ago

... The imaginary part is an artifact because it can be used on State objects (which doesn't necessarily needs to be an eigenstate.

Let me correct myself... The above is wrong for non-colinear calculations.

Ok, so in the fatbands does it make sense to show the absolute value?

What about when summing orbitals? Do we take abs first and then sum the absolute values or sum first and then take the abs value of the sum?

  1. When dealing with non-orthogonal basis sets, one can encounter projections which are negative. This relates to the bonding nature of the orbital. The diagonal contributions are always positive, yet, the non-orthogonality can produce finite negative contributions for individual orbitals. So there is nothing weird going on here. It just depends heavily on the basis sets because this is basically a Mulliken analysis.
  2. Generally when doing fatbands, it should correspond to the equivalent of a PDOS calculation for that band (without energy considerations).
  3. Generally for colinear calculations, the projected imaginary part of each orbital will be minimal and should have no physical meaning (other than the basis set choice). So for colinear you'll take the real part of the projection, and discard the imaginary value (but negative values are still meaningful because it tells something about the orbitals physical placement).
  4. And here, you are doing something which is spin-orbit. So then the imaginary part contains information on the spin-projections. What you should do depends on what you want to show? If it is the charge of the projection (the default understanding of fatbands), then you would sum the $\Re(\uparrow + \downarrow)$ contributions, which amounts to something like norm[0, :2].real.sum(1) == band-weight for the first orbital. Again, this can be negative, but may still contain some physical understanding of the bonding nature due to the basis set (same as above). So probably a raised warning for negative values, and printing the absolute value (of only the real part) seems like it could work.