[src] Return InterpModel in read_w90_tb

Note due to bugs in Julia v1.8.3, the tests fail.
JuliaLang/julia#47685

The tests work fine with a master branch julia
JuliaLang/julia@c8ea33d

Project.toml

@@ -27,6 +27,7 @@ WannierIO = "cb1bc77f-5443-4951-af9f-05b616a3e422"
 Bravais = "0.1.8"
 Brillouin = "0.5.8"
 Comonicon = "1.0.0"
+DelimitedFiles = "1.9.0"
 LazyGrids = "0.4.0"
 NLSolversBase = "7.8.2"
 NearestNeighbors = "0.4.11"

docs/src/api/interpolation.md

@@ -109,7 +109,7 @@ while [`Model`](@ref Model) works only for Wannierization.
 
 ```@docs
 InterpModel
-InterpModel(kRvectors, H, kpath)
+InterpModel(kRvectors, kpath, H)
 InterpModel(model::Model; mdrs::Bool=true)
 ```
 

src/cli/fermisurf.jl

@@ -28,15 +28,15 @@ Interpolate Fermi surface.
     ws::Bool=false,
 )
     if tb
-        Rvecs, H, r = read_w90_tb(seedname)
+        model = read_w90_tb(seedname)
     else
         if isempty(amn)
             model = read_w90_post(seedname)
         else
             model = read_w90_post(seedname; chk=false, amn=amn)
         end
-        Rvecs = model.kRvectors.Rvectors
     end
+    Rvecs = model.kRvectors.Rvectors
     _print_type(Rvecs)
 
     # If you really want, you can use WS interpolation.
@@ -45,7 +45,7 @@ Interpolate Fermi surface.
     if ws && (Rvecs isa RVectorsMDRS)
         Rvecs = Rvecs.Rvectors
     end
-    kpoints, E = fermi_surface(Rvecs, H; n_k=nk)
+    kpoints, E = fermi_surface(Rvecs, model.H; n_k=nk)
 
     origin = zeros(Float64, 3)
     recip_latt = get_recip_lattice(Rvecs.lattice)

src/interp/fermisurf.jl

@@ -32,18 +32,10 @@ function fermi_surface(
 ) where {
     T<:Real,RV<:Union{RVectors{T},RVectorsMDRS{T}},KT<:Union{AbstractVector{Int},Integer}
 }
-    length(n_k) in [1, 3] || error("n_k must be an integer or a vector of length 3")
-
     n_wann, _, n_rvecs = size(H)
     n_rvecs == Rvectors.n_rvecs || error("n_rvecs of H != Rvectors.n_rvecs")
 
-    # bxsf need general grid, i.e., the last kpoint is the periodic image of the first one
-    # so I increase the n_k by 1
-    if length(n_k) == 3
-        n_kx, n_ky, n_kz = [n + 1 for n in n_k]
-    else
-        n_kx = n_ky = n_kz = n_k + 1
-    end
+    n_kx, n_ky, n_kz = _expand_nk(n_k)
     # kpoints are in fractional coordinates
     kpoints = get_kpoints([n_kx, n_ky, n_kz]; endpoint=true)
     n_kpts = n_kx * n_ky * n_kz
@@ -81,3 +73,31 @@ function fermi_surface(
 
     return kpoints, E
 end
+
+"""
+    fermi_surface(model; n_k)
+
+Interpolate Fermi surface.
+
+# Arguments
+- `model`: [`InterpModel`](@ref)
+
+See also [`fermi_surface`](@ref fermi_surface(Rvectors::RV, H::AbstractArray{Complex{T},3}; n_k::KT)).
+"""
+function fermi_surface(model::InterpModel; n_k)
+    return fermi_surface(model.kRvectors.Rvectors, model.H; n_k=n_k)
+end
+
+function _expand_nk(n_k::T) where {T<:Union{AbstractVector{Int},Integer}}
+    length(n_k) in [1, 3] || error("n_k must be an integer or a vector of length 3")
+
+    # bxsf need general grid, i.e., the last kpoint is the periodic image of the first one
+    # so I increase the n_k by 1
+    if length(n_k) == 3
+        n_kx, n_ky, n_kz = [n + 1 for n in n_k]
+    else
+        n_kx = n_ky = n_kz = n_k + 1
+    end
+
+    return n_kx, n_ky, n_kz
+end

src/interp/interp_model.jl

@@ -10,19 +10,28 @@ Usually after Wannierization of a [`Model`](@ref Model), we construct this
 
 # Fields
 - `kRvectors`: the kpoints and R vectors
-- `H`: `n_wann * n_wann * n_rvecs`, the Hamiltonian in real space
 - `kpath`: the kpoint path for band structure
+- `H`: `n_wann * n_wann * n_rvecs`, the Hamiltonian in real space
+- `r`: `n_wann * n_wann * n_rvecs`, the position operator in real space
 - `S`: `n_wann * n_wann * n_rvecs * 3`, optional, the spin operator
+
+!!! note
+    For MDRS interpolation, the `n_rvecs` in this docstring is actually the
+    number of ``\tilde{\bm{R}}`` vectors, i.e., `KRVectors.Rvectors.n_r̃vecs`.
+    For WS interpolation, it is just `KRVectors.Rvectors.n_rvecs`.
 """
 struct InterpModel{T<:Real}
     # R vectors for Fourier transform
     kRvectors::KRVectors{T}
 
+    # kpoint path for band structure
+    kpath::KPath
+
     # Hamiltonian H(R), n_wann * n_wann * n_rvecs
     H::Array{Complex{T},3}
 
-    # kpoint path for band structure
-    kpath::KPath
+    # position operator r(R), n_wann * n_wann * n_rvecs * 3
+    r::Array{Complex{T},4}
 
     # spin matrix S(R), n_wann * n_wann * n_rvecs * 3
     S::Array{Complex{T},4}
@@ -57,33 +66,62 @@ function Base.show(io::IO, model::InterpModel)
     return nothing
 end
 
-function InterpModel(kRvectors, H, kpath, S)
+function InterpModel(
+    kRvectors::KRVectors, kpath::KPath, H::Array{T,3}, r::Array{T,4}, S::Array{T,4}
+) where {T<:Complex}
     n_wann = size(H, 1)
-    return InterpModel(kRvectors, H, kpath, S, n_wann)
+    return InterpModel(kRvectors, kpath, H, r, S, n_wann)
+end
+
+function _get_nrvecs(kRvectors::KRVectors)
+    Rvecs = kRvectors.Rvectors
+    if Rvecs isa RVectorsMDRS
+        n_rvecs = Rvecs.n_r̃vecs
+    elseif Rvecs isa RVectors
+        n_rvecs = Rvecs.n_rvecs
+    else
+        error("Unknown Rvectors type: $(typeof(Rvecs))")
+    end
+    return n_rvecs
 end
 
 """
-    InterpModel(kRvectors, H, kpath)
+    InterpModel(kRvectors, H, kpath, H, r)
 
 A `InterpModel` constructor ignoring spin operator matrices.
 
 # Arguments
 - `kRvectors`: the kpoint and R vectors
+- `kpath`: the kpoint path for band structure
 - `H`: `n_wann * n_wann * n_rvecs`, the Hamiltonian in real space
+- `r`: `n_wann * n_wann * n_rvecs * 3`, the position operator in real space
+"""
+function InterpModel(
+    kRvectors::KRVectors, kpath::KPath, H::Array{T,3}, r::Array{T,4}
+) where {T<:Complex}
+    n_rvecs = _get_nrvecs(kRvectors)
+    n_wann = size(H, 1)
+    S = Array{T,4}(undef, n_wann, n_wann, n_rvecs, 3)
+
+    return InterpModel(kRvectors, kpath, H, r, S)
+end
+
+"""
+    InterpModel(kRvectors, H, kpath, H)
+
+A `InterpModel` constructor ignoring position and spin operator matrices.
+
+# Arguments
+- `kRvectors`: the kpoint and R vectors
 - `kpath`: the kpoint path for band structure
+- `H`: `n_wann * n_wann * n_rvecs`, the Hamiltonian in real space
 """
-function InterpModel(kRvectors, H, kpath)
-    Rvecs = kRvectors.Rvectors
-    if typeof(Rvecs) <: RVectorsMDRS
-        n_rvecs = Rvecs.n_r̃vecs
-    elseif typeof(Rvecs) <: RVectors
-        n_rvecs = Rvecs.n_rvecs
-    else
-        error("Unknown Rvectors type: $(typeof(Rvecs))")
-    end
+function InterpModel(kRvectors::KRVectors, kpath::KPath, H::Array{T,3}) where {T<:Complex}
+    n_rvecs = _get_nrvecs(kRvectors)
     n_wann = size(H, 1)
-    S = Array{eltype(H),4}(undef, n_wann, n_wann, n_rvecs, 3)
-    return InterpModel(kRvectors, H, kpath, S)
+    r = Array{T,4}(undef, n_wann, n_wann, n_rvecs, 3)
+    S = Array{T,4}(undef, n_wann, n_wann, n_rvecs, 3)
+    return InterpModel(kRvectors, kpath, H, r, S)
 end
 
 """
@@ -116,5 +154,5 @@ function InterpModel(model::Model; mdrs::Bool=true)
     Hᵏ = get_Hk(model.E, model.A)
     Hᴿ = fourier(kRvecs, Hᵏ)
 
-    return InterpModel(kRvecs, Hᴿ, kpath)
+    return InterpModel(kRvecs, kpath, Hᴿ)
 end

src/interp/krvector.jl

@@ -118,6 +118,27 @@ function KRVectors(
     return KRVectors(lattice, kgrid, kpoints, k_xyz, xyz_k, Rvectors)
 end
 
+"""
+    KRVectors(Rvectors)
+
+Construct a `KRVectors` from `RVectors` or `RVectorsMDRS`.
+
+The kpoints part are just empty.
+
+!!! note
+
+    If we only need to Wannier interpolate operators, e.g., interpolate band structure
+    from `tb.dat` file, then we don't need info on the kpoint grid of Wannierization,
+    so we can just leave it empty.
+"""
+function KRVectors(Rvectors::Union{RVectors{T},RVectorsMDRS{T}}) where {T<:Real}
+    kgrid = Vec3{Int}([0, 0, 0])  # no info on kgrid
+    k_xyz = Vector{Vec3{Int}}(undef, 0)
+    xyz_k = Array{Int,3}(undef, 0, 0, 0)
+    kpoints = zeros(T, 3, 0)
+    return KRVectors(Rvectors.lattice, kgrid, kpoints, k_xyz, xyz_k, Rvectors)
+end
+
 function Base.show(io::IO, x::KRVectors)
     @printf(io, "recip_lattice: Å⁻¹\n")
     for i in 1:3

src/io/w90/model.jl

@@ -169,7 +169,7 @@ function read_w90_post(
     Hᵏ = get_Hk(model.E, model.A)
     Hᴿ = fourier(kRvecs, Hᵏ)
 
-    return InterpModel(kRvecs, Hᴿ, kpath)
+    return InterpModel(kRvecs, kpath, Hᴿ)
 end
 
 """

src/io/w90/tb.jl

@@ -1,7 +1,7 @@
 export read_w90_tb
 
 """
-    read_w90_tb(seedname::AbstractString)
+    _read_w90_tb(seedname::AbstractString)
 
 Read `seedname_tb.dat` and `seedname_wsvec.dat`.
 
@@ -10,7 +10,7 @@ Read `seedname_tb.dat` and `seedname_wsvec.dat`.
 - Hamiltonian
 - position operator
 """
-function read_w90_tb(seedname::AbstractString)
+function _read_w90_tb(seedname::AbstractString)
     mdrs, wsvec = read_w90_wsvec(seedname * "_wsvec.dat")
     R = wsvec.R
     if mdrs
@@ -33,3 +33,65 @@ function read_w90_tb(seedname::AbstractString)
     Rvecs_mdrs = RVectorsMDRS(Rvecs, T, Nᵀ)
     return Rvecs_mdrs, tbdat.H, tbdat.r
 end
+
+"""
+    read_w90_tb(seedname; kpoints, atom_positions, atom_labels)
+
+Read `seedname_tb.dat` and `seedname_wsvec.dat`, return an [`InterpModel`](@ref).
+
+# Arguments
+- `seedname`: the seedname of `seedname_tb.dat` and `seedname_wsvec.dat`
+
+# Keyword Arguments
+- `kpoints`: the kpoints used in Wannierization, each column is a fractional coordinate.
+- `atom_positions`: columns are the fractional coordinates of atoms
+- `atom_labels`: labels of atoms
+
+# Return
+- An [`InterpModel`](@ref).
+
+!!! note
+
+    Usually, if you only need to interpolate operators (e.g., the band structure),
+    that is only inverse Fourier transform [`invfourier`](@ref) is needed,
+    then you don't need to provide `kpoints`.
+    If in some cases, you want to generate Wannier gauge operators, then passing
+    the `kpoints` allows you to construct a complete [`KRVectors`](@ref),
+    and you can subsequently call [`fourier`](@ref) on Bloch gauge operators
+    to get Wannier gauge operators.
+
+!!! warning
+
+    Since no atomic positions and labels are provided, the auto-generated `KPath`
+    is probably wrong.
+    It is recommended that you pass the `atom_positions` and `atom_labels` arguments
+    so that this function can auto generate the correct [`KPath`](@ref).
+"""
+function read_w90_tb(
+    seedname::AbstractString;
+    kpoints::Union{Nothing,AbstractMatrix}=nothing,
+    atom_positions::Union{Nothing,AbstractMatrix}=nothing,
+    atom_labels::Union{Nothing,AbstractVector{String}}=nothing,
+)
+    Rvecs, H, r = _read_w90_tb(seedname)
+    if kpoints === nothing
+        kRvecs = KRVectors(Rvecs)
+    else
+        size(kpoints, 1) == 3 || error("kpoints should be 3 * n_kpts")
+        kgrid = Vec3{Int}(get_kgrid(kpoints))
+        kRvecs = KRVectors(Rvecs.lattice, kgrid, kpoints, Rvecs)
+    end
+
+    if atom_positions === nothing || atom_labels === nothing
+        # generate a fake atom
+        atom_positions = zeros(3, 1)
+        atom_labels = ["H"]
+    end
+    kpath = get_kpath(Rvecs.lattice, atom_positions, atom_labels)
+
+    n_wann = size(H, 1)
+    n_rvecs = size(H, 3)
+    S = Array{ComplexF64,4}(undef, n_wann, n_wann, n_rvecs, 3)
+
+    return InterpModel(kRvecs, kpath, H, r, S)
+end

test/interp/fermisurf.jl

@@ -1,10 +1,10 @@
 
 @testset "Fermi surface" begin
-    Rvecs, H, r = read_w90_tb(joinpath(FIXTURE_PATH, "valence/band/mdrs/silicon"))
+    model = read_w90_tb(joinpath(FIXTURE_PATH, "valence/band/mdrs/silicon"))
 
     ref_bxsf = WannierIO.read_bxsf(joinpath(FIXTURE_PATH, "valence/band/mdrs/wjl.bxsf"))
 
-    kpoints, E = Wannier.fermi_surface(Rvecs, H; n_k=2)
+    kpoints, E = Wannier.fermi_surface(model; n_k=2)
 
     @test all(isapprox.(E, ref_bxsf.E; atol=1e-7))
 end

test/interp/fourier.jl

@@ -1,34 +1,36 @@
 model = read_w90(joinpath(FIXTURE_PATH, "valence/band/silicon"); amn=false)
 model.A .= get_A(read_chk(joinpath(FIXTURE_PATH, "valence/band/silicon.chk.fmt")))
-model_ws = read_w90_post(joinpath(FIXTURE_PATH, "valence/band/silicon"); mdrs=false)
-model_mdrs = read_w90_post(joinpath(FIXTURE_PATH, "valence/band/silicon"); mdrs=true)
-tb_ws = read_w90_tb(joinpath(FIXTURE_PATH, "valence/band/ws/silicon"))
-tb_mdrs = read_w90_tb(joinpath(FIXTURE_PATH, "valence/band/mdrs/silicon"))
+model_ws = read_w90_tb(
+    joinpath(FIXTURE_PATH, "valence/band/ws/silicon"); kpoints=model.kpoints
+)
+model_mdrs = read_w90_tb(
+    joinpath(FIXTURE_PATH, "valence/band/mdrs/silicon"); kpoints=model.kpoints
+)
 
 @testset "fourier WS" begin
     Hᵏ = Wannier.get_Hk(model.E, model.A)
     Hᴿ = Wannier.fourier(model_ws.kRvectors, Hᵏ)
-    Rvecs, ref_Hᴿ, positions = tb_ws
+    ref_Hᴿ = model_ws.H
     @test all(isapprox.(Hᴿ, ref_Hᴿ; atol=1e-7))
 end
 
 @testset "invfourier WS" begin
     ref_Hᵏ = Wannier.get_Hk(model.E, model.A)
-    Rvecs, Hᴿ, positions = tb_ws
+    Hᴿ = model_ws.H
     Hᵏ = Wannier.invfourier(model_ws.kRvectors, Hᴿ, model.kpoints)
     @test all(isapprox.(Hᵏ, ref_Hᵏ; atol=1e-7))
 end
 
 @testset "fourier MDRS v1" begin
     Hᵏ = Wannier.get_Hk(model.E, model.A)
     Hᴿ = Wannier.fourier(model_mdrs.kRvectors, Hᵏ; version=:v1)
-    Rvecs, ref_Hᴿ, positions = tb_mdrs
+    ref_Hᴿ = model_mdrs.H
     @test all(isapprox.(Hᴿ, ref_Hᴿ; atol=1e-7))
 end
 
 @testset "invfourier MDRS v1" begin
     ref_Hᵏ = Wannier.get_Hk(model.E, model.A)
-    Rvecs, Hᴿ, positions = tb_mdrs
+    Hᴿ = model_mdrs.H
     Hᵏ = Wannier.invfourier(model_mdrs.kRvectors, Hᴿ, model.kpoints; version=:v1)
     @test all(isapprox.(Hᵏ, ref_Hᵏ; atol=1e-7))
 end