From 7f142d2241a80c7418a1b3e11fdccb7fc363ca00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Mon, 11 Dec 2023 09:22:14 +0100 Subject: [PATCH 01/22] correction for collinear spins --- src/occupation.jl | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index fbfb32ca92..4beb708dc7 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -146,23 +146,33 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT "occupation $filled_occ. Typically this indicates that you need to put " * "a temperature or switch to a calculation with collinear spin polarization.") end - n_fill = div(n_electrons, n_spin * filled_occ, RoundUp) + n_fill = div(n_electrons, filled_occ, RoundUp) + + # For a given kpoint, all the eigenvalues corresponding to different spins must be treated together. + kpt_coordinates = unique([basis.kpoints[i].coordinate for i in 1:length(basis.kpoints)]) + if length(kpt_coordinates) < length(basis.kpoints) + merged_spin_eigenvalues = [sort(vcat(eigenvalues[findall(j->basis.kpoints[j].coordinate==coordinate, 1:length(basis.kpoints))]...)) + for coordinate in kpt_coordinates] + else + merged_spin_eigenvalues = eigenvalues + end # For zero temperature, two cases arise: either there are as many bands # as electrons, in which case we set εF to the highest energy level # reached, or there are unoccupied conduction bands and we take # εF as the midpoint between valence and conduction bands. + if n_fill == length(eigenvalues[1]) - εF = maximum(maximum, eigenvalues) + 1 + εF = maximum(maximum, merged_spin_eigenvalues) + 1 εF = mpi_max(εF, basis.comm_kpts) else # highest occupied energy level - HOMO = maximum([εk[n_fill] for εk in eigenvalues]) + HOMO = maximum([εk[n_fill] for εk in merged_spin_eigenvalues]) HOMO = mpi_max(HOMO, basis.comm_kpts) # lowest unoccupied energy level, be careful that not all k-points # might have at least n_fill+1 energy levels so we have to take care # of that by specifying init to minimum - LUMO = minimum(minimum.([εk[n_fill+1:end] for εk in eigenvalues]; init=T(Inf))) + LUMO = minimum(minimum.([εk[n_fill+1:end] for εk in merged_spin_eigenvalues]; init=T(Inf))) LUMO = mpi_min(LUMO, basis.comm_kpts) εF = (HOMO + LUMO) / 2 end From e6c9f2a4b83e86e66accb2e7a6e8f6be402a7dc7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 12 Dec 2023 11:36:47 +0100 Subject: [PATCH 02/22] treat all kpoints together --- src/occupation.jl | 49 +++++++++++++++++++++-------------------------- 1 file changed, 22 insertions(+), 27 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 4beb708dc7..9a338fee0e 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -137,44 +137,39 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT temperature, smearing, tol_n_elec) where {T} filled_occ = filled_occupation(basis.model) n_electrons = basis.model.n_electrons - n_spin = basis.model.n_spin_components @assert iszero(temperature) # Sanity check that we can indeed fill the appropriate number of states - if n_electrons % (n_spin * filled_occ) != 0 + if n_electrons % (filled_occ) != 0 error("$n_electrons electrons cannot be attained by filling states with " * "occupation $filled_occ. Typically this indicates that you need to put " * "a temperature or switch to a calculation with collinear spin polarization.") end - n_fill = div(n_electrons, filled_occ, RoundUp) - # For a given kpoint, all the eigenvalues corresponding to different spins must be treated together. - kpt_coordinates = unique([basis.kpoints[i].coordinate for i in 1:length(basis.kpoints)]) - if length(kpt_coordinates) < length(basis.kpoints) - merged_spin_eigenvalues = [sort(vcat(eigenvalues[findall(j->basis.kpoints[j].coordinate==coordinate, 1:length(basis.kpoints))]...)) - for coordinate in kpt_coordinates] + all_eigenvalues = sort(vcat(eigenvalues...)) + + i_min = 1, i_max = length(all_eigenvalues) + if excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) == 0 + εF = all_eigenvalues[i_max] + elseif excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) * + excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_min]; temperature, smearing) < 0 + while i_max - i_min > 1 + i = div(i_min+i_max, 2) + εF = all_eigenvalues[i] + if excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) <= 0 + i_min = i + else + i_max = i + end + end + εF = 1/2*(all_eigenvalues[i_min]+all_eigenvalues[i_max]) else - merged_spin_eigenvalues = eigenvalues + error("") end - # For zero temperature, two cases arise: either there are as many bands - # as electrons, in which case we set εF to the highest energy level - # reached, or there are unoccupied conduction bands and we take - # εF as the midpoint between valence and conduction bands. - - if n_fill == length(eigenvalues[1]) - εF = maximum(maximum, merged_spin_eigenvalues) + 1 - εF = mpi_max(εF, basis.comm_kpts) - else - # highest occupied energy level - HOMO = maximum([εk[n_fill] for εk in merged_spin_eigenvalues]) - HOMO = mpi_max(HOMO, basis.comm_kpts) - # lowest unoccupied energy level, be careful that not all k-points - # might have at least n_fill+1 energy levels so we have to take care - # of that by specifying init to minimum - LUMO = minimum(minimum.([εk[n_fill+1:end] for εk in merged_spin_eigenvalues]; init=T(Inf))) - LUMO = mpi_min(LUMO, basis.comm_kpts) - εF = (HOMO + LUMO) / 2 + if not allequal(compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation) + @warn("All kpoints don't have the same occupations, this could indicate "* + "that a metalic system is being treated with zero temperature.") end excess(εF) = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) From 16d281816cb84dc62d0ef847c28e801513331fdd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 12 Dec 2023 11:40:30 +0100 Subject: [PATCH 03/22] wip --- src/occupation.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/occupation.jl b/src/occupation.jl index 9a338fee0e..fd564b1711 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -148,7 +148,8 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT all_eigenvalues = sort(vcat(eigenvalues...)) - i_min = 1, i_max = length(all_eigenvalues) + i_min = 1 + i_max = length(all_eigenvalues) if excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) == 0 εF = all_eigenvalues[i_max] elseif excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) * From 862eaa00376be9177f8e9b6391f8bd86f781dc8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 12 Dec 2023 13:56:05 +0100 Subject: [PATCH 04/22] fix --- src/occupation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/occupation.jl b/src/occupation.jl index fd564b1711..b8d8ae97a1 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -168,7 +168,7 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT error("") end - if not allequal(compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation) + if !allequal(compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation) @warn("All kpoints don't have the same occupations, this could indicate "* "that a metalic system is being treated with zero temperature.") end From 3fccd531b77d830d65eed205d99f1653f5bc3e15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 12 Dec 2023 15:36:36 +0100 Subject: [PATCH 05/22] to debug --- src/occupation.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index b8d8ae97a1..3d18d3fd98 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -163,14 +163,15 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT i_max = i end end + println((i_min, i_max)) εF = 1/2*(all_eigenvalues[i_min]+all_eigenvalues[i_max]) else - error("") + error("Not enough eigenvalues computed") end if !allequal(compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation) @warn("All kpoints don't have the same occupations, this could indicate "* - "that a metalic system is being treated with zero temperature.") + "that a metallic system is being treated with zero temperature.") end excess(εF) = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) From f6deab01f429bf8cca19023836974dbffea2421f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 12 Dec 2023 15:43:46 +0100 Subject: [PATCH 06/22] to debug --- src/occupation.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/occupation.jl b/src/occupation.jl index 3d18d3fd98..a507e4ef1a 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -155,6 +155,9 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT elseif excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) * excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_min]; temperature, smearing) < 0 while i_max - i_min > 1 + println((i_min, i_max)) + println((all_eigenvalues[i_min], all_eigenvalues[i_max])) + println((excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_min]; temperature, smearing), excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing))) i = div(i_min+i_max, 2) εF = all_eigenvalues[i] if excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) <= 0 @@ -163,7 +166,6 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT i_max = i end end - println((i_min, i_max)) εF = 1/2*(all_eigenvalues[i_min]+all_eigenvalues[i_max]) else error("Not enough eigenvalues computed") From 20824f2bb470fd9521ae7b9297e3fb7f349d3e6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 12 Dec 2023 16:32:12 +0100 Subject: [PATCH 07/22] remove debug --- src/occupation.jl | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index a507e4ef1a..537d5ef662 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -148,16 +148,14 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT all_eigenvalues = sort(vcat(eigenvalues...)) + # Bissection method to find the index of the eigenvalue such that excess_n_electrons = 0 i_min = 1 i_max = length(all_eigenvalues) + if excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) == 0 εF = all_eigenvalues[i_max] - elseif excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) * - excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_min]; temperature, smearing) < 0 + else while i_max - i_min > 1 - println((i_min, i_max)) - println((all_eigenvalues[i_min], all_eigenvalues[i_max])) - println((excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_min]; temperature, smearing), excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing))) i = div(i_min+i_max, 2) εF = all_eigenvalues[i] if excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) <= 0 @@ -167,13 +165,11 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT end end εF = 1/2*(all_eigenvalues[i_min]+all_eigenvalues[i_max]) - else - error("Not enough eigenvalues computed") end if !allequal(compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation) - @warn("All kpoints don't have the same occupations, this could indicate "* - "that a metallic system is being treated with zero temperature.") + @warn("Not all kpoints have the same number of occupied states, which could mean "* + "that a metallic system is treated at zero temperature.") end excess(εF) = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) From 4b158db276b0427b27b3af6427ab6cb2294b95f7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 19 Dec 2023 09:41:03 +0100 Subject: [PATCH 08/22] Warning correction and typo --- src/occupation.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 537d5ef662..05a778c24a 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -1,4 +1,5 @@ import Roots +#import MPI # Functions for finding the Fermi level and occupation numbers for bands # The goal is to find εF so that @@ -148,7 +149,7 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT all_eigenvalues = sort(vcat(eigenvalues...)) - # Bissection method to find the index of the eigenvalue such that excess_n_electrons = 0 + # Bisection method to find the index of the eigenvalue such that excess_n_electrons = 0 i_min = 1 i_max = length(all_eigenvalues) @@ -167,7 +168,9 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT εF = 1/2*(all_eigenvalues[i_min]+all_eigenvalues[i_max]) end - if !allequal(compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation) + occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation + merged_spin_occupations = sum([occ[krange_spin(basis, i)] for i in 1:basis.model.n_spin_components]) + if !allequal(merged_spin_occupations) @warn("Not all kpoints have the same number of occupied states, which could mean "* "that a metallic system is treated at zero temperature.") end From 6fe210aeaca3fd426bda6a02f9e867b3c519bed7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Wed, 20 Dec 2023 16:35:23 +0100 Subject: [PATCH 09/22] correction for MPI --- src/occupation.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 05a778c24a..f3ad5770f9 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -147,9 +147,13 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT "a temperature or switch to a calculation with collinear spin polarization.") end - all_eigenvalues = sort(vcat(eigenvalues...)) - + # Gathering the eigenvalues distributed over the kpoints in MPI + n_bands = length(eigenvalues[1]) # assuming that the same number of bands are computed for each kpoint + counts = [ Int32(n_bands*length(basis.krange_allprocs[rank])) for rank in 1:MPI.Comm_size(basis.comm_kpts)] + all_eigenvalues = MPI.Allgatherv(vcat(eigenvalues...), counts, basis.comm_kpts) + # Bisection method to find the index of the eigenvalue such that excess_n_electrons = 0 + all_eigenvalues = sort(all_eigenvalues) i_min = 1 i_max = length(all_eigenvalues) From 004eadd5a012fe32c1e022175c90723d9eec8698 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Thu, 16 May 2024 15:01:01 +0200 Subject: [PATCH 10/22] Bissection correction --- src/occupation.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index f3ad5770f9..8f63025b8e 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -150,26 +150,29 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT # Gathering the eigenvalues distributed over the kpoints in MPI n_bands = length(eigenvalues[1]) # assuming that the same number of bands are computed for each kpoint counts = [ Int32(n_bands*length(basis.krange_allprocs[rank])) for rank in 1:MPI.Comm_size(basis.comm_kpts)] - all_eigenvalues = MPI.Allgatherv(vcat(eigenvalues...), counts, basis.comm_kpts) - + all_eigenvalues = MPI.Allgatherv(reduce(vcat, eigenvalues), counts, basis.comm_kpts) + # mapreduce to sort better # Bisection method to find the index of the eigenvalue such that excess_n_electrons = 0 all_eigenvalues = sort(all_eigenvalues) i_min = 1 i_max = length(all_eigenvalues) - if excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) == 0 εF = all_eigenvalues[i_max] else while i_max - i_min > 1 i = div(i_min+i_max, 2) εF = all_eigenvalues[i] - if excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) <= 0 + excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) + if excess < 0 i_min = i + elseif excess > 0 + i_max = i else + i_min = i i_max = i end end - εF = 1/2*(all_eigenvalues[i_min]+all_eigenvalues[i_max]) + εF = all_eigenvalues[i_min]+1/2*(all_eigenvalues[i_max]-all_eigenvalues[i_min]) end occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation @@ -178,9 +181,8 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT @warn("Not all kpoints have the same number of occupied states, which could mean "* "that a metallic system is treated at zero temperature.") end - - excess(εF) = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - if abs(excess(εF)) > tol_n_elec + excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) + if abs(excess) > tol_n_elec error("Unable to find non-fractional occupations that have the " * "correct number of electrons. You should add a temperature.") end From d34d51a2349297aff6c9a85e4ead2c1e1bfd167a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Fri, 24 May 2024 10:18:38 +0200 Subject: [PATCH 11/22] Working bissection - Allgather TODO --- src/occupation.jl | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 8f63025b8e..56d57fbaf7 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -148,31 +148,35 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT end # Gathering the eigenvalues distributed over the kpoints in MPI - n_bands = length(eigenvalues[1]) # assuming that the same number of bands are computed for each kpoint - counts = [ Int32(n_bands*length(basis.krange_allprocs[rank])) for rank in 1:MPI.Comm_size(basis.comm_kpts)] - all_eigenvalues = MPI.Allgatherv(reduce(vcat, eigenvalues), counts, basis.comm_kpts) - # mapreduce to sort better + #n_bands = length(eigenvalues[1]) # assuming that the same number of bands are computed for each kpoint + #counts = [ Int32(n_bands*length(basis.krange_allprocs[rank])) for rank in 1:MPI.Comm_size(basis.comm_kpts)] + #all_eigenvalues = similar(reduce(vcat, eigenvalues)) + #all_eigenvalues_vbuf = MPI.VBuffer(all_eigenvalues, counts) + #MPI.Allgatherv!(reduce(vcat, eigenvalues), all_eigenvalues_vbuf, basis.comm_kpts) + #all_eigenvalues = MPI.Allgatherv(reduce(vcat, eigenvalues), counts, basis.comm_kpts) + + all_eigenvalues = reduce(vcat, eigenvalues) + # Bisection method to find the index of the eigenvalue such that excess_n_electrons = 0 - all_eigenvalues = sort(all_eigenvalues) + unique!(sort!((all_eigenvalues))) i_min = 1 i_max = length(all_eigenvalues) - if excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max]; temperature, smearing) == 0 - εF = all_eigenvalues[i_max] + if abs(excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max] + T(1); temperature, smearing)) < tol_n_elec + εF = all_eigenvalues[i_max] + T(1) else - while i_max - i_min > 1 + while i_max - i_min > 0 i = div(i_min+i_max, 2) - εF = all_eigenvalues[i] + εF = all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - if excess < 0 + if excess < -tol_n_elec i_min = i - elseif excess > 0 + elseif excess > tol_n_elec i_max = i else i_min = i i_max = i end end - εF = all_eigenvalues[i_min]+1/2*(all_eigenvalues[i_max]-all_eigenvalues[i_min]) end occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation @@ -186,6 +190,5 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT error("Unable to find non-fractional occupations that have the " * "correct number of electrons. You should add a temperature.") end - εF end From 921ce2f42554df33be6d3933fbd1eeeb5f229a29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Fri, 14 Jun 2024 15:20:10 +0200 Subject: [PATCH 12/22] Avoid infinite loop --- src/occupation.jl | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 56d57fbaf7..c16045432c 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -161,13 +161,28 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT unique!(sort!((all_eigenvalues))) i_min = 1 i_max = length(all_eigenvalues) - if abs(excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max] + T(1); temperature, smearing)) < tol_n_elec + excess_min = excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_min] + T(1/2) + *(all_eigenvalues[i_min+1]-all_eigenvalues[i_min]); + temperature, smearing) + excess_max = excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max] + T(1); + temperature, smearing) + + @show length(all_eigenvalues) + @show excess_min, excess_max + + if abs(excess_max) < tol_n_elec # Try to fill all the bands εF = all_eigenvalues[i_max] + T(1) - else - while i_max - i_min > 0 + elseif abs(excess_min) < tol_n_elec # Try to fill only the smallest eigenvalues + εF = all_eigenvalues[i_min] + T(1/2) + *(all_eigenvalues[i_min+1]-all_eigenvalues[i_min]) + + elseif excess_min < -tol_n_elec && excess_max > tol_n_elec + # If excess_min < 0 and excess_max > 0 we can start the bisection method + while i_max - i_min > 1 i = div(i_min+i_max, 2) εF = all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) + if excess < -tol_n_elec i_min = i elseif excess > tol_n_elec @@ -176,11 +191,22 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT i_min = i i_max = i end + + end + + else + if excess_max < -tol_n_elec + error("Could not obtain required number of electrons by filling every state. " * + "Increase n_bands.") + else # excess_min > tol_n_elec + error("Unable to find non-fractional occupations that have the " * + "correct number of electrons. You should add a temperature.") end end occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation - merged_spin_occupations = sum([occ[krange_spin(basis, i)] for i in 1:basis.model.n_spin_components]) + merged_spin_occupations = sum([ occ[krange_spin(basis, i)] + for i in 1:basis.model.n_spin_components ]) if !allequal(merged_spin_occupations) @warn("Not all kpoints have the same number of occupied states, which could mean "* "that a metallic system is treated at zero temperature.") From af5f835c850b5597ff5dd12053fd0fdef3dcf756 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Fri, 14 Jun 2024 16:35:19 +0200 Subject: [PATCH 13/22] MPI correction --- src/occupation.jl | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index c16045432c..3622e0494a 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -148,14 +148,13 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT end # Gathering the eigenvalues distributed over the kpoints in MPI - #n_bands = length(eigenvalues[1]) # assuming that the same number of bands are computed for each kpoint - #counts = [ Int32(n_bands*length(basis.krange_allprocs[rank])) for rank in 1:MPI.Comm_size(basis.comm_kpts)] - #all_eigenvalues = similar(reduce(vcat, eigenvalues)) - #all_eigenvalues_vbuf = MPI.VBuffer(all_eigenvalues, counts) - #MPI.Allgatherv!(reduce(vcat, eigenvalues), all_eigenvalues_vbuf, basis.comm_kpts) - #all_eigenvalues = MPI.Allgatherv(reduce(vcat, eigenvalues), counts, basis.comm_kpts) - - all_eigenvalues = reduce(vcat, eigenvalues) + n_bands = length(eigenvalues[1]) # assuming that the same number of bands + # are computed for each kpoint + counts = [ Int32(n_bands*length(basis.krange_allprocs[rank][])) + for rank in 1:MPI.Comm_size(basis.comm_kpts) ] + all_eigenvalues = Array{T}(undef, sum(counts)) + all_eigenvalues_vbuf = MPI.VBuffer(all_eigenvalues, counts) + MPI.Allgatherv!(reduce(vcat, eigenvalues), all_eigenvalues_vbuf, basis.comm_kpts) # Bisection method to find the index of the eigenvalue such that excess_n_electrons = 0 unique!(sort!((all_eigenvalues))) @@ -167,9 +166,6 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT excess_max = excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max] + T(1); temperature, smearing) - @show length(all_eigenvalues) - @show excess_min, excess_max - if abs(excess_max) < tol_n_elec # Try to fill all the bands εF = all_eigenvalues[i_max] + T(1) elseif abs(excess_min) < tol_n_elec # Try to fill only the smallest eigenvalues From fa6629dee2f284f8600c0402553dbf0625b1942f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Wed, 3 Jul 2024 09:46:04 +0200 Subject: [PATCH 14/22] Correct count with spin --- src/occupation.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 3622e0494a..780848358c 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -150,12 +150,12 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT # Gathering the eigenvalues distributed over the kpoints in MPI n_bands = length(eigenvalues[1]) # assuming that the same number of bands # are computed for each kpoint - counts = [ Int32(n_bands*length(basis.krange_allprocs[rank][])) + counts = [ Int32(n_bands*sum(length.(basis.krange_allprocs[rank]))) for rank in 1:MPI.Comm_size(basis.comm_kpts) ] all_eigenvalues = Array{T}(undef, sum(counts)) all_eigenvalues_vbuf = MPI.VBuffer(all_eigenvalues, counts) MPI.Allgatherv!(reduce(vcat, eigenvalues), all_eigenvalues_vbuf, basis.comm_kpts) - + @show counts # Bisection method to find the index of the eigenvalue such that excess_n_electrons = 0 unique!(sort!((all_eigenvalues))) i_min = 1 @@ -178,7 +178,6 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT i = div(i_min+i_max, 2) εF = all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - if excess < -tol_n_elec i_min = i elseif excess > tol_n_elec From 88326c1306f7eefa1afc7218144a80f6a95376e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Wed, 3 Jul 2024 14:54:19 +0200 Subject: [PATCH 15/22] using searchsorted --- src/occupation.jl | 50 ++++++++++++----------------------------------- 1 file changed, 12 insertions(+), 38 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 780848358c..f2446596c4 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -155,48 +155,22 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT all_eigenvalues = Array{T}(undef, sum(counts)) all_eigenvalues_vbuf = MPI.VBuffer(all_eigenvalues, counts) MPI.Allgatherv!(reduce(vcat, eigenvalues), all_eigenvalues_vbuf, basis.comm_kpts) - @show counts - # Bisection method to find the index of the eigenvalue such that excess_n_electrons = 0 + + # Searching for the Fermi level in between the eigenvalues unique!(sort!((all_eigenvalues))) - i_min = 1 - i_max = length(all_eigenvalues) - excess_min = excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_min] + T(1/2) - *(all_eigenvalues[i_min+1]-all_eigenvalues[i_min]); - temperature, smearing) - excess_max = excess_n_electrons(basis, eigenvalues, all_eigenvalues[i_max] + T(1); - temperature, smearing) - - if abs(excess_max) < tol_n_elec # Try to fill all the bands - εF = all_eigenvalues[i_max] + T(1) - elseif abs(excess_min) < tol_n_elec # Try to fill only the smallest eigenvalues - εF = all_eigenvalues[i_min] + T(1/2) - *(all_eigenvalues[i_min+1]-all_eigenvalues[i_min]) - - elseif excess_min < -tol_n_elec && excess_max > tol_n_elec - # If excess_min < 0 and excess_max > 0 we can start the bisection method - while i_max - i_min > 1 - i = div(i_min+i_max, 2) - εF = all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) - excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - if excess < -tol_n_elec - i_min = i - elseif excess > tol_n_elec - i_max = i - else - i_min = i - i_max = i - end - - end - - else - if excess_max < -tol_n_elec + εFs = [ all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) + for i=1:length(all_eigenvalues)-1 ] + i = searchsortedfirst( map(εF -> excess_n_electrons(basis, eigenvalues, εF; + temperature, smearing), εFs), 0) + if i > length(εFs) + εF = last(all_eigenvalues) + T(1) + excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) + if excess < 0 error("Could not obtain required number of electrons by filling every state. " * "Increase n_bands.") - else # excess_min > tol_n_elec - error("Unable to find non-fractional occupations that have the " * - "correct number of electrons. You should add a temperature.") end + else + εF = εFs[i] end occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation From 7c9b6c9ec6f5e0f1b827de050305a378020b3e94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Wed, 3 Jul 2024 15:40:59 +0200 Subject: [PATCH 16/22] adding tol_n_elec --- src/occupation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index f2446596c4..c54f0ccbb9 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -161,11 +161,11 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT εFs = [ all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) for i=1:length(all_eigenvalues)-1 ] i = searchsortedfirst( map(εF -> excess_n_electrons(basis, eigenvalues, εF; - temperature, smearing), εFs), 0) + temperature, smearing), εFs), -tol_n_elec) if i > length(εFs) εF = last(all_eigenvalues) + T(1) excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - if excess < 0 + if excess < -tol_n_elec error("Could not obtain required number of electrons by filling every state. " * "Increase n_bands.") end From d34be28f59231e66ce5f381307d483d3a0b5327a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Thu, 4 Jul 2024 15:35:52 +0200 Subject: [PATCH 17/22] machine precision in unique and less expensive searchsortedfirst --- src/occupation.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index c54f0ccbb9..680aa35fab 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -146,8 +146,7 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT "occupation $filled_occ. Typically this indicates that you need to put " * "a temperature or switch to a calculation with collinear spin polarization.") end - - # Gathering the eigenvalues distributed over the kpoints in MPI + # Gather the eigenvalues distributed over the kpoints in MPI n_bands = length(eigenvalues[1]) # assuming that the same number of bands # are computed for each kpoint counts = [ Int32(n_bands*sum(length.(basis.krange_allprocs[rank]))) @@ -156,12 +155,18 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT all_eigenvalues_vbuf = MPI.VBuffer(all_eigenvalues, counts) MPI.Allgatherv!(reduce(vcat, eigenvalues), all_eigenvalues_vbuf, basis.comm_kpts) - # Searching for the Fermi level in between the eigenvalues - unique!(sort!((all_eigenvalues))) + # Search for the Fermi level in between the eigenvalues + sort!((all_eigenvalues)) εFs = [ all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) for i=1:length(all_eigenvalues)-1 ] - i = searchsortedfirst( map(εF -> excess_n_electrons(basis, eigenvalues, εF; - temperature, smearing), εFs), -tol_n_elec) + # Remove candidate Fermi levels that are between two identical eigenvalues + # (at machine precision) + εFs = εFs[ [abs(all_eigenvalues[i+1]-all_eigenvalues[i]) > all_eigenvalues[i]*2*eps(T) + for i=1:length(all_eigenvalues)-1] ] + excess_εFs_only(εF) = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) + excess_εFs_only(x::NamedTuple) = x.val # To apply "by" only to searched array εFs ... + # First index such that excess(εFs[i]) > -tol_n_elec + i = searchsortedfirst( εFs, (;val=-tol_n_elec); by=excess_εFs_only) if i > length(εFs) εF = last(all_eigenvalues) + T(1) excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) From 728b4d83e6ab12ec2ba96c211e81954258a62e82 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Fri, 5 Jul 2024 10:36:50 +0200 Subject: [PATCH 18/22] Bissection method by hand --- src/occupation.jl | 38 ++++++++++++++++++++++++++++---------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 680aa35fab..7cdbb06d60 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -146,6 +146,7 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT "occupation $filled_occ. Typically this indicates that you need to put " * "a temperature or switch to a calculation with collinear spin polarization.") end + # Gather the eigenvalues distributed over the kpoints in MPI n_bands = length(eigenvalues[1]) # assuming that the same number of bands # are computed for each kpoint @@ -163,21 +164,38 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT # (at machine precision) εFs = εFs[ [abs(all_eigenvalues[i+1]-all_eigenvalues[i]) > all_eigenvalues[i]*2*eps(T) for i=1:length(all_eigenvalues)-1] ] - excess_εFs_only(εF) = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - excess_εFs_only(x::NamedTuple) = x.val # To apply "by" only to searched array εFs ... - # First index such that excess(εFs[i]) > -tol_n_elec - i = searchsortedfirst( εFs, (;val=-tol_n_elec); by=excess_εFs_only) - if i > length(εFs) - εF = last(all_eigenvalues) + T(1) - excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - if excess < -tol_n_elec + push!(εFs, last(all_eigenvalues) + T(1)) + + # Bisection method to find εF such that abs(excess) < tol_n_elec + i_min = 1 + i_max = length(εFs) + excess_min = excess_n_electrons(basis, eigenvalues, εFs[1]; temperature, smearing) + excess_max = excess_n_electrons(basis, eigenvalues, last(εFs); temperature, smearing) + if excess_max =< tol_n_elec # Try to fill all the bands + εF = last(εFs) + if excess_max < -tol_n_elec error("Could not obtain required number of electrons by filling every state. " * "Increase n_bands.") end + elseif excess_min >= -tol_n_elec # Try to fill only the smallest band(s) + εF = εFs[1] else - εF = εFs[i] - end + while i_max - i_min > 1 + i = div(i_min+i_max, 2) + εF = εFs[i] + excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) + if excess < -tol_n_elec + i_min = i + elseif excess > tol_n_elec + i_max = i + else + i_min = i + i_max = i + end + end + end + occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation merged_spin_occupations = sum([ occ[krange_spin(basis, i)] for i in 1:basis.model.n_spin_components ]) From b49c88933dc777a27e13b61d856154186746b10f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Fri, 5 Jul 2024 10:55:23 +0200 Subject: [PATCH 19/22] typo --- src/occupation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/occupation.jl b/src/occupation.jl index 7cdbb06d60..bc998e8460 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -171,7 +171,7 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT i_max = length(εFs) excess_min = excess_n_electrons(basis, eigenvalues, εFs[1]; temperature, smearing) excess_max = excess_n_electrons(basis, eigenvalues, last(εFs); temperature, smearing) - if excess_max =< tol_n_elec # Try to fill all the bands + if excess_max <= tol_n_elec # Try to fill all the bands εF = last(εFs) if excess_max < -tol_n_elec error("Could not obtain required number of electrons by filling every state. " * From a4fb489f8a309c62edd745ade67c1c963da2c0fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Mon, 16 Sep 2024 09:41:09 +0200 Subject: [PATCH 20/22] formatting --- src/occupation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index bc998e8460..820a5aa5cf 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -171,15 +171,15 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT i_max = length(εFs) excess_min = excess_n_electrons(basis, eigenvalues, εFs[1]; temperature, smearing) excess_max = excess_n_electrons(basis, eigenvalues, last(εFs); temperature, smearing) - if excess_max <= tol_n_elec # Try to fill all the bands + if excess_max <= tol_n_elec # Try to fill all the bands εF = last(εFs) if excess_max < -tol_n_elec error("Could not obtain required number of electrons by filling every state. " * "Increase n_bands.") end - elseif excess_min >= -tol_n_elec # Try to fill only the smallest band(s) - εF = εFs[1] - else + elseif excess_min >= -tol_n_elec # Try to fill only the smallest band(s) + εF = first(εFs) + else # Bissection to find first(εFs) < εF < last(εFs) while i_max - i_min > 1 i = div(i_min+i_max, 2) εF = εFs[i] From 487761058862308508b293e2e65af478db0c5e45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Tue, 5 Nov 2024 19:45:48 +0100 Subject: [PATCH 21/22] Bissection function --- src/occupation.jl | 66 ++++++++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 29 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 820a5aa5cf..199fe3f961 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -130,6 +130,41 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiTwoSt Roots.find_zero(excess, εF, Roots.Secant(), Roots.Bisection(); atol=eps(T)) end +""" +Bisection method to find εF in the array εFs such that |excess_funcion(εF)| < tol_n_elec +where excess_function is an increasing function and εFs is a sorted array (increasing). +""" +function bisection(εFs, excess_function, tol_n_elec) + i_min = 1 + i_max = length(εFs) + excess_min = excess_function(εFs[1]) + excess_max = excess_function(last(εFs)) + if excess_max <= tol_n_elec # Try to fill all the bands + εF = last(εFs) + if excess_max < -tol_n_elec + error("Could not obtain required number of electrons by filling every state. " * + "Increase n_bands.") + end + elseif excess_min >= -tol_n_elec # Try to fill only the smallest band(s) + εF = first(εFs) + else # Bissection to find first(εFs) < εF < last(εFs) + while i_max - i_min > 1 + i = div(i_min+i_max, 2) + εF = εFs[i] + excess = excess_function(εF) + if excess < -tol_n_elec + i_min = i + elseif excess > tol_n_elec + i_max = i + else + i_min = i + i_max = i + end + + end + end + εF +end # Note: This is not exported, but only called by the above algorithms for # the zero-temperature case. @@ -166,35 +201,8 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT for i=1:length(all_eigenvalues)-1] ] push!(εFs, last(all_eigenvalues) + T(1)) - # Bisection method to find εF such that abs(excess) < tol_n_elec - i_min = 1 - i_max = length(εFs) - excess_min = excess_n_electrons(basis, eigenvalues, εFs[1]; temperature, smearing) - excess_max = excess_n_electrons(basis, eigenvalues, last(εFs); temperature, smearing) - if excess_max <= tol_n_elec # Try to fill all the bands - εF = last(εFs) - if excess_max < -tol_n_elec - error("Could not obtain required number of electrons by filling every state. " * - "Increase n_bands.") - end - elseif excess_min >= -tol_n_elec # Try to fill only the smallest band(s) - εF = first(εFs) - else # Bissection to find first(εFs) < εF < last(εFs) - while i_max - i_min > 1 - i = div(i_min+i_max, 2) - εF = εFs[i] - excess = excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - if excess < -tol_n_elec - i_min = i - elseif excess > tol_n_elec - i_max = i - else - i_min = i - i_max = i - end - - end - end + excess_function = εF->excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) + εF = bisection(εFs, excess_function, tol_n_elec) occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation merged_spin_occupations = sum([ occ[krange_spin(basis, i)] From 551c83ccb8c96dcd48af09fa2f88e6810eea560d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9mentine=20Barat?= Date: Sun, 1 Dec 2024 17:12:55 +0100 Subject: [PATCH 22/22] cleaning --- src/occupation.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/occupation.jl b/src/occupation.jl index 199fe3f961..4c1c47ac74 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -131,10 +131,11 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiTwoSt end """ -Bisection method to find εF in the array εFs such that |excess_funcion(εF)| < tol_n_elec +Discrete bisection method to find εF in the array εFs such that +|excess_funcion(εF)| < tol_n_elec where excess_function is an increasing function and εFs is a sorted array (increasing). """ -function bisection(εFs, excess_function, tol_n_elec) +function discrete_find_zero(excess_function, εFs, tol_n_elec) i_min = 1 i_max = length(εFs) excess_min = excess_function(εFs[1]) @@ -160,7 +161,6 @@ function bisection(εFs, excess_function, tol_n_elec) i_min = i i_max = i end - end end εF @@ -193,20 +193,18 @@ function compute_fermi_level(basis::PlaneWaveBasis{T}, eigenvalues, ::FermiZeroT # Search for the Fermi level in between the eigenvalues sort!((all_eigenvalues)) - εFs = [ all_eigenvalues[i] + T(1/2)*(all_eigenvalues[i+1]-all_eigenvalues[i]) - for i=1:length(all_eigenvalues)-1 ] + εFs = all_eigenvalues[1:end-1] .+ T(1/2)*diff(all_eigenvalues) # Candidates Fermi-levels # Remove candidate Fermi levels that are between two identical eigenvalues # (at machine precision) - εFs = εFs[ [abs(all_eigenvalues[i+1]-all_eigenvalues[i]) > all_eigenvalues[i]*2*eps(T) - for i=1:length(all_eigenvalues)-1] ] + εFs = εFs[ diff(all_eigenvalues) .> 2*eps(T)*all_eigenvalues[1:end-1] ] push!(εFs, last(all_eigenvalues) + T(1)) excess_function = εF->excess_n_electrons(basis, eigenvalues, εF; temperature, smearing) - εF = bisection(εFs, excess_function, tol_n_elec) + εF = discrete_find_zero(excess_function, εFs, tol_n_elec) occ = compute_occupation(basis, eigenvalues, εF; temperature, smearing).occupation - merged_spin_occupations = sum([ occ[krange_spin(basis, i)] - for i in 1:basis.model.n_spin_components ]) + merged_spin_occupations = sum( occ[krange_spin(basis, i)] + for i=1:basis.model.n_spin_components ) if !allequal(merged_spin_occupations) @warn("Not all kpoints have the same number of occupied states, which could mean "* "that a metallic system is treated at zero temperature.")