Change of dimension in last kpoint when using non-converging scf

[jonas-pueschel]

Hi all,

I currently have a problem using non converging SCF. I'm using ~1-3 steps of SCF in order to improve an initial guess for an iterative method. The problem here is, that the dimensions of ψ1 and ψ2 sometimes differ, making the method fail. This happens around ~30% of the time even for fixed ψ1.

Here is a minimal working example.

using DFTK
using LinearAlgebra

function  graphene_setup(;Ecut = 15)
    L = 20  # height of the simulation box
    kgrid = [6, 6, 1]


    # Define the geometry and pseudopotential
    a = 4.66  # lattice constant
    a1 = a*[1/2,-sqrt(3)/2, 0]
    a2 = a*[1/2, sqrt(3)/2, 0]
    a3 = L*[0  , 0        , 1]
    lattice = [a1 a2 a3]
    C1 = [1/3,-1/3,0.0]  # in reduced coordinates
    C2 = -C1
    positions = [C1, C2]
    C = ElementPsp(:C; psp=load_psp("hgh/pbe/c-q4"))
    atoms = [C, C]

    # Run SCF
    model = model_PBE(lattice, atoms, positions)
    basis = PlaneWaveBasis(model; Ecut, kgrid)

    return [model, basis]
end

model, basis = graphene_setup(Ecut = 15);

# Convergence we desire in the density
tol = 1e-6

filled_occ = DFTK.filled_occupation(model)
n_spin = model.n_spin_components
n_bands = div(model.n_electrons, n_spin * filled_occ, RoundUp)


ψ1 = [DFTK.random_orbitals(basis, kpt, n_bands) for kpt in basis.kpoints];

#scf messes up dimensions
scfres_start = self_consistent_field(basis; ψ = ψ1 , maxiter = 1);
ψ2 = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ;

# different dims in the last kpoint!
println(length(ψ1[length(ψ1)])) #4248
println(length(ψ2[length(ψ2)])) #sometimes 4248, sometimes 5310

I've consulted with @gkemlin and he came to the following conclusion: In this case, there is an additional occupied state in the last kpoint of ψ2. In scfres_start.occupation, we have the following :

7-element Vector{Vector{Float64}}:
 [2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0]
 [2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0]
 [2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0]
 [2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0]
 [2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0]
 [2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0]
 [2.0, 2.0, 2.0, 0.0, 2.0, 0.0, 0.0] 

Here we can see, that 4th state is 0 and the 5th state is filled. This is because the two eigenvalues in scfres.eigenvalues[7][4:5] are very close, as seen here:

2-element Vector{Float64}:
 -0.015228616549833577
 -0.015228616549801455 

Now for the randomness of this occuring, if you do just one iteration of the SCF, the eigenvalues are sometimes not sorted.

Is this intended behaviour? And if it is, is there any way to avoid this in order to always get a "working" ψ2? As the method I'm using is a Riemannian scheme, it is unfortunately not possible to introduce a very small numerical temperature to allow for fractional occupation number, as this destroys the manifold structure of the problem (at least as far as I know).

Thanks for your help and best regards,
Jonas

[antoine-levitt]

Good catch! With small differences the eigenvalues might not be sorted because they are recomputed after the Rayleigh-Ritz. Does https://github.com/JuliaMolSim/DFTK.jl/pull/964/files fix it?

[gkemlin]

It seems to work yes! I was not sure if sorting directly here would break something elsewhere 😄

[antoine-levitt]

I don't see what it would break? This is a bit of a hack, a better solution would be to not recompute the eigenvalues but get them from the Rayleigh Ritz

[jonas-pueschel]

Good catch! With small differences the eigenvalues might not be sorted because they are recomputed after the Rayleigh-Ritz. Does https://github.com/JuliaMolSim/DFTK.jl/pull/964/files fix it?

Yes, this fixes my problem. Thank you very much :)