diff --git a/src/DFTK.jl b/src/DFTK.jl index a812986bbb..e9c491ba45 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -52,6 +52,7 @@ export PspUpf include("pseudo/NormConservingPsp.jl") include("pseudo/PspHgh.jl") include("pseudo/PspUpf.jl") +include("pseudo/PspLinComb.jl") export ElementPsp export ElementCohenBergstresser @@ -60,6 +61,7 @@ export ElementGaussian export charge_nuclear, charge_ionic export n_elec_valence, n_elec_core export element_symbol, mass, species # Note: Re-exported from AtomsBase +export virtual_crystal_element include("elements.jl") export SymOp diff --git a/src/elements.jl b/src/elements.jl index 68b43c63fc..15c4f07634 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -52,10 +52,8 @@ function core_charge_density_fourier(::Element, ::T)::T where {T <: Real} error("Abstract elements do not necesesarily provide core charge density.") end -# Fallback print function: Base.show(io::IO, el::Element) = print(io, "$(typeof(el))(:$(species(el)))") - # # ElementCoulomb # @@ -97,15 +95,19 @@ struct ElementPsp{P} <: Element family::Union{PseudoFamily,Nothing} # PseudoFamily if known, else Nothing mass # Atomic mass end -function Base.show(io::IO, el::ElementPsp) +function Base.show(io::IO, el::ElementPsp{P}) where {P} print(io, "ElementPsp(:$(el.species), ") if !isnothing(el.family) - print(io, el.family, ")") + print(io, el.family) elseif isempty(el.psp.identifier) - print(io, "custom psp)") + print(io, "custom psp") else - print(io, "\"$(el.psp.identifier)\")") + print(io, "\"$(el.psp.identifier)\"") + end + if P <: PspLinComb + print(io, ", ", el.psp.description) end + print(io, ")") end pseudofamily(el::ElementPsp) = el.family @@ -275,3 +277,28 @@ function local_potential_fourier(el::ElementGaussian, p::Real) -el.α * exp(- (p * el.L)^2 / 2) # = ∫_ℝ³ V(x) exp(-ix⋅p) dx end # TODO Strictly speaking needs a eval_psp_energy_correction + +# +# Helper functions +# + +# TODO Better name ? +# TODO docstring +function virtual_crystal_element(coefficients::Vector{<:Number}, elements::Vector{<:ElementPsp{<:NormConservingPsp}}; + species=ChemicalSpecies(0)) + length(coefficients) == length(elements) || throw( + ArgumentError("Expect coefficients and elements to have equal length.")) + sum(coefficients) ≈ one(eltype(coefficients)) || throw( + ArgumentError("Expect coefficients to sum to 1")) + family = allequal(p -> p.family, elements) ? first(elements).family : nothing + + description = "lin. comb. of " + description *= join([(@sprintf "%.2f*%s" c "$(el.species)") for (c, el) in zip(coefficients, elements)], " ") + lincomb_psp = PspLinComb(coefficients, [el.psp for el in elements]; description) + lincomb_mass = sum(c * mass(el) for (c, el) in zip(coefficients, elements)) + ElementPsp(species, lincomb_psp, family, lincomb_mass) +end +function virtual_crystal_element(α::Number, αpsp::ElementPsp{<:NormConservingPsp}, + β::Number, βpsp::ElementPsp{<:NormConservingPsp}; kwargs...) + virtual_crystal_element([α, β], [αpsp, βpsp]; kwargs...) +end diff --git a/src/pseudo/PspLinComb.jl b/src/pseudo/PspLinComb.jl new file mode 100644 index 0000000000..126b5341aa --- /dev/null +++ b/src/pseudo/PspLinComb.jl @@ -0,0 +1,83 @@ +# Linear combination of two pseudos +# Note that this data structure is deliberately not yet exported for now +# as we may want to find a different solution for representing virtual crystal elements +# in the future; users should use virtual_crystal_element to implicitly construct +# these objects. +# +struct PspLinComb{T, P <: NormConservingPsp} <: NormConservingPsp + lmax::Int + h::Vector{Matrix{T}} + hidx_to_psp_proj::Vector{Vector{Tuple{Int,Int}}} # map per angular momentum and projector index of this object to the tuple (psp index, projector index) + coefficients::Vector{T} + pseudos::Vector{P} + identifier::String # String identifying the PSP + description::String # Descriptive string +end + +function PspLinComb(coefficients::Vector{T}, pseudos::Vector{<:NormConservingPsp}; + description="psp linear combination") where {T} + @assert length(coefficients) == length(pseudos) + @assert sum(coefficients) ≈ one(eltype(coefficients)) + @assert !any(p -> p isa PspLinComb, pseudos) + + # These assumptions we could lift, but we make tham to make our lifes easier + @assert allequal(charge_ionic, pseudos) + @assert allequal(has_valence_density, pseudos) + @assert allequal(has_core_density, pseudos) + @assert allequal(p -> p.lmax, pseudos) + + TT = promote_type(T, eltype(eltype(first(pseudos).h))) + lmax = first(pseudos).lmax + h = Matrix{TT}[] + hidx_to_psp_proj = Vector{Tuple{Int,Int}}[] # (ipsp, iproj) + for l in 0:lmax + proj_l_total = sum(p -> count_n_proj_radial(p, l), pseudos) + hl = zeros(T, proj_l_total, proj_l_total) + splits = Tuple{Int,Int}[] # (ipsp, iproj) + + for (ipsp, p) in enumerate(pseudos) + iproj_offset = isempty(splits) ? 0 : last(splits)[2] + rnge = iproj_offset .+ (1:count_n_proj_radial(p, l)) + # TODO Is this the correct scaling of the coefficient in h ? + hl[rnge, rnge] = p.h[l+1] * coefficients[ipsp] + append!(splits, [(ipsp, iproj) for iproj in 1:count_n_proj_radial(p, l)]) + end + push!(h, hl) + push!(hidx_to_psp_proj, splits) + end + @assert length(h) == length(hidx_to_psp_proj) == lmax+1 + + identifier = "" + PspLinComb(lmax, h, hidx_to_psp_proj, coefficients, pseudos, identifier, description) +end + +charge_ionic(psp::PspLinComb) = charge_ionic(first(psp.pseudos)) +has_valence_density(psp::PspLinComb) = all(has_valence_density, psp.pseudos) +has_core_density(psp::PspLinComb) = any(has_core_density, psp.pseudos) + +function eval_psp_projector_real(psp::PspLinComb, i, l, r::Real) + (ipsp, iproj) = psp.hidx_to_psp_proj[l+1][i] + eval_psp_projector_real(psp.pseudos[ipsp], iproj, l, r) +end +function eval_psp_projector_fourier(psp, i, l, p::Real) + (ipsp, iproj) = psp.hidx_to_psp_proj[l+1][i] + eval_psp_projector_fourier(psp.pseudos[ipsp], iproj, l, p) +end + +function eval_psp_energy_correction(T, psp::PspLinComb) + sum(c * eval_psp_energy_correction(T, p) for (c, p) in zip(psp.coefficients, psp.pseudos)) +end + +macro make_psplincomb_call(fn) + quote + function $fn(psp::PspLinComb, arg::Real) + sum(c * $fn(pp, arg) for (c, pp) in zip(psp.coefficients, psp.pseudos)) + end + end +end +@make_psplincomb_call DFTK.eval_psp_local_real +@make_psplincomb_call DFTK.eval_psp_local_fourier +@make_psplincomb_call DFTK.eval_psp_density_valence_real +@make_psplincomb_call DFTK.eval_psp_density_valence_fourier +@make_psplincomb_call DFTK.eval_psp_density_core_real +@make_psplincomb_call DFTK.eval_psp_density_core_fourier diff --git a/src/pseudo/pseudopotential_data.jl b/src/pseudo/pseudopotential_data.jl index e998eb728e..9430c16ef3 100644 --- a/src/pseudo/pseudopotential_data.jl +++ b/src/pseudo/pseudopotential_data.jl @@ -29,7 +29,7 @@ Return the recommended kinetic energy cutoff, supersampling and density cutoff f Values may be `missing` if the data cannot be determined. """ function PseudoPotentialData.recommended_cutoff(el::Element) - if isnothing(pseudofamily(el)) + if isnothing(pseudofamily(el)) || :X == element_symbol(el) return (; Ecut=missing, supersampling=missing, Ecut_density=missing) else return recommended_cutoff(pseudofamily(el), element_symbol(el)) @@ -45,7 +45,7 @@ Effectively this returns `pseudometa(pseudofamily(element), element_symbol(eleme """ function PseudoPotentialData.pseudometa(el::Element) family = pseudofamily(el) - if isnothing(family) + if isnothing(family) || :X == element_symbol(el) return Dict{String,Any}() else return pseudometa(family, element_symbol(el)) diff --git a/src/terms/operators.jl b/src/terms/operators.jl index 673d91c394..020e8b72cb 100644 --- a/src/terms/operators.jl +++ b/src/terms/operators.jl @@ -1,6 +1,8 @@ ### Linear operators operating on quantities in real or fourier space # This is the optimized low-level interface. Use the functions in Hamiltonian for high-level usage +import Base: Matrix, Array + """ Linear operators that act on tuples (real, fourier) The main entry point is `apply!(out, op, in)` which performs the operation `out += op*in` @@ -41,7 +43,7 @@ struct NoopOperator{T <: Real} <: RealFourierOperator kpoint::Kpoint{T} end apply!(Hψ, op::NoopOperator, ψ) = nothing -function Base.Matrix(op::NoopOperator) +function Matrix(op::NoopOperator) n_Gk = length(G_vectors(op.basis, op.kpoint)) zeros_like(G_vectors(op.basis), eltype(op.basis), n_Gk, n_Gk) end @@ -60,7 +62,7 @@ end function apply!(Hψ, op::RealSpaceMultiplication, ψ) Hψ.real .+= op.potential .* ψ.real end -function Base.Matrix(op::RealSpaceMultiplication) +function Matrix(op::RealSpaceMultiplication) # V(G, G') = = 1/sqrt(Ω) pot_fourier = fft(op.basis, op.potential) n_G = length(G_vectors(op.basis, op.kpoint)) diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl index c9295085f9..1174a72bf9 100644 --- a/src/workarounds/forwarddiff_rules.jl +++ b/src/workarounds/forwarddiff_rules.jl @@ -202,6 +202,16 @@ function construct_value(psp::PspUpf{T,I}) where {T <: AbstractFloat, I <: Abstr # but does not yet permit response derivatives w.r.t. UPF parameters. psp end +construct_value(psp::PspLinComb) = psp +function construct_value(psp::PspLinComb{<: Dual}) + PspLinComb(psp.lmax, + [ForwardDiff.value.(hl) for hl in psp.h], + psp.hidx_to_psp_proj, + ForwardDiff.value.(psp.coefficients), + construct_value.(psp.pseudos), + psp.identifier, + psp.description) +end function construct_value(basis::PlaneWaveBasis{T}) where {T <: Dual} # NOTE: This is a pretty slow function as it *recomputes* basically diff --git a/test/PspLinComb.jl b/test/PspLinComb.jl new file mode 100644 index 0000000000..e9e660210a --- /dev/null +++ b/test/PspLinComb.jl @@ -0,0 +1,6 @@ + + +# Test whether doing a 0.5 of one element with itself yields the identical answer as not doing the lincomb at all + + +# Test whether the linear combination pseudo does not break any important properities (like normalisation, being norm-conserving etc.)