/
dos.jl
68 lines (61 loc) · 2.56 KB
/
dos.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# Compute densities of states
# IDOS (integrated density of states)
# N(ε) = sum_n f_n = sum_n f((εn-ε)/temperature)
# DOS (density of states)
# D(ε) = N'(ε)
# LDOS (local density of states)
# LD = sum_n f_n |ψn|^2
"""
Total density of states at energy ε
"""
function compute_dos(ε, basis, eigenvalues; smearing=basis.model.smearing,
temperature=basis.model.temperature)
if (temperature == 0) || smearing isa Smearing.None
error("compute_dos only supports finite temperature")
end
filled_occ = filled_occupation(basis.model)
D = zeros(typeof(ε), basis.model.n_spin_components)
for σ in 1:basis.model.n_spin_components, ik = krange_spin(basis, σ)
for (iband, εnk) in enumerate(eigenvalues[ik])
enred = (εnk - ε) / temperature
D[σ] -= (filled_occ * basis.kweights[ik] / temperature
* Smearing.occupation_derivative(smearing, enred))
end
end
D = mpi_sum(D, basis.comm_kpts)
end
function compute_dos(scfres::NamedTuple; ε=scfres.εF, kwargs...)
compute_dos(ε, scfres.basis, scfres.eigenvalues; kwargs...)
end
"""
Local density of states, in real space. `weight_threshold` is a threshold
to screen away small contributions to the LDOS.
"""
function compute_ldos(ε, basis::PlaneWaveBasis{T}, eigenvalues, ψ;
smearing=basis.model.smearing,
temperature=basis.model.temperature,
weight_threshold=eps(T)) where {T}
if (temperature == 0) || smearing isa Smearing.None
error("compute_ldos only supports finite temperature")
end
filled_occ = filled_occupation(basis.model)
weights = deepcopy(eigenvalues)
for (ik, εk) in enumerate(eigenvalues)
for (iband, εnk) in enumerate(εk)
enred = (εnk - ε) / temperature
weights[ik][iband] = (-filled_occ / temperature
* Smearing.occupation_derivative(smearing, enred))
end
end
# Use compute_density routine to compute LDOS, using just the modified
# weights (as "occupations") at each k-point. Note, that this automatically puts in the
# required symmetrization with respect to kpoints and BZ symmetry
compute_density(basis, ψ, weights; occupation_threshold=weight_threshold)
end
function compute_ldos(scfres::NamedTuple; ε=scfres.εF, kwargs...)
compute_ldos(ε, scfres.basis, scfres.eigenvalues, scfres.ψ; kwargs...)
end
"""
Plot the density of states over a reasonable range. Requires to load `Plots.jl` beforehand.
"""
function plot_dos end