Fix compute_fermi_level for FermiZeroTemperature and collinear spins #928
Conversation
|
I think you can just brutally concatenate everything and get rid of n_spin, no? Also should be able to get rid of some of the rest of that function? I think that function should just look like the bisection in FermiBisection but doing the bisection on the index of the eigenvalues rather than on the eigenvalues themselves, ie find the index i such that excess(all_eigs[i]) == 0 and set the fermi level to all_eigs[i] + all_eigs[i+1]. Doesn't look like there's a integer bisection already coded, so just do it by hand I guess. @gkemlin this came up because we had trouble with SCF for magnetic systems (without temperature), there's no reason to want to have the same number of occupied states in every kpoint, is there? |
mfherbst
changed the title
|
The krange_spin function with loop over model.n_spin_components might also be helpful. For zero temperature there might be hidden assumptions about how many bands are actually converged in the diagonalisation or sth. like that. |
|
We discussed the semantics here, and it seems reasonable to allow a variable number of occupied states per kpoint, but emit a warning. It's helpful when eg something is an insulator but the SCF goes through a phase where it thinks it's a metal. Also doing a metal at zero temperature is relatively sensible too. I'm not sure what the other codes do here. |
|
Abinit says by default |
|
I agree, sounds reasonable and useful to have (referring to the different number of occupied bands per k-Point and spin) |
|
Hi, that's a good question. From my very sparse experience with spins (I think it only concerns what we did with Fe2MnAl), I remember that
However, this is was with positive temperature but I guess that at zero temperature there is still no reason to enforce same occupancy for every kpt and spin so that's something reasonable to have, yes. Still, just to make sure I get it, at zero temperature (e.g. for ferromagnetic insulators), if you look at every kpts with both spin channels taken together, the number of occupied states should be the same everywhere right ? Otherwise I can't see how the Fermi level can en up in a spectral gap. |
Yes, if the thing is an insulator. However, it can possibly happen that during SCF iterations the system gets lost and temporarily believes it's a metal (I think I encountered this issue before). Also, there's no particular reason we shouldn't be able to do metals at zero temperature (it might blow up, but it might not, and we should be able to try). |
src/occupation.jl
Outdated
|
|
||
| all_eigenvalues = sort(vcat(eigenvalues...)) | ||
|
|
||
| # Bissection method to find the index of the eigenvalue such that excess_n_electrons = 0 |
There was a problem hiding this comment.
One s at bisection, and there's too much space between find and the
src/occupation.jl
Outdated
| i_max = length(all_eigenvalues) | ||
|
|
||
| if excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) == 0 | ||
| εF = all_eigenvalues[i_max] |
There was a problem hiding this comment.
We used +1 before but I'm not sure why, this is fine.
src/occupation.jl
Outdated
| εF = maximum(maximum, eigenvalues) + 1 | ||
| εF = mpi_max(εF, basis.comm_kpts) | ||
|
|
||
| all_eigenvalues = sort(vcat(eigenvalues...)) |
There was a problem hiding this comment.
Careful that this will not work with MPI (hopefully the tests should catch this)
There was a problem hiding this comment.
You need to wrap and call MPI.Allgatherv here, so that each process does a bisection on the full array. (You cannot do a Gatherv because excess_n_electrons calls mpi_sum and will hang if not all processes participate.)
|
|
||
| if !allequal(compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation) | ||
| @warn("Not all kpoints have the same number of occupied states, which could mean "* | ||
| "that a metallic system is treated at zero temperature.") |
There was a problem hiding this comment.
We want this to be fine for an insulating spin-polarized system (ie where occupations are different between different spins, but not different kpoints)
…eBarat/DFTK.jl into compute_fermi_level_correction
This PR changes compute_fermi_level to have the two spin components of a same kpoint be treated together in the case of collinear spins. The eigenvalues of the two spin components are merged before searching for the HOMO and LUMO.