Ensure contiguous occupations

[abussy]

This PR adds checks to ensure that occupations are contiguous. In various places in the code, it is simply assumed that occupied and empty orbitals are clearly separated (e.g. here or here). This is currently not a safe assumption, as it is not enforced.

This PR adds a check on eigenvalues monotonocity when computing the occupations. Additionally, a helper function providing masks in the form of ranges (useful for taking views into GPU arrays) was added. Finally, the creation of supercells now enforces consistent reordering of the eigenvalues and wave-functions.

Note: range masks are not used in densities.jl, as the compute_density() function is sometime hijacked in other contexts, where the assumption on occupations does not hold (e.g. here)

This PR is a necessary fix for #1187 tests to pass.

[mfherbst]

LGTM, thanks @abussy

[antoine-levitt]

This is wrong when using methfessel paxton smearing

[abussy]

I can see that this could fail with Methfessel Paxton smearing:

DFTK.jl/src/occupation.jl

Lines 221 to 223 in 3551ad3

	n = count(occ_i -> abs(occ_i) > occupation_threshold, occ)
	# Check that all occupied elements are contiguous, otherwise range is wrong
	@assert all(occ[1:n] .> occupation_threshold)

Looking at the rest of the code, it is still assumed that there is a clear distinction between occupied and empty orbitals (see links in initial post). There, the last occupation with abs(occ) > threshold is considered to be the boundary.

I think that it is still healthy to assert that eigenvalues are monotonic down to epsilon, and to enforce it when building supercells. However, the occupied_empty_masks() function can be re-written to reflect the same convention on occupation boundary.

[mfherbst]

Be careful @abussy. I think @antoine-levitt is right (and I unfortunately missed it). The important difference is

 n = count(occ_i -> abs(occ_i) > occupation_threshold, occ)

versus something like

N = [something(findlast(x -> abs(x) > threshold, occk), 0) for occk in occupation]

which has been used so far. Whereas the former fails for Methfessel-Paxton and Marzari-Vanderbilt, the latter will still work.

BTW we should have a test for this kind of edge case. At least not the first time I make this mistake.

[mfherbst]

Edit: In reading your post again, I think you say the same thing, all right then !

Indeed I think it's fine to make the definition:

• occupied: Occupations can be smaller than threshold in magnitude, but are usually larger
• non-occupied: Occupations definitely smaller than threshold in magnitude

[abussy]

Yes, we are in agreement here.

I'll open a new PR with a fix. I'll also try and find a new test that breaks the current state.