@@ -68,65 +68,196 @@ function compute_ldos(scfres::NamedTuple; ε=scfres.εF, kwargs...)
6868end
6969
7070"""
71- Projected density of states at energy ε for an atom with given i and l.
71+ compute_pdos(εs, basis, ψ, eigenvalues; [positions, smearing, temperature])
72+
73+ Compute the projected density of states (PDOS) for all atoms and orbitals.
74+
75+ Overview of inputs:
76+ - `εs`: Vector of energies at which the PDOS will be computed
77+
78+ Overview of outputs:
79+ - `pdos`: 3D array of PDOS, pdos[iε_idx, iproj, σ] = PDOS at energy εs[iε_idx] for projector iproj and spin σ
80+ - `projector_labels` : Vector of tuples (iatom, species, label, n, l, m) for each projector,
81+ that maps the iproj index to the corresponding atomic orbital.
82+ For details see the documentation for `atomic_orbital_projectors`.
83+
84+ Notes:
85+ - The atomic orbital projectors are taken from the pseudopotential files used to build the atoms.
86+ In doubt, consult the pseudopotential file for the list of available atomic orbitals.
7287"""
73- function compute_pdos (ε, basis, eigenvalues, ψ, i, l, psp, position;
74- smearing= basis. model. smearing,
75- temperature= basis. model. temperature)
88+ function compute_pdos (εs, basis:: PlaneWaveBasis{T} , ψ, eigenvalues;
89+ positions= basis. model. positions,
90+ smearing= basis. model. smearing,
91+ temperature= basis. model. temperature) where {T}
7692 if (temperature == 0 ) || smearing isa Smearing. None
7793 error (" compute_pdos only supports finite temperature"
7894 end
7995 filled_occ = filled_occupation (basis. model)
96+ projections, projector_labels = atomic_orbital_projections (basis, ψ; positions)
97+ nprojs = length (projector_labels)
8098
81- projs = compute_pdos_projs (basis, eigenvalues, ψ, i, l, psp, position)
82-
83- D = zeros (typeof (ε[1 ]), length (ε), 2 l+ 1 )
84- for (i, iε) in enumerate (ε)
85- for (ik, projk) in enumerate (projs)
99+ D = zeros (typeof (εs[1 ]), length (εs), nprojs, basis. model. n_spin_components)
100+ for (iε, ε) in enumerate (εs)
101+ for σ in 1 : basis. model. n_spin_components, ik = krange_spin (basis, σ)
102+ projsk = projections[ik]
86103 @views for (iband, εnk) in enumerate (eigenvalues[ik])
87- enred = (εnk - iε) / temperature
88- D[i, :] .- = (filled_occ .* basis. kweights[ik] .* projk[iband, :]
89- ./ temperature
90- .* Smearing. occupation_derivative (smearing, enred))
104+ enred = (εnk - ε) / temperature
105+ for iproj in 1 : size (projsk, 2 )
106+ D[iε, iproj, σ] -= (filled_occ * basis. kweights[ik] * projsk[iband, iproj]
107+ / temperature
108+ * Smearing. occupation_derivative (smearing, enred))
109+ end
91110 end
92111 end
93112 end
94- D = mpi_sum (D, basis. comm_kpts)
113+ pdos = mpi_sum (D, basis. comm_kpts)
114+
115+ return (; pdos, projector_labels, εs)
95116end
96117
97- function compute_pdos (scfres:: NamedTuple , iatom, i, l; ε= scfres. εF, kwargs... )
98- psp = scfres. basis. model. atoms[iatom]. psp
99- position = scfres. basis. model. positions[iatom]
100- # TODO do the symmetrization instead of unfolding
101- scfres = unfold_bz (scfres)
102- compute_pdos (ε, scfres. basis, scfres. eigenvalues, scfres. ψ, i, l, psp, position; kwargs... )
118+ function compute_pdos (εs, bands; kwargs... )
119+ compute_pdos (εs, bands. basis, bands. ψ, bands. eigenvalues; kwargs... )
103120end
104121
105- # Build atomic orbitals projections projs[ik][iband,m] = |<ψnk, ϕ>|^2 for a single atom.
106- # TODO optimization ? accept a range of positions, in case we want to compute the PDOS
107- # for all atoms of the same type (and reuse the computation of the atomic orbitals)
108- function compute_pdos_projs (basis, eigenvalues, ψ, i, l, psp:: NormConservingPsp , position)
109- # Precompute the form factors on all kpoints at once so we can better use the cache (memory-compute tradeoff).
110- # Revisit this (pass the cache around explicitly) if RAM becomes an issue.
122+ """
123+ Structure for manifold choice and projectors extraction.
124+
125+ Overview of fields:
126+ - `iatom`: Atom position in the atoms array.
127+ - `species`: Chemical Element as in ElementPsp.
128+ - `label`: Orbital name, e.g.: "3S".
129+
130+ All fields are optional, only the given ones will be used for selection.
131+ Can be called with an orbital NamedTuple and returns a boolean
132+ stating whether the orbital belongs to the manifold.
133+ """
134+ @kwdef struct OrbitalManifold
135+ iatom :: Union{Int64, Nothing} = nothing
136+ species :: Union{Symbol, AtomsBase.ChemicalSpecies, Nothing} = nothing
137+ label :: Union{String, Nothing} = nothing
138+ end
139+ function (s:: OrbitalManifold )(orb)
140+ iatom_match = isnothing (s. iatom) || (s. iatom == orb. iatom)
141+ # See JuliaMolSim/AtomsBase.jl#139 why both species equalities are needed
142+ species_match = isnothing (s. species) || (s. species == orb. species) || (orb. species == s. species)
143+ label_match = isnothing (s. label) || (s. label == orb. label)
144+ iatom_match && species_match && label_match
145+ end
146+
147+ """
148+ atomic_orbital_projectors(basis; [isonmanifold, positions])
149+
150+ Build the matrices of projectors onto the pseudoatomic orbitals.
151+
152+ projector[ik][iG, iproj] = 1/√Ω FT{ ϕperₙₗₘ(. - Rᵢ) }(k+G) + orthogonalization
153+
154+ where Ω is the unit cell volume, ϕperₙₗₘ(. - Rᵢ) is the periodized pseudoatomic orbital (n, l, m) centered at Rᵢ
155+ and iproj is the index corresponding to atom i and the quantum numbers (n, l, m).
156+ This correspondance is recorded in `labels`.
157+
158+ The projectors are computed by decomposition into a form factor multiplied by a structure factor:
159+ FT{ ϕperₙₗₘ(. - Rᵢ) }(k+G) = Fourier transform of periodized atomic orbital ϕₙₗₘ (form factor)
160+ * structure factor for atom center exp(-i<Rᵢ,k+G>)
161+
162+ Overview of inputs:
163+ - `positions` : Positions of the atoms in the unit cell. Default is model.positions.
164+ - `isonmanifold` (opt) : (see notes below) OrbitalManifold struct to select only a subset of orbitals
165+ for the computation.
166+
167+ Overview of outputs:
168+ - `projectors`: Vector of matrices of projectors
169+ - `labels`: Vector of NamedTuples containing iatom, species, n, l, m and orbital name for each projector
170+
171+ Notes:
172+ - The orbitals used for the projectors are all orthogonalized against each other.
173+ This corresponds to ortho-atomic projectors in Quantum Espresso.
174+ - 'n' in labels is not exactly the principal quantum number, but rather the index of the radial function
175+ in the pseudopotential. As an example, if the only S orbitals in the pseudopotential are 3S and 4S,
176+ then those are indexed as n=1, l=0 and n=2, l=0 respectively.
177+ - Use 'isonmanifold' kwarg with caution, since the resulting projectors would be
178+ orthonormalized only against the manifold basis.
179+ Most applications require the whole projectors basis to be orthonormal instead.
180+ """
181+ function atomic_orbital_projectors (basis:: PlaneWaveBasis{T} ;
182+ isonmanifold = l -> true ,
183+ positions = basis. model. positions) where {T}
184+
111185 G_plus_k_all = [Gplusk_vectors (basis, basis. kpoints[ik])
112186 for ik = 1 : length (basis. kpoints)]
113187 G_plus_k_all_cart = [map (recip_vector_red_to_cart (basis. model), gpk)
114188 for gpk in G_plus_k_all]
115189
116- # Build form factors of pseudo-wavefunctions centered at 0.
117- fun (p) = eval_psp_pswfc_fourier (psp, i, l, p)
118- # form_factors_all[ik][p,m]
119- form_factors_all = build_form_factors (fun, l, G_plus_k_all_cart)
120-
121- projs = Vector {Matrix} (undef, length (basis. kpoints))
122- for (ik, ψk) in enumerate (ψ)
123- structure_factor = [cis2pi (- dot (position, p)) for p in G_plus_k_all[ik]]
124- # TODO orthogonalize pseudo-atomic wave functions?
125- proj_vectors = structure_factor .* form_factors_all[ik] ./ sqrt (basis. model. unit_cell_volume)
126- projs[ik] = abs2 .(ψk' * proj_vectors) # contract on p -> projs[ik][iband,m]
190+ projectors = [Matrix {Complex{T}} (undef, length (G_plus_k), 0 ) for G_plus_k in G_plus_k_all]
191+ labels = []
192+ for (iatom, atom) in enumerate (basis. model. atoms)
193+ psp = atom. psp
194+ for l in 0 : psp. lmax, n in 1 : DFTK. count_n_pswfc_radial (psp, l)
195+ label = DFTK. pswfc_label (psp, n, l)
196+ if ! isonmanifold ((; iatom, atom. species, label))
197+ continue
198+ end
199+ fun (p) = eval_psp_pswfc_fourier (psp, n, l, p)
200+ form_factors_l = build_form_factors (fun, l, G_plus_k_all_cart)
201+ for ik in 1 : length (G_plus_k_all_cart)
202+ structure_factor = [cis2pi (- dot (positions[iatom], p)) for p in G_plus_k_all[ik]]
203+ projectors[ik] = hcat (projectors[ik],
204+ form_factors_l[ik] .* structure_factor
205+ ./ sqrt (basis. model. unit_cell_volume))
206+ end
207+ for m in - l: l
208+ push! (labels, (; iatom, atom. species, n, l, m, label))
209+ end
210+ end
211+ end
212+
213+ projectors = ortho_lowdin .(projectors)
214+
215+ return (; projectors, labels)
216+ end
217+
218+ """
219+ atomic_orbital_projections(basis, ψ; [isonmanifold, positions])
220+
221+ Build the projection matrices of ψ onto each pseudo-atomic orbital.
222+
223+ projection[ik][iband, iproj] = ‖ ψ'[ik][:, iband] * projector[ik][:, iproj] ‖²
224+
225+ For more details, see documentation for [`atomic_orbital_projectors`](@ref).
226+ """
227+ function atomic_orbital_projections (basis:: PlaneWaveBasis{T} , ψ;
228+ isonmanifold = l -> true ,
229+ positions = basis. model. positions) where {T}
230+ projectors, labels = atomic_orbital_projectors (basis; isonmanifold, positions)
231+ projections = map (ψ, projectors) do ψk, projectorsk
232+ abs2 .(ψk' * projectorsk)
127233 end
128234
129- projs
235+ return (; projections, labels)
236+ end
237+
238+ """
239+ sum_pdos(pdos_res, manifolds)
240+
241+ This function extracts and sums up all the PDOSes, directly from the output of the `compute_pdos` function,
242+ that match any of the manifolds.
243+
244+ Overview of inputs:
245+ - `pdos_res`: Whole output from compute_pdos.
246+ - `manifolds`: Vector of OrbitalManifolds to select the desired projectors pdos.
247+
248+ Overview of outputs:
249+ - `pdos`: Vector containing the pdos(ε).
250+ """
251+ function sum_pdos (pdos_res, manifolds:: AbstractVector )
252+ pdos = zeros (Float64, length (pdos_res. εs), size (pdos_res. pdos, 3 ))
253+ for σ in 1 : size (pdos_res. pdos, 3 )
254+ for (j, orb) in enumerate (pdos_res. projector_labels)
255+ if any (manifold (orb) for manifold in manifolds)
256+ pdos[:, σ] += pdos_res. pdos[:, j, σ]
257+ end
258+ end
259+ end
260+ return pdos
130261end
131262
132263"""
0 commit comments