Update conversion with shconv and manual sparsity

Author: vitartas

Project.toml

@@ -40,9 +40,11 @@ HDF5 = "0.17"
 LoggingExtras = "1.1"
 MPI = "0.20"
 NeighbourLists = "0.5"
+OffsetArrays = "1.17"
 Reexport = "1.2"
 SafeTestsets = "0.1"
 StaticArrays = "1.9"
+StructArrays = "0.7"
 Unitful = "1.25"
 UnitfulAtomic = "1"
 julia = "1.9"

bin/quoll.jl

@@ -62,6 +62,24 @@ for idir in my_idirs
 
     # TODO: in the future could allow for <other canonical formats/shortcut conversions>
     @info "Converting operators into canonical format"
-    # operators = Dict(Quoll.RealBSparseOperator, operators)
-    
-end
+
+    operators = Dict([
+        operatorkind => Quoll.RealBSparseOperator(operator, params.input.radii)
+        for (operatorkind, operator) in operators
+    ])
+
+    # TODO: include a function here which would initialise the k-point grid if its required
+    # so that we don't have to do logic everywhere else multiple times
+
+    # Basis projection
+    # - requires k-point grid
+    # - incompatible with error metrics (TODO: probably should throw an error but I could do that at parse time)
+
+    # Postprocessing
+    # - most of them require k-point grid (except band structure)
+    # - should store eigenvalues, fermi level and pass them to error metrics if required
+
+    # Error metrics
+    # - some require k-point grid and eigenvalues (e.g. eigenvalue error)
+
+end

examples/SiC/input_file.toml

@@ -5,10 +5,10 @@
 format = "FHI-aims"
 directories = ["SiC_FHIaims"]
 operators = ["H", "S"]
-radii = [
-    "C 12.0 ang",
-    "Si 12.0 ang",
-]
+# radii = [
+#     "C 12.0 ang",
+#     "Si 12.0 ang",
+# ]
 
 # Convert
 # - Output format

src/conversions/fhiaims/fhiaims_csc.jl

@@ -1,19 +1,25 @@
 
 # TODO: might need to write a more specialized function which doesn't do conversions for
 # each operator separately assuming the metadata is the same
-function RealBSparseOperator(in_operator::FHIaimsCSCOperator)
-    ## Metadata
-    # Convert sparsity
-    out_sparsity = RealBlockSparsity(in_operator.metadata.sparsity, in_operator.metadata.basisset)
+# TODO: Allow to compute sparsity from radii
+function RealBSparseOperator(in_operator::FHIaimsCSCOperator, radii)
+
+    # Obtain sparsity
+    if isnothing(radii)
+        out_sparsity = RealBlockSparsity(in_operator.metadata.sparsity, in_operator.metadata.basisset)
+    else
+        @info "Computing sparsity manually from radii"
+        out_sparsity = RealBlockSparsity(in_operator.metadata.atoms, radii, hermitian = true)
+    end
 
     # Compute z1z2_ij2interval
     z1z2_ij2interval = compute_z1z2_ij2interval(in_operator.metadata.atoms, out_sparsity)
 
+    # Construct metadata
     out_metadata = RealBSparseMetadata(
         in_operator.metadata.atoms, out_sparsity, in_operator.metadata.basisset, in_operator.metadata.spins, z1z2_ij2interval
     )
 
-    ## Data
     # Initialize out_operator with zeros
     out_operator = RealBSparseOperator(in_operator.kind, out_metadata)
 
@@ -41,9 +47,11 @@ function populate!(out_keydata, out_sparsity, out_basisset, in_data, in_sparsity
 
     atom2basis_offset = get_offsets(out_basisset)
     iglobal2ilocal = get_iglobal2ilocal(out_sparsity)
+    shconv = SHConversion(out_type) ∘ inv(SHConversion(in_type))
 
-    # shconv = SHConversion(out_type) ∘ inv(SHConversion(in_type))
-    # shifts, phases = precompute_shiftphases(out_basisset, shconv)
+    # Use inv(shconv) for shifts and phases because matrix is being reordered
+    # using `type2` reordering, for more information see tests/unit/shconversion.jl
+    shifts, phases = precompute_shiftphases(out_basisset, inv(shconv))
 
     for i_cell in axes(in_sparsity.colcellptr, 2)
         image = in_sparsity.images[i_cell]
@@ -65,24 +73,31 @@ function populate!(out_keydata, out_sparsity, out_basisset, in_data, in_sparsity
                 i_basis_row_local = i_basis_row - atom2basis_offset[i_atom_row]
                 i_basis_col_local = i_basis_col - atom2basis_offset[i_atom_col]
 
-                out_keydata[(i_atom_row, i_atom_col)][i_basis_row_local, i_basis_col_local, i_cell_local] = in_data[i_index]
+                out_keydata[(i_atom_row, i_atom_col)][
+                    i_basis_row_local + shifts[out_basisset.atom2species[i_atom_row]][i_basis_row_local],
+                    i_basis_col_local + shifts[out_basisset.atom2species[i_atom_col]][i_basis_col_local],
+                    i_cell_local
+                ] = (
+                    in_data[i_index]
+                    * phases[out_basisset.atom2species[i_atom_row]][i_basis_row_local]
+                    * phases[out_basisset.atom2species[i_atom_col]][i_basis_col_local]
+                )
             end
         end
     end
-end
 
-# struct FHIaimsCSCMetadata{A<:AbstractSystem, E} <: AbstractFHIaimsMetadata
-#     atoms::A
-#     sparsity::RealCSCSparsity
-#     basisset::BasisSetMetadata{E}
-#     spinset::Union{SpinsMetadata, Nothing}
-#     # TODO: Is Union{SpinsMetadata, Nothing} best approach here?
-#     # Making FHIaimsCSCMetadata a parametric type wrt typeof(spins)
-#     # or making SpinsMetadata a parametric type might be an overkill
-# end
-
-# struct FHIaimsCSCOperator{O<:AbstractOperatorKind, T<:AbstractFloat, A<:AbstractSystem, E} <: AbstractFHIaimsOperator
-#     kind::O
-#     data::Vector{T}
-#     metadata::FHIaimsCSCMetadata{A, E}
-# end
+    # Fill in the 'L' part of on-site blocks that, due of Hermitian CSC sparsity,
+    # were not filled in the loop above
+    n_atoms = length(atom2basis_offset)
+    for i_atom in 1:n_atoms
+        for k in axes(out_keydata[(i_atom, i_atom)], 3)
+            for j in axes(out_keydata[(i_atom, i_atom)], 2)
+                @inbounds for i in axes(out_keydata[(i_atom, i_atom)], 1)
+                    j > i || continue
+                    out_keydata[(i_atom, i_atom)][j, i, k] = out_keydata[(i_atom, i_atom)][i, j, k]
+                end
+            end
+        end
+    end
+
+end

src/parser/input_output.jl

@@ -95,7 +95,7 @@ function parse_radius(radius_string)
     r = parse(Float64, radius_match.captures[2])
     unit = radius_match.captures[3]
 
-    if unit === nothing
+    if isnothing(unit)
         @warn "Units to radii field not supplied, assuming Å"
         r = r * u"Å"
     else

src/sparsity.jl

@@ -104,6 +104,7 @@ function RealBlockSparsity(colcellptr, rowval, images, basis2atom)
     return RealBlockSparsity(ij2images, images, true)
 end
 
+# A map from image indices in `images` to local image indices in `ij2images` for a given ij
 function get_iglobal2ilocal(sparsity::RealBlockSparsity)
     return Dictionary(
         keys(sparsity.ij2images),

test/unit/shconversion.jl

@@ -3,7 +3,7 @@ using AtomsBase
 
 using Test
 
-function reorder_matrix(A1, basis, shconv)
+function reorder_matrix_type1(A1, basis, shconv)
     A2 = zero(A1)
 
     # Precompute shifts and phases
@@ -30,6 +30,34 @@ function reorder_matrix(A1, basis, shconv)
     return A2
 end
 
+function reorder_matrix_type2(A1, basis, inv_shconv)
+    A2 = zero(A1)
+
+    # Precompute shifts and phases
+    n = size(A1, 2)
+    shifts = Vector{Int}(undef, n)
+    phases = Vector{Int}(undef, n)
+
+    for j in 1:n
+        shifts[j], phases[j] = Quoll.get_shiftphase(basis[j], inv_shconv)
+    end
+
+    for j in 1:n
+        shift_j = shifts[j]
+        phase_j = phases[j]
+
+        for i in 1:n
+            shift_i = shifts[i]
+            phase_i = phases[i]
+
+            # Note shifts on A2 instead of A1 as in type1 conversion
+            A2[i + shift_i, j + shift_j] = A1[i, j] * phase_i * phase_j
+        end
+    end
+
+    return A2
+end
+
 
 @testset "SHConversion" begin
     orders = [[1], [3,  1,  2], [1,  4,  5,  3,  2]]
@@ -89,7 +117,7 @@ end
     shconv2 = Quoll.SHConversion(orders2, phases2)
 
     @testset "Basic" begin
-        A2 = reorder_matrix(A1, basis, shconv1)
+        A2 = reorder_matrix_type1(A1, basis, shconv1)
         @test A2 == [
             0  2 -1  3;
             2  4 -3  5;
@@ -106,19 +134,23 @@ end
             2 3 5 4;
         ]
 
-        A2 = reorder_matrix(A1, basis, shconv1)
-        A3 = reorder_matrix(A2, basis, shconv2)
+        A2 = reorder_matrix_type1(A1, basis, shconv1)
+        A3 = reorder_matrix_type1(A2, basis, shconv2)
         @test A3 == ref
 
-        A3_direct = reorder_matrix(A1, basis, shconv2 ∘ shconv1)
+        A3_direct = reorder_matrix_type1(A1, basis, shconv2 ∘ shconv1)
         @test A3_direct == ref
     end
 
     @testset "Inverse" begin
-        A2 = reorder_matrix(A1, basis, shconv1)
-        A3 = reorder_matrix(A2, basis, inv(shconv1))
+        A2 = reorder_matrix_type1(A1, basis, shconv1)
+        A3 = reorder_matrix_type1(A2, basis, inv(shconv1))
         @test A1 == A3
 
+        A2type1 = reorder_matrix_type1(A1, basis, shconv1)
+        A2type2 = reorder_matrix_type2(A1, basis, inv(shconv1))
+        @test A2type1 == A2type2
+
         identity_shconv = inv(shconv1) ∘ shconv1
 
         l = 0